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