Back to index

salome-med  6.5.0
MEDCouplingRemapper.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 #include "MEDCouplingRemapper.hxx"
00021 #include "MEDCouplingMemArray.hxx"
00022 #include "MEDCouplingFieldDouble.hxx"
00023 #include "MEDCouplingFieldTemplate.hxx"
00024 #include "MEDCouplingFieldDiscretization.hxx"
00025 #include "MEDCouplingExtrudedMesh.hxx"
00026 #include "MEDCouplingNormalizedUnstructuredMesh.txx"
00027 
00028 #include "Interpolation1D.txx"
00029 #include "Interpolation2DCurve.hxx"
00030 #include "Interpolation2D.txx"
00031 #include "Interpolation3D.txx"
00032 #include "Interpolation3DSurf.hxx"
00033 #include "Interpolation2D1D.txx"
00034 #include "Interpolation3D2D.txx"
00035 
00036 using namespace ParaMEDMEM;
00037 
00038 MEDCouplingRemapper::MEDCouplingRemapper():_src_mesh(0),_target_mesh(0),_nature_of_deno(NoNature),_time_deno_update(0)
00039 {
00040 }
00041 
00042 MEDCouplingRemapper::~MEDCouplingRemapper()
00043 {
00044   releaseData(false);
00045 }
00046 
00047 int MEDCouplingRemapper::prepare(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const char *method) throw(INTERP_KERNEL::Exception)
00048 {
00049   releaseData(true);
00050   _src_mesh=const_cast<MEDCouplingMesh *>(srcMesh); _target_mesh=const_cast<MEDCouplingMesh *>(targetMesh);
00051   _src_mesh->incrRef(); _target_mesh->incrRef();
00052   int meshInterpType=((int)_src_mesh->getType()*16)+(int)_target_mesh->getType();
00053   switch(meshInterpType)
00054     {
00055     case 85://Unstructured-Unstructured
00056       return prepareUU(method);
00057     case 136://Extruded-Extruded
00058       return prepareEE(method);
00059     default:
00060       throw INTERP_KERNEL::Exception("Not managed type of meshes !");
00061     }
00062 }
00063 
00064 int MEDCouplingRemapper::prepareEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target) throw(INTERP_KERNEL::Exception)
00065 {
00066   std::string meth(src->getDiscretization()->getStringRepr());
00067   meth+=target->getDiscretization()->getStringRepr();
00068   return prepare(src->getMesh(),target->getMesh(),meth.c_str());
00069 }
00070 
00078 void MEDCouplingRemapper::transfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, double dftValue) throw(INTERP_KERNEL::Exception)
00079 {
00080   transferUnderground(srcField,targetField,true,dftValue);
00081 }
00082 
00092 void MEDCouplingRemapper::partialTransfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField) throw(INTERP_KERNEL::Exception)
00093 {
00094   transferUnderground(srcField,targetField,false,std::numeric_limits<double>::max());
00095 }
00096 
00097 void MEDCouplingRemapper::reverseTransfer(MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *targetField, double dftValue) throw(INTERP_KERNEL::Exception)
00098 {
00099   if(_src_method!=srcField->getDiscretization()->getStringRepr())
00100     throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
00101   if(_target_method!=targetField->getDiscretization()->getStringRepr())
00102     throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
00103   if(srcField->getNature()!=targetField->getNature())
00104     throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
00105   DataArrayDouble *array=srcField->getArray();
00106   int trgNbOfCompo=targetField->getNumberOfComponents();
00107   if(array)
00108     {
00109       if(trgNbOfCompo!=srcField->getNumberOfComponents())
00110         throw INTERP_KERNEL::Exception("Number of components mismatch !");
00111     }
00112   else
00113     {
00114       array=DataArrayDouble::New();
00115       array->alloc(srcField->getNumberOfTuples(),trgNbOfCompo);
00116       srcField->setArray(array);
00117       array->decrRef();
00118     }
00119   computeDeno(srcField->getNature(),srcField,targetField);
00120   double *resPointer=array->getPointer();
00121   const double *inputPointer=targetField->getArray()->getConstPointer();
00122   computeReverseProduct(inputPointer,trgNbOfCompo,dftValue,resPointer);
00123 }
00124 
00125 MEDCouplingFieldDouble *MEDCouplingRemapper::transferField(const MEDCouplingFieldDouble *srcField, double dftValue) throw(INTERP_KERNEL::Exception)
00126 {
00127   if(_src_method!=srcField->getDiscretization()->getStringRepr())
00128     throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
00129   MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(MEDCouplingFieldDiscretization::getTypeOfFieldFromStringRepr(_target_method.c_str()),srcField->getTimeDiscretization());
00130   ret->copyTinyAttrFrom(srcField);
00131   ret->setNature(srcField->getNature());
00132   ret->setMesh(_target_mesh);
00133   transfer(srcField,ret,dftValue);
00134   return ret;
00135 }
00136 
00137 MEDCouplingFieldDouble *MEDCouplingRemapper::reverseTransferField(const MEDCouplingFieldDouble *targetField, double dftValue) throw(INTERP_KERNEL::Exception)
00138 {
00139   if(_target_method!=targetField->getDiscretization()->getStringRepr())
00140     throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
00141   MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(MEDCouplingFieldDiscretization::getTypeOfFieldFromStringRepr(_target_method.c_str()),targetField->getTimeDiscretization());
00142   ret->copyTinyAttrFrom(targetField);
00143   ret->setNature(targetField->getNature());
00144   ret->setMesh(_src_mesh);
00145   reverseTransfer(ret,targetField,dftValue);
00146   return ret;
00147 }
00148 
00149 bool MEDCouplingRemapper::setOptionInt(const std::string& key, int value)
00150 {
00151   return INTERP_KERNEL::InterpolationOptions::setOptionInt(key,value);
00152 }
00153 
00154 bool MEDCouplingRemapper::setOptionDouble(const std::string& key, double value)
00155 {
00156   return INTERP_KERNEL::InterpolationOptions::setOptionDouble(key,value);
00157 }
00158 
00159 bool MEDCouplingRemapper::setOptionString(const std::string& key, const std::string& value)
00160 {
00161   return INTERP_KERNEL::InterpolationOptions::setOptionString(key,value);
00162 }
00163 
00164 int MEDCouplingRemapper::prepareUU(const char *method) throw(INTERP_KERNEL::Exception)
00165 {
00166   MEDCouplingUMesh *src_mesh=(MEDCouplingUMesh *)_src_mesh;
00167   MEDCouplingUMesh *target_mesh=(MEDCouplingUMesh *)_target_mesh;
00168   INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::checkAndSplitInterpolationMethod(method,_src_method,_target_method);
00169   const int srcMeshDim=src_mesh->getMeshDimension();
00170   int srcSpaceDim=-1;
00171   if(srcMeshDim!=-1)
00172     srcSpaceDim=src_mesh->getSpaceDimension();
00173   const int trgMeshDim=target_mesh->getMeshDimension();
00174   int trgSpaceDim=-1;
00175   if(trgMeshDim!=-1)
00176     trgSpaceDim=target_mesh->getSpaceDimension();
00177   if(trgSpaceDim!=srcSpaceDim)
00178     if(trgSpaceDim!=-1 && srcSpaceDim!=-1)
00179       throw INTERP_KERNEL::Exception("Incoherent space dimension detected between target and source.");
00180   int nbCols;
00181   if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==1)
00182     {
00183       MEDCouplingNormalizedUnstructuredMesh<1,1> source_mesh_wrapper(src_mesh);
00184       MEDCouplingNormalizedUnstructuredMesh<1,1> target_mesh_wrapper(target_mesh);
00185       INTERP_KERNEL::Interpolation1D interpolation(*this);
00186       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
00187     }
00188   else if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==2)
00189     {
00190       MEDCouplingNormalizedUnstructuredMesh<2,1> source_mesh_wrapper(src_mesh);
00191       MEDCouplingNormalizedUnstructuredMesh<2,1> target_mesh_wrapper(target_mesh);
00192       INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
00193       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
00194     }
00195   else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==2)
00196     {
00197       MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
00198       MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
00199       INTERP_KERNEL::Interpolation2D interpolation(*this);
00200       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
00201     }
00202   else if(srcMeshDim==3 && trgMeshDim==3 && srcSpaceDim==3)
00203     {
00204       MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
00205       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
00206       INTERP_KERNEL::Interpolation3D interpolation(*this);
00207       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
00208     }
00209   else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==3)
00210     {
00211       MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
00212       MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh);
00213       INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
00214       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
00215     }
00216   else if(srcMeshDim==3 && trgMeshDim==1 && srcSpaceDim==3)
00217     {
00218       if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
00219         throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
00220       MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
00221       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
00222       INTERP_KERNEL::Interpolation3D interpolation(*this);
00223       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
00224     }
00225   else if(srcMeshDim==1 && trgMeshDim==3 && srcSpaceDim==3)
00226     {
00227       if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
00228         throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
00229       MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
00230       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
00231       INTERP_KERNEL::Interpolation3D interpolation(*this);
00232       std::vector<std::map<int,double> > matrixTmp;
00233       nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
00234       reverseMatrix(matrixTmp,nbCols,_matrix);
00235       nbCols=matrixTmp.size();
00236     }
00237   else if(srcMeshDim==2 && trgMeshDim==1 && srcSpaceDim==2)
00238     {
00239       if(getIntersectionType()==INTERP_KERNEL::PointLocator)
00240         {
00241           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
00242           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
00243           INTERP_KERNEL::Interpolation2D interpolation(*this);
00244           nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
00245         }
00246       else
00247         {
00248           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
00249           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
00250           INTERP_KERNEL::Interpolation2D1D interpolation(*this);
00251           std::vector<std::map<int,double> > matrixTmp;
00252           nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
00253           reverseMatrix(matrixTmp,nbCols,_matrix);
00254           nbCols=matrixTmp.size();
00255           INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
00256           if(!duplicateFaces.empty())
00257             {
00258               std::ostringstream oss; oss << "An unexpected situation happend ! For the following 1D Cells are part of edges shared by 2D cells :\n";
00259               for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
00260                 {
00261                   oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
00262                   std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
00263                   oss << std::endl;
00264                 }
00265             }
00266         }
00267     }
00268   else if(srcMeshDim==1 && trgMeshDim==2 && srcSpaceDim==2)
00269     {
00270       if(getIntersectionType()==INTERP_KERNEL::PointLocator)
00271         {
00272           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
00273           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
00274           INTERP_KERNEL::Interpolation2D interpolation(*this);
00275           std::vector<std::map<int,double> > matrixTmp;
00276           nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
00277           reverseMatrix(matrixTmp,nbCols,_matrix);
00278           nbCols=matrixTmp.size();
00279         }
00280       else
00281         {
00282           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
00283           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
00284           INTERP_KERNEL::Interpolation2D1D interpolation(*this);
00285           nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
00286           INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
00287           if(!duplicateFaces.empty())
00288             {
00289               std::ostringstream oss; oss << "An unexpected situation happend ! For the following 1D Cells are part of edges shared by 2D cells :\n";
00290               for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
00291                 {
00292                   oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
00293                   std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
00294                   oss << std::endl;
00295                 }
00296             }
00297         }
00298     }
00299   else if(srcMeshDim==2 && trgMeshDim==3 && srcSpaceDim==3)
00300     {
00301       MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
00302       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
00303       INTERP_KERNEL::Interpolation3D2D interpolation(*this);
00304       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
00305       INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
00306       if(!duplicateFaces.empty())
00307         {
00308           std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
00309           for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
00310             {
00311               oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
00312               std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
00313               oss << std::endl;
00314             }
00315         }
00316     }
00317   else if(srcMeshDim==3 && trgMeshDim==2 && srcSpaceDim==3)
00318     {
00319       MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
00320       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
00321       INTERP_KERNEL::Interpolation3D2D interpolation(*this);
00322       std::vector<std::map<int,double> > matrixTmp;
00323       nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
00324       reverseMatrix(matrixTmp,nbCols,_matrix);
00325       nbCols=matrixTmp.size();
00326       INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
00327       if(!duplicateFaces.empty())
00328         {
00329           std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
00330           for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
00331             {
00332               oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
00333               std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
00334               oss << std::endl;
00335             }
00336         }
00337     }
00338   else if(trgMeshDim==-1)
00339     {
00340       if(srcMeshDim==2 && srcSpaceDim==2)
00341         {
00342           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
00343           INTERP_KERNEL::Interpolation2D interpolation(*this);
00344           nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,_src_method.c_str());
00345         }
00346       else if(srcMeshDim==3 && srcSpaceDim==3)
00347         {
00348           MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
00349           INTERP_KERNEL::Interpolation3D interpolation(*this);
00350           nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,_src_method.c_str());
00351         }
00352       else if(srcMeshDim==2 && srcSpaceDim==3)
00353         {
00354           MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
00355           INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
00356           nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,_src_method.c_str());
00357         }
00358       else
00359         throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
00360     }
00361   else if(srcMeshDim==-1)
00362     {
00363       if(trgMeshDim==2 && trgSpaceDim==2)
00364         {
00365           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(target_mesh);
00366           INTERP_KERNEL::Interpolation2D interpolation(*this);
00367           nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,_target_method.c_str());
00368         }
00369       else if(trgMeshDim==3 && trgSpaceDim==3)
00370         {
00371           MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(target_mesh);
00372           INTERP_KERNEL::Interpolation3D interpolation(*this);
00373           nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,_target_method.c_str());
00374         }
00375       else if(trgMeshDim==2 && trgSpaceDim==3)
00376         {
00377           MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(target_mesh);
00378           INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
00379           nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,_target_method.c_str());
00380         }
00381       else
00382         throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh from -1D sourceMesh");
00383     }
00384   else
00385     throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension");
00386   _deno_multiply.clear();
00387   _deno_multiply.resize(_matrix.size());
00388   _deno_reverse_multiply.clear();
00389   _deno_reverse_multiply.resize(nbCols);
00390   declareAsNew();
00391   return 1;
00392 }
00393 
00394 int MEDCouplingRemapper::prepareEE(const char *method) throw(INTERP_KERNEL::Exception)
00395 {
00396   MEDCouplingExtrudedMesh *src_mesh=(MEDCouplingExtrudedMesh *)_src_mesh;
00397   MEDCouplingExtrudedMesh *target_mesh=(MEDCouplingExtrudedMesh *)_target_mesh;
00398   std::string methC(method);
00399   if(methC!="P0P0")
00400     throw INTERP_KERNEL::Exception("Only P0P0 method implemented for Extruded/Extruded meshes !");
00401   INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::checkAndSplitInterpolationMethod(method,_src_method,_target_method);
00402   MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh->getMesh2D());
00403   MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh->getMesh2D());
00404   INTERP_KERNEL::Interpolation3DSurf interpolation2D(*this);
00405   std::vector<std::map<int,double> > matrix2D;
00406   int nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,method);
00407   MEDCouplingUMesh *s1D,*t1D;
00408   double v[3];
00409   MEDCouplingExtrudedMesh::Project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v);
00410   MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D);
00411   MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
00412   std::vector<std::map<int,double> > matrix1D;
00413   INTERP_KERNEL::Interpolation1D interpolation1D(*this);
00414   int nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,method);
00415   s1D->decrRef();
00416   t1D->decrRef();
00417   buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
00418                                              target_mesh->getMesh3DIds()->getConstPointer());
00419   //
00420   _deno_multiply.clear();
00421   _deno_multiply.resize(_matrix.size());
00422   _deno_reverse_multiply.clear();
00423   _deno_reverse_multiply.resize(nbCols2D*nbCols1D);
00424   declareAsNew();
00425   return 1;
00426 }
00427 
00428 void MEDCouplingRemapper::updateTime() const
00429 {
00430 }
00431 
00432 void MEDCouplingRemapper::releaseData(bool matrixSuppression)
00433 {
00434   if(_src_mesh)
00435     _src_mesh->decrRef();
00436   if(_target_mesh)
00437     _target_mesh->decrRef();
00438   _src_mesh=0;
00439   _target_mesh=0;
00440   if(matrixSuppression)
00441     {
00442       _matrix.clear();
00443       _deno_multiply.clear();
00444       _deno_reverse_multiply.clear();
00445     }
00446 }
00447 
00448 void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue) throw(INTERP_KERNEL::Exception)
00449 {
00450   if(_src_method!=srcField->getDiscretization()->getStringRepr())
00451     throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
00452   if(_target_method!=targetField->getDiscretization()->getStringRepr())
00453     throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
00454   if(srcField->getNature()!=targetField->getNature())
00455     throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
00456   DataArrayDouble *array=targetField->getArray();
00457   int srcNbOfCompo=srcField->getNumberOfComponents();
00458   if(array)
00459     {
00460       if(srcNbOfCompo!=targetField->getNumberOfComponents())
00461         throw INTERP_KERNEL::Exception("Number of components mismatch !");
00462     }
00463   else
00464     {
00465       if(!isDftVal)
00466         throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !");
00467       array=DataArrayDouble::New();
00468       array->alloc(targetField->getNumberOfTuples(),srcNbOfCompo);
00469       targetField->setArray(array);
00470       array->decrRef();
00471     }
00472   computeDeno(srcField->getNature(),srcField,targetField);
00473   double *resPointer=array->getPointer();
00474   const double *inputPointer=srcField->getArray()->getConstPointer();
00475   computeProduct(inputPointer,srcNbOfCompo,isDftVal,dftValue,resPointer);
00476 }
00477 
00478 void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
00479 {
00480   if(nat==NoNature)
00481     return computeDenoFromScratch(nat,srcField,trgField);
00482   else if(nat!=_nature_of_deno)
00483    return computeDenoFromScratch(nat,srcField,trgField);
00484   else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis())
00485     return computeDenoFromScratch(nat,srcField,trgField);
00486 }
00487 
00488 void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField) throw(INTERP_KERNEL::Exception)
00489 {
00490   _nature_of_deno=nat;
00491   _time_deno_update=getTimeOfThis();
00492   switch(_nature_of_deno)
00493     {
00494     case ConservativeVolumic:
00495       {
00496         computeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
00497         break;
00498       }
00499     case Integral:
00500       {
00501         MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),true);
00502         MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),true);
00503         const double *denoPtr=deno->getArray()->getConstPointer();
00504         const double *denoRPtr=denoR->getArray()->getConstPointer();
00505         if(trgField->getMesh()->getMeshDimension()==-1)
00506           {
00507             double *denoRPtr2=denoR->getArray()->getPointer();
00508             denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
00509           }
00510         if(srcField->getMesh()->getMeshDimension()==-1)
00511           {
00512             double *denoPtr2=deno->getArray()->getPointer();
00513             denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
00514           }
00515         int idx=0;
00516         for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
00517           for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
00518             {
00519               _deno_multiply[idx][(*iter2).first]=denoPtr[(*iter2).first];
00520               _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[idx];
00521             }
00522         deno->decrRef();
00523         denoR->decrRef();
00524         break;
00525       }
00526     case IntegralGlobConstraint:
00527       {
00528         computeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
00529         break;
00530       }
00531     case RevIntegral:
00532       {
00533         MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),true);
00534         MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),true);
00535         const double *denoPtr=deno->getArray()->getConstPointer();
00536         const double *denoRPtr=denoR->getArray()->getConstPointer();
00537         if(trgField->getMesh()->getMeshDimension()==-1)
00538           {
00539             double *denoRPtr2=denoR->getArray()->getPointer();
00540             denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
00541           }
00542         if(srcField->getMesh()->getMeshDimension()==-1)
00543           {
00544             double *denoPtr2=deno->getArray()->getPointer();
00545             denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
00546           }
00547         int idx=0;
00548         for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
00549           for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
00550             {
00551               _deno_multiply[idx][(*iter2).first]=denoPtr[idx];
00552               _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[(*iter2).first];
00553             }
00554         deno->decrRef();
00555         denoR->decrRef();
00556         break;
00557       }
00558     case NoNature:
00559       throw INTERP_KERNEL::Exception("No nature specified ! Select one !");
00560     }
00561 }
00562 
00563 void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)
00564 {
00565   int idx=0;
00566   double *tmp=new double[inputNbOfCompo];
00567   for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
00568     {
00569       if((*iter1).empty())
00570         {
00571           if(isDftVal)
00572             std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
00573           continue;
00574         }
00575       else
00576         std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,0.);
00577       std::map<int,double>::const_iterator iter3=_deno_multiply[idx].begin();
00578       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter3++)
00579         {
00580           std::transform(inputPointer+(*iter2).first*inputNbOfCompo,inputPointer+((*iter2).first+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/(*iter3).second));
00581           std::transform(tmp,tmp+inputNbOfCompo,resPointer+idx*inputNbOfCompo,resPointer+idx*inputNbOfCompo,std::plus<double>());
00582         }
00583     }
00584   delete [] tmp;
00585 }
00586 
00587 void MEDCouplingRemapper::computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer)
00588 {
00589   std::vector<bool> isReached(_deno_reverse_multiply.size(),false);
00590   int idx=0;
00591   double *tmp=new double[inputNbOfCompo];
00592   std::fill(resPointer,resPointer+inputNbOfCompo*_deno_reverse_multiply.size(),0.);
00593   for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
00594     {
00595       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
00596         {
00597           isReached[(*iter2).first]=true;
00598           std::transform(inputPointer+idx*inputNbOfCompo,inputPointer+(idx+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/_deno_reverse_multiply[(*iter2).first][idx]));
00599           std::transform(tmp,tmp+inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,std::plus<double>());
00600         }
00601     }
00602   delete [] tmp;
00603   idx=0;
00604   for(std::vector<bool>::const_iterator iter3=isReached.begin();iter3!=isReached.end();iter3++,idx++)
00605     if(!*iter3)
00606       std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
00607 }
00608 
00609 void MEDCouplingRemapper::reverseMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
00610 {
00611   matOut.resize(nbColsMatIn);
00612   int id=0;
00613   for(std::vector<std::map<int,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
00614     for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
00615       matOut[(*iter2).first][id]=(*iter2).second;
00616 }
00617 
00618 void MEDCouplingRemapper::computeRowSumAndColSum(const std::vector<std::map<int,double> >& matrixDeno,
00619                                                  std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
00620 {
00621   std::map<int,double> values;
00622   int idx=0;
00623   for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
00624     {
00625       double sum=0.;
00626       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
00627         {
00628           sum+=(*iter2).second;
00629           values[(*iter2).first]+=(*iter2).second;
00630         }
00631       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
00632         deno[idx][(*iter2).first]=sum;
00633     }
00634   idx=0;
00635   for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
00636     {
00637       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
00638         denoReverse[(*iter2).first][idx]=values[(*iter2).first];
00639     }
00640 }
00641 
00642 void MEDCouplingRemapper::computeColSumAndRowSum(const std::vector<std::map<int,double> >& matrixDeno,
00643                                                  std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
00644 {
00645   std::map<int,double> values;
00646   int idx=0;
00647   for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
00648     {
00649       double sum=0.;
00650       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
00651         {
00652           sum+=(*iter2).second;
00653           values[(*iter2).first]+=(*iter2).second;
00654         }
00655       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
00656         denoReverse[(*iter2).first][idx]=sum;
00657     }
00658   idx=0;
00659   for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
00660     {
00661       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
00662         deno[idx][(*iter2).first]=values[(*iter2).first];
00663     }
00664 }
00665 
00666 void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<int,double> >& m1D,
00667                                                                      const std::vector< std::map<int,double> >& m2D,
00668                                                                      const int *corrCellIdSrc, int nbOf2DCellsSrc, int nbOf1DCellsSrc,
00669                                                                      const int *corrCellIdTrg)
00670 {
00671   int nbOf2DCellsTrg=m2D.size();
00672   int nbOf1DCellsTrg=m1D.size();
00673   int nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
00674   _matrix.resize(nbOf3DCellsTrg);
00675   int id2R=0;
00676   for(std::vector< std::map<int,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
00677     {
00678       for(std::map<int,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
00679         {
00680           int id1R=0;
00681           for(std::vector< std::map<int,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
00682             {
00683               for(std::map<int,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
00684                 {
00685                   _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
00686                 }
00687             }
00688         }
00689     }
00690 }
00691 
00692 void MEDCouplingRemapper::PrintMatrix(const std::vector<std::map<int,double> >& m)
00693 {
00694   int id=0;
00695   for(std::vector<std::map<int,double> >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++)
00696     {
00697       std::cout << "Target Cell # " << id << " : ";
00698       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
00699         std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), ";
00700       std::cout << std::endl;
00701     }
00702 }
00703 
00704 const std::vector<std::map<int,double> >& MEDCouplingRemapper::getCrudeMatrix() const
00705 {
00706   return _matrix;
00707 }