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   if(InterpolationOptions::getIntersectionType()!=INTERP_KERNEL::Triangulation)
 681       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC: only 'Triangulation' intersection type supported!");
 682   const MEDCouplingUMesh *src_mesh=static_cast<const MEDCouplingUMesh *>(_src_ft->getMesh());
 683   const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
 684   const int srcMeshDim=src_mesh->getMeshDimension();
 685   const int srcSpceDim=src_mesh->getSpaceDimension();
 686   const int trgMeshDim=target_mesh->getMeshDimension();
 687   if(srcMeshDim!=srcSpceDim || srcMeshDim!=trgMeshDim)
 688     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC: space dimension of unstructured source mesh should be equal to mesh dimension of unstructured source mesh, and should also be equal to target cartesian dimension!");
 689   std::vector<std::map<int,double> > res;
 690   switch(srcMeshDim)
 691   {
 692     case 1:
 693       {
 694         MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
 695         MEDCouplingNormalizedUnstructuredMesh<1,1> sourceWrapper(src_mesh);
 696         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 697         myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
 698         break;
 699       }
 700     case 2:
 701       {
 702         MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
 703         MEDCouplingNormalizedUnstructuredMesh<2,2> sourceWrapper(src_mesh);
 704         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 705         myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
 706         break;
 707       }
 708     case 3:
 709       {
 710         MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
 711         MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(src_mesh);
 712         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 713         myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
 714         break;
 715       }
 716     default:
 717       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only dimension 1 2 or 3 supported !");
 718   }
 719   ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix);
 720   nullifiedTinyCoeffInCrudeMatrixAbs(0.);
 721   //
 722   _deno_multiply.clear();
 723   _deno_multiply.resize(_matrix.size());
 724   _deno_reverse_multiply.clear();
 725   _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
 726   declareAsNew();
 727   return 1;
 728 }
 729
 730 int MEDCouplingRemapper::prepareInterpKernelOnlyCU()
 731 {
 732   std::string srcMeth,trgMeth;
 733   std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
 734   if(methodCpp!="P0P0")
 735     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only P0P0 interpolation supported for the moment !");
 736   if(InterpolationOptions::getIntersectionType()!=INTERP_KERNEL::Triangulation)
 737     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU: only 'Triangulation' intersection type supported!");
 738   const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
 739   const MEDCouplingUMesh *target_mesh=static_cast<const MEDCouplingUMesh *>(_target_ft->getMesh());
 740   const int srcMeshDim=src_mesh->getMeshDimension();
 741   const int trgMeshDim=target_mesh->getMeshDimension();
 742   const int trgSpceDim=target_mesh->getSpaceDimension();
 743   if(trgMeshDim!=trgSpceDim || trgMeshDim!=srcMeshDim)
 744     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC: space dimension of unstructured target mesh should be equal to mesh dimension of unstructured target mesh, and should also be equal to source cartesian dimension!");
 745   switch(srcMeshDim)
 746   {
 747     case 1:
 748       {
 749         MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
 750         MEDCouplingNormalizedUnstructuredMesh<1,1> targetWrapper(target_mesh);
 751         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 752         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 753         break;
 754       }
 755     case 2:
 756       {
 757         MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
 758         MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(target_mesh);
 759         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 760         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 761         break;
 762       }
 763     case 3:
 764       {
 765         MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
 766         MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(target_mesh);
 767         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 768         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 769         break;
 770       }
 771     default:
 772       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only dimension 1 2 or 3 supported !");
 773   }
 774   nullifiedTinyCoeffInCrudeMatrixAbs(0.);
 775   //
 776   _deno_multiply.clear();
 777   _deno_multiply.resize(_matrix.size());
 778   _deno_reverse_multiply.clear();
 779   _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
 780   declareAsNew();
 781   return 1;
 782 }
 783
 784 int MEDCouplingRemapper::prepareInterpKernelOnlyCC()
 785 {
 786   std::string srcMeth,trgMeth;
 787   std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
 788   if(methodCpp!="P0P0")
 789     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only P0P0 interpolation supported for the moment !");
 790   if(InterpolationOptions::getIntersectionType()!=INTERP_KERNEL::Triangulation)
 791     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC: only 'Triangulation' intersection type supported!");
 792   const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
 793   const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
 794   const int srcMeshDim=src_mesh->getMeshDimension();
 795   const int trgMeshDim=target_mesh->getMeshDimension();
 796   if(trgMeshDim!=srcMeshDim)
 797     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : dimension of target cartesian mesh should be equal to dimension of source cartesian mesh !");
 798   switch(srcMeshDim)
 799   {
 800     case 1:
 801       {
 802         MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
 803         MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
 804         INTERP_KERNEL::InterpolationCC myInterpolator(*this);
 805         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 806         break;
 807       }
 808     case 2:
 809       {
 810         MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
 811         MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
 812         INTERP_KERNEL::InterpolationCC myInterpolator(*this);
 813         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 814         break;
 815       }
 816     case 3:
 817       {
 818         MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
 819         MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
 820         INTERP_KERNEL::InterpolationCC myInterpolator(*this);
 821         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 822         break;
 823       }
 824     default:
 825       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only dimension 1 2 or 3 supported !");
 826   }
 827   nullifiedTinyCoeffInCrudeMatrixAbs(0.);
 828   //
 829   _deno_multiply.clear();
 830   _deno_multiply.resize(_matrix.size());
 831   _deno_reverse_multiply.clear();
 832   _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
 833   declareAsNew();
 834   return 1;
 835 }
 836
 837 int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss()
 838 {
 839   if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
 840     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 !");
 841   MCAuto<DataArrayDouble> trgLoc=_target_ft->getLocalizationOfDiscr();
 842   const double *trgLocPtr=trgLoc->begin();
 843   int trgSpaceDim=trgLoc->getNumberOfComponents();
 844   MCAuto<DataArrayInt> srcOffsetArr=_src_ft->getDiscretization()->getOffsetArr(_src_ft->getMesh());
 845   if(trgSpaceDim!=_src_ft->getMesh()->getSpaceDimension())
 846     {
 847       std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss : space dimensions mismatch between source and target !";
 848       oss << " Target discretization localization has dimension " << trgSpaceDim << ", whereas the space dimension of source is equal to ";
 849       oss << _src_ft->getMesh()->getSpaceDimension() << " !";
 850       throw INTERP_KERNEL::Exception(oss.str().c_str());
 851     }
 852   const int *srcOffsetArrPtr=srcOffsetArr->begin();
 853   MCAuto<DataArrayDouble> srcLoc=_src_ft->getLocalizationOfDiscr();
 854   const double *srcLocPtr=srcLoc->begin();
 855   MCAuto<DataArrayInt> eltsArr,eltsIndexArr;
 856   int trgNbOfGaussPts=trgLoc->getNumberOfTuples();
 857   _matrix.resize(trgNbOfGaussPts);
 858   _src_ft->getMesh()->getCellsContainingPoints(trgLoc->begin(),trgNbOfGaussPts,getPrecision(),eltsArr,eltsIndexArr);
 859   const int *elts(eltsArr->begin()),*eltsIndex(eltsIndexArr->begin());
 860   MCAuto<DataArrayInt> nbOfSrcCellsShTrgPts(eltsIndexArr->deltaShiftIndex());
 861   MCAuto<DataArrayInt> ids0=nbOfSrcCellsShTrgPts->findIdsNotEqual(0);
 862   for(const int *trgId=ids0->begin();trgId!=ids0->end();trgId++)
 863     {
 864       const double *ptTrg=trgLocPtr+trgSpaceDim*(*trgId);
 865       int srcCellId=elts[eltsIndex[*trgId]];
 866       double dist=std::numeric_limits<double>::max();
 867       int srcEntry=-1;
 868       for(int srcId=srcOffsetArrPtr[srcCellId];srcId<srcOffsetArrPtr[srcCellId+1];srcId++)
 869         {
 870           const double *ptSrc=srcLocPtr+trgSpaceDim*srcId;
 871           double tmp=0.;
 872           for(int i=0;i<trgSpaceDim;i++)
 873             tmp+=(ptTrg[i]-ptSrc[i])*(ptTrg[i]-ptSrc[i]);
 874           if(tmp<dist)
 875             { dist=tmp; srcEntry=srcId; }
 876         }
 877       _matrix[*trgId][srcEntry]=1.;
 878     }
 879   if(ids0->getNumberOfTuples()!=trgNbOfGaussPts)
 880     {
 881       MCAuto<DataArrayInt> orphanTrgIds=nbOfSrcCellsShTrgPts->findIdsEqual(0);
 882       MCAuto<DataArrayDouble> orphanTrg=trgLoc->selectByTupleId(orphanTrgIds->begin(),orphanTrgIds->end());
 883       MCAuto<DataArrayInt> srcIdPerTrg=srcLoc->findClosestTupleId(orphanTrg);
 884       const int *srcIdPerTrgPtr=srcIdPerTrg->begin();
 885       for(const int *orphanTrgId=orphanTrgIds->begin();orphanTrgId!=orphanTrgIds->end();orphanTrgId++,srcIdPerTrgPtr++)
 886         _matrix[*orphanTrgId][*srcIdPerTrgPtr]=2.;
 887     }
 888   _deno_multiply.clear();
 889   _deno_multiply.resize(_matrix.size());
 890   _deno_reverse_multiply.clear();
 891   _deno_reverse_multiply.resize(srcLoc->getNumberOfTuples());
 892   declareAsNew();
 893   return 1;
 894 }
 895
 896 /*!
 897  * This method checks that the input interpolation \a method is managed by not INTERP_KERNEL only methods.
 898  * If no an INTERP_KERNEL::Exception will be thrown. If yes, a magic number will be returned to switch in the MEDCouplingRemapper::prepareNotInterpKernelOnly method.
 899  */
 900 int MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel(const std::string& method)
 901 {
 902   if(method=="GAUSSGAUSS")
 903     return 0;
 904   std::ostringstream oss; oss << "MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel : ";
 905   oss << "The method \"" << method << "\" is not manageable by not INTERP_KERNEL only method.";
 906   oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS !";
 907   throw INTERP_KERNEL::Exception(oss.str().c_str());
 908 }
 909
 910 /*!
 911  * This method determines regarding \c _interp_matrix_pol attribute ( set by MEDCouplingRemapper::setInterpolationMatrixPolicy and by default equal
 912  * 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
 913  * only INTERP_KERNEL method should be applied.
 914  */
 915 bool MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly() const
 916 {
 917   std::string srcm,trgm,method;
 918   method=checkAndGiveInterpolationMethodStr(srcm,trgm);
 919   switch(_interp_matrix_pol)
 920   {
 921     case IK_ONLY_PREFERED:
 922       {
 923         try
 924         {
 925             std::string tmp1,tmp2;
 926             INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,tmp1,tmp2);
 927             return true;
 928         }
 929         catch(INTERP_KERNEL::Exception& /*e*/)
 930         {
 931             return false;
 932         }
 933       }
 934     case NOT_IK_ONLY_PREFERED:
 935       {
 936         try
 937         {
 938             CheckInterpolationMethodManageableByNotOnlyInterpKernel(method);
 939             return false;
 940         }
 941         catch(INTERP_KERNEL::Exception& /*e*/)
 942         {
 943             return true;
 944         }
 945       }
 946     case IK_ONLY_FORCED:
 947       return true;
 948     case NOT_IK_ONLY_FORCED:
 949       return false;
 950     default:
 951       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly : internal error ! The interpolation matrix policy is not managed ! Try to change it using MEDCouplingRemapper::setInterpolationMatrixPolicy !");
 952   }
 953 }
 954
 955 void MEDCouplingRemapper::updateTime() const
 956 {
 957 }
 958
 959 void MEDCouplingRemapper::checkPrepare() const
 960 {
 961   const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
 962   if(!s || !t)
 963     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that MEDCouplingRemapper::prepare(Ex) has not been called !");
 964   if(!s->getMesh() || !t->getMesh())
 965     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that no all field templates have their mesh set !");
 966 }
 967
 968 /*!
 969  * This method builds a code considering already set field discretization int \a this : \a _src_ft and \a _target_ft.
 970  * This method returns 3 information (2 in output parameters and 1 in return).
 971  *
 972  * \param [out] srcMeth the string code of the discretization of source field template
 973  * \param [out] trgMeth the string code of the discretization of target field template
 974  * \return the standardized string code (compatible with INTERP_KERNEL) for matrix of numerators (in \a _matrix)
 975  */
 976 std::string MEDCouplingRemapper::checkAndGiveInterpolationMethodStr(std::string& srcMeth, std::string& trgMeth) const
 977 {
 978   const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
 979   if(!s || !t)
 980     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have been set !");
 981   if(!s->getMesh() || !t->getMesh())
 982     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have their mesh set !");
 983   srcMeth=_src_ft->getDiscretization()->getRepr();
 984   trgMeth=_target_ft->getDiscretization()->getRepr();
 985   return BuildMethodFrom(srcMeth,trgMeth);
 986 }
 987
 988 std::string MEDCouplingRemapper::BuildMethodFrom(const std::string& meth1, const std::string& meth2)
 989 {
 990   std::string method(meth1); method+=meth2;
 991   return method;
 992 }
 993
 994 void MEDCouplingRemapper::releaseData(bool matrixSuppression)
 995 {
 996   _src_ft=0;
 997   _target_ft=0;
 998   if(matrixSuppression)
 999     {
1000       _matrix.clear();
1001       _deno_multiply.clear();
1002       _deno_reverse_multiply.clear();
1003     }
1004 }
1005
1006 void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue)
1007 {
1008   if(!srcField || !targetField)
1009     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferUnderground : srcField or targetField is NULL !");
1010   srcField->checkConsistencyLight();
1011   checkPrepare();
1012   if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
1013     throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
1014   if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
1015     throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
1016   if(srcField->getNature()!=targetField->getNature())
1017     throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
1018   if(srcField->getNumberOfTuplesExpected()!=_src_ft->getNumberOfTuplesExpected())
1019     {
1020       std::ostringstream oss;
1021       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();
1022       oss << " ! It appears that the source support is not the same between the prepare and the transfer !";
1023       throw INTERP_KERNEL::Exception(oss.str().c_str());
1024     }
1025   DataArrayDouble *array(targetField->getArray());
1026   int srcNbOfCompo(srcField->getNumberOfComponents());
1027   if(array)
1028     {
1029       targetField->checkConsistencyLight();
1030       if(srcNbOfCompo!=targetField->getNumberOfComponents())
1031         throw INTERP_KERNEL::Exception("Number of components mismatch !");
1032     }
1033   else
1034     {
1035       if(!isDftVal)
1036         throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !");
1037       MCAuto<DataArrayDouble> tmp(DataArrayDouble::New());
1038       tmp->alloc(targetField->getNumberOfTuples(),srcNbOfCompo);
1039       targetField->setArray(tmp);
1040     }
1041   computeDeno(srcField->getNature(),srcField,targetField);
1042   double *resPointer(targetField->getArray()->getPointer());
1043   const double *inputPointer(srcField->getArray()->getConstPointer());
1044   computeProduct(inputPointer,srcNbOfCompo,isDftVal,dftValue,resPointer);
1045 }
1046
1047 void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1048 {
1049   if(nat==NoNature)
1050     return computeDenoFromScratch(nat,srcField,trgField);
1051   else if(nat!=_nature_of_deno)
1052     return computeDenoFromScratch(nat,srcField,trgField);
1053   else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis())
1054     return computeDenoFromScratch(nat,srcField,trgField);
1055 }
1056
1057 void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1058 {
1059   _nature_of_deno=nat;
1060   _time_deno_update=getTimeOfThis();
1061   switch(_nature_of_deno)
1062   {
1063     case IntensiveMaximum:
1064       {
1065         ComputeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1066         break;
1067       }
1068     case ExtensiveMaximum:
1069       {
1070         MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1071         MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1072         const double *denoPtr=deno->getArray()->getConstPointer();
1073         const double *denoRPtr=denoR->getArray()->getConstPointer();
1074         if(trgField->getMesh()->getMeshDimension()==-1)
1075           {
1076             double *denoRPtr2=denoR->getArray()->getPointer();
1077             denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1078           }
1079         if(srcField->getMesh()->getMeshDimension()==-1)
1080           {
1081             double *denoPtr2=deno->getArray()->getPointer();
1082             denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1083           }
1084         int idx=0;
1085         for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1086           for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1087             {
1088               _deno_multiply[idx][(*iter2).first]=denoPtr[(*iter2).first];
1089               _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[idx];
1090             }
1091         deno->decrRef();
1092         denoR->decrRef();
1093         break;
1094       }
1095     case ExtensiveConservation:
1096       {
1097         ComputeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1098         break;
1099       }
1100     case IntensiveConservation:
1101       {
1102         MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1103         MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1104         const double *denoPtr=deno->getArray()->getConstPointer();
1105         const double *denoRPtr=denoR->getArray()->getConstPointer();
1106         if(trgField->getMesh()->getMeshDimension()==-1)
1107           {
1108             double *denoRPtr2=denoR->getArray()->getPointer();
1109             denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1110           }
1111         if(srcField->getMesh()->getMeshDimension()==-1)
1112           {
1113             double *denoPtr2=deno->getArray()->getPointer();
1114             denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1115           }
1116         int idx=0;
1117         for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1118           for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1119             {
1120               _deno_multiply[idx][(*iter2).first]=denoPtr[idx];
1121               _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[(*iter2).first];
1122             }
1123         deno->decrRef();
1124         denoR->decrRef();
1125         break;
1126       }
1127     case NoNature:
1128       throw INTERP_KERNEL::Exception("No nature specified ! Select one !");
1129   }
1130 }
1131
1132 void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)
1133 {
1134   int idx=0;
1135   double *tmp=new double[inputNbOfCompo];
1136   for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1137     {
1138       if((*iter1).empty())
1139         {
1140           if(isDftVal)
1141             std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1142           continue;
1143         }
1144       else
1145         std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,0.);
1146       std::map<int,double>::const_iterator iter3=_deno_multiply[idx].begin();
1147       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter3++)
1148         {
1149           std::transform(inputPointer+(*iter2).first*inputNbOfCompo,inputPointer+((*iter2).first+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/(*iter3).second));
1150           std::transform(tmp,tmp+inputNbOfCompo,resPointer+idx*inputNbOfCompo,resPointer+idx*inputNbOfCompo,std::plus<double>());
1151         }
1152     }
1153   delete [] tmp;
1154 }
1155
1156 void MEDCouplingRemapper::computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer)
1157 {
1158   std::vector<bool> isReached(_deno_reverse_multiply.size(),false);
1159   int idx=0;
1160   double *tmp=new double[inputNbOfCompo];
1161   std::fill(resPointer,resPointer+inputNbOfCompo*_deno_reverse_multiply.size(),0.);
1162   for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1163     {
1164       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1165         {
1166           isReached[(*iter2).first]=true;
1167           std::transform(inputPointer+idx*inputNbOfCompo,inputPointer+(idx+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/_deno_reverse_multiply[(*iter2).first][idx]));
1168           std::transform(tmp,tmp+inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,std::plus<double>());
1169         }
1170     }
1171   delete [] tmp;
1172   idx=0;
1173   for(std::vector<bool>::const_iterator iter3=isReached.begin();iter3!=isReached.end();iter3++,idx++)
1174     if(!*iter3)
1175       std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1176 }
1177
1178 void MEDCouplingRemapper::ReverseMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
1179 {
1180   matOut.resize(nbColsMatIn);
1181   int id=0;
1182   for(std::vector<std::map<int,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
1183     for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1184       matOut[(*iter2).first][id]=(*iter2).second;
1185 }
1186
1187 void MEDCouplingRemapper::ComputeRowSumAndColSum(const std::vector<std::map<int,double> >& matrixDeno,
1188                                                  std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1189 {
1190   std::map<int,double> values;
1191   int idx=0;
1192   for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1193     {
1194       double sum=0.;
1195       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1196         {
1197           sum+=(*iter2).second;
1198           values[(*iter2).first]+=(*iter2).second;
1199         }
1200       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1201         deno[idx][(*iter2).first]=sum;
1202     }
1203   idx=0;
1204   for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1205     {
1206       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1207         denoReverse[(*iter2).first][idx]=values[(*iter2).first];
1208     }
1209 }
1210
1211 void MEDCouplingRemapper::ComputeColSumAndRowSum(const std::vector<std::map<int,double> >& matrixDeno,
1212                                                  std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1213 {
1214   std::map<int,double> values;
1215   int idx=0;
1216   for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1217     {
1218       double sum=0.;
1219       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1220         {
1221           sum+=(*iter2).second;
1222           values[(*iter2).first]+=(*iter2).second;
1223         }
1224       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1225         denoReverse[(*iter2).first][idx]=sum;
1226     }
1227   idx=0;
1228   for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1229     {
1230       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1231         deno[idx][(*iter2).first]=values[(*iter2).first];
1232     }
1233 }
1234
1235 void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<int,double> >& m1D,
1236                                                                      const std::vector< std::map<int,double> >& m2D,
1237                                                                      const int *corrCellIdSrc, int nbOf2DCellsSrc, int nbOf1DCellsSrc,
1238                                                                      const int *corrCellIdTrg)
1239 {
1240   int nbOf2DCellsTrg=m2D.size();
1241   int nbOf1DCellsTrg=m1D.size();
1242   int nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
1243   _matrix.resize(nbOf3DCellsTrg);
1244   int id2R=0;
1245   for(std::vector< std::map<int,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
1246     {
1247       for(std::map<int,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
1248         {
1249           int id1R=0;
1250           for(std::vector< std::map<int,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
1251             {
1252               for(std::map<int,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
1253                 {
1254                   _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
1255                 }
1256             }
1257         }
1258     }
1259 }
1260
1261 void MEDCouplingRemapper::PrintMatrix(const std::vector<std::map<int,double> >& m)
1262 {
1263   int id=0;
1264   for(std::vector<std::map<int,double> >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++)
1265     {
1266       std::cout << "Target Cell # " << id << " : ";
1267       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1268         std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), ";
1269       std::cout << std::endl;
1270     }
1271 }
1272
1273 const std::vector<std::map<int,double> >& MEDCouplingRemapper::getCrudeMatrix() const
1274 {
1275   return _matrix;
1276 }
1277
1278 /*!
1279  * Returns the number of columns of matrix returned by MEDCouplingRemapper::getCrudeMatrix method.
1280  */
1281 int MEDCouplingRemapper::getNumberOfColsOfMatrix() const
1282 {
1283   return (int)_deno_reverse_multiply.size();
1284 }
1285
1286 /*!
1287  * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1288  * If not the behaviour is unpredictable.
1289  * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than \a maxValAbs this coefficient is
1290  * 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.
1291  * 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
1292  * can occur.
1293  *
1294  * \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.
1295  * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged.
1296  * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix
1297  */
1298 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs(double maxValAbs)
1299 {
1300   int ret=0;
1301   std::vector<std::map<int,double> > matrixNew(_matrix.size());
1302   int i=0;
1303   for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++,i++)
1304     {
1305       std::map<int,double>& rowNew=matrixNew[i];
1306       for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1307         {
1308           if(fabs((*it2).second)>maxValAbs)
1309             rowNew[(*it2).first]=(*it2).second;
1310           else
1311             ret++;
1312         }
1313     }
1314   if(ret>0)
1315     _matrix=matrixNew;
1316   return ret;
1317 }
1318
1319 /*!
1320  * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1321  * If not the behaviour is unpredictable.
1322  * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than delta multiplied by \a scaleFactor this coefficient is
1323  * 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.
1324  * delta is the value returned by MEDCouplingRemapper::getMaxValueInCrudeMatrix method.
1325  * 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
1326  * can occur.
1327  *
1328  * \param [in] scaleFactor is the scale factor from which coefficients lower than \a scaleFactor times range width of coefficients are set to zero.
1329  * \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
1330  *         that all coefficients are null.
1331  * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs
1332  */
1333 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix(double scaleFactor)
1334 {
1335   double maxVal=getMaxValueInCrudeMatrix();
1336   if(maxVal==0.)
1337     return -1;
1338   return nullifiedTinyCoeffInCrudeMatrixAbs(scaleFactor*maxVal);
1339 }
1340
1341 /*!
1342  * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1343  * If not the behaviour is unpredictable.
1344  * This method returns the maximum of the absolute values of coefficients into the sparse crude matrix.
1345  * The returned value is positive.
1346  */
1347 double MEDCouplingRemapper::getMaxValueInCrudeMatrix() const
1348 {
1349   double ret=0.;
1350   for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++)
1351     for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1352       if(fabs((*it2).second)>ret)
1353         ret=fabs((*it2).second);
1354   return ret;
1355 }