//
// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//
+// Author : Anthony Geay (CEA/DEN)
#include "MEDCouplingRemapper.hxx"
#include "MEDCouplingMemArray.hxx"
#include "MEDCouplingFieldTemplate.hxx"
#include "MEDCouplingFieldDiscretization.hxx"
#include "MEDCouplingExtrudedMesh.hxx"
+#include "MEDCouplingCMesh.hxx"
#include "MEDCouplingNormalizedUnstructuredMesh.txx"
+#include "MEDCouplingNormalizedCartesianMesh.txx"
#include "Interpolation1D.txx"
#include "Interpolation2DCurve.hxx"
#include "Interpolation3DSurf.hxx"
#include "Interpolation2D1D.txx"
#include "Interpolation3D2D.txx"
+#include "InterpolationCU.txx"
+#include "InterpolationCC.txx"
using namespace ParaMEDMEM;
{
case 85://Unstructured-Unstructured
return prepareUU(method);
+ case 87://Unstructured-Cartesian
+ return prepareUC(method);
+ case 117://Cartesian-Unstructured
+ return prepareCU(method);
+ case 119://Cartesian-Cartesian
+ return prepareCC(method);
case 136://Extruded-Extruded
return prepareEE(method);
default:
INTERP_KERNEL::Interpolation3D interpolation(*this);
std::vector<std::map<int,double> > matrixTmp;
nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
- reverseMatrix(matrixTmp,nbCols,_matrix);
+ ReverseMatrix(matrixTmp,nbCols,_matrix);
nbCols=matrixTmp.size();
}
else if(srcMeshDim==2 && trgMeshDim==1 && srcSpaceDim==2)
INTERP_KERNEL::Interpolation2D1D interpolation(*this);
std::vector<std::map<int,double> > matrixTmp;
nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
- reverseMatrix(matrixTmp,nbCols,_matrix);
+ ReverseMatrix(matrixTmp,nbCols,_matrix);
nbCols=matrixTmp.size();
INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
if(!duplicateFaces.empty())
INTERP_KERNEL::Interpolation2D interpolation(*this);
std::vector<std::map<int,double> > matrixTmp;
nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
- reverseMatrix(matrixTmp,nbCols,_matrix);
+ ReverseMatrix(matrixTmp,nbCols,_matrix);
nbCols=matrixTmp.size();
}
else
INTERP_KERNEL::Interpolation3D2D interpolation(*this);
std::vector<std::map<int,double> > matrixTmp;
nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
- reverseMatrix(matrixTmp,nbCols,_matrix);
+ ReverseMatrix(matrixTmp,nbCols,_matrix);
nbCols=matrixTmp.size();
INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
if(!duplicateFaces.empty())
return 1;
}
+int MEDCouplingRemapper::prepareUC(const char *method) throw(INTERP_KERNEL::Exception)
+{
+ std::string methodCpp(method);
+ if(methodCpp!="P0P0")
+ throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareUC : only P0P0 interpolation supported for the moment !");
+ INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::checkAndSplitInterpolationMethod(method,_src_method,_target_method);
+ MEDCouplingUMesh *src_mesh=static_cast<MEDCouplingUMesh *>(_src_mesh);
+ MEDCouplingCMesh *target_mesh=static_cast<MEDCouplingCMesh *>(_target_mesh);
+ const int srcMeshDim=src_mesh->getMeshDimension();
+ const int srcSpceDim=src_mesh->getSpaceDimension();
+ const int trgMeshDim=target_mesh->getMeshDimension();
+ if(srcMeshDim!=srcSpceDim || srcMeshDim!=trgMeshDim)
+ throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareUC : 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);
+ MEDCouplingNormalizedUnstructuredMesh<1,1> sourceWrapper(src_mesh);
+ INTERP_KERNEL::InterpolationCU myInterpolator(*this);
+ myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
+ break;
+ }
+ case 2:
+ {
+ MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
+ MEDCouplingNormalizedUnstructuredMesh<2,2> sourceWrapper(src_mesh);
+ INTERP_KERNEL::InterpolationCU myInterpolator(*this);
+ myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
+ break;
+ }
+ case 3:
+ {
+ MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
+ MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(src_mesh);
+ INTERP_KERNEL::InterpolationCU myInterpolator(*this);
+ myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
+ break;
+ }
+ default:
+ throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareUC : only dimension 1 2 or 3 supported !");
+ }
+ ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix);
+ nullifiedTinyCoeffInCrudeMatrixAbs(0.);
+ //
+ _deno_multiply.clear();
+ _deno_multiply.resize(_matrix.size());
+ _deno_reverse_multiply.clear();
+ _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
+ declareAsNew();
+ return 1;
+}
+
+int MEDCouplingRemapper::prepareCU(const char *method) throw(INTERP_KERNEL::Exception)
+{
+ std::string methodCpp(method);
+ if(methodCpp!="P0P0")
+ throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareCU : only P0P0 interpolation supported for the moment !");
+ INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::checkAndSplitInterpolationMethod(method,_src_method,_target_method);
+ MEDCouplingCMesh *src_mesh=static_cast<MEDCouplingCMesh *>(_src_mesh);
+ MEDCouplingUMesh *target_mesh=static_cast<MEDCouplingUMesh *>(_target_mesh);
+ const int srcMeshDim=src_mesh->getMeshDimension();
+ const int trgMeshDim=target_mesh->getMeshDimension();
+ const int trgSpceDim=target_mesh->getSpaceDimension();
+ if(trgMeshDim!=trgSpceDim || trgMeshDim!=srcMeshDim)
+ throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareCU : 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);
+ MEDCouplingNormalizedUnstructuredMesh<1,1> targetWrapper(target_mesh);
+ INTERP_KERNEL::InterpolationCU myInterpolator(*this);
+ myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
+ break;
+ }
+ case 2:
+ {
+ MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
+ MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(target_mesh);
+ INTERP_KERNEL::InterpolationCU myInterpolator(*this);
+ myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
+ break;
+ }
+ case 3:
+ {
+ MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
+ MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(target_mesh);
+ INTERP_KERNEL::InterpolationCU myInterpolator(*this);
+ myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
+ break;
+ }
+ default:
+ throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareCU : only dimension 1 2 or 3 supported !");
+ }
+ nullifiedTinyCoeffInCrudeMatrixAbs(0.);
+ //
+ _deno_multiply.clear();
+ _deno_multiply.resize(_matrix.size());
+ _deno_reverse_multiply.clear();
+ _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
+ declareAsNew();
+ return 1;
+}
+
+int MEDCouplingRemapper::prepareCC(const char *method) throw(INTERP_KERNEL::Exception)
+{
+ std::string methodCpp(method);
+ if(methodCpp!="P0P0")
+ throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareCC : only P0P0 interpolation supported for the moment !");
+ INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::checkAndSplitInterpolationMethod(method,_src_method,_target_method);
+ MEDCouplingCMesh *src_mesh=static_cast<MEDCouplingCMesh *>(_src_mesh);
+ MEDCouplingCMesh *target_mesh=static_cast<MEDCouplingCMesh *>(_target_mesh);
+ const int srcMeshDim=src_mesh->getMeshDimension();
+ const int trgMeshDim=target_mesh->getMeshDimension();
+ if(trgMeshDim!=srcMeshDim)
+ throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareCC : dim of target cartesian should be equal to dim of source cartesian dimension !");
+ switch(srcMeshDim)
+ {
+ case 1:
+ {
+ MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
+ MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
+ INTERP_KERNEL::InterpolationCC myInterpolator(*this);
+ myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
+ break;
+ }
+ case 2:
+ {
+ MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
+ MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
+ INTERP_KERNEL::InterpolationCC myInterpolator(*this);
+ myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
+ break;
+ }
+ case 3:
+ {
+ MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
+ MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
+ INTERP_KERNEL::InterpolationCC myInterpolator(*this);
+ myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
+ break;
+ }
+ default:
+ throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareCC : only dimension 1 2 or 3 supported !");
+ }
+ nullifiedTinyCoeffInCrudeMatrixAbs(0.);
+ //
+ _deno_multiply.clear();
+ _deno_multiply.resize(_matrix.size());
+ _deno_reverse_multiply.clear();
+ _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
+ declareAsNew();
+ return 1;
+}
+
void MEDCouplingRemapper::updateTime() const
{
}
{
case ConservativeVolumic:
{
- computeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
+ ComputeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
break;
}
case Integral:
{
- MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),true);
- MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),true);
+ MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
+ MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
const double *denoPtr=deno->getArray()->getConstPointer();
const double *denoRPtr=denoR->getArray()->getConstPointer();
if(trgField->getMesh()->getMeshDimension()==-1)
}
case IntegralGlobConstraint:
{
- computeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
+ ComputeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
break;
}
case RevIntegral:
{
- MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),true);
- MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),true);
+ MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
+ MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
const double *denoPtr=deno->getArray()->getConstPointer();
const double *denoRPtr=denoR->getArray()->getConstPointer();
if(trgField->getMesh()->getMeshDimension()==-1)
std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
}
-void MEDCouplingRemapper::reverseMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
+void MEDCouplingRemapper::ReverseMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
{
matOut.resize(nbColsMatIn);
int id=0;
matOut[(*iter2).first][id]=(*iter2).second;
}
-void MEDCouplingRemapper::computeRowSumAndColSum(const std::vector<std::map<int,double> >& matrixDeno,
+void MEDCouplingRemapper::ComputeRowSumAndColSum(const std::vector<std::map<int,double> >& matrixDeno,
std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
{
std::map<int,double> values;
}
}
-void MEDCouplingRemapper::computeColSumAndRowSum(const std::vector<std::map<int,double> >& matrixDeno,
+void MEDCouplingRemapper::ComputeColSumAndRowSum(const std::vector<std::map<int,double> >& matrixDeno,
std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
{
std::map<int,double> values;
{
return _matrix;
}
+
+/*!
+ * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
+ * If not the behaviour is unpredictable.
+ * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than \a maxValAbs this coefficient is
+ * set to 0. That is to say that its entry disappear from the map storing the corresponding row in the data storage of sparse crude matrix.
+ * This method is useful to correct at a high level some problems linked to precision. Indeed, with some \ref NatureOfField "natures of field" some threshold effect
+ * can occur.
+ *
+ * \param [in] maxValAbs is a limit behind which a coefficient is set to 0. \a maxValAbs is expected to be positive, if not this method do nothing.
+ * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged.
+ * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix
+ */
+int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs(double maxValAbs) throw(INTERP_KERNEL::Exception)
+{
+ int ret=0;
+ std::vector<std::map<int,double> > matrixNew(_matrix.size());
+ int i=0;
+ for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++,i++)
+ {
+ std::map<int,double>& rowNew=matrixNew[i];
+ for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
+ {
+ if(fabs((*it2).second)>maxValAbs)
+ rowNew[(*it2).first]=(*it2).second;
+ else
+ ret++;
+ }
+ }
+ if(ret>0)
+ _matrix=matrixNew;
+ return ret;
+}
+
+/*!
+ * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
+ * If not the behaviour is unpredictable.
+ * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than delta multiplied by \a scaleFactor this coefficient is
+ * set to 0. That is to say that its entry disappear from the map storing the corresponding row in the data storage of sparse crude matrix.
+ * delta is the value returned by MEDCouplingRemapper::getMaxValueInCrudeMatrix method.
+ * This method is useful to correct at a high level some problems linked to precision. Indeed, with some \ref NatureOfField "natures of field" some threshold effect
+ * can occur.
+ *
+ * \param [in] scaleFactor is the scale factor from which coefficients lower than \a scaleFactor times range width of coefficients are set to zero.
+ * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged. If -1 is returned it means
+ * that all coefficients are null.
+ * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs
+ */
+int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix(double scaleFactor) throw(INTERP_KERNEL::Exception)
+{
+ double maxVal=getMaxValueInCrudeMatrix();
+ if(maxVal==0.)
+ return -1;
+ return nullifiedTinyCoeffInCrudeMatrixAbs(scaleFactor*maxVal);
+}
+
+/*!
+ * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
+ * If not the behaviour is unpredictable.
+ * This method returns the maximum of the absolute values of coefficients into the sparse crude matrix.
+ * The returned value is positive.
+ */
+double MEDCouplingRemapper::getMaxValueInCrudeMatrix() const throw(INTERP_KERNEL::Exception)
+{
+ double ret=0.;
+ for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++)
+ for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
+ if(fabs((*it2).second)>ret)
+ ret=fabs((*it2).second);
+ return ret;
+}
