src/MEDCoupling/MEDCouplingRemapper.cxx

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