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