X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FMEDCoupling%2FMEDCouplingRemapper.cxx;h=aaffd1d14e5364d087498eaaa5e804fe5d045230;hb=378cb2ebe08f8f4543ef632b2bd5f77fe180f978;hp=8700012b043616331f1f8c996bef749acc242ce2;hpb=8763c12d01e33d6845dd53be65b001514d00bd42;p=tools%2Fmedcoupling.git diff --git a/src/MEDCoupling/MEDCouplingRemapper.cxx b/src/MEDCoupling/MEDCouplingRemapper.cxx old mode 100644 new mode 100755 index 8700012b0..aaffd1d14 --- a/src/MEDCoupling/MEDCouplingRemapper.cxx +++ b/src/MEDCoupling/MEDCouplingRemapper.cxx @@ -81,7 +81,7 @@ int MEDCouplingRemapper::prepareInterpKernelOnly() { int meshInterpType=((int)_src_ft->getMesh()->getType()*16)+(int)_target_ft->getMesh()->getType(); switch(meshInterpType) - { + { case 90: case 91: case 165: @@ -106,7 +106,7 @@ int MEDCouplingRemapper::prepareInterpKernelOnly() return prepareInterpKernelOnlyEE(); default: throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnly : Not managed type of meshes ! Dealt meshes type are : Unstructured<->Unstructured, Unstructured<->Cartesian, Cartesian<->Cartesian, Extruded<->Extruded !"); - } + } } int MEDCouplingRemapper::prepareNotInterpKernelOnly() @@ -114,7 +114,7 @@ int MEDCouplingRemapper::prepareNotInterpKernelOnly() std::string srcm,trgm,method; method=checkAndGiveInterpolationMethodStr(srcm,trgm); switch(CheckInterpolationMethodManageableByNotOnlyInterpKernel(method)) - { + { case 0: return prepareNotInterpKernelOnlyGaussGauss(); default: @@ -122,7 +122,7 @@ int MEDCouplingRemapper::prepareNotInterpKernelOnly() std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnly : INTERNAL ERROR ! the method \"" << method << "\" declared as managed bu not implemented !"; throw INTERP_KERNEL::Exception(oss.str().c_str()); } - } + } } /*! @@ -134,6 +134,8 @@ int MEDCouplingRemapper::prepareNotInterpKernelOnly() */ void MEDCouplingRemapper::transfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, double dftValue) { + if(!srcField || !targetField) + throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transfer : input field must be both not NULL !"); transferUnderground(srcField,targetField,true,dftValue); } @@ -148,34 +150,46 @@ void MEDCouplingRemapper::transfer(const MEDCouplingFieldDouble *srcField, MEDCo */ void MEDCouplingRemapper::partialTransfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField) { + if(!srcField || !targetField) + throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : input field must be both not NULL !"); transferUnderground(srcField,targetField,false,std::numeric_limits::max()); } void MEDCouplingRemapper::reverseTransfer(MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *targetField, double dftValue) { + if(!srcField || !targetField) + throw INTERP_KERNEL::Exception("MEDCouplingRemapper::reverseTransfer : input fields must be both not NULL !"); checkPrepare(); + targetField->checkCoherency(); if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr()) throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field"); if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr()) throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field"); if(srcField->getNature()!=targetField->getNature()) throw INTERP_KERNEL::Exception("Natures of fields mismatch !"); - DataArrayDouble *array=srcField->getArray(); + if(targetField->getNumberOfTuplesExpected()!=_target_ft->getNumberOfTuplesExpected()) + { + std::ostringstream oss; + oss << "MEDCouplingRemapper::reverseTransfer : in given source field the number of tuples required is " << _target_ft->getNumberOfTuplesExpected() << " (on prepare) and number of tuples in given target field is " << targetField->getNumberOfTuplesExpected(); + oss << " ! It appears that the target support is not the same between the prepare and the transfer !"; + throw INTERP_KERNEL::Exception(oss.str().c_str()); + } + DataArrayDouble *array(srcField->getArray()); int trgNbOfCompo=targetField->getNumberOfComponents(); if(array) { - if(trgNbOfCompo!=srcField->getNumberOfComponents()) + srcField->checkCoherency(); + if(trgNbOfCompo!=srcField->getNumberOfTuplesExpected()) throw INTERP_KERNEL::Exception("Number of components mismatch !"); } else { - array=DataArrayDouble::New(); - array->alloc(srcField->getNumberOfTuples(),trgNbOfCompo); - srcField->setArray(array); - array->decrRef(); + MEDCouplingAutoRefCountObjectPtr tmp(DataArrayDouble::New()); + tmp->alloc(srcField->getNumberOfTuplesExpected(),trgNbOfCompo); + srcField->setArray(tmp); } computeDeno(srcField->getNature(),srcField,targetField); - double *resPointer=array->getPointer(); + double *resPointer(srcField->getArray()->getPointer()); const double *inputPointer=targetField->getArray()->getConstPointer(); computeReverseProduct(inputPointer,trgNbOfCompo,dftValue,resPointer); } @@ -183,6 +197,9 @@ void MEDCouplingRemapper::reverseTransfer(MEDCouplingFieldDouble *srcField, cons MEDCouplingFieldDouble *MEDCouplingRemapper::transferField(const MEDCouplingFieldDouble *srcField, double dftValue) { checkPrepare(); + if(!srcField) + throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferField : input srcField is NULL !"); + srcField->checkCoherency(); if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr()) throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field"); MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(*_target_ft,srcField->getTimeDiscretization()); @@ -194,6 +211,9 @@ MEDCouplingFieldDouble *MEDCouplingRemapper::transferField(const MEDCouplingFiel MEDCouplingFieldDouble *MEDCouplingRemapper::reverseTransferField(const MEDCouplingFieldDouble *targetField, double dftValue) { + if(!targetField) + throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferField : input targetField is NULL !"); + targetField->checkCoherency(); checkPrepare(); if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr()) throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field"); @@ -281,7 +301,7 @@ int MEDCouplingRemapper::getInterpolationMatrixPolicy() const void MEDCouplingRemapper::setInterpolationMatrixPolicy(int newInterpMatPol) { switch(newInterpMatPol) - { + { case 0: _interp_matrix_pol=IK_ONLY_PREFERED; break; @@ -296,7 +316,7 @@ void MEDCouplingRemapper::setInterpolationMatrixPolicy(int newInterpMatPol) break; default: throw INTERP_KERNEL::Exception("MEDCouplingRemapper::setInterpolationMatrixPolicy : invalid input integer value ! Should be in [0 (IK_PREFERED) , 1 (NOT_IK_PREFERED), 2 (IK_ONLY_FORCED), 3 (NOT_IK_ONLY_FORCED)] ! For information, the default is IK_PREFERED=0 !"); - } + } } int MEDCouplingRemapper::prepareInterpKernelOnlyUU() @@ -304,7 +324,7 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUU() const MEDCouplingPointSet *src_mesh=static_cast(_src_ft->getMesh()); const MEDCouplingPointSet *target_mesh=static_cast(_target_ft->getMesh()); std::string srcMeth,trgMeth; - std::string method=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth); + std::string method(checkAndGiveInterpolationMethodStr(srcMeth,trgMeth)); const int srcMeshDim=src_mesh->getMeshDimension(); int srcSpaceDim=-1; if(srcMeshDim!=-1) @@ -369,7 +389,8 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUU() MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh); INTERP_KERNEL::Interpolation3D interpolation(*this); std::vector > matrixTmp; - nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method); + std::string revMethod(BuildMethodFrom(trgMeth,srcMeth)); + nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod); ReverseMatrix(matrixTmp,nbCols,_matrix); nbCols=matrixTmp.size(); } @@ -388,7 +409,8 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUU() MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh); INTERP_KERNEL::Interpolation2D1D interpolation(*this); std::vector > matrixTmp; - nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method); + std::string revMethod(BuildMethodFrom(trgMeth,srcMeth)); + nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod); ReverseMatrix(matrixTmp,nbCols,_matrix); nbCols=matrixTmp.size(); INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces(); @@ -412,7 +434,8 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUU() MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh); INTERP_KERNEL::Interpolation2D interpolation(*this); std::vector > matrixTmp; - nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method); + std::string revMethod(BuildMethodFrom(trgMeth,srcMeth)); + nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod); ReverseMatrix(matrixTmp,nbCols,_matrix); nbCols=matrixTmp.size(); } @@ -459,7 +482,8 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUU() MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh); INTERP_KERNEL::Interpolation3D2D interpolation(*this); std::vector > matrixTmp; - nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method); + std::string revMethod(BuildMethodFrom(trgMeth,srcMeth)); + nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod); ReverseMatrix(matrixTmp,nbCols,_matrix); nbCols=matrixTmp.size(); INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces(); @@ -579,7 +603,7 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUC() throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : space dim of src unstructured should be equal to mesh dim of src unstructured and should be equal also equal to trg cartesian dimension !"); std::vector > res; switch(srcMeshDim) - { + { case 1: { MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh); @@ -606,7 +630,7 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUC() } default: throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only dimension 1 2 or 3 supported !"); - } + } ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix); nullifiedTinyCoeffInCrudeMatrixAbs(0.); // @@ -632,7 +656,7 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyCU() if(trgMeshDim!=trgSpceDim || trgMeshDim!=srcMeshDim) throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : space dim of target unstructured should be equal to mesh dim of target unstructured and should be equal also equal to source cartesian dimension !"); switch(srcMeshDim) - { + { case 1: { MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh); @@ -659,7 +683,7 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyCU() } default: throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only dimension 1 2 or 3 supported !"); - } + } nullifiedTinyCoeffInCrudeMatrixAbs(0.); // _deno_multiply.clear(); @@ -683,7 +707,7 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyCC() if(trgMeshDim!=srcMeshDim) throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : dim of target cartesian should be equal to dim of source cartesian dimension !"); switch(srcMeshDim) - { + { case 1: { MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh); @@ -710,7 +734,7 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyCC() } default: throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only dimension 1 2 or 3 supported !"); - } + } nullifiedTinyCoeffInCrudeMatrixAbs(0.); // _deno_multiply.clear(); @@ -724,7 +748,7 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyCC() int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss() { if(getIntersectionType()!=INTERP_KERNEL::PointLocator) - throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss : The intersection type is not supported ! Only PointLocator is supported for Gauss->Gauss interpolation ! Please invoke setIntersectionType(PointLocator) on the MEDCouplingRemapper instance !"); + throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss : The intersection type is not supported ! Only PointLocator is supported for Gauss->Gauss interpolation ! Please invoke setIntersectionType(PointLocator) on the MEDCouplingRemapper instance !"); MEDCouplingAutoRefCountObjectPtr trgLoc=_target_ft->getLocalizationOfDiscr(); const double *trgLocPtr=trgLoc->begin(); int trgSpaceDim=trgLoc->getNumberOfComponents(); @@ -804,31 +828,31 @@ bool MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly() const std::string srcm,trgm,method; method=checkAndGiveInterpolationMethodStr(srcm,trgm); switch(_interp_matrix_pol) - { + { case IK_ONLY_PREFERED: { try - { + { std::string tmp1,tmp2; INTERP_KERNEL::Interpolation::CheckAndSplitInterpolationMethod(method,tmp1,tmp2); return true; - } + } catch(INTERP_KERNEL::Exception& /*e*/) - { + { return false; - } + } } case NOT_IK_ONLY_PREFERED: { try - { + { CheckInterpolationMethodManageableByNotOnlyInterpKernel(method); return false; - } + } catch(INTERP_KERNEL::Exception& /*e*/) - { + { return true; - } + } } case IK_ONLY_FORCED: return true; @@ -836,7 +860,7 @@ bool MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly() const return false; default: throw INTERP_KERNEL::Exception("MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly : internal error ! The interpolation matrix policy is not managed ! Try to change it using MEDCouplingRemapper::setInterpolationMatrixPolicy !"); - } + } } void MEDCouplingRemapper::updateTime() const @@ -869,7 +893,12 @@ std::string MEDCouplingRemapper::checkAndGiveInterpolationMethodStr(std::string& throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have their mesh set !"); srcMeth=_src_ft->getDiscretization()->getRepr(); trgMeth=_target_ft->getDiscretization()->getRepr(); - std::string method(srcMeth); method+=trgMeth; + return BuildMethodFrom(srcMeth,trgMeth); +} + +std::string MEDCouplingRemapper::BuildMethodFrom(const std::string& meth1, const std::string& meth2) +{ + std::string method(meth1); method+=meth2; return method; } @@ -887,6 +916,9 @@ void MEDCouplingRemapper::releaseData(bool matrixSuppression) void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue) { + if(!srcField || !targetField) + throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferUnderground : srcField or targetField is NULL !"); + srcField->checkCoherency(); checkPrepare(); if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr()) throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field"); @@ -894,10 +926,18 @@ void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcF throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field"); if(srcField->getNature()!=targetField->getNature()) throw INTERP_KERNEL::Exception("Natures of fields mismatch !"); - DataArrayDouble *array=targetField->getArray(); - int srcNbOfCompo=srcField->getNumberOfComponents(); + if(srcField->getNumberOfTuplesExpected()!=_src_ft->getNumberOfTuplesExpected()) + { + std::ostringstream oss; + oss << "MEDCouplingRemapper::transferUnderground : in given source field the number of tuples required is " << _src_ft->getNumberOfTuplesExpected() << " (on prepare) and number of tuples in given source field is " << srcField->getNumberOfTuplesExpected(); + oss << " ! It appears that the source support is not the same between the prepare and the transfer !"; + throw INTERP_KERNEL::Exception(oss.str().c_str()); + } + DataArrayDouble *array(targetField->getArray()); + int srcNbOfCompo(srcField->getNumberOfComponents()); if(array) { + targetField->checkCoherency(); if(srcNbOfCompo!=targetField->getNumberOfComponents()) throw INTERP_KERNEL::Exception("Number of components mismatch !"); } @@ -905,14 +945,13 @@ void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcF { if(!isDftVal) throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !"); - array=DataArrayDouble::New(); - array->alloc(targetField->getNumberOfTuples(),srcNbOfCompo); - targetField->setArray(array); - array->decrRef(); + MEDCouplingAutoRefCountObjectPtr tmp(DataArrayDouble::New()); + tmp->alloc(targetField->getNumberOfTuples(),srcNbOfCompo); + targetField->setArray(tmp); } computeDeno(srcField->getNature(),srcField,targetField); - double *resPointer=array->getPointer(); - const double *inputPointer=srcField->getArray()->getConstPointer(); + double *resPointer(targetField->getArray()->getPointer()); + const double *inputPointer(srcField->getArray()->getConstPointer()); computeProduct(inputPointer,srcNbOfCompo,isDftVal,dftValue,resPointer); } @@ -921,7 +960,7 @@ void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldD if(nat==NoNature) return computeDenoFromScratch(nat,srcField,trgField); else if(nat!=_nature_of_deno) - return computeDenoFromScratch(nat,srcField,trgField); + return computeDenoFromScratch(nat,srcField,trgField); else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis()) return computeDenoFromScratch(nat,srcField,trgField); } @@ -931,7 +970,7 @@ void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCou _nature_of_deno=nat; _time_deno_update=getTimeOfThis(); switch(_nature_of_deno) - { + { case ConservativeVolumic: { ComputeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply); @@ -998,7 +1037,7 @@ void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCou } case NoNature: throw INTERP_KERNEL::Exception("No nature specified ! Select one !"); - } + } } void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)