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