-// Copyright (C) 2007-2013 CEA/DEN, EDF R&D
+// Copyright (C) 2007-2015 CEA/DEN, EDF R&D
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
-// version 2.1 of the License.
+// version 2.1 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
releaseData(false);
}
+/*!
+ * This method is the second step of the remapping process. The remapping process works in three phases :
+ *
+ * - Set remapping options appropriately
+ * - The computation of remapping matrix
+ * - Apply the matrix vector multiply to obtain the result of the remapping
+ *
+ * This method performs the second step (computation of remapping matrix) which may be CPU-time consuming. This phase is also the most critical (where the most tricky algorithm) in the remapping process.
+ * Strictly speaking to perform the computation of the remapping matrix the field templates source-side and target-side is required (which is the case of MEDCouplingRemapper::prepareEx).
+ * So this method is less precise but a more user friendly way to compute a remapping matrix.
+ *
+ * \param [in] srcMesh the source mesh
+ * \param [in] targetMesh the target mesh
+ * \param [in] method A string obtained by aggregation of the spatial discretisation string representation of source field and target field. The string representation is those returned by MEDCouplingFieldDiscretization::getStringRepr.
+ * Example : "P0" is for cell discretization. "P1" is for node discretization. So "P0P1" for \a method parameter means from a source cell field (lying on \a srcMesh) to a target node field (lying on \a targetMesh).
+ *
+ * \sa MEDCouplingRemapper::prepareEx
+ */
int MEDCouplingRemapper::prepare(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method)
{
if(!srcMesh || !targetMesh)
return prepareEx(src,target);
}
+/*!
+ * This method is the generalization of MEDCouplingRemapper::prepare. Indeed, MEDCouplingFieldTemplate instances gives all required information to compute the remapping matrix.
+ * This method must be used instead of MEDCouplingRemapper::prepare if Gauss point to Gauss point must be applied.
+ *
+ * \param [in] src is the field template source side.
+ * \param [in] target is the field template target side.
+ *
+ * \sa MEDCouplingRemapper::prepare
+ */
int MEDCouplingRemapper::prepareEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target)
{
if(!src || !target)
{
int meshInterpType=((int)_src_ft->getMesh()->getType()*16)+(int)_target_ft->getMesh()->getType();
switch(meshInterpType)
- {
+ {
case 90:
case 91:
case 165:
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()
std::string srcm,trgm,method;
method=checkAndGiveInterpolationMethodStr(srcm,trgm);
switch(CheckInterpolationMethodManageableByNotOnlyInterpKernel(method))
- {
+ {
case 0:
return prepareNotInterpKernelOnlyGaussGauss();
default:
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());
}
- }
+ }
}
/*!
* If meshes of \b srcField and \b targetField do not match exactly those given into \ref ParaMEDMEM::MEDCouplingRemapper::prepare "prepare method" an exception will be thrown.
*
* \param [in] srcField is the source field from which the interpolation will be done. The mesh into \b srcField should be the same than those specified on ParaMEDMEM::MEDCouplingRemapper::prepare.
- * \param [out] targetField the destination field with the allocated array in which all tuples will be overwritten.
+ * \param [in/out] targetField the destination field with the allocated array in which all tuples will be overwritten.
+ * \param [in] dftValue is the value that will be assigned in the targetField to each entity of target mesh (entity depending on the method selected on prepare invocation) that is not intercepted by any entity of source mesh.
+ * For example in "P0P0" case (cell-cell) if a cell in target mesh is not overlapped by any source cell the \a dftValue value will be attached on that cell in the returned \a targetField. In some cases a target
+ * cell not intercepted by any source cell is a bug so in this case it is advised to set a huge value (1e300 for example) to \a dftValue to quickly point to the problem. But for users doing parallelism a target cell can
+ * be intercepted by a source cell on a different process. In this case 0. assigned to \a dftValue is more appropriate.
+ *
+ * \sa transferField
*/
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);
}
*/
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<double>::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<DataArrayDouble > 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);
}
+/*!
+ * This method performs the operation source to target using matrix computed in ParaMEDMEM::MEDCouplingRemapper::prepare method.
+ * If mesh of \b srcField does not match exactly those given into \ref ParaMEDMEM::MEDCouplingRemapper::prepare "prepare method" an exception will be thrown.
+ *
+ * \param [in] srcField is the source field from which the interpolation will be done. The mesh into \b srcField should be the same than those specified on ParaMEDMEM::MEDCouplingRemapper::prepare.
+ * \param [in] dftValue is the value that will be assigned in the targetField to each entity of target mesh (entity depending on the method selected on prepare invocation) that is not intercepted by any entity of source mesh.
+ * For example in "P0P0" case (cell-cell) if a cell in target mesh is not overlapped by any source cell the \a dftValue value will be attached on that cell in the returned \a targetField. In some cases a target
+ * cell not intercepted by any source cell is a bug so in this case it is advised to set a huge value (1e300 for example) to \a dftValue to quickly point to the problem. But for users doing parallelism a target cell can
+ * be intercepted by a source cell on a different process. In this case 0. assigned to \a dftValue is more appropriate.
+ * \return destination field to be deallocated by the caller.
+ *
+ * \sa transfer
+ */
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());
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");
void MEDCouplingRemapper::setInterpolationMatrixPolicy(int newInterpMatPol)
{
switch(newInterpMatPol)
- {
+ {
case 0:
_interp_matrix_pol=IK_ONLY_PREFERED;
break;
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()
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<std::map<int,double> > res;
switch(srcMeshDim)
- {
+ {
case 1:
{
MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
}
default:
throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only dimension 1 2 or 3 supported !");
- }
+ }
ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix);
nullifiedTinyCoeffInCrudeMatrixAbs(0.);
//
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);
}
default:
throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only dimension 1 2 or 3 supported !");
- }
+ }
nullifiedTinyCoeffInCrudeMatrixAbs(0.);
//
_deno_multiply.clear();
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);
}
default:
throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only dimension 1 2 or 3 supported !");
- }
+ }
nullifiedTinyCoeffInCrudeMatrixAbs(0.);
//
_deno_multiply.clear();
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<DataArrayDouble> trgLoc=_target_ft->getLocalizationOfDiscr();
const double *trgLocPtr=trgLoc->begin();
int trgSpaceDim=trgLoc->getNumberOfComponents();
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<INTERP_KERNEL::Interpolation3D>::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;
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
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");
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 !");
}
{
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<DataArrayDouble> 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);
}
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);
}
_nature_of_deno=nat;
_time_deno_update=getTimeOfThis();
switch(_nature_of_deno)
- {
+ {
case ConservativeVolumic:
{
ComputeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
}
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)
