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 "MEDCouplingNormalizedUnstructuredMesh.txx"
  28
  29 #include "Interpolation1D.txx"
  30 #include "Interpolation2DCurve.hxx"
  31 #include "Interpolation2D.txx"
  32 #include "Interpolation3D.txx"
  33 #include "Interpolation3DSurf.hxx"
  34 #include "Interpolation2D1D.txx"
  35 #include "Interpolation3D2D.txx"
  36
  37 using namespace ParaMEDMEM;
  38
  39 MEDCouplingRemapper::MEDCouplingRemapper():_src_mesh(0),_target_mesh(0),_nature_of_deno(NoNature),_time_deno_update(0)
  40 {
  41 }
  42
  43 MEDCouplingRemapper::~MEDCouplingRemapper()
  44 {
  45   releaseData(false);
  46 }
  47
  48 int MEDCouplingRemapper::prepare(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const char *method) throw(INTERP_KERNEL::Exception)
  49 {
  50   releaseData(true);
  51   _src_mesh=const_cast<MEDCouplingMesh *>(srcMesh); _target_mesh=const_cast<MEDCouplingMesh *>(targetMesh);
  52   _src_mesh->incrRef(); _target_mesh->incrRef();
  53   int meshInterpType=((int)_src_mesh->getType()*16)+(int)_target_mesh->getType();
  54   switch(meshInterpType)
  55     {
  56     case 85://Unstructured-Unstructured
  57       return prepareUU(method);
  58     case 136://Extruded-Extruded
  59       return prepareEE(method);
  60     default:
  61       throw INTERP_KERNEL::Exception("Not managed type of meshes !");
  62     }
  63 }
  64
  65 int MEDCouplingRemapper::prepareEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target) throw(INTERP_KERNEL::Exception)
  66 {
  67   std::string meth(src->getDiscretization()->getStringRepr());
  68   meth+=target->getDiscretization()->getStringRepr();
  69   return prepare(src->getMesh(),target->getMesh(),meth.c_str());
  70 }
  71
  72 /*!
  73  * This method performs the operation source to target using matrix computed in ParaMEDMEM::MEDCouplingRemapper::prepare method.
  74  * 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.
  75  *
  76  * \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.
  77  * \param [out] targetField the destination field with the allocated array in which all tuples will be overwritten.
  78  */
  79 void MEDCouplingRemapper::transfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, double dftValue) throw(INTERP_KERNEL::Exception)
  80 {
  81   transferUnderground(srcField,targetField,true,dftValue);
  82 }
  83
  84 /*!
  85  * This method is equivalent to ParaMEDMEM::MEDCouplingRemapper::transfer except that here \b targetField is a in/out parameter.
  86  * 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
  87  * let unchanged.
  88  * This method requires that \b targetField was fully defined and allocated. If the array is not allocated an exception will be thrown.
  89  *
  90  * \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.
  91  * \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.
  92  */
  93 void MEDCouplingRemapper::partialTransfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField) throw(INTERP_KERNEL::Exception)
  94 {
  95   transferUnderground(srcField,targetField,false,std::numeric_limits<double>::max());
  96 }
  97
  98 void MEDCouplingRemapper::reverseTransfer(MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *targetField, double dftValue) throw(INTERP_KERNEL::Exception)
  99 {
 100   if(_src_method!=srcField->getDiscretization()->getStringRepr())
 101     throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
 102   if(_target_method!=targetField->getDiscretization()->getStringRepr())
 103     throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
 104   if(srcField->getNature()!=targetField->getNature())
 105     throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
 106   DataArrayDouble *array=srcField->getArray();
 107   int trgNbOfCompo=targetField->getNumberOfComponents();
 108   if(array)
 109     {
 110       if(trgNbOfCompo!=srcField->getNumberOfComponents())
 111         throw INTERP_KERNEL::Exception("Number of components mismatch !");
 112     }
 113   else
 114     {
 115       array=DataArrayDouble::New();
 116       array->alloc(srcField->getNumberOfTuples(),trgNbOfCompo);
 117       srcField->setArray(array);
 118       array->decrRef();
 119     }
 120   computeDeno(srcField->getNature(),srcField,targetField);
 121   double *resPointer=array->getPointer();
 122   const double *inputPointer=targetField->getArray()->getConstPointer();
 123   computeReverseProduct(inputPointer,trgNbOfCompo,dftValue,resPointer);
 124 }
 125
 126 MEDCouplingFieldDouble *MEDCouplingRemapper::transferField(const MEDCouplingFieldDouble *srcField, double dftValue) throw(INTERP_KERNEL::Exception)
 127 {
 128   if(_src_method!=srcField->getDiscretization()->getStringRepr())
 129     throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
 130   MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(MEDCouplingFieldDiscretization::getTypeOfFieldFromStringRepr(_target_method.c_str()),srcField->getTimeDiscretization());
 131   ret->copyTinyAttrFrom(srcField);
 132   ret->setNature(srcField->getNature());
 133   ret->setMesh(_target_mesh);
 134   transfer(srcField,ret,dftValue);
 135   return ret;
 136 }
 137
 138 MEDCouplingFieldDouble *MEDCouplingRemapper::reverseTransferField(const MEDCouplingFieldDouble *targetField, double dftValue) throw(INTERP_KERNEL::Exception)
 139 {
 140   if(_target_method!=targetField->getDiscretization()->getStringRepr())
 141     throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
 142   MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(MEDCouplingFieldDiscretization::getTypeOfFieldFromStringRepr(_target_method.c_str()),targetField->getTimeDiscretization());
 143   ret->copyTinyAttrFrom(targetField);
 144   ret->setNature(targetField->getNature());
 145   ret->setMesh(_src_mesh);
 146   reverseTransfer(ret,targetField,dftValue);
 147   return ret;
 148 }
 149
 150 bool MEDCouplingRemapper::setOptionInt(const std::string& key, int value)
 151 {
 152   return INTERP_KERNEL::InterpolationOptions::setOptionInt(key,value);
 153 }
 154
 155 bool MEDCouplingRemapper::setOptionDouble(const std::string& key, double value)
 156 {
 157   return INTERP_KERNEL::InterpolationOptions::setOptionDouble(key,value);
 158 }
 159
 160 bool MEDCouplingRemapper::setOptionString(const std::string& key, const std::string& value)
 161 {
 162   return INTERP_KERNEL::InterpolationOptions::setOptionString(key,value);
 163 }
 164
 165 int MEDCouplingRemapper::prepareUU(const char *method) throw(INTERP_KERNEL::Exception)
 166 {
 167   MEDCouplingUMesh *src_mesh=(MEDCouplingUMesh *)_src_mesh;
 168   MEDCouplingUMesh *target_mesh=(MEDCouplingUMesh *)_target_mesh;
 169   INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::checkAndSplitInterpolationMethod(method,_src_method,_target_method);
 170   const int srcMeshDim=src_mesh->getMeshDimension();
 171   int srcSpaceDim=-1;
 172   if(srcMeshDim!=-1)
 173     srcSpaceDim=src_mesh->getSpaceDimension();
 174   const int trgMeshDim=target_mesh->getMeshDimension();
 175   int trgSpaceDim=-1;
 176   if(trgMeshDim!=-1)
 177     trgSpaceDim=target_mesh->getSpaceDimension();
 178   if(trgSpaceDim!=srcSpaceDim)
 179     if(trgSpaceDim!=-1 && srcSpaceDim!=-1)
 180       throw INTERP_KERNEL::Exception("Incoherent space dimension detected between target and source.");
 181   int nbCols;
 182   if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==1)
 183     {
 184       MEDCouplingNormalizedUnstructuredMesh<1,1> source_mesh_wrapper(src_mesh);
 185       MEDCouplingNormalizedUnstructuredMesh<1,1> target_mesh_wrapper(target_mesh);
 186       INTERP_KERNEL::Interpolation1D interpolation(*this);
 187       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 188     }
 189   else if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==2)
 190     {
 191       MEDCouplingNormalizedUnstructuredMesh<2,1> source_mesh_wrapper(src_mesh);
 192       MEDCouplingNormalizedUnstructuredMesh<2,1> target_mesh_wrapper(target_mesh);
 193       INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
 194       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 195     }
 196   else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==2)
 197     {
 198       MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 199       MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
 200       INTERP_KERNEL::Interpolation2D interpolation(*this);
 201       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 202     }
 203   else if(srcMeshDim==3 && trgMeshDim==3 && srcSpaceDim==3)
 204     {
 205       MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 206       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
 207       INTERP_KERNEL::Interpolation3D interpolation(*this);
 208       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 209     }
 210   else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==3)
 211     {
 212       MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
 213       MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh);
 214       INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
 215       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 216     }
 217   else if(srcMeshDim==3 && trgMeshDim==1 && srcSpaceDim==3)
 218     {
 219       if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
 220         throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
 221       MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 222       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
 223       INTERP_KERNEL::Interpolation3D interpolation(*this);
 224       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 225     }
 226   else if(srcMeshDim==1 && trgMeshDim==3 && srcSpaceDim==3)
 227     {
 228       if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
 229         throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
 230       MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 231       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
 232       INTERP_KERNEL::Interpolation3D interpolation(*this);
 233       std::vector<std::map<int,double> > matrixTmp;
 234       nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
 235       reverseMatrix(matrixTmp,nbCols,_matrix);
 236       nbCols=matrixTmp.size();
 237     }
 238   else if(srcMeshDim==2 && trgMeshDim==1 && srcSpaceDim==2)
 239     {
 240       if(getIntersectionType()==INTERP_KERNEL::PointLocator)
 241         {
 242           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 243           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
 244           INTERP_KERNEL::Interpolation2D interpolation(*this);
 245           nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 246         }
 247       else
 248         {
 249           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 250           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
 251           INTERP_KERNEL::Interpolation2D1D interpolation(*this);
 252           std::vector<std::map<int,double> > matrixTmp;
 253           nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
 254           reverseMatrix(matrixTmp,nbCols,_matrix);
 255           nbCols=matrixTmp.size();
 256           INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
 257           if(!duplicateFaces.empty())
 258             {
 259               std::ostringstream oss; oss << "An unexpected situation happend ! For the following 1D Cells are part of edges shared by 2D cells :\n";
 260               for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
 261                 {
 262                   oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
 263                   std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
 264                   oss << std::endl;
 265                 }
 266             }
 267         }
 268     }
 269   else if(srcMeshDim==1 && trgMeshDim==2 && srcSpaceDim==2)
 270     {
 271       if(getIntersectionType()==INTERP_KERNEL::PointLocator)
 272         {
 273           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 274           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
 275           INTERP_KERNEL::Interpolation2D interpolation(*this);
 276           std::vector<std::map<int,double> > matrixTmp;
 277           nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
 278           reverseMatrix(matrixTmp,nbCols,_matrix);
 279           nbCols=matrixTmp.size();
 280         }
 281       else
 282         {
 283           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 284           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
 285           INTERP_KERNEL::Interpolation2D1D interpolation(*this);
 286           nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 287           INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
 288           if(!duplicateFaces.empty())
 289             {
 290               std::ostringstream oss; oss << "An unexpected situation happend ! For the following 1D Cells are part of edges shared by 2D cells :\n";
 291               for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
 292                 {
 293                   oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
 294                   std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
 295                   oss << std::endl;
 296                 }
 297             }
 298         }
 299     }
 300   else if(srcMeshDim==2 && trgMeshDim==3 && srcSpaceDim==3)
 301     {
 302       MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 303       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
 304       INTERP_KERNEL::Interpolation3D2D interpolation(*this);
 305       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 306       INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
 307       if(!duplicateFaces.empty())
 308         {
 309           std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
 310           for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
 311             {
 312               oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
 313               std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
 314               oss << std::endl;
 315             }
 316         }
 317     }
 318   else if(srcMeshDim==3 && trgMeshDim==2 && srcSpaceDim==3)
 319     {
 320       MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 321       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
 322       INTERP_KERNEL::Interpolation3D2D interpolation(*this);
 323       std::vector<std::map<int,double> > matrixTmp;
 324       nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
 325       reverseMatrix(matrixTmp,nbCols,_matrix);
 326       nbCols=matrixTmp.size();
 327       INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
 328       if(!duplicateFaces.empty())
 329         {
 330           std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
 331           for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
 332             {
 333               oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
 334               std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
 335               oss << std::endl;
 336             }
 337         }
 338     }
 339   else if(trgMeshDim==-1)
 340     {
 341       if(srcMeshDim==2 && srcSpaceDim==2)
 342         {
 343           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 344           INTERP_KERNEL::Interpolation2D interpolation(*this);
 345           nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,_src_method.c_str());
 346         }
 347       else if(srcMeshDim==3 && srcSpaceDim==3)
 348         {
 349           MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 350           INTERP_KERNEL::Interpolation3D interpolation(*this);
 351           nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,_src_method.c_str());
 352         }
 353       else if(srcMeshDim==2 && srcSpaceDim==3)
 354         {
 355           MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
 356           INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
 357           nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,_src_method.c_str());
 358         }
 359       else
 360         throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
 361     }
 362   else if(srcMeshDim==-1)
 363     {
 364       if(trgMeshDim==2 && trgSpaceDim==2)
 365         {
 366           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(target_mesh);
 367           INTERP_KERNEL::Interpolation2D interpolation(*this);
 368           nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,_target_method.c_str());
 369         }
 370       else if(trgMeshDim==3 && trgSpaceDim==3)
 371         {
 372           MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(target_mesh);
 373           INTERP_KERNEL::Interpolation3D interpolation(*this);
 374           nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,_target_method.c_str());
 375         }
 376       else if(trgMeshDim==2 && trgSpaceDim==3)
 377         {
 378           MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(target_mesh);
 379           INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
 380           nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,_target_method.c_str());
 381         }
 382       else
 383         throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh from -1D sourceMesh");
 384     }
 385   else
 386     throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension");
 387   _deno_multiply.clear();
 388   _deno_multiply.resize(_matrix.size());
 389   _deno_reverse_multiply.clear();
 390   _deno_reverse_multiply.resize(nbCols);
 391   declareAsNew();
 392   return 1;
 393 }
 394
 395 int MEDCouplingRemapper::prepareEE(const char *method) throw(INTERP_KERNEL::Exception)
 396 {
 397   MEDCouplingExtrudedMesh *src_mesh=(MEDCouplingExtrudedMesh *)_src_mesh;
 398   MEDCouplingExtrudedMesh *target_mesh=(MEDCouplingExtrudedMesh *)_target_mesh;
 399   std::string methC(method);
 400   if(methC!="P0P0")
 401     throw INTERP_KERNEL::Exception("Only P0P0 method implemented for Extruded/Extruded meshes !");
 402   INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::checkAndSplitInterpolationMethod(method,_src_method,_target_method);
 403   MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh->getMesh2D());
 404   MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh->getMesh2D());
 405   INTERP_KERNEL::Interpolation3DSurf interpolation2D(*this);
 406   std::vector<std::map<int,double> > matrix2D;
 407   int nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,method);
 408   MEDCouplingUMesh *s1D,*t1D;
 409   double v[3];
 410   MEDCouplingExtrudedMesh::Project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v);
 411   MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D);
 412   MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
 413   std::vector<std::map<int,double> > matrix1D;
 414   INTERP_KERNEL::Interpolation1D interpolation1D(*this);
 415   int nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,method);
 416   s1D->decrRef();
 417   t1D->decrRef();
 418   buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
 419                                              target_mesh->getMesh3DIds()->getConstPointer());
 420   //
 421   _deno_multiply.clear();
 422   _deno_multiply.resize(_matrix.size());
 423   _deno_reverse_multiply.clear();
 424   _deno_reverse_multiply.resize(nbCols2D*nbCols1D);
 425   declareAsNew();
 426   return 1;
 427 }
 428
 429 void MEDCouplingRemapper::updateTime() const
 430 {
 431 }
 432
 433 void MEDCouplingRemapper::releaseData(bool matrixSuppression)
 434 {
 435   if(_src_mesh)
 436     _src_mesh->decrRef();
 437   if(_target_mesh)
 438     _target_mesh->decrRef();
 439   _src_mesh=0;
 440   _target_mesh=0;
 441   if(matrixSuppression)
 442     {
 443       _matrix.clear();
 444       _deno_multiply.clear();
 445       _deno_reverse_multiply.clear();
 446     }
 447 }
 448
 449 void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue) throw(INTERP_KERNEL::Exception)
 450 {
 451   if(_src_method!=srcField->getDiscretization()->getStringRepr())
 452     throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
 453   if(_target_method!=targetField->getDiscretization()->getStringRepr())
 454     throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
 455   if(srcField->getNature()!=targetField->getNature())
 456     throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
 457   DataArrayDouble *array=targetField->getArray();
 458   int srcNbOfCompo=srcField->getNumberOfComponents();
 459   if(array)
 460     {
 461       if(srcNbOfCompo!=targetField->getNumberOfComponents())
 462         throw INTERP_KERNEL::Exception("Number of components mismatch !");
 463     }
 464   else
 465     {
 466       if(!isDftVal)
 467         throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !");
 468       array=DataArrayDouble::New();
 469       array->alloc(targetField->getNumberOfTuples(),srcNbOfCompo);
 470       targetField->setArray(array);
 471       array->decrRef();
 472     }
 473   computeDeno(srcField->getNature(),srcField,targetField);
 474   double *resPointer=array->getPointer();
 475   const double *inputPointer=srcField->getArray()->getConstPointer();
 476   computeProduct(inputPointer,srcNbOfCompo,isDftVal,dftValue,resPointer);
 477 }
 478
 479 void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
 480 {
 481   if(nat==NoNature)
 482     return computeDenoFromScratch(nat,srcField,trgField);
 483   else if(nat!=_nature_of_deno)
 484    return computeDenoFromScratch(nat,srcField,trgField);
 485   else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis())
 486     return computeDenoFromScratch(nat,srcField,trgField);
 487 }
 488
 489 void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField) throw(INTERP_KERNEL::Exception)
 490 {
 491   _nature_of_deno=nat;
 492   _time_deno_update=getTimeOfThis();
 493   switch(_nature_of_deno)
 494     {
 495     case ConservativeVolumic:
 496       {
 497         computeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
 498         break;
 499       }
 500     case Integral:
 501       {
 502         MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),true);
 503         MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),true);
 504         const double *denoPtr=deno->getArray()->getConstPointer();
 505         const double *denoRPtr=denoR->getArray()->getConstPointer();
 506         if(trgField->getMesh()->getMeshDimension()==-1)
 507           {
 508             double *denoRPtr2=denoR->getArray()->getPointer();
 509             denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
 510           }
 511         if(srcField->getMesh()->getMeshDimension()==-1)
 512           {
 513             double *denoPtr2=deno->getArray()->getPointer();
 514             denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
 515           }
 516         int idx=0;
 517         for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
 518           for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 519             {
 520               _deno_multiply[idx][(*iter2).first]=denoPtr[(*iter2).first];
 521               _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[idx];
 522             }
 523         deno->decrRef();
 524         denoR->decrRef();
 525         break;
 526       }
 527     case IntegralGlobConstraint:
 528       {
 529         computeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
 530         break;
 531       }
 532     case RevIntegral:
 533       {
 534         MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),true);
 535         MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),true);
 536         const double *denoPtr=deno->getArray()->getConstPointer();
 537         const double *denoRPtr=denoR->getArray()->getConstPointer();
 538         if(trgField->getMesh()->getMeshDimension()==-1)
 539           {
 540             double *denoRPtr2=denoR->getArray()->getPointer();
 541             denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
 542           }
 543         if(srcField->getMesh()->getMeshDimension()==-1)
 544           {
 545             double *denoPtr2=deno->getArray()->getPointer();
 546             denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
 547           }
 548         int idx=0;
 549         for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
 550           for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 551             {
 552               _deno_multiply[idx][(*iter2).first]=denoPtr[idx];
 553               _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[(*iter2).first];
 554             }
 555         deno->decrRef();
 556         denoR->decrRef();
 557         break;
 558       }
 559     case NoNature:
 560       throw INTERP_KERNEL::Exception("No nature specified ! Select one !");
 561     }
 562 }
 563
 564 void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)
 565 {
 566   int idx=0;
 567   double *tmp=new double[inputNbOfCompo];
 568   for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
 569     {
 570       if((*iter1).empty())
 571         {
 572           if(isDftVal)
 573             std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
 574           continue;
 575         }
 576       else
 577         std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,0.);
 578       std::map<int,double>::const_iterator iter3=_deno_multiply[idx].begin();
 579       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter3++)
 580         {
 581           std::transform(inputPointer+(*iter2).first*inputNbOfCompo,inputPointer+((*iter2).first+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/(*iter3).second));
 582           std::transform(tmp,tmp+inputNbOfCompo,resPointer+idx*inputNbOfCompo,resPointer+idx*inputNbOfCompo,std::plus<double>());
 583         }
 584     }
 585   delete [] tmp;
 586 }
 587
 588 void MEDCouplingRemapper::computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer)
 589 {
 590   std::vector<bool> isReached(_deno_reverse_multiply.size(),false);
 591   int idx=0;
 592   double *tmp=new double[inputNbOfCompo];
 593   std::fill(resPointer,resPointer+inputNbOfCompo*_deno_reverse_multiply.size(),0.);
 594   for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
 595     {
 596       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 597         {
 598           isReached[(*iter2).first]=true;
 599           std::transform(inputPointer+idx*inputNbOfCompo,inputPointer+(idx+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/_deno_reverse_multiply[(*iter2).first][idx]));
 600           std::transform(tmp,tmp+inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,std::plus<double>());
 601         }
 602     }
 603   delete [] tmp;
 604   idx=0;
 605   for(std::vector<bool>::const_iterator iter3=isReached.begin();iter3!=isReached.end();iter3++,idx++)
 606     if(!*iter3)
 607       std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
 608 }
 609
 610 void MEDCouplingRemapper::reverseMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
 611 {
 612   matOut.resize(nbColsMatIn);
 613   int id=0;
 614   for(std::vector<std::map<int,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
 615     for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 616       matOut[(*iter2).first][id]=(*iter2).second;
 617 }
 618
 619 void MEDCouplingRemapper::computeRowSumAndColSum(const std::vector<std::map<int,double> >& matrixDeno,
 620                                                  std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
 621 {
 622   std::map<int,double> values;
 623   int idx=0;
 624   for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
 625     {
 626       double sum=0.;
 627       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 628         {
 629           sum+=(*iter2).second;
 630           values[(*iter2).first]+=(*iter2).second;
 631         }
 632       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 633         deno[idx][(*iter2).first]=sum;
 634     }
 635   idx=0;
 636   for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
 637     {
 638       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 639         denoReverse[(*iter2).first][idx]=values[(*iter2).first];
 640     }
 641 }
 642
 643 void MEDCouplingRemapper::computeColSumAndRowSum(const std::vector<std::map<int,double> >& matrixDeno,
 644                                                  std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
 645 {
 646   std::map<int,double> values;
 647   int idx=0;
 648   for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
 649     {
 650       double sum=0.;
 651       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 652         {
 653           sum+=(*iter2).second;
 654           values[(*iter2).first]+=(*iter2).second;
 655         }
 656       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 657         denoReverse[(*iter2).first][idx]=sum;
 658     }
 659   idx=0;
 660   for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
 661     {
 662       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 663         deno[idx][(*iter2).first]=values[(*iter2).first];
 664     }
 665 }
 666
 667 void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<int,double> >& m1D,
 668                                                                      const std::vector< std::map<int,double> >& m2D,
 669                                                                      const int *corrCellIdSrc, int nbOf2DCellsSrc, int nbOf1DCellsSrc,
 670                                                                      const int *corrCellIdTrg)
 671 {
 672   int nbOf2DCellsTrg=m2D.size();
 673   int nbOf1DCellsTrg=m1D.size();
 674   int nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
 675   _matrix.resize(nbOf3DCellsTrg);
 676   int id2R=0;
 677   for(std::vector< std::map<int,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
 678     {
 679       for(std::map<int,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
 680         {
 681           int id1R=0;
 682           for(std::vector< std::map<int,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
 683             {
 684               for(std::map<int,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
 685                 {
 686                   _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
 687                 }
 688             }
 689         }
 690     }
 691 }
 692
 693 void MEDCouplingRemapper::PrintMatrix(const std::vector<std::map<int,double> >& m)
 694 {
 695   int id=0;
 696   for(std::vector<std::map<int,double> >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++)
 697     {
 698       std::cout << "Target Cell # " << id << " : ";
 699       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 700         std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), ";
 701       std::cout << std::endl;
 702     }
 703 }
 704
 705 const std::vector<std::map<int,double> >& MEDCouplingRemapper::getCrudeMatrix() const
 706 {
 707   return _matrix;
 708 }
 709
 710 /*!
 711  * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
 712  * If not the behaviour is unpredictable.
 713  * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than \a maxValAbs this coefficient is
 714  * 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.
 715  * 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
 716  * can occur.
 717  *
 718  * \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.
 719  * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged.
 720  * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix
 721  */
 722 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs(double maxValAbs) throw(INTERP_KERNEL::Exception)
 723 {
 724   int ret=0;
 725   std::vector<std::map<int,double> > matrixNew(_matrix.size());
 726   int i=0;
 727   for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++,i++)
 728     {
 729       std::map<int,double>& rowNew=matrixNew[i];
 730       for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 731         {
 732           if(fabs((*it2).second)>maxValAbs)
 733             rowNew[(*it2).first]=(*it2).second;
 734           else
 735             ret++;
 736         }
 737     }
 738   if(ret>0)
 739     _matrix=matrixNew;
 740   return ret;
 741 }
 742
 743 /*!
 744  * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
 745  * If not the behaviour is unpredictable.
 746  * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than delta multiplied by \a scaleFactor this coefficient is
 747  * 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.
 748  * delta is the value returned by MEDCouplingRemapper::getMaxValueInCrudeMatrix method.
 749  * 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
 750  * can occur.
 751  *
 752  * \param [in] scaleFactor is the scale factor from which coefficients lower than \a scaleFactor times range width of coefficients are set to zero.
 753  * \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
 754  *         that all coefficients are null.
 755  * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs
 756  */
 757 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix(double scaleFactor) throw(INTERP_KERNEL::Exception)
 758 {
 759   double maxVal=getMaxValueInCrudeMatrix();
 760   if(maxVal==0.)
 761     return -1;
 762   return nullifiedTinyCoeffInCrudeMatrixAbs(scaleFactor*maxVal);
 763 }
 764
 765 /*!
 766  * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
 767  * If not the behaviour is unpredictable.
 768  * This method returns the maximum of the absolute values of coefficients into the sparse crude matrix.
 769  * The returned value is positive.
 770  */
 771 double MEDCouplingRemapper::getMaxValueInCrudeMatrix() const throw(INTERP_KERNEL::Exception)
 772 {
 773   double ret=0.;
 774   for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++)
 775     for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
 776       if(fabs((*it2).second)>ret)
 777         ret=fabs((*it2).second);
 778   return ret;
 779 }