src/MEDCoupling/MEDCouplingRemapper.cxx

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