src/MEDCoupling/MEDCouplingRemapper.cxx

   1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (CEA/DEN)
  20
  21 #include "MEDCouplingRemapper.hxx"
  22 #include "MEDCouplingMemArray.hxx"
  23 #include "MEDCouplingFieldDouble.hxx"
  24 #include "MEDCouplingFieldTemplate.hxx"
  25 #include "MEDCouplingFieldDiscretization.hxx"
  26 #include "MEDCouplingExtrudedMesh.hxx"
  27 #include "MEDCouplingCMesh.hxx"
  28 #include "MEDCouplingNormalizedUnstructuredMesh.txx"
  29 #include "MEDCouplingNormalizedCartesianMesh.txx"
  30
  31 #include "Interpolation1D.txx"
  32 #include "Interpolation2DCurve.hxx"
  33 #include "Interpolation2D.txx"
  34 #include "Interpolation3D.txx"
  35 #include "Interpolation3DSurf.hxx"
  36 #include "Interpolation2D1D.txx"
  37 #include "Interpolation3D2D.txx"
  38 #include "InterpolationCU.txx"
  39 #include "InterpolationCC.txx"
  40
  41 using namespace ParaMEDMEM;
  42
  43 MEDCouplingRemapper::MEDCouplingRemapper():_src_ft(0),_target_ft(0),_interp_matrix_pol(IK_ONLY_PREFERED),_nature_of_deno(NoNature),_time_deno_update(0)
  44 {
  45 }
  46
  47 MEDCouplingRemapper::~MEDCouplingRemapper()
  48 {
  49   releaseData(false);
  50 }
  51
  52 int MEDCouplingRemapper::prepare(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const char *method) throw(INTERP_KERNEL::Exception)
  53 {
  54   if(!srcMesh || !targetMesh)
  55     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepare : presence of NULL input pointer !");
  56   std::string srcMethod,targetMethod;
  57   INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,srcMethod,targetMethod);
  58   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldTemplate> src=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(srcMethod.c_str()));
  59   src->setMesh(srcMesh);
  60   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldTemplate> target=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(targetMethod.c_str()));
  61   target->setMesh(targetMesh);
  62   return prepareEx(src,target);
  63 }
  64
  65 int MEDCouplingRemapper::prepareEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target) throw(INTERP_KERNEL::Exception)
  66 {
  67   if(!src || !target)
  68     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareEx : presence of NULL input pointer !");
  69   if(!src->getMesh() || !target->getMesh())
  70     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareEx : presence of NULL mesh pointer in given field template !");
  71   releaseData(true);
  72   _src_ft=const_cast<MEDCouplingFieldTemplate *>(src); _src_ft->incrRef();
  73   _target_ft=const_cast<MEDCouplingFieldTemplate *>(target); _target_ft->incrRef();
  74   if(isInterpKernelOnlyOrNotOnly())
  75     return prepareInterpKernelOnly();
  76   else
  77     return prepareNotInterpKernelOnly();
  78 }
  79
  80 int MEDCouplingRemapper::prepareInterpKernelOnly() throw(INTERP_KERNEL::Exception)
  81 {
  82   int meshInterpType=((int)_src_ft->getMesh()->getType()*16)+(int)_target_ft->getMesh()->getType();
  83   switch(meshInterpType)
  84     {
  85     case 85://Unstructured-Unstructured
  86       return prepareInterpKernelOnlyUU();
  87     case 87://Unstructured-Cartesian
  88       return prepareInterpKernelOnlyUC();
  89     case 117://Cartesian-Unstructured
  90       return prepareInterpKernelOnlyCU();
  91     case 119://Cartesian-Cartesian
  92       return prepareInterpKernelOnlyCC();
  93     case 136://Extruded-Extruded
  94       return prepareInterpKernelOnlyEE();
  95     default:
  96       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnly : Not managed type of meshes ! Dealt meshes type are : Unstructured<->Unstructured, Unstructured<->Cartesian, Cartesian<->Cartesian, Extruded<->Extruded !");
  97     }
  98 }
  99
 100 int MEDCouplingRemapper::prepareNotInterpKernelOnly() throw(INTERP_KERNEL::Exception)
 101 {
 102   std::string srcm,trgm,method;
 103   method=checkAndGiveInterpolationMethodStr(srcm,trgm);
 104   switch(CheckInterpolationMethodManageableByNotOnlyInterpKernel(method))
 105     {
 106     case 0:
 107       return prepareNotInterpKernelOnlyGaussGauss();
 108     default:
 109       {
 110         std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnly : INTERNAL ERROR ! the method \"" << method << "\" declared as managed bu not implemented !";
 111         throw INTERP_KERNEL::Exception(oss.str().c_str());
 112       }
 113     }
 114 }
 115
 116 /*!
 117  * This method performs the operation source to target using matrix computed in ParaMEDMEM::MEDCouplingRemapper::prepare method.
 118  * 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.
 119  *
 120  * \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.
 121  * \param [out] targetField the destination field with the allocated array in which all tuples will be overwritten.
 122  */
 123 void MEDCouplingRemapper::transfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, double dftValue) throw(INTERP_KERNEL::Exception)
 124 {
 125   transferUnderground(srcField,targetField,true,dftValue);
 126 }
 127
 128 /*!
 129  * This method is equivalent to ParaMEDMEM::MEDCouplingRemapper::transfer except that here \b targetField is a in/out parameter.
 130  * 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
 131  * let unchanged.
 132  * This method requires that \b targetField was fully defined and allocated. If the array is not allocated an exception will be thrown.
 133  *
 134  * \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.
 135  * \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.
 136  */
 137 void MEDCouplingRemapper::partialTransfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField) throw(INTERP_KERNEL::Exception)
 138 {
 139   transferUnderground(srcField,targetField,false,std::numeric_limits<double>::max());
 140 }
 141
 142 void MEDCouplingRemapper::reverseTransfer(MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *targetField, double dftValue) throw(INTERP_KERNEL::Exception)
 143 {
 144   checkPrepare();
 145   if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
 146     throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
 147   if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
 148     throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
 149   if(srcField->getNature()!=targetField->getNature())
 150     throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
 151   DataArrayDouble *array=srcField->getArray();
 152   int trgNbOfCompo=targetField->getNumberOfComponents();
 153   if(array)
 154     {
 155       if(trgNbOfCompo!=srcField->getNumberOfComponents())
 156         throw INTERP_KERNEL::Exception("Number of components mismatch !");
 157     }
 158   else
 159     {
 160       array=DataArrayDouble::New();
 161       array->alloc(srcField->getNumberOfTuples(),trgNbOfCompo);
 162       srcField->setArray(array);
 163       array->decrRef();
 164     }
 165   computeDeno(srcField->getNature(),srcField,targetField);
 166   double *resPointer=array->getPointer();
 167   const double *inputPointer=targetField->getArray()->getConstPointer();
 168   computeReverseProduct(inputPointer,trgNbOfCompo,dftValue,resPointer);
 169 }
 170
 171 MEDCouplingFieldDouble *MEDCouplingRemapper::transferField(const MEDCouplingFieldDouble *srcField, double dftValue) throw(INTERP_KERNEL::Exception)
 172 {
 173   checkPrepare();
 174   if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
 175     throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
 176   MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(*_target_ft,srcField->getTimeDiscretization());
 177   ret->setNature(srcField->getNature());
 178   transfer(srcField,ret,dftValue);
 179   ret->copyAllTinyAttrFrom(srcField);//perform copy of tiny strings after and not before transfer because the array will be created on transfer
 180   return ret;
 181 }
 182
 183 MEDCouplingFieldDouble *MEDCouplingRemapper::reverseTransferField(const MEDCouplingFieldDouble *targetField, double dftValue) throw(INTERP_KERNEL::Exception)
 184 {
 185   checkPrepare();
 186   if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
 187     throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
 188   MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(*_src_ft,targetField->getTimeDiscretization());
 189   ret->setNature(targetField->getNature());
 190   reverseTransfer(ret,targetField,dftValue);
 191   ret->copyAllTinyAttrFrom(targetField);//perform copy of tiny strings after and not before reverseTransfer because the array will be created on reverseTransfer
 192   return ret;
 193 }
 194
 195 /*!
 196  * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
 197  * is here only for automatic CORBA generators.
 198  */
 199 bool MEDCouplingRemapper::setOptionInt(const std::string& key, int value)
 200 {
 201   return INTERP_KERNEL::InterpolationOptions::setOptionInt(key,value);
 202 }
 203
 204 /*!
 205  * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
 206  * is here only for automatic CORBA generators.
 207  */
 208 bool MEDCouplingRemapper::setOptionDouble(const std::string& key, double value)
 209 {
 210   return INTERP_KERNEL::InterpolationOptions::setOptionDouble(key,value);
 211 }
 212
 213 /*!
 214  * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
 215  * is here only for automatic CORBA generators.
 216  */
 217 bool MEDCouplingRemapper::setOptionString(const std::string& key, const std::string& value)
 218 {
 219   return INTERP_KERNEL::InterpolationOptions::setOptionString(key,value);
 220 }
 221
 222 /*!
 223  * This method returns the interpolation matrix policy. This policy specifies which interpolation matrix method to keep or prefered.
 224  * If interpolation matrix policy is :
 225  *
 226  * - set to IK_ONLY_PREFERED (0) (the default) : the INTERP_KERNEL only method is prefered. That is to say, if it is possible to treat the case
 227  *   regarding spatial discretization of source and target with INTERP_KERNEL only method, INTERP_KERNEL only method will be performed.
 228  *   If not, the \b not only INTERP_KERNEL method will be attempt.
 229  *
 230  * - set to NOT_IK_ONLY_PREFERED (1) : the \b NOT only INTERP_KERNEL method is prefered. That is to say, if it is possible to treat the case
 231  *   regarding spatial discretization of source and target with \b NOT only INTERP_KERNEL method, \b NOT only INTERP_KERNEL method, will be performed.
 232  *   If not, the INTERP_KERNEL only method will be attempt.
 233  *
 234  * - IK_ONLY_FORCED (2) : Only INTERP_KERNEL only method will be launched.
 235  *
 236  * - NOT_IK_ONLY_FORCED (3) : Only \b NOT INTERP_KERNEL only method will be launched.
 237  *
 238  * \sa MEDCouplingRemapper::setInterpolationMatrixPolicy
 239  */
 240 int MEDCouplingRemapper::getInterpolationMatrixPolicy() const
 241 {
 242   return _interp_matrix_pol;
 243 }
 244
 245 /*!
 246  * This method sets a new interpolation matrix policy. The default one is IK_PREFERED (0). The input is of type \c int to be dealt by standard Salome
 247  * CORBA component generators. This method throws an INTERP_KERNEL::Exception if a the input integer is not in the available possibilities, that is to say not in
 248  * [0 (IK_PREFERED) , 1 (NOT_IK_PREFERED), 2 (IK_ONLY_FORCED), 3 (NOT_IK_ONLY_FORCED)].
 249  *
 250  * If interpolation matrix policy is :
 251  *
 252  * - set to IK_ONLY_PREFERED (0) (the default) : the INTERP_KERNEL only method is prefered. That is to say, if it is possible to treat the case
 253  *   regarding spatial discretization of source and target with INTERP_KERNEL only method, INTERP_KERNEL only method will be performed.
 254  *   If not, the \b not only INTERP_KERNEL method will be attempt.
 255  *
 256  * - set to NOT_IK_ONLY_PREFERED (1) : the \b NOT only INTERP_KERNEL method is prefered. That is to say, if it is possible to treat the case
 257  *   regarding spatial discretization of source and target with \b NOT only INTERP_KERNEL method, \b NOT only INTERP_KERNEL method, will be performed.
 258  *   If not, the INTERP_KERNEL only method will be attempt.
 259  *
 260  * - IK_ONLY_FORCED (2) : Only INTERP_KERNEL only method will be launched.
 261  *
 262  * - NOT_IK_ONLY_FORCED (3) : Only \b NOT INTERP_KERNEL only method will be launched.
 263  *
 264  * \input newInterpMatPol the new interpolation matrix method policy. This parameter is of type \c int and not of type \c ParaMEDMEM::InterpolationMatrixPolicy
 265  *                        for automatic generation of CORBA component.
 266  *
 267  * \sa MEDCouplingRemapper::getInterpolationMatrixPolicy
 268  */
 269 void MEDCouplingRemapper::setInterpolationMatrixPolicy(int newInterpMatPol) throw(INTERP_KERNEL::Exception)
 270 {
 271   switch(newInterpMatPol)
 272     {
 273     case 0:
 274       _interp_matrix_pol=IK_ONLY_PREFERED;
 275       break;
 276     case 1:
 277       _interp_matrix_pol=NOT_IK_ONLY_PREFERED;
 278       break;
 279     case 2:
 280       _interp_matrix_pol=IK_ONLY_FORCED;
 281       break;
 282     case 3:
 283       _interp_matrix_pol=NOT_IK_ONLY_FORCED;
 284       break;
 285     default:
 286       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::setInterpolationMatrixPolicy : invalid input integer value ! Should be in [0 (IK_PREFERED) , 1 (NOT_IK_PREFERED), 2 (IK_ONLY_FORCED), 3 (NOT_IK_ONLY_FORCED)] ! For information, the default is IK_PREFERED=0 !");
 287     }
 288 }
 289
 290 int MEDCouplingRemapper::prepareInterpKernelOnlyUU() throw(INTERP_KERNEL::Exception)
 291 {
 292   const MEDCouplingUMesh *src_mesh=static_cast<const MEDCouplingUMesh *>(_src_ft->getMesh());
 293   const MEDCouplingUMesh *target_mesh=static_cast<const MEDCouplingUMesh *>(_target_ft->getMesh());
 294   std::string srcMeth,trgMeth;
 295   std::string method=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
 296   const int srcMeshDim=src_mesh->getMeshDimension();
 297   int srcSpaceDim=-1;
 298   if(srcMeshDim!=-1)
 299     srcSpaceDim=src_mesh->getSpaceDimension();
 300   const int trgMeshDim=target_mesh->getMeshDimension();
 301   int trgSpaceDim=-1;
 302   if(trgMeshDim!=-1)
 303     trgSpaceDim=target_mesh->getSpaceDimension();
 304   if(trgSpaceDim!=srcSpaceDim)
 305     if(trgSpaceDim!=-1 && srcSpaceDim!=-1)
 306       throw INTERP_KERNEL::Exception("Incoherent space dimension detected between target and source.");
 307   int nbCols;
 308   if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==1)
 309     {
 310       MEDCouplingNormalizedUnstructuredMesh<1,1> source_mesh_wrapper(src_mesh);
 311       MEDCouplingNormalizedUnstructuredMesh<1,1> target_mesh_wrapper(target_mesh);
 312       INTERP_KERNEL::Interpolation1D interpolation(*this);
 313       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method.c_str());
 314     }
 315   else if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==2)
 316     {
 317       MEDCouplingNormalizedUnstructuredMesh<2,1> source_mesh_wrapper(src_mesh);
 318       MEDCouplingNormalizedUnstructuredMesh<2,1> target_mesh_wrapper(target_mesh);
 319       INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
 320       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method.c_str());
 321     }
 322   else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==2)
 323     {
 324       MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 325       MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
 326       INTERP_KERNEL::Interpolation2D interpolation(*this);
 327       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method.c_str());
 328     }
 329   else if(srcMeshDim==3 && trgMeshDim==3 && srcSpaceDim==3)
 330     {
 331       MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 332       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
 333       INTERP_KERNEL::Interpolation3D interpolation(*this);
 334       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method.c_str());
 335     }
 336   else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==3)
 337     {
 338       MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
 339       MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh);
 340       INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
 341       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method.c_str());
 342     }
 343   else if(srcMeshDim==3 && trgMeshDim==1 && srcSpaceDim==3)
 344     {
 345       if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
 346         throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
 347       MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 348       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
 349       INTERP_KERNEL::Interpolation3D interpolation(*this);
 350       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method.c_str());
 351     }
 352   else if(srcMeshDim==1 && trgMeshDim==3 && srcSpaceDim==3)
 353     {
 354       if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
 355         throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
 356       MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 357       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
 358       INTERP_KERNEL::Interpolation3D interpolation(*this);
 359       std::vector<std::map<int,double> > matrixTmp;
 360       nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method.c_str());
 361       ReverseMatrix(matrixTmp,nbCols,_matrix);
 362       nbCols=matrixTmp.size();
 363     }
 364   else if(srcMeshDim==2 && trgMeshDim==1 && srcSpaceDim==2)
 365     {
 366       if(getIntersectionType()==INTERP_KERNEL::PointLocator)
 367         {
 368           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 369           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
 370           INTERP_KERNEL::Interpolation2D interpolation(*this);
 371           nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method.c_str());
 372         }
 373       else
 374         {
 375           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 376           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
 377           INTERP_KERNEL::Interpolation2D1D interpolation(*this);
 378           std::vector<std::map<int,double> > matrixTmp;
 379           nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method.c_str());
 380           ReverseMatrix(matrixTmp,nbCols,_matrix);
 381           nbCols=matrixTmp.size();
 382           INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
 383           if(!duplicateFaces.empty())
 384             {
 385               std::ostringstream oss; oss << "An unexpected situation happend ! For the following 1D Cells are part of edges shared by 2D cells :\n";
 386               for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
 387                 {
 388                   oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
 389                   std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
 390                   oss << std::endl;
 391                 }
 392             }
 393         }
 394     }
 395   else if(srcMeshDim==1 && trgMeshDim==2 && srcSpaceDim==2)
 396     {
 397       if(getIntersectionType()==INTERP_KERNEL::PointLocator)
 398         {
 399           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 400           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
 401           INTERP_KERNEL::Interpolation2D interpolation(*this);
 402           std::vector<std::map<int,double> > matrixTmp;
 403           nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method.c_str());
 404           ReverseMatrix(matrixTmp,nbCols,_matrix);
 405           nbCols=matrixTmp.size();
 406         }
 407       else
 408         {
 409           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 410           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
 411           INTERP_KERNEL::Interpolation2D1D interpolation(*this);
 412           nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method.c_str());
 413           INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
 414           if(!duplicateFaces.empty())
 415             {
 416               std::ostringstream oss; oss << "An unexpected situation happend ! For the following 1D Cells are part of edges shared by 2D cells :\n";
 417               for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
 418                 {
 419                   oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
 420                   std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
 421                   oss << std::endl;
 422                 }
 423             }
 424         }
 425     }
 426   else if(srcMeshDim==2 && trgMeshDim==3 && srcSpaceDim==3)
 427     {
 428       MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 429       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
 430       INTERP_KERNEL::Interpolation3D2D interpolation(*this);
 431       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method.c_str());
 432       INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
 433       if(!duplicateFaces.empty())
 434         {
 435           std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
 436           for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
 437             {
 438               oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
 439               std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
 440               oss << std::endl;
 441             }
 442         }
 443     }
 444   else if(srcMeshDim==3 && trgMeshDim==2 && srcSpaceDim==3)
 445     {
 446       MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 447       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
 448       INTERP_KERNEL::Interpolation3D2D interpolation(*this);
 449       std::vector<std::map<int,double> > matrixTmp;
 450       nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method.c_str());
 451       ReverseMatrix(matrixTmp,nbCols,_matrix);
 452       nbCols=matrixTmp.size();
 453       INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
 454       if(!duplicateFaces.empty())
 455         {
 456           std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
 457           for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
 458             {
 459               oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
 460               std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
 461               oss << std::endl;
 462             }
 463         }
 464     }
 465   else if(trgMeshDim==-1)
 466     {
 467       if(srcMeshDim==2 && srcSpaceDim==2)
 468         {
 469           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 470           INTERP_KERNEL::Interpolation2D interpolation(*this);
 471           nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth.c_str());
 472         }
 473       else if(srcMeshDim==3 && srcSpaceDim==3)
 474         {
 475           MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 476           INTERP_KERNEL::Interpolation3D interpolation(*this);
 477           nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth.c_str());
 478         }
 479       else if(srcMeshDim==2 && srcSpaceDim==3)
 480         {
 481           MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
 482           INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
 483           nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth.c_str());
 484         }
 485       else
 486         throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
 487     }
 488   else if(srcMeshDim==-1)
 489     {
 490       if(trgMeshDim==2 && trgSpaceDim==2)
 491         {
 492           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(target_mesh);
 493           INTERP_KERNEL::Interpolation2D interpolation(*this);
 494           nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth.c_str());
 495         }
 496       else if(trgMeshDim==3 && trgSpaceDim==3)
 497         {
 498           MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(target_mesh);
 499           INTERP_KERNEL::Interpolation3D interpolation(*this);
 500           nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth.c_str());
 501         }
 502       else if(trgMeshDim==2 && trgSpaceDim==3)
 503         {
 504           MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(target_mesh);
 505           INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
 506           nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth.c_str());
 507         }
 508       else
 509         throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh from -1D sourceMesh");
 510     }
 511   else
 512     throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension");
 513   _deno_multiply.clear();
 514   _deno_multiply.resize(_matrix.size());
 515   _deno_reverse_multiply.clear();
 516   _deno_reverse_multiply.resize(nbCols);
 517   declareAsNew();
 518   return 1;
 519 }
 520
 521 int MEDCouplingRemapper::prepareInterpKernelOnlyEE() throw(INTERP_KERNEL::Exception)
 522 {
 523   std::string srcMeth,trgMeth;
 524   std::string methC=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
 525   const MEDCouplingExtrudedMesh *src_mesh=static_cast<const MEDCouplingExtrudedMesh *>(_src_ft->getMesh());
 526   const MEDCouplingExtrudedMesh *target_mesh=static_cast<const MEDCouplingExtrudedMesh *>(_target_ft->getMesh());
 527   if(methC!="P0P0")
 528     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyEE : Only P0P0 method implemented for Extruded/Extruded meshes !");
 529   MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh->getMesh2D());
 530   MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh->getMesh2D());
 531   INTERP_KERNEL::Interpolation3DSurf interpolation2D(*this);
 532   std::vector<std::map<int,double> > matrix2D;
 533   int nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,methC.c_str());
 534   MEDCouplingUMesh *s1D,*t1D;
 535   double v[3];
 536   MEDCouplingExtrudedMesh::Project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v);
 537   MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D);
 538   MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
 539   std::vector<std::map<int,double> > matrix1D;
 540   INTERP_KERNEL::Interpolation1D interpolation1D(*this);
 541   int nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,methC.c_str());
 542   s1D->decrRef();
 543   t1D->decrRef();
 544   buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
 545                                              target_mesh->getMesh3DIds()->getConstPointer());
 546   //
 547   _deno_multiply.clear();
 548   _deno_multiply.resize(_matrix.size());
 549   _deno_reverse_multiply.clear();
 550   _deno_reverse_multiply.resize(nbCols2D*nbCols1D);
 551   declareAsNew();
 552   return 1;
 553 }
 554
 555 int MEDCouplingRemapper::prepareInterpKernelOnlyUC() throw(INTERP_KERNEL::Exception)
 556 {
 557   std::string srcMeth,trgMeth;
 558   std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
 559   if(methodCpp!="P0P0")
 560     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only P0P0 interpolation supported for the moment !");
 561   const MEDCouplingUMesh *src_mesh=static_cast<const MEDCouplingUMesh *>(_src_ft->getMesh());
 562   const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
 563   const int srcMeshDim=src_mesh->getMeshDimension();
 564   const int srcSpceDim=src_mesh->getSpaceDimension();
 565   const int trgMeshDim=target_mesh->getMeshDimension();
 566   if(srcMeshDim!=srcSpceDim || srcMeshDim!=trgMeshDim)
 567     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : space dim of src unstructured should be equal to mesh dim of src unstructured and should be equal also equal to trg cartesian dimension !");
 568   std::vector<std::map<int,double> > res;
 569   switch(srcMeshDim)
 570     {
 571     case 1:
 572       {
 573         MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
 574         MEDCouplingNormalizedUnstructuredMesh<1,1> sourceWrapper(src_mesh);
 575         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 576         myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
 577         break;
 578       }
 579     case 2:
 580       {
 581         MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
 582         MEDCouplingNormalizedUnstructuredMesh<2,2> sourceWrapper(src_mesh);
 583         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 584         myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
 585         break;
 586       }
 587     case 3:
 588       {
 589         MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
 590         MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(src_mesh);
 591         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 592         myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
 593         break;
 594       }
 595     default:
 596       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only dimension 1 2 or 3 supported !");
 597     }
 598   ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix);
 599   nullifiedTinyCoeffInCrudeMatrixAbs(0.);
 600   //
 601   _deno_multiply.clear();
 602   _deno_multiply.resize(_matrix.size());
 603   _deno_reverse_multiply.clear();
 604   _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
 605   declareAsNew();
 606   return 1;
 607 }
 608
 609 int MEDCouplingRemapper::prepareInterpKernelOnlyCU() throw(INTERP_KERNEL::Exception)
 610 {
 611   std::string srcMeth,trgMeth;
 612   std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
 613   if(methodCpp!="P0P0")
 614     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only P0P0 interpolation supported for the moment !");
 615   const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
 616   const MEDCouplingUMesh *target_mesh=static_cast<const MEDCouplingUMesh *>(_target_ft->getMesh());
 617   const int srcMeshDim=src_mesh->getMeshDimension();
 618   const int trgMeshDim=target_mesh->getMeshDimension();
 619   const int trgSpceDim=target_mesh->getSpaceDimension();
 620   if(trgMeshDim!=trgSpceDim || trgMeshDim!=srcMeshDim)
 621     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : space dim of target unstructured should be equal to mesh dim of target unstructured and should be equal also equal to source cartesian dimension !");
 622   switch(srcMeshDim)
 623     {
 624     case 1:
 625       {
 626         MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
 627         MEDCouplingNormalizedUnstructuredMesh<1,1> targetWrapper(target_mesh);
 628         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 629         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 630         break;
 631       }
 632     case 2:
 633       {
 634         MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
 635         MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(target_mesh);
 636         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 637         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 638         break;
 639       }
 640     case 3:
 641       {
 642         MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
 643         MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(target_mesh);
 644         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 645         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 646         break;
 647       }
 648     default:
 649       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only dimension 1 2 or 3 supported !");
 650     }
 651   nullifiedTinyCoeffInCrudeMatrixAbs(0.);
 652   //
 653   _deno_multiply.clear();
 654   _deno_multiply.resize(_matrix.size());
 655   _deno_reverse_multiply.clear();
 656   _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
 657   declareAsNew();
 658   return 1;
 659 }
 660
 661 int MEDCouplingRemapper::prepareInterpKernelOnlyCC() throw(INTERP_KERNEL::Exception)
 662 {
 663   std::string srcMeth,trgMeth;
 664   std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
 665   if(methodCpp!="P0P0")
 666     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only P0P0 interpolation supported for the moment !");
 667   const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
 668   const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
 669   const int srcMeshDim=src_mesh->getMeshDimension();
 670   const int trgMeshDim=target_mesh->getMeshDimension();
 671   if(trgMeshDim!=srcMeshDim)
 672     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : dim of target cartesian should be equal to dim of source cartesian dimension !");
 673   switch(srcMeshDim)
 674     {
 675     case 1:
 676       {
 677         MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
 678         MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
 679         INTERP_KERNEL::InterpolationCC myInterpolator(*this);
 680         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 681         break;
 682       }
 683     case 2:
 684       {
 685         MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
 686         MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
 687         INTERP_KERNEL::InterpolationCC myInterpolator(*this);
 688         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 689         break;
 690       }
 691     case 3:
 692       {
 693         MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
 694         MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
 695         INTERP_KERNEL::InterpolationCC myInterpolator(*this);
 696         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 697         break;
 698       }
 699     default:
 700       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only dimension 1 2 or 3 supported !");
 701     }
 702   nullifiedTinyCoeffInCrudeMatrixAbs(0.);
 703   //
 704   _deno_multiply.clear();
 705   _deno_multiply.resize(_matrix.size());
 706   _deno_reverse_multiply.clear();
 707   _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
 708   declareAsNew();
 709   return 1;
 710 }
 711
 712 int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss() throw(INTERP_KERNEL::Exception)
 713 {
 714   if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
 715       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss : The intersection type is not supported ! Only PointLocator is supported for Gauss->Gauss interpolation ! Please invoke setIntersectionType(PointLocator) on the MEDCouplingRemapper instance !");
 716   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> trgLoc=_target_ft->getLocalizationOfDiscr();
 717   const double *trgLocPtr=trgLoc->begin();
 718   int trgSpaceDim=trgLoc->getNumberOfComponents();
 719   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> srcOffsetArr=_src_ft->getDiscretization()->getOffsetArr(_src_ft->getMesh());
 720   if(trgSpaceDim!=_src_ft->getMesh()->getSpaceDimension())
 721     {
 722       std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss : space dimensions mismatch between source and target !";
 723       oss << " Target discretization localization has dimension " << trgSpaceDim << ", whereas the space dimension of source is equal to ";
 724       oss << _src_ft->getMesh()->getSpaceDimension() << " !";
 725       throw INTERP_KERNEL::Exception(oss.str().c_str());
 726     }
 727   const int *srcOffsetArrPtr=srcOffsetArr->begin();
 728   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> srcLoc=_src_ft->getLocalizationOfDiscr();
 729   const double *srcLocPtr=srcLoc->begin();
 730   std::vector<int> elts,eltsIndex;
 731   int trgNbOfGaussPts=trgLoc->getNumberOfTuples();
 732   _matrix.resize(trgNbOfGaussPts);
 733   _src_ft->getMesh()->getCellsContainingPoints(trgLoc->begin(),trgNbOfGaussPts,getPrecision(),elts,eltsIndex);
 734   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> eltsIndex2=DataArrayInt::New(); eltsIndex2->useArray(&eltsIndex[0],false,CPP_DEALLOC,(int)eltsIndex.size(),1);
 735   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nbOfSrcCellsShTrgPts=eltsIndex2->deltaShiftIndex();
 736   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ids0=nbOfSrcCellsShTrgPts->getIdsNotEqual(0);
 737   for(const int *trgId=ids0->begin();trgId!=ids0->end();trgId++)
 738     {
 739       const double *ptTrg=trgLocPtr+trgSpaceDim*(*trgId);
 740       int srcCellId=elts[eltsIndex[*trgId]];
 741       double dist=std::numeric_limits<double>::max();
 742       int srcEntry=-1;
 743       for(int srcId=srcOffsetArrPtr[srcCellId];srcId<srcOffsetArrPtr[srcCellId+1];srcId++)
 744         {
 745           const double *ptSrc=srcLocPtr+trgSpaceDim*srcId;
 746           double tmp=0.;
 747           for(int i=0;i<trgSpaceDim;i++)
 748             tmp+=(ptTrg[i]-ptSrc[i])*(ptTrg[i]-ptSrc[i]);
 749           if(tmp<dist)
 750             { dist=tmp; srcEntry=srcId; }
 751         }
 752       _matrix[*trgId][srcEntry]=1.;
 753     }
 754   if(ids0->getNumberOfTuples()!=trgNbOfGaussPts)
 755     {
 756       MEDCouplingAutoRefCountObjectPtr<DataArrayInt> orphanTrgIds=nbOfSrcCellsShTrgPts->getIdsEqual(0);
 757       MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> orphanTrg=trgLoc->selectByTupleId(orphanTrgIds->begin(),orphanTrgIds->end());
 758       MEDCouplingAutoRefCountObjectPtr<DataArrayInt> srcIdPerTrg=srcLoc->findClosestTupleId(orphanTrg);
 759       const int *srcIdPerTrgPtr=srcIdPerTrg->begin();
 760       for(const int *orphanTrgId=orphanTrgIds->begin();orphanTrgId!=orphanTrgIds->end();orphanTrgId++,srcIdPerTrgPtr++)
 761         _matrix[*orphanTrgId][*srcIdPerTrgPtr]=2.;
 762     }
 763   _deno_multiply.clear();
 764   _deno_multiply.resize(_matrix.size());
 765   _deno_reverse_multiply.clear();
 766   _deno_reverse_multiply.resize(srcLoc->getNumberOfTuples());
 767   declareAsNew();
 768   return 1;
 769 }
 770
 771 /*!
 772  * This method checks that the input interpolation \a method is managed by not INTERP_KERNEL only methods.
 773  * If no an INTERP_KERNEL::Exception will be thrown. If yes, a magic number will be returned to switch in the MEDCouplingRemapper::prepareNotInterpKernelOnly method.
 774  */
 775 int MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel(const std::string& method) throw(INTERP_KERNEL::Exception)
 776 {
 777   if(method=="GAUSSGAUSS")
 778     return 0;
 779   std::ostringstream oss; oss << "MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel : ";
 780   oss << "The method \"" << method << "\" is not manageable by not INTERP_KERNEL only method.";
 781   oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS !";
 782   throw INTERP_KERNEL::Exception(oss.str().c_str());
 783 }
 784
 785 /*!
 786  * This method determines regarding \c _interp_matrix_pol attribute ( set by MEDCouplingRemapper::setInterpolationMatrixPolicy and by default equal
 787  * to IK_ONLY_PREFERED = 0 ) , which method will be applied. If \c true is returned the INTERP_KERNEL only method should be applied to \c false the \b not
 788  * only INTERP_KERNEL method should be applied.
 789  */
 790 bool MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly() const throw(INTERP_KERNEL::Exception)
 791 {
 792   std::string srcm,trgm,method;
 793   method=checkAndGiveInterpolationMethodStr(srcm,trgm);
 794   switch(_interp_matrix_pol)
 795     {
 796     case IK_ONLY_PREFERED:
 797       {
 798         try
 799           {
 800             std::string tmp1,tmp2;
 801             INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method.c_str(),tmp1,tmp2);
 802             return true;
 803           }
 804         catch(INTERP_KERNEL::Exception& e)
 805           {
 806             return false;
 807           }
 808       }
 809     case NOT_IK_ONLY_PREFERED:
 810       {
 811         try
 812           {
 813             CheckInterpolationMethodManageableByNotOnlyInterpKernel(method);
 814             return false;
 815           }
 816         catch(INTERP_KERNEL::Exception& e)
 817           {
 818             return true;
 819           }
 820       }
 821     case IK_ONLY_FORCED:
 822       return true;
 823     case NOT_IK_ONLY_FORCED:
 824       return false;
 825     default:
 826       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly : internal error ! The interpolation matrix policy is not managed ! Try to change it using MEDCouplingRemapper::setInterpolationMatrixPolicy !");
 827     }
 828 }
 829
 830 void MEDCouplingRemapper::updateTime() const
 831 {
 832 }
 833
 834 void MEDCouplingRemapper::checkPrepare() const throw(INTERP_KERNEL::Exception)
 835 {
 836   const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
 837   if(!s || !t)
 838     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that MEDCouplingRemapper::prepare(Ex) has not been called !");
 839   if(!s->getMesh() || !t->getMesh())
 840     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that no all field templates have their mesh set !");
 841 }
 842
 843 /*!
 844  * This method builds a code considering already set field discretization int \a this : \a _src_ft and \a _target_ft.
 845  * This method returns 3 informations (2 in ouput parameters and 1 in return).
 846  *
 847  * \param [out] srcMeth the string code of the discretization of source field template
 848  * \param [out] trgMeth the string code of the discretization of target field template
 849  * \return the standardized string code (compatible with INTERP_KERNEL) for matrix of numerators (in \a _matrix)
 850  */
 851 std::string MEDCouplingRemapper::checkAndGiveInterpolationMethodStr(std::string& srcMeth, std::string& trgMeth) const throw(INTERP_KERNEL::Exception)
 852 {
 853   const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
 854   if(!s || !t)
 855     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have been set !");
 856   if(!s->getMesh() || !t->getMesh())
 857     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have their mesh set !");
 858   srcMeth=_src_ft->getDiscretization()->getRepr();
 859   trgMeth=_target_ft->getDiscretization()->getRepr();
 860   std::string method(srcMeth); method+=trgMeth;
 861   return method;
 862 }
 863
 864 void MEDCouplingRemapper::releaseData(bool matrixSuppression)
 865 {
 866   _src_ft=0;
 867   _target_ft=0;
 868   if(matrixSuppression)
 869     {
 870       _matrix.clear();
 871       _deno_multiply.clear();
 872       _deno_reverse_multiply.clear();
 873     }
 874 }
 875
 876 void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue) throw(INTERP_KERNEL::Exception)
 877 {
 878   checkPrepare();
 879   if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
 880     throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
 881   if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
 882     throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
 883   if(srcField->getNature()!=targetField->getNature())
 884     throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
 885   DataArrayDouble *array=targetField->getArray();
 886   int srcNbOfCompo=srcField->getNumberOfComponents();
 887   if(array)
 888     {
 889       if(srcNbOfCompo!=targetField->getNumberOfComponents())
 890         throw INTERP_KERNEL::Exception("Number of components mismatch !");
 891     }
 892   else
 893     {
 894       if(!isDftVal)
 895         throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !");
 896       array=DataArrayDouble::New();
 897       array->alloc(targetField->getNumberOfTuples(),srcNbOfCompo);
 898       targetField->setArray(array);
 899       array->decrRef();
 900     }
 901   computeDeno(srcField->getNature(),srcField,targetField);
 902   double *resPointer=array->getPointer();
 903   const double *inputPointer=srcField->getArray()->getConstPointer();
 904   computeProduct(inputPointer,srcNbOfCompo,isDftVal,dftValue,resPointer);
 905 }
 906
 907 void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
 908 {
 909   if(nat==NoNature)
 910     return computeDenoFromScratch(nat,srcField,trgField);
 911   else if(nat!=_nature_of_deno)
 912    return computeDenoFromScratch(nat,srcField,trgField);
 913   else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis())
 914     return computeDenoFromScratch(nat,srcField,trgField);
 915 }
 916
 917 void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField) throw(INTERP_KERNEL::Exception)
 918 {
 919   _nature_of_deno=nat;
 920   _time_deno_update=getTimeOfThis();
 921   switch(_nature_of_deno)
 922     {
 923     case ConservativeVolumic:
 924       {
 925         ComputeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
 926         break;
 927       }
 928     case Integral:
 929       {
 930         MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
 931         MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
 932         const double *denoPtr=deno->getArray()->getConstPointer();
 933         const double *denoRPtr=denoR->getArray()->getConstPointer();
 934         if(trgField->getMesh()->getMeshDimension()==-1)
 935           {
 936             double *denoRPtr2=denoR->getArray()->getPointer();
 937             denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
 938           }
 939         if(srcField->getMesh()->getMeshDimension()==-1)
 940           {
 941             double *denoPtr2=deno->getArray()->getPointer();
 942             denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
 943           }
 944         int idx=0;
 945         for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
 946           for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 947             {
 948               _deno_multiply[idx][(*iter2).first]=denoPtr[(*iter2).first];
 949               _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[idx];
 950             }
 951         deno->decrRef();
 952         denoR->decrRef();
 953         break;
 954       }
 955     case IntegralGlobConstraint:
 956       {
 957         ComputeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
 958         break;
 959       }
 960     case RevIntegral:
 961       {
 962         MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
 963         MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
 964         const double *denoPtr=deno->getArray()->getConstPointer();
 965         const double *denoRPtr=denoR->getArray()->getConstPointer();
 966         if(trgField->getMesh()->getMeshDimension()==-1)
 967           {
 968             double *denoRPtr2=denoR->getArray()->getPointer();
 969             denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
 970           }
 971         if(srcField->getMesh()->getMeshDimension()==-1)
 972           {
 973             double *denoPtr2=deno->getArray()->getPointer();
 974             denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
 975           }
 976         int idx=0;
 977         for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
 978           for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 979             {
 980               _deno_multiply[idx][(*iter2).first]=denoPtr[idx];
 981               _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[(*iter2).first];
 982             }
 983         deno->decrRef();
 984         denoR->decrRef();
 985         break;
 986       }
 987     case NoNature:
 988       throw INTERP_KERNEL::Exception("No nature specified ! Select one !");
 989     }
 990 }
 991
 992 void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)
 993 {
 994   int idx=0;
 995   double *tmp=new double[inputNbOfCompo];
 996   for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
 997     {
 998       if((*iter1).empty())
 999         {
1000           if(isDftVal)
1001             std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1002           continue;
1003         }
1004       else
1005         std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,0.);
1006       std::map<int,double>::const_iterator iter3=_deno_multiply[idx].begin();
1007       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter3++)
1008         {
1009           std::transform(inputPointer+(*iter2).first*inputNbOfCompo,inputPointer+((*iter2).first+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/(*iter3).second));
1010           std::transform(tmp,tmp+inputNbOfCompo,resPointer+idx*inputNbOfCompo,resPointer+idx*inputNbOfCompo,std::plus<double>());
1011         }
1012     }
1013   delete [] tmp;
1014 }
1015
1016 void MEDCouplingRemapper::computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer)
1017 {
1018   std::vector<bool> isReached(_deno_reverse_multiply.size(),false);
1019   int idx=0;
1020   double *tmp=new double[inputNbOfCompo];
1021   std::fill(resPointer,resPointer+inputNbOfCompo*_deno_reverse_multiply.size(),0.);
1022   for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1023     {
1024       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1025         {
1026           isReached[(*iter2).first]=true;
1027           std::transform(inputPointer+idx*inputNbOfCompo,inputPointer+(idx+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/_deno_reverse_multiply[(*iter2).first][idx]));
1028           std::transform(tmp,tmp+inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,std::plus<double>());
1029         }
1030     }
1031   delete [] tmp;
1032   idx=0;
1033   for(std::vector<bool>::const_iterator iter3=isReached.begin();iter3!=isReached.end();iter3++,idx++)
1034     if(!*iter3)
1035       std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1036 }
1037
1038 void MEDCouplingRemapper::ReverseMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
1039 {
1040   matOut.resize(nbColsMatIn);
1041   int id=0;
1042   for(std::vector<std::map<int,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
1043     for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1044       matOut[(*iter2).first][id]=(*iter2).second;
1045 }
1046
1047 void MEDCouplingRemapper::ComputeRowSumAndColSum(const std::vector<std::map<int,double> >& matrixDeno,
1048                                                  std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1049 {
1050   std::map<int,double> values;
1051   int idx=0;
1052   for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1053     {
1054       double sum=0.;
1055       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1056         {
1057           sum+=(*iter2).second;
1058           values[(*iter2).first]+=(*iter2).second;
1059         }
1060       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1061         deno[idx][(*iter2).first]=sum;
1062     }
1063   idx=0;
1064   for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1065     {
1066       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1067         denoReverse[(*iter2).first][idx]=values[(*iter2).first];
1068     }
1069 }
1070
1071 void MEDCouplingRemapper::ComputeColSumAndRowSum(const std::vector<std::map<int,double> >& matrixDeno,
1072                                                  std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1073 {
1074   std::map<int,double> values;
1075   int idx=0;
1076   for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1077     {
1078       double sum=0.;
1079       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1080         {
1081           sum+=(*iter2).second;
1082           values[(*iter2).first]+=(*iter2).second;
1083         }
1084       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1085         denoReverse[(*iter2).first][idx]=sum;
1086     }
1087   idx=0;
1088   for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1089     {
1090       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1091         deno[idx][(*iter2).first]=values[(*iter2).first];
1092     }
1093 }
1094
1095 void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<int,double> >& m1D,
1096                                                                      const std::vector< std::map<int,double> >& m2D,
1097                                                                      const int *corrCellIdSrc, int nbOf2DCellsSrc, int nbOf1DCellsSrc,
1098                                                                      const int *corrCellIdTrg)
1099 {
1100   int nbOf2DCellsTrg=m2D.size();
1101   int nbOf1DCellsTrg=m1D.size();
1102   int nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
1103   _matrix.resize(nbOf3DCellsTrg);
1104   int id2R=0;
1105   for(std::vector< std::map<int,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
1106     {
1107       for(std::map<int,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
1108         {
1109           int id1R=0;
1110           for(std::vector< std::map<int,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
1111             {
1112               for(std::map<int,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
1113                 {
1114                   _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
1115                 }
1116             }
1117         }
1118     }
1119 }
1120
1121 void MEDCouplingRemapper::PrintMatrix(const std::vector<std::map<int,double> >& m)
1122 {
1123   int id=0;
1124   for(std::vector<std::map<int,double> >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++)
1125     {
1126       std::cout << "Target Cell # " << id << " : ";
1127       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1128         std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), ";
1129       std::cout << std::endl;
1130     }
1131 }
1132
1133 const std::vector<std::map<int,double> >& MEDCouplingRemapper::getCrudeMatrix() const
1134 {
1135   return _matrix;
1136 }
1137
1138 /*!
1139  * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1140  * If not the behaviour is unpredictable.
1141  * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than \a maxValAbs this coefficient is
1142  * 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.
1143  * 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
1144  * can occur.
1145  *
1146  * \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.
1147  * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged.
1148  * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix
1149  */
1150 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs(double maxValAbs) throw(INTERP_KERNEL::Exception)
1151 {
1152   int ret=0;
1153   std::vector<std::map<int,double> > matrixNew(_matrix.size());
1154   int i=0;
1155   for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++,i++)
1156     {
1157       std::map<int,double>& rowNew=matrixNew[i];
1158       for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1159         {
1160           if(fabs((*it2).second)>maxValAbs)
1161             rowNew[(*it2).first]=(*it2).second;
1162           else
1163             ret++;
1164         }
1165     }
1166   if(ret>0)
1167     _matrix=matrixNew;
1168   return ret;
1169 }
1170
1171 /*!
1172  * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1173  * If not the behaviour is unpredictable.
1174  * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than delta multiplied by \a scaleFactor this coefficient is
1175  * 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.
1176  * delta is the value returned by MEDCouplingRemapper::getMaxValueInCrudeMatrix method.
1177  * 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
1178  * can occur.
1179  *
1180  * \param [in] scaleFactor is the scale factor from which coefficients lower than \a scaleFactor times range width of coefficients are set to zero.
1181  * \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
1182  *         that all coefficients are null.
1183  * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs
1184  */
1185 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix(double scaleFactor) throw(INTERP_KERNEL::Exception)
1186 {
1187   double maxVal=getMaxValueInCrudeMatrix();
1188   if(maxVal==0.)
1189     return -1;
1190   return nullifiedTinyCoeffInCrudeMatrixAbs(scaleFactor*maxVal);
1191 }
1192
1193 /*!
1194  * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1195  * If not the behaviour is unpredictable.
1196  * This method returns the maximum of the absolute values of coefficients into the sparse crude matrix.
1197  * The returned value is positive.
1198  */
1199 double MEDCouplingRemapper::getMaxValueInCrudeMatrix() const throw(INTERP_KERNEL::Exception)
1200 {
1201   double ret=0.;
1202   for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++)
1203     for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1204       if(fabs((*it2).second)>ret)
1205         ret=fabs((*it2).second);
1206   return ret;
1207 }