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