src/MEDCoupling/MEDCouplingRemapper.cxx

   1 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (CEA/DEN)
  20
  21 #include "MEDCouplingRemapper.hxx"
  22 #include "MEDCouplingMemArray.hxx"
  23 #include "MEDCouplingFieldDouble.hxx"
  24 #include "MEDCouplingFieldTemplate.hxx"
  25 #include "MEDCouplingFieldDiscretization.hxx"
  26 #include "MEDCouplingExtrudedMesh.hxx"
  27 #include "MEDCouplingCMesh.hxx"
  28 #include "MEDCouplingNormalizedUnstructuredMesh.txx"
  29 #include "MEDCouplingNormalizedCartesianMesh.txx"
  30
  31 #include "Interpolation1D.txx"
  32 #include "Interpolation2DCurve.hxx"
  33 #include "Interpolation2D.txx"
  34 #include "Interpolation3D.txx"
  35 #include "Interpolation3DSurf.hxx"
  36 #include "Interpolation2D1D.txx"
  37 #include "Interpolation3D2D.txx"
  38 #include "InterpolationCU.txx"
  39 #include "InterpolationCC.txx"
  40
  41 using namespace ParaMEDMEM;
  42
  43 MEDCouplingRemapper::MEDCouplingRemapper():_src_mesh(0),_target_mesh(0),_nature_of_deno(NoNature),_time_deno_update(0)
  44 {
  45 }
  46
  47 MEDCouplingRemapper::~MEDCouplingRemapper()
  48 {
  49   releaseData(false);
  50 }
  51
  52 int MEDCouplingRemapper::prepare(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const char *method) throw(INTERP_KERNEL::Exception)
  53 {
  54   releaseData(true);
  55   _src_mesh=const_cast<MEDCouplingMesh *>(srcMesh); _target_mesh=const_cast<MEDCouplingMesh *>(targetMesh);
  56   _src_mesh->incrRef(); _target_mesh->incrRef();
  57   int meshInterpType=((int)_src_mesh->getType()*16)+(int)_target_mesh->getType();
  58   switch(meshInterpType)
  59     {
  60     case 85://Unstructured-Unstructured
  61       return prepareUU(method);
  62     case 87://Unstructured-Cartesian
  63       return prepareUC(method);
  64     case 117://Cartesian-Unstructured
  65       return prepareCU(method);
  66     case 119://Cartesian-Cartesian
  67       return prepareCC(method);
  68     case 136://Extruded-Extruded
  69       return prepareEE(method);
  70     default:
  71       throw INTERP_KERNEL::Exception("Not managed type of meshes !");
  72     }
  73 }
  74
  75 int MEDCouplingRemapper::prepareEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target) throw(INTERP_KERNEL::Exception)
  76 {
  77   std::string meth(src->getDiscretization()->getStringRepr());
  78   meth+=target->getDiscretization()->getStringRepr();
  79   return prepare(src->getMesh(),target->getMesh(),meth.c_str());
  80 }
  81
  82 /*!
  83  * This method performs the operation source to target using matrix computed in ParaMEDMEM::MEDCouplingRemapper::prepare method.
  84  * 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.
  85  *
  86  * \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.
  87  * \param [out] targetField the destination field with the allocated array in which all tuples will be overwritten.
  88  */
  89 void MEDCouplingRemapper::transfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, double dftValue) throw(INTERP_KERNEL::Exception)
  90 {
  91   transferUnderground(srcField,targetField,true,dftValue);
  92 }
  93
  94 /*!
  95  * This method is equivalent to ParaMEDMEM::MEDCouplingRemapper::transfer except that here \b targetField is a in/out parameter.
  96  * If an entity (cell for example) in targetField is not fetched by any entity (cell for example) of \b srcField, the value in targetField is
  97  * let unchanged.
  98  * This method requires that \b targetField was fully defined and allocated. If the array is not allocated an exception will be thrown.
  99  *
 100  * \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.
 101  * \param [in,out] targetField the destination field with the allocated array in which only tuples whose entities are fetched by interpolation will be overwritten only.
 102  */
 103 void MEDCouplingRemapper::partialTransfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField) throw(INTERP_KERNEL::Exception)
 104 {
 105   transferUnderground(srcField,targetField,false,std::numeric_limits<double>::max());
 106 }
 107
 108 void MEDCouplingRemapper::reverseTransfer(MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *targetField, double dftValue) throw(INTERP_KERNEL::Exception)
 109 {
 110   if(_src_method!=srcField->getDiscretization()->getStringRepr())
 111     throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
 112   if(_target_method!=targetField->getDiscretization()->getStringRepr())
 113     throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
 114   if(srcField->getNature()!=targetField->getNature())
 115     throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
 116   DataArrayDouble *array=srcField->getArray();
 117   int trgNbOfCompo=targetField->getNumberOfComponents();
 118   if(array)
 119     {
 120       if(trgNbOfCompo!=srcField->getNumberOfComponents())
 121         throw INTERP_KERNEL::Exception("Number of components mismatch !");
 122     }
 123   else
 124     {
 125       array=DataArrayDouble::New();
 126       array->alloc(srcField->getNumberOfTuples(),trgNbOfCompo);
 127       srcField->setArray(array);
 128       array->decrRef();
 129     }
 130   computeDeno(srcField->getNature(),srcField,targetField);
 131   double *resPointer=array->getPointer();
 132   const double *inputPointer=targetField->getArray()->getConstPointer();
 133   computeReverseProduct(inputPointer,trgNbOfCompo,dftValue,resPointer);
 134 }
 135
 136 MEDCouplingFieldDouble *MEDCouplingRemapper::transferField(const MEDCouplingFieldDouble *srcField, double dftValue) throw(INTERP_KERNEL::Exception)
 137 {
 138   if(_src_method!=srcField->getDiscretization()->getStringRepr())
 139     throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
 140   MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(MEDCouplingFieldDiscretization::getTypeOfFieldFromStringRepr(_target_method.c_str()),srcField->getTimeDiscretization());
 141   ret->copyAllTinyAttrFrom(srcField);
 142   ret->setNature(srcField->getNature());
 143   ret->setMesh(_target_mesh);
 144   transfer(srcField,ret,dftValue);
 145   return ret;
 146 }
 147
 148 MEDCouplingFieldDouble *MEDCouplingRemapper::reverseTransferField(const MEDCouplingFieldDouble *targetField, double dftValue) throw(INTERP_KERNEL::Exception)
 149 {
 150   if(_target_method!=targetField->getDiscretization()->getStringRepr())
 151     throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
 152   MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(MEDCouplingFieldDiscretization::getTypeOfFieldFromStringRepr(_src_method.c_str()),targetField->getTimeDiscretization());
 153   ret->copyAllTinyAttrFrom(targetField);
 154   ret->setNature(targetField->getNature());
 155   ret->setMesh(_src_mesh);
 156   reverseTransfer(ret,targetField,dftValue);
 157   return ret;
 158 }
 159
 160 bool MEDCouplingRemapper::setOptionInt(const std::string& key, int value)
 161 {
 162   return INTERP_KERNEL::InterpolationOptions::setOptionInt(key,value);
 163 }
 164
 165 bool MEDCouplingRemapper::setOptionDouble(const std::string& key, double value)
 166 {
 167   return INTERP_KERNEL::InterpolationOptions::setOptionDouble(key,value);
 168 }
 169
 170 bool MEDCouplingRemapper::setOptionString(const std::string& key, const std::string& value)
 171 {
 172   return INTERP_KERNEL::InterpolationOptions::setOptionString(key,value);
 173 }
 174
 175 int MEDCouplingRemapper::prepareUU(const char *method) throw(INTERP_KERNEL::Exception)
 176 {
 177   MEDCouplingUMesh *src_mesh=(MEDCouplingUMesh *)_src_mesh;
 178   MEDCouplingUMesh *target_mesh=(MEDCouplingUMesh *)_target_mesh;
 179   INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::checkAndSplitInterpolationMethod(method,_src_method,_target_method);
 180   const int srcMeshDim=src_mesh->getMeshDimension();
 181   int srcSpaceDim=-1;
 182   if(srcMeshDim!=-1)
 183     srcSpaceDim=src_mesh->getSpaceDimension();
 184   const int trgMeshDim=target_mesh->getMeshDimension();
 185   int trgSpaceDim=-1;
 186   if(trgMeshDim!=-1)
 187     trgSpaceDim=target_mesh->getSpaceDimension();
 188   if(trgSpaceDim!=srcSpaceDim)
 189     if(trgSpaceDim!=-1 && srcSpaceDim!=-1)
 190       throw INTERP_KERNEL::Exception("Incoherent space dimension detected between target and source.");
 191   int nbCols;
 192   if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==1)
 193     {
 194       MEDCouplingNormalizedUnstructuredMesh<1,1> source_mesh_wrapper(src_mesh);
 195       MEDCouplingNormalizedUnstructuredMesh<1,1> target_mesh_wrapper(target_mesh);
 196       INTERP_KERNEL::Interpolation1D interpolation(*this);
 197       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 198     }
 199   else if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==2)
 200     {
 201       MEDCouplingNormalizedUnstructuredMesh<2,1> source_mesh_wrapper(src_mesh);
 202       MEDCouplingNormalizedUnstructuredMesh<2,1> target_mesh_wrapper(target_mesh);
 203       INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
 204       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 205     }
 206   else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==2)
 207     {
 208       MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 209       MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
 210       INTERP_KERNEL::Interpolation2D interpolation(*this);
 211       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 212     }
 213   else if(srcMeshDim==3 && trgMeshDim==3 && srcSpaceDim==3)
 214     {
 215       MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 216       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
 217       INTERP_KERNEL::Interpolation3D interpolation(*this);
 218       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 219     }
 220   else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==3)
 221     {
 222       MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
 223       MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh);
 224       INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
 225       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 226     }
 227   else if(srcMeshDim==3 && trgMeshDim==1 && srcSpaceDim==3)
 228     {
 229       if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
 230         throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
 231       MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 232       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
 233       INTERP_KERNEL::Interpolation3D interpolation(*this);
 234       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 235     }
 236   else if(srcMeshDim==1 && trgMeshDim==3 && srcSpaceDim==3)
 237     {
 238       if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
 239         throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
 240       MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 241       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
 242       INTERP_KERNEL::Interpolation3D interpolation(*this);
 243       std::vector<std::map<int,double> > matrixTmp;
 244       nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
 245       ReverseMatrix(matrixTmp,nbCols,_matrix);
 246       nbCols=matrixTmp.size();
 247     }
 248   else if(srcMeshDim==2 && trgMeshDim==1 && srcSpaceDim==2)
 249     {
 250       if(getIntersectionType()==INTERP_KERNEL::PointLocator)
 251         {
 252           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 253           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
 254           INTERP_KERNEL::Interpolation2D interpolation(*this);
 255           nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 256         }
 257       else
 258         {
 259           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 260           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
 261           INTERP_KERNEL::Interpolation2D1D interpolation(*this);
 262           std::vector<std::map<int,double> > matrixTmp;
 263           nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
 264           ReverseMatrix(matrixTmp,nbCols,_matrix);
 265           nbCols=matrixTmp.size();
 266           INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
 267           if(!duplicateFaces.empty())
 268             {
 269               std::ostringstream oss; oss << "An unexpected situation happend ! For the following 1D Cells are part of edges shared by 2D cells :\n";
 270               for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
 271                 {
 272                   oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
 273                   std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
 274                   oss << std::endl;
 275                 }
 276             }
 277         }
 278     }
 279   else if(srcMeshDim==1 && trgMeshDim==2 && srcSpaceDim==2)
 280     {
 281       if(getIntersectionType()==INTERP_KERNEL::PointLocator)
 282         {
 283           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 284           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
 285           INTERP_KERNEL::Interpolation2D interpolation(*this);
 286           std::vector<std::map<int,double> > matrixTmp;
 287           nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
 288           ReverseMatrix(matrixTmp,nbCols,_matrix);
 289           nbCols=matrixTmp.size();
 290         }
 291       else
 292         {
 293           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 294           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
 295           INTERP_KERNEL::Interpolation2D1D interpolation(*this);
 296           nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 297           INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
 298           if(!duplicateFaces.empty())
 299             {
 300               std::ostringstream oss; oss << "An unexpected situation happend ! For the following 1D Cells are part of edges shared by 2D cells :\n";
 301               for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
 302                 {
 303                   oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
 304                   std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
 305                   oss << std::endl;
 306                 }
 307             }
 308         }
 309     }
 310   else if(srcMeshDim==2 && trgMeshDim==3 && srcSpaceDim==3)
 311     {
 312       MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 313       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
 314       INTERP_KERNEL::Interpolation3D2D interpolation(*this);
 315       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 316       INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
 317       if(!duplicateFaces.empty())
 318         {
 319           std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
 320           for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
 321             {
 322               oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
 323               std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
 324               oss << std::endl;
 325             }
 326         }
 327     }
 328   else if(srcMeshDim==3 && trgMeshDim==2 && srcSpaceDim==3)
 329     {
 330       MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 331       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
 332       INTERP_KERNEL::Interpolation3D2D interpolation(*this);
 333       std::vector<std::map<int,double> > matrixTmp;
 334       nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
 335       ReverseMatrix(matrixTmp,nbCols,_matrix);
 336       nbCols=matrixTmp.size();
 337       INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
 338       if(!duplicateFaces.empty())
 339         {
 340           std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
 341           for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
 342             {
 343               oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
 344               std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
 345               oss << std::endl;
 346             }
 347         }
 348     }
 349   else if(trgMeshDim==-1)
 350     {
 351       if(srcMeshDim==2 && srcSpaceDim==2)
 352         {
 353           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 354           INTERP_KERNEL::Interpolation2D interpolation(*this);
 355           nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,_src_method.c_str());
 356         }
 357       else if(srcMeshDim==3 && srcSpaceDim==3)
 358         {
 359           MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 360           INTERP_KERNEL::Interpolation3D interpolation(*this);
 361           nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,_src_method.c_str());
 362         }
 363       else if(srcMeshDim==2 && srcSpaceDim==3)
 364         {
 365           MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
 366           INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
 367           nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,_src_method.c_str());
 368         }
 369       else
 370         throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
 371     }
 372   else if(srcMeshDim==-1)
 373     {
 374       if(trgMeshDim==2 && trgSpaceDim==2)
 375         {
 376           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(target_mesh);
 377           INTERP_KERNEL::Interpolation2D interpolation(*this);
 378           nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,_target_method.c_str());
 379         }
 380       else if(trgMeshDim==3 && trgSpaceDim==3)
 381         {
 382           MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(target_mesh);
 383           INTERP_KERNEL::Interpolation3D interpolation(*this);
 384           nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,_target_method.c_str());
 385         }
 386       else if(trgMeshDim==2 && trgSpaceDim==3)
 387         {
 388           MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(target_mesh);
 389           INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
 390           nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,_target_method.c_str());
 391         }
 392       else
 393         throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh from -1D sourceMesh");
 394     }
 395   else
 396     throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension");
 397   _deno_multiply.clear();
 398   _deno_multiply.resize(_matrix.size());
 399   _deno_reverse_multiply.clear();
 400   _deno_reverse_multiply.resize(nbCols);
 401   declareAsNew();
 402   return 1;
 403 }
 404
 405 int MEDCouplingRemapper::prepareEE(const char *method) throw(INTERP_KERNEL::Exception)
 406 {
 407   MEDCouplingExtrudedMesh *src_mesh=(MEDCouplingExtrudedMesh *)_src_mesh;
 408   MEDCouplingExtrudedMesh *target_mesh=(MEDCouplingExtrudedMesh *)_target_mesh;
 409   std::string methC(method);
 410   if(methC!="P0P0")
 411     throw INTERP_KERNEL::Exception("Only P0P0 method implemented for Extruded/Extruded meshes !");
 412   INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::checkAndSplitInterpolationMethod(method,_src_method,_target_method);
 413   MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh->getMesh2D());
 414   MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh->getMesh2D());
 415   INTERP_KERNEL::Interpolation3DSurf interpolation2D(*this);
 416   std::vector<std::map<int,double> > matrix2D;
 417   int nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,method);
 418   MEDCouplingUMesh *s1D,*t1D;
 419   double v[3];
 420   MEDCouplingExtrudedMesh::Project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v);
 421   MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D);
 422   MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
 423   std::vector<std::map<int,double> > matrix1D;
 424   INTERP_KERNEL::Interpolation1D interpolation1D(*this);
 425   int nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,method);
 426   s1D->decrRef();
 427   t1D->decrRef();
 428   buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
 429                                              target_mesh->getMesh3DIds()->getConstPointer());
 430   //
 431   _deno_multiply.clear();
 432   _deno_multiply.resize(_matrix.size());
 433   _deno_reverse_multiply.clear();
 434   _deno_reverse_multiply.resize(nbCols2D*nbCols1D);
 435   declareAsNew();
 436   return 1;
 437 }
 438
 439 int MEDCouplingRemapper::prepareUC(const char *method) throw(INTERP_KERNEL::Exception)
 440 {
 441   std::string methodCpp(method);
 442   if(methodCpp!="P0P0")
 443     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareUC : only P0P0 interpolation supported for the moment !");
 444   INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::checkAndSplitInterpolationMethod(method,_src_method,_target_method);
 445   MEDCouplingUMesh *src_mesh=static_cast<MEDCouplingUMesh *>(_src_mesh);
 446   MEDCouplingCMesh *target_mesh=static_cast<MEDCouplingCMesh *>(_target_mesh);
 447   const int srcMeshDim=src_mesh->getMeshDimension();
 448   const int srcSpceDim=src_mesh->getSpaceDimension();
 449   const int trgMeshDim=target_mesh->getMeshDimension();
 450   if(srcMeshDim!=srcSpceDim || srcMeshDim!=trgMeshDim)
 451     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 !");
 452   std::vector<std::map<int,double> > res;
 453   switch(srcMeshDim)
 454     {
 455     case 1:
 456       {
 457         MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
 458         MEDCouplingNormalizedUnstructuredMesh<1,1> sourceWrapper(src_mesh);
 459         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 460         myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
 461         break;
 462       }
 463     case 2:
 464       {
 465         MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
 466         MEDCouplingNormalizedUnstructuredMesh<2,2> sourceWrapper(src_mesh);
 467         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 468         myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
 469         break;
 470       }
 471     case 3:
 472       {
 473         MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
 474         MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(src_mesh);
 475         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 476         myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
 477         break;
 478       }
 479     default:
 480       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareUC : only dimension 1 2 or 3 supported !");
 481     }
 482   ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix);
 483   nullifiedTinyCoeffInCrudeMatrixAbs(0.);
 484   //
 485   _deno_multiply.clear();
 486   _deno_multiply.resize(_matrix.size());
 487   _deno_reverse_multiply.clear();
 488   _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
 489   declareAsNew();
 490   return 1;
 491 }
 492
 493 int MEDCouplingRemapper::prepareCU(const char *method) throw(INTERP_KERNEL::Exception)
 494 {
 495   std::string methodCpp(method);
 496   if(methodCpp!="P0P0")
 497     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareCU : only P0P0 interpolation supported for the moment !");
 498   INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::checkAndSplitInterpolationMethod(method,_src_method,_target_method);
 499   MEDCouplingCMesh *src_mesh=static_cast<MEDCouplingCMesh *>(_src_mesh);
 500   MEDCouplingUMesh *target_mesh=static_cast<MEDCouplingUMesh *>(_target_mesh);
 501   const int srcMeshDim=src_mesh->getMeshDimension();
 502   const int trgMeshDim=target_mesh->getMeshDimension();
 503   const int trgSpceDim=target_mesh->getSpaceDimension();
 504   if(trgMeshDim!=trgSpceDim || trgMeshDim!=srcMeshDim)
 505     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 !");
 506   switch(srcMeshDim)
 507     {
 508     case 1:
 509       {
 510         MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
 511         MEDCouplingNormalizedUnstructuredMesh<1,1> targetWrapper(target_mesh);
 512         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 513         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 514         break;
 515       }
 516     case 2:
 517       {
 518         MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
 519         MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(target_mesh);
 520         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 521         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 522         break;
 523       }
 524     case 3:
 525       {
 526         MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
 527         MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(target_mesh);
 528         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 529         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 530         break;
 531       }
 532     default:
 533       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareCU : only dimension 1 2 or 3 supported !");
 534     }
 535   nullifiedTinyCoeffInCrudeMatrixAbs(0.);
 536   //
 537   _deno_multiply.clear();
 538   _deno_multiply.resize(_matrix.size());
 539   _deno_reverse_multiply.clear();
 540   _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
 541   declareAsNew();
 542   return 1;
 543 }
 544
 545 int MEDCouplingRemapper::prepareCC(const char *method) throw(INTERP_KERNEL::Exception)
 546 {
 547   std::string methodCpp(method);
 548   if(methodCpp!="P0P0")
 549     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareCC : only P0P0 interpolation supported for the moment !");
 550   INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::checkAndSplitInterpolationMethod(method,_src_method,_target_method);
 551   MEDCouplingCMesh *src_mesh=static_cast<MEDCouplingCMesh *>(_src_mesh);
 552   MEDCouplingCMesh *target_mesh=static_cast<MEDCouplingCMesh *>(_target_mesh);
 553   const int srcMeshDim=src_mesh->getMeshDimension();
 554   const int trgMeshDim=target_mesh->getMeshDimension();
 555   if(trgMeshDim!=srcMeshDim)
 556     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareCC : dim of target cartesian should be equal to dim of source cartesian dimension !");
 557   switch(srcMeshDim)
 558     {
 559     case 1:
 560       {
 561         MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
 562         MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
 563         INTERP_KERNEL::InterpolationCC myInterpolator(*this);
 564         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 565         break;
 566       }
 567     case 2:
 568       {
 569         MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
 570         MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
 571         INTERP_KERNEL::InterpolationCC myInterpolator(*this);
 572         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 573         break;
 574       }
 575     case 3:
 576       {
 577         MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
 578         MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
 579         INTERP_KERNEL::InterpolationCC myInterpolator(*this);
 580         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 581         break;
 582       }
 583     default:
 584       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareCC : only dimension 1 2 or 3 supported !");
 585     }
 586   nullifiedTinyCoeffInCrudeMatrixAbs(0.);
 587   //
 588   _deno_multiply.clear();
 589   _deno_multiply.resize(_matrix.size());
 590   _deno_reverse_multiply.clear();
 591   _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
 592   declareAsNew();
 593   return 1;
 594 }
 595
 596 void MEDCouplingRemapper::updateTime() const
 597 {
 598 }
 599
 600 void MEDCouplingRemapper::releaseData(bool matrixSuppression)
 601 {
 602   if(_src_mesh)
 603     _src_mesh->decrRef();
 604   if(_target_mesh)
 605     _target_mesh->decrRef();
 606   _src_mesh=0;
 607   _target_mesh=0;
 608   if(matrixSuppression)
 609     {
 610       _matrix.clear();
 611       _deno_multiply.clear();
 612       _deno_reverse_multiply.clear();
 613     }
 614 }
 615
 616 void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue) throw(INTERP_KERNEL::Exception)
 617 {
 618   if(_src_method!=srcField->getDiscretization()->getStringRepr())
 619     throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
 620   if(_target_method!=targetField->getDiscretization()->getStringRepr())
 621     throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
 622   if(srcField->getNature()!=targetField->getNature())
 623     throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
 624   DataArrayDouble *array=targetField->getArray();
 625   int srcNbOfCompo=srcField->getNumberOfComponents();
 626   if(array)
 627     {
 628       if(srcNbOfCompo!=targetField->getNumberOfComponents())
 629         throw INTERP_KERNEL::Exception("Number of components mismatch !");
 630     }
 631   else
 632     {
 633       if(!isDftVal)
 634         throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !");
 635       array=DataArrayDouble::New();
 636       array->alloc(targetField->getNumberOfTuples(),srcNbOfCompo);
 637       targetField->setArray(array);
 638       array->decrRef();
 639     }
 640   computeDeno(srcField->getNature(),srcField,targetField);
 641   double *resPointer=array->getPointer();
 642   const double *inputPointer=srcField->getArray()->getConstPointer();
 643   computeProduct(inputPointer,srcNbOfCompo,isDftVal,dftValue,resPointer);
 644 }
 645
 646 void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
 647 {
 648   if(nat==NoNature)
 649     return computeDenoFromScratch(nat,srcField,trgField);
 650   else if(nat!=_nature_of_deno)
 651    return computeDenoFromScratch(nat,srcField,trgField);
 652   else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis())
 653     return computeDenoFromScratch(nat,srcField,trgField);
 654 }
 655
 656 void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField) throw(INTERP_KERNEL::Exception)
 657 {
 658   _nature_of_deno=nat;
 659   _time_deno_update=getTimeOfThis();
 660   switch(_nature_of_deno)
 661     {
 662     case ConservativeVolumic:
 663       {
 664         ComputeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
 665         break;
 666       }
 667     case Integral:
 668       {
 669         MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
 670         MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
 671         const double *denoPtr=deno->getArray()->getConstPointer();
 672         const double *denoRPtr=denoR->getArray()->getConstPointer();
 673         if(trgField->getMesh()->getMeshDimension()==-1)
 674           {
 675             double *denoRPtr2=denoR->getArray()->getPointer();
 676             denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
 677           }
 678         if(srcField->getMesh()->getMeshDimension()==-1)
 679           {
 680             double *denoPtr2=deno->getArray()->getPointer();
 681             denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
 682           }
 683         int idx=0;
 684         for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
 685           for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 686             {
 687               _deno_multiply[idx][(*iter2).first]=denoPtr[(*iter2).first];
 688               _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[idx];
 689             }
 690         deno->decrRef();
 691         denoR->decrRef();
 692         break;
 693       }
 694     case IntegralGlobConstraint:
 695       {
 696         ComputeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
 697         break;
 698       }
 699     case RevIntegral:
 700       {
 701         MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
 702         MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
 703         const double *denoPtr=deno->getArray()->getConstPointer();
 704         const double *denoRPtr=denoR->getArray()->getConstPointer();
 705         if(trgField->getMesh()->getMeshDimension()==-1)
 706           {
 707             double *denoRPtr2=denoR->getArray()->getPointer();
 708             denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
 709           }
 710         if(srcField->getMesh()->getMeshDimension()==-1)
 711           {
 712             double *denoPtr2=deno->getArray()->getPointer();
 713             denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
 714           }
 715         int idx=0;
 716         for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
 717           for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 718             {
 719               _deno_multiply[idx][(*iter2).first]=denoPtr[idx];
 720               _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[(*iter2).first];
 721             }
 722         deno->decrRef();
 723         denoR->decrRef();
 724         break;
 725       }
 726     case NoNature:
 727       throw INTERP_KERNEL::Exception("No nature specified ! Select one !");
 728     }
 729 }
 730
 731 void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)
 732 {
 733   int idx=0;
 734   double *tmp=new double[inputNbOfCompo];
 735   for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
 736     {
 737       if((*iter1).empty())
 738         {
 739           if(isDftVal)
 740             std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
 741           continue;
 742         }
 743       else
 744         std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,0.);
 745       std::map<int,double>::const_iterator iter3=_deno_multiply[idx].begin();
 746       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter3++)
 747         {
 748           std::transform(inputPointer+(*iter2).first*inputNbOfCompo,inputPointer+((*iter2).first+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/(*iter3).second));
 749           std::transform(tmp,tmp+inputNbOfCompo,resPointer+idx*inputNbOfCompo,resPointer+idx*inputNbOfCompo,std::plus<double>());
 750         }
 751     }
 752   delete [] tmp;
 753 }
 754
 755 void MEDCouplingRemapper::computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer)
 756 {
 757   std::vector<bool> isReached(_deno_reverse_multiply.size(),false);
 758   int idx=0;
 759   double *tmp=new double[inputNbOfCompo];
 760   std::fill(resPointer,resPointer+inputNbOfCompo*_deno_reverse_multiply.size(),0.);
 761   for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
 762     {
 763       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 764         {
 765           isReached[(*iter2).first]=true;
 766           std::transform(inputPointer+idx*inputNbOfCompo,inputPointer+(idx+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/_deno_reverse_multiply[(*iter2).first][idx]));
 767           std::transform(tmp,tmp+inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,std::plus<double>());
 768         }
 769     }
 770   delete [] tmp;
 771   idx=0;
 772   for(std::vector<bool>::const_iterator iter3=isReached.begin();iter3!=isReached.end();iter3++,idx++)
 773     if(!*iter3)
 774       std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
 775 }
 776
 777 void MEDCouplingRemapper::ReverseMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
 778 {
 779   matOut.resize(nbColsMatIn);
 780   int id=0;
 781   for(std::vector<std::map<int,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
 782     for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 783       matOut[(*iter2).first][id]=(*iter2).second;
 784 }
 785
 786 void MEDCouplingRemapper::ComputeRowSumAndColSum(const std::vector<std::map<int,double> >& matrixDeno,
 787                                                  std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
 788 {
 789   std::map<int,double> values;
 790   int idx=0;
 791   for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
 792     {
 793       double sum=0.;
 794       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 795         {
 796           sum+=(*iter2).second;
 797           values[(*iter2).first]+=(*iter2).second;
 798         }
 799       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 800         deno[idx][(*iter2).first]=sum;
 801     }
 802   idx=0;
 803   for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
 804     {
 805       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 806         denoReverse[(*iter2).first][idx]=values[(*iter2).first];
 807     }
 808 }
 809
 810 void MEDCouplingRemapper::ComputeColSumAndRowSum(const std::vector<std::map<int,double> >& matrixDeno,
 811                                                  std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
 812 {
 813   std::map<int,double> values;
 814   int idx=0;
 815   for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
 816     {
 817       double sum=0.;
 818       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 819         {
 820           sum+=(*iter2).second;
 821           values[(*iter2).first]+=(*iter2).second;
 822         }
 823       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 824         denoReverse[(*iter2).first][idx]=sum;
 825     }
 826   idx=0;
 827   for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
 828     {
 829       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 830         deno[idx][(*iter2).first]=values[(*iter2).first];
 831     }
 832 }
 833
 834 void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<int,double> >& m1D,
 835                                                                      const std::vector< std::map<int,double> >& m2D,
 836                                                                      const int *corrCellIdSrc, int nbOf2DCellsSrc, int nbOf1DCellsSrc,
 837                                                                      const int *corrCellIdTrg)
 838 {
 839   int nbOf2DCellsTrg=m2D.size();
 840   int nbOf1DCellsTrg=m1D.size();
 841   int nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
 842   _matrix.resize(nbOf3DCellsTrg);
 843   int id2R=0;
 844   for(std::vector< std::map<int,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
 845     {
 846       for(std::map<int,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
 847         {
 848           int id1R=0;
 849           for(std::vector< std::map<int,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
 850             {
 851               for(std::map<int,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
 852                 {
 853                   _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
 854                 }
 855             }
 856         }
 857     }
 858 }
 859
 860 void MEDCouplingRemapper::PrintMatrix(const std::vector<std::map<int,double> >& m)
 861 {
 862   int id=0;
 863   for(std::vector<std::map<int,double> >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++)
 864     {
 865       std::cout << "Target Cell # " << id << " : ";
 866       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 867         std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), ";
 868       std::cout << std::endl;
 869     }
 870 }
 871
 872 const std::vector<std::map<int,double> >& MEDCouplingRemapper::getCrudeMatrix() const
 873 {
 874   return _matrix;
 875 }
 876
 877 /*!
 878  * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
 879  * If not the behaviour is unpredictable.
 880  * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than \a maxValAbs this coefficient is
 881  * 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.
 882  * 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
 883  * can occur.
 884  *
 885  * \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.
 886  * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged.
 887  * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix
 888  */
 889 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs(double maxValAbs) throw(INTERP_KERNEL::Exception)
 890 {
 891   int ret=0;
 892   std::vector<std::map<int,double> > matrixNew(_matrix.size());
 893   int i=0;
 894   for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++,i++)
 895     {
 896       std::map<int,double>& rowNew=matrixNew[i];
 897       for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 898         {
 899           if(fabs((*it2).second)>maxValAbs)
 900             rowNew[(*it2).first]=(*it2).second;
 901           else
 902             ret++;
 903         }
 904     }
 905   if(ret>0)
 906     _matrix=matrixNew;
 907   return ret;
 908 }
 909
 910 /*!
 911  * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
 912  * If not the behaviour is unpredictable.
 913  * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than delta multiplied by \a scaleFactor this coefficient is
 914  * 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.
 915  * delta is the value returned by MEDCouplingRemapper::getMaxValueInCrudeMatrix method.
 916  * 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
 917  * can occur.
 918  *
 919  * \param [in] scaleFactor is the scale factor from which coefficients lower than \a scaleFactor times range width of coefficients are set to zero.
 920  * \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
 921  *         that all coefficients are null.
 922  * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs
 923  */
 924 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix(double scaleFactor) throw(INTERP_KERNEL::Exception)
 925 {
 926   double maxVal=getMaxValueInCrudeMatrix();
 927   if(maxVal==0.)
 928     return -1;
 929   return nullifiedTinyCoeffInCrudeMatrixAbs(scaleFactor*maxVal);
 930 }
 931
 932 /*!
 933  * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
 934  * If not the behaviour is unpredictable.
 935  * This method returns the maximum of the absolute values of coefficients into the sparse crude matrix.
 936  * The returned value is positive.
 937  */
 938 double MEDCouplingRemapper::getMaxValueInCrudeMatrix() const throw(INTERP_KERNEL::Exception)
 939 {
 940   double ret=0.;
 941   for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++)
 942     for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 943       if(fabs((*it2).second)>ret)
 944         ret=fabs((*it2).second);
 945   return ret;
 946 }