Salome HOME
Spread and condense on cell fields on images meshes with ghost management in 1D and 2D.
[tools/medcoupling.git] / src / MEDCoupling / MEDCouplingRemapper.cxx
index 1c5e17a917958fe933db1eea4040d2759d160b58..33a0debffedb8116c402d9dd858d2f0fc0ae9518 100644 (file)
@@ -1,9 +1,9 @@
-// Copyright (C) 2007-2013  CEA/DEN, EDF R&D
+// Copyright (C) 2007-2014  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
@@ -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());
       }
-    }
+  }
 }
 
 /*!
@@ -281,7 +281,7 @@ int MEDCouplingRemapper::getInterpolationMatrixPolicy() const
 void MEDCouplingRemapper::setInterpolationMatrixPolicy(int newInterpMatPol)
 {
   switch(newInterpMatPol)
-    {
+  {
     case 0:
       _interp_matrix_pol=IK_ONLY_PREFERED;
       break;
@@ -296,7 +296,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 +304,7 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUU()
   const MEDCouplingPointSet *src_mesh=static_cast<const MEDCouplingPointSet *>(_src_ft->getMesh());
   const MEDCouplingPointSet *target_mesh=static_cast<const MEDCouplingPointSet *>(_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 +369,8 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUU()
       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
       INTERP_KERNEL::Interpolation3D interpolation(*this);
       std::vector<std::map<int,double> > 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 +389,8 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUU()
           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
           INTERP_KERNEL::Interpolation2D1D interpolation(*this);
           std::vector<std::map<int,double> > 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 +414,8 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUU()
           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
           INTERP_KERNEL::Interpolation2D interpolation(*this);
           std::vector<std::map<int,double> > 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 +462,8 @@ int MEDCouplingRemapper::prepareInterpKernelOnlyUU()
       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
       INTERP_KERNEL::Interpolation3D2D interpolation(*this);
       std::vector<std::map<int,double> > 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 +583,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<std::map<int,double> > res;
   switch(srcMeshDim)
-    {
+  {
     case 1:
       {
         MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
@@ -606,7 +610,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 +636,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 +663,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 +687,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 +714,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 +728,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<DataArrayDouble> trgLoc=_target_ft->getLocalizationOfDiscr();
   const double *trgLocPtr=trgLoc->begin();
   int trgSpaceDim=trgLoc->getNumberOfComponents();
@@ -804,31 +808,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<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;
@@ -836,7 +840,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 +873,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;
 }
 
@@ -921,7 +930,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 +940,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 +1007,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)