src/MEDCoupling/MEDCouplingRemapper.cxx

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