src/MEDCoupling/MEDCouplingRemapper.cxx

   1 // Copyright (C) 2007-2019  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<mcIdType,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<mcIdType,double> >& m)
 105 {
 106   restartUsing(src,target);
 107   if(ToIdType(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   std::size_t trgNbOfCompo=targetField->getNumberOfComponents();
 249   if(array)
 250     {
 251       srcField->checkConsistencyLight();
 252       if(ToIdType(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,(int)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   mcIdType 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==3 && trgMeshDim==0 && srcSpaceDim==3)
 468     {
 469       if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
 470         throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 0D ! 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::Interpolation3D1D interpolation(*this);//Not a bug : 3D1D deal with 3D0D
 474       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 475     }
 476   else if(srcMeshDim==1 && trgMeshDim==0 && srcSpaceDim==3)
 477     {
 478       if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
 479         throw INTERP_KERNEL::Exception("Invalid interpolation requested between 1D and 0D into 3D space ! 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::Interpolation1D0D interpolation(*this);
 483       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 484     }
 485   else if(srcMeshDim==1 && trgMeshDim==3 && srcSpaceDim==3)
 486     {
 487       if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
 488         throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
 489       MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 490       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
 491       INTERP_KERNEL::Interpolation3D1D interpolation(*this);
 492       std::vector<std::map<mcIdType,double> > matrixTmp;
 493       std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
 494       nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
 495       ReverseMatrix(matrixTmp,nbCols,_matrix);
 496       nbCols=ToIdType(matrixTmp.size());
 497     }
 498   else if(srcMeshDim==2 && trgMeshDim==1 && srcSpaceDim==2)
 499     {
 500       if(getIntersectionType()==INTERP_KERNEL::PointLocator)
 501         {
 502           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 503           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
 504           INTERP_KERNEL::Interpolation2D interpolation(*this);
 505           nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 506         }
 507       else
 508         {
 509           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 510           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
 511           INTERP_KERNEL::Interpolation2D1D interpolation(*this);
 512           std::vector<std::map<mcIdType,double> > matrixTmp;
 513           std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
 514           nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
 515           ReverseMatrix(matrixTmp,nbCols,_matrix);
 516           nbCols=ToIdType(matrixTmp.size());
 517           INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
 518           if(!duplicateFaces.empty())
 519             {
 520               std::ostringstream oss; oss << "An unexpected situation happened ! For the following 1D Cells are part of edges shared by 2D cells :\n";
 521               for(std::map<mcIdType,std::set<mcIdType> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
 522                 {
 523                   oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
 524                   std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<mcIdType>(oss," "));
 525                   oss << std::endl;
 526                 }
 527             }
 528         }
 529     }
 530   else if(srcMeshDim==2 && trgMeshDim==0 && srcSpaceDim==2)
 531     {
 532       if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
 533         throw INTERP_KERNEL::Exception("Invalid interpolation requested between 2D and 0D ! Select PointLocator as intersection type !");
 534       MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 535       MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
 536       INTERP_KERNEL::Interpolation2D interpolation(*this);
 537       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 538     }
 539   else if(srcMeshDim==1 && trgMeshDim==2 && srcSpaceDim==2)
 540     {
 541       if(getIntersectionType()==INTERP_KERNEL::PointLocator)
 542         {
 543           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 544           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
 545           INTERP_KERNEL::Interpolation2D interpolation(*this);
 546           std::vector<std::map<mcIdType,double> > matrixTmp;
 547           std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
 548           nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
 549           ReverseMatrix(matrixTmp,nbCols,_matrix);
 550           nbCols=ToIdType(matrixTmp.size());
 551         }
 552       else
 553         {
 554           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 555           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
 556           INTERP_KERNEL::Interpolation2D1D interpolation(*this);
 557           nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 558           INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
 559           if(!duplicateFaces.empty())
 560             {
 561               std::ostringstream oss; oss << "An unexpected situation happened ! For the following 1D Cells are part of edges shared by 2D cells :\n";
 562               for(std::map<mcIdType,std::set<mcIdType> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
 563                 {
 564                   oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
 565                   std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<mcIdType>(oss," "));
 566                   oss << std::endl;
 567                 }
 568             }
 569         }
 570     }
 571   else if(srcMeshDim==2 && trgMeshDim==3 && srcSpaceDim==3)
 572     {
 573       if(getIntersectionType()==INTERP_KERNEL::PointLocator)
 574         throw INTERP_KERNEL::Exception("Using point locator to transfer a mesh of dim 2 to a mesh of dim 3 does not make sense: 3D centers of mass can not be localized in a mesh having mesh-dim=2, space-dim=3!!");
 575       else
 576         {
 577           MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 578           MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
 579           INTERP_KERNEL::Interpolation2D3D interpolation(*this);
 580           nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 581           INTERP_KERNEL::Interpolation2D3D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
 582           if(!duplicateFaces.empty())
 583             {
 584               std::ostringstream oss; oss << "An unexpected situation happened ! For the following 2D Cells are part of edges shared by 3D cells :\n";
 585               for(std::map<mcIdType,std::set<mcIdType> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
 586                 {
 587                   oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
 588                   std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<mcIdType>(oss," "));
 589                   oss << std::endl;
 590                 }
 591             }
 592         }
 593     }
 594   else if(srcMeshDim==3 && trgMeshDim==2 && srcSpaceDim==3)
 595     {
 596       if(getIntersectionType()==INTERP_KERNEL::PointLocator)
 597         {
 598           MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 599           MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
 600           INTERP_KERNEL::Interpolation3D interpolation(*this);
 601           nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 602         }
 603       else
 604         {
 605           MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 606           MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
 607           INTERP_KERNEL::Interpolation2D3D interpolation(*this);
 608           std::vector<std::map<mcIdType,double> > matrixTmp;
 609           std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
 610           nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
 611           ReverseMatrix(matrixTmp,nbCols,_matrix);
 612           nbCols=ToIdType(matrixTmp.size());
 613           INTERP_KERNEL::Interpolation2D3D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
 614           if(!duplicateFaces.empty())
 615             {
 616               std::ostringstream oss; oss << "An unexpected situation happened ! For the following 2D Cells are part of edges shared by 3D cells :\n";
 617               for(std::map<mcIdType,std::set<mcIdType> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
 618                 {
 619                   oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
 620                   std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<mcIdType>(oss," "));
 621                   oss << std::endl;
 622                 }
 623             }
 624         }
 625     }
 626   else if(trgMeshDim==-1)
 627     {
 628       if(srcMeshDim==2 && srcSpaceDim==2)
 629         {
 630           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 631           INTERP_KERNEL::Interpolation2D interpolation(*this);
 632           nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
 633         }
 634       else if(srcMeshDim==3 && srcSpaceDim==3)
 635         {
 636           MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 637           INTERP_KERNEL::Interpolation3D interpolation(*this);
 638           nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
 639         }
 640       else if(srcMeshDim==2 && srcSpaceDim==3)
 641         {
 642           MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
 643           INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
 644           nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
 645         }
 646       else
 647         throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
 648     }
 649   else if(srcMeshDim==-1)
 650     {
 651       if(trgMeshDim==2 && trgSpaceDim==2)
 652         {
 653           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(target_mesh);
 654           INTERP_KERNEL::Interpolation2D interpolation(*this);
 655           nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
 656         }
 657       else if(trgMeshDim==3 && trgSpaceDim==3)
 658         {
 659           MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(target_mesh);
 660           INTERP_KERNEL::Interpolation3D interpolation(*this);
 661           nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
 662         }
 663       else if(trgMeshDim==2 && trgSpaceDim==3)
 664         {
 665           MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(target_mesh);
 666           INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
 667           nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
 668         }
 669       else
 670         throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh from -1D sourceMesh");
 671     }
 672   else
 673     throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension");
 674   _deno_multiply.clear();
 675   _deno_multiply.resize(_matrix.size());
 676   _deno_reverse_multiply.clear();
 677   _deno_reverse_multiply.resize(nbCols);
 678   declareAsNew();
 679   return 1;
 680 }
 681
 682 int MEDCouplingRemapper::prepareInterpKernelOnlyEE()
 683 {
 684   std::string srcMeth,trgMeth;
 685   std::string methC=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
 686   const MEDCouplingMappedExtrudedMesh *src_mesh=static_cast<const MEDCouplingMappedExtrudedMesh *>(_src_ft->getMesh());
 687   const MEDCouplingMappedExtrudedMesh *target_mesh=static_cast<const MEDCouplingMappedExtrudedMesh *>(_target_ft->getMesh());
 688   if(methC!="P0P0")
 689     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyEE : Only P0P0 method implemented for Extruded/Extruded meshes !");
 690   MCAuto<MEDCouplingUMesh> src2D(src_mesh->getMesh2D()->clone(false)); src2D->changeSpaceDimension(2,0.);
 691   MCAuto<MEDCouplingUMesh> trg2D(target_mesh->getMesh2D()->clone(false)); trg2D->changeSpaceDimension(2,0.);
 692   MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src2D);
 693   MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(trg2D);
 694   INTERP_KERNEL::Interpolation2D interpolation2D(*this);
 695   std::vector<std::map<mcIdType,double> > matrix2D;
 696   mcIdType nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,methC);
 697   MEDCouplingUMesh *s1D,*t1D;
 698   double v[3];
 699   MEDCouplingMappedExtrudedMesh::Project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v);
 700   MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D);
 701   MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
 702   std::vector<std::map<mcIdType,double> > matrix1D;
 703   INTERP_KERNEL::Interpolation1D interpolation1D(*this);
 704   if(interpolation1D.getIntersectionType()==INTERP_KERNEL::Geometric2D)// For intersection type of 1D, Geometric2D do not deal with it ! -> make interpolation1D not inherite from this
 705     interpolation1D.setIntersectionType(INTERP_KERNEL::Triangulation);//
 706   mcIdType nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,methC);
 707   s1D->decrRef();
 708   t1D->decrRef();
 709   buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
 710                                              target_mesh->getMesh3DIds()->getConstPointer());
 711   //
 712   _deno_multiply.clear();
 713   _deno_multiply.resize(_matrix.size());
 714   _deno_reverse_multiply.clear();
 715   _deno_reverse_multiply.resize(nbCols2D*nbCols1D);
 716   declareAsNew();
 717   return 1;
 718 }
 719
 720 int MEDCouplingRemapper::prepareInterpKernelOnlyUC()
 721 {
 722   std::string srcMeth,trgMeth;
 723   std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
 724   if(methodCpp!="P0P0")
 725     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC: only P0P0 interpolation supported for the moment !");
 726   if(InterpolationOptions::getIntersectionType()!=INTERP_KERNEL::Triangulation)
 727       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC: only 'Triangulation' intersection type supported!");
 728   const MEDCouplingUMesh *src_mesh=static_cast<const MEDCouplingUMesh *>(_src_ft->getMesh());
 729   const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
 730   const int srcMeshDim=src_mesh->getMeshDimension();
 731   const int srcSpceDim=src_mesh->getSpaceDimension();
 732   const int trgMeshDim=target_mesh->getMeshDimension();
 733   if(srcMeshDim!=srcSpceDim || srcMeshDim!=trgMeshDim)
 734     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!");
 735   std::vector<std::map<mcIdType,double> > res;
 736   switch(srcMeshDim)
 737   {
 738     case 1:
 739       {
 740         MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
 741         MEDCouplingNormalizedUnstructuredMesh<1,1> sourceWrapper(src_mesh);
 742         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 743         myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
 744         break;
 745       }
 746     case 2:
 747       {
 748         MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
 749         MEDCouplingNormalizedUnstructuredMesh<2,2> sourceWrapper(src_mesh);
 750         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 751         myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
 752         break;
 753       }
 754     case 3:
 755       {
 756         MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
 757         MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(src_mesh);
 758         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 759         myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
 760         break;
 761       }
 762     default:
 763       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only dimension 1 2 or 3 supported !");
 764   }
 765   ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix);
 766   nullifiedTinyCoeffInCrudeMatrixAbs(0.);
 767   //
 768   _deno_multiply.clear();
 769   _deno_multiply.resize(_matrix.size());
 770   _deno_reverse_multiply.clear();
 771   _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
 772   declareAsNew();
 773   return 1;
 774 }
 775
 776 int MEDCouplingRemapper::prepareInterpKernelOnlyCU()
 777 {
 778   std::string srcMeth,trgMeth;
 779   std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
 780   if(methodCpp!="P0P0")
 781     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only P0P0 interpolation supported for the moment !");
 782   if(InterpolationOptions::getIntersectionType()!=INTERP_KERNEL::Triangulation)
 783     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU: only 'Triangulation' intersection type supported!");
 784   const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
 785   const MEDCouplingUMesh *target_mesh=static_cast<const MEDCouplingUMesh *>(_target_ft->getMesh());
 786   const int srcMeshDim=src_mesh->getMeshDimension();
 787   const int trgMeshDim=target_mesh->getMeshDimension();
 788   const int trgSpceDim=target_mesh->getSpaceDimension();
 789   if(trgMeshDim!=trgSpceDim || trgMeshDim!=srcMeshDim)
 790     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!");
 791   switch(srcMeshDim)
 792   {
 793     case 1:
 794       {
 795         MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
 796         MEDCouplingNormalizedUnstructuredMesh<1,1> targetWrapper(target_mesh);
 797         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 798         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 799         break;
 800       }
 801     case 2:
 802       {
 803         MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
 804         MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(target_mesh);
 805         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 806         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 807         break;
 808       }
 809     case 3:
 810       {
 811         MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
 812         MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(target_mesh);
 813         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 814         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 815         break;
 816       }
 817     default:
 818       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only dimension 1 2 or 3 supported !");
 819   }
 820   nullifiedTinyCoeffInCrudeMatrixAbs(0.);
 821   //
 822   _deno_multiply.clear();
 823   _deno_multiply.resize(_matrix.size());
 824   _deno_reverse_multiply.clear();
 825   _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
 826   declareAsNew();
 827   return 1;
 828 }
 829
 830 int MEDCouplingRemapper::prepareInterpKernelOnlyCC()
 831 {
 832   std::string srcMeth,trgMeth;
 833   std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
 834   if(methodCpp!="P0P0")
 835     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only P0P0 interpolation supported for the moment !");
 836   if(InterpolationOptions::getIntersectionType()!=INTERP_KERNEL::Triangulation)
 837     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC: only 'Triangulation' intersection type supported!");
 838   const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
 839   const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
 840   const int srcMeshDim=src_mesh->getMeshDimension();
 841   const int trgMeshDim=target_mesh->getMeshDimension();
 842   if(trgMeshDim!=srcMeshDim)
 843     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : dimension of target cartesian mesh should be equal to dimension of source cartesian mesh !");
 844   switch(srcMeshDim)
 845   {
 846     case 1:
 847       {
 848         MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
 849         MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
 850         INTERP_KERNEL::InterpolationCC myInterpolator(*this);
 851         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 852         break;
 853       }
 854     case 2:
 855       {
 856         MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
 857         MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
 858         INTERP_KERNEL::InterpolationCC myInterpolator(*this);
 859         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 860         break;
 861       }
 862     case 3:
 863       {
 864         MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
 865         MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
 866         INTERP_KERNEL::InterpolationCC myInterpolator(*this);
 867         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 868         break;
 869       }
 870     default:
 871       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only dimension 1 2 or 3 supported !");
 872   }
 873   nullifiedTinyCoeffInCrudeMatrixAbs(0.);
 874   //
 875   _deno_multiply.clear();
 876   _deno_multiply.resize(_matrix.size());
 877   _deno_reverse_multiply.clear();
 878   _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
 879   declareAsNew();
 880   return 1;
 881 }
 882
 883 int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss()
 884 {
 885   if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
 886     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 !");
 887   MCAuto<DataArrayDouble> trgLoc=_target_ft->getLocalizationOfDiscr();
 888   const double *trgLocPtr=trgLoc->begin();
 889   mcIdType trgSpaceDim=ToIdType(trgLoc->getNumberOfComponents());
 890   MCAuto<DataArrayIdType> srcOffsetArr=_src_ft->getDiscretization()->getOffsetArr(_src_ft->getMesh());
 891   if(trgSpaceDim!=_src_ft->getMesh()->getSpaceDimension())
 892     {
 893       std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss : space dimensions mismatch between source and target !";
 894       oss << " Target discretization localization has dimension " << trgSpaceDim << ", whereas the space dimension of source is equal to ";
 895       oss << _src_ft->getMesh()->getSpaceDimension() << " !";
 896       throw INTERP_KERNEL::Exception(oss.str().c_str());
 897     }
 898   const mcIdType *srcOffsetArrPtr=srcOffsetArr->begin();
 899   MCAuto<DataArrayDouble> srcLoc=_src_ft->getLocalizationOfDiscr();
 900   const double *srcLocPtr=srcLoc->begin();
 901   MCAuto<DataArrayIdType> eltsArr,eltsIndexArr;
 902   mcIdType trgNbOfGaussPts=trgLoc->getNumberOfTuples();
 903   _matrix.resize(trgNbOfGaussPts);
 904   _src_ft->getMesh()->getCellsContainingPointsLinearPartOnlyOnNonDynType(trgLoc->begin(),trgNbOfGaussPts,getPrecision(),eltsArr,eltsIndexArr);
 905   const mcIdType *elts(eltsArr->begin()),*eltsIndex(eltsIndexArr->begin());
 906   MCAuto<DataArrayIdType> nbOfSrcCellsShTrgPts(eltsIndexArr->deltaShiftIndex());
 907   MCAuto<DataArrayIdType> ids0=nbOfSrcCellsShTrgPts->findIdsNotEqual(0);
 908   for(const mcIdType *trgId=ids0->begin();trgId!=ids0->end();trgId++)
 909     {
 910       const double *ptTrg=trgLocPtr+trgSpaceDim*(*trgId);
 911       mcIdType srcCellId=elts[eltsIndex[*trgId]];
 912       double dist=std::numeric_limits<double>::max();
 913       mcIdType srcEntry=-1;
 914       for(mcIdType srcId=srcOffsetArrPtr[srcCellId];srcId<srcOffsetArrPtr[srcCellId+1];srcId++)
 915         {
 916           const double *ptSrc=srcLocPtr+trgSpaceDim*srcId;
 917           double tmp=0.;
 918           for(mcIdType i=0;i<trgSpaceDim;i++)
 919             tmp+=(ptTrg[i]-ptSrc[i])*(ptTrg[i]-ptSrc[i]);
 920           if(tmp<dist)
 921             { dist=tmp; srcEntry=srcId; }
 922         }
 923       _matrix[*trgId][srcEntry]=1.;
 924     }
 925   if(ids0->getNumberOfTuples()!=trgNbOfGaussPts)
 926     {
 927       MCAuto<DataArrayIdType> orphanTrgIds=nbOfSrcCellsShTrgPts->findIdsEqual(0);
 928       MCAuto<DataArrayDouble> orphanTrg=trgLoc->selectByTupleId(orphanTrgIds->begin(),orphanTrgIds->end());
 929       MCAuto<DataArrayIdType> srcIdPerTrg=srcLoc->findClosestTupleId(orphanTrg);
 930       const mcIdType *srcIdPerTrgPtr=srcIdPerTrg->begin();
 931       for(const mcIdType *orphanTrgId=orphanTrgIds->begin();orphanTrgId!=orphanTrgIds->end();orphanTrgId++,srcIdPerTrgPtr++)
 932         _matrix[*orphanTrgId][*srcIdPerTrgPtr]=2.;
 933     }
 934   _deno_multiply.clear();
 935   _deno_multiply.resize(_matrix.size());
 936   _deno_reverse_multiply.clear();
 937   _deno_reverse_multiply.resize(srcLoc->getNumberOfTuples());
 938   declareAsNew();
 939   return 1;
 940 }
 941
 942 /*!
 943  * This method checks that the input interpolation \a method is managed by not INTERP_KERNEL only methods.
 944  * If no an INTERP_KERNEL::Exception will be thrown. If yes, a magic number will be returned to switch in the MEDCouplingRemapper::prepareNotInterpKernelOnly method.
 945  */
 946 int MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel(const std::string& method)
 947 {
 948   if(method=="GAUSSGAUSS")
 949     return 0;
 950   std::ostringstream oss; oss << "MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel : ";
 951   oss << "The method \"" << method << "\" is not manageable by not INTERP_KERNEL only method.";
 952   oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS !";
 953   throw INTERP_KERNEL::Exception(oss.str().c_str());
 954 }
 955
 956 /*!
 957  * This method determines regarding \c _interp_matrix_pol attribute ( set by MEDCouplingRemapper::setInterpolationMatrixPolicy and by default equal
 958  * 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
 959  * only INTERP_KERNEL method should be applied.
 960  */
 961 bool MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly() const
 962 {
 963   std::string srcm,trgm,method;
 964   method=checkAndGiveInterpolationMethodStr(srcm,trgm);
 965   switch(_interp_matrix_pol)
 966   {
 967     case IK_ONLY_PREFERED:
 968       {
 969         try
 970         {
 971             std::string tmp1,tmp2;
 972             INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,tmp1,tmp2);
 973             return true;
 974         }
 975         catch(INTERP_KERNEL::Exception& /*e*/)
 976         {
 977             return false;
 978         }
 979       }
 980     case NOT_IK_ONLY_PREFERED:
 981       {
 982         try
 983         {
 984             CheckInterpolationMethodManageableByNotOnlyInterpKernel(method);
 985             return false;
 986         }
 987         catch(INTERP_KERNEL::Exception& /*e*/)
 988         {
 989             return true;
 990         }
 991       }
 992     case IK_ONLY_FORCED:
 993       return true;
 994     case NOT_IK_ONLY_FORCED:
 995       return false;
 996     default:
 997       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly : internal error ! The interpolation matrix policy is not managed ! Try to change it using MEDCouplingRemapper::setInterpolationMatrixPolicy !");
 998   }
 999 }
1000
1001 void MEDCouplingRemapper::updateTime() const
1002 {
1003 }
1004
1005 void MEDCouplingRemapper::checkPrepare() const
1006 {
1007   const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
1008   if(!s || !t)
1009     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that MEDCouplingRemapper::prepare(Ex) has not been called !");
1010   if(!s->getMesh() || !t->getMesh())
1011     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that no all field templates have their mesh set !");
1012 }
1013
1014 /*!
1015  * This method builds a code considering already set field discretization int \a this : \a _src_ft and \a _target_ft.
1016  * This method returns 3 information (2 in output parameters and 1 in return).
1017  *
1018  * \param [out] srcMeth the string code of the discretization of source field template
1019  * \param [out] trgMeth the string code of the discretization of target field template
1020  * \return the standardized string code (compatible with INTERP_KERNEL) for matrix of numerators (in \a _matrix)
1021  */
1022 std::string MEDCouplingRemapper::checkAndGiveInterpolationMethodStr(std::string& srcMeth, std::string& trgMeth) const
1023 {
1024   const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
1025   if(!s || !t)
1026     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have been set !");
1027   if(!s->getMesh() || !t->getMesh())
1028     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have their mesh set !");
1029   srcMeth=_src_ft->getDiscretization()->getRepr();
1030   trgMeth=_target_ft->getDiscretization()->getRepr();
1031   return BuildMethodFrom(srcMeth,trgMeth);
1032 }
1033
1034 std::string MEDCouplingRemapper::BuildMethodFrom(const std::string& meth1, const std::string& meth2)
1035 {
1036   std::string method(meth1); method+=meth2;
1037   return method;
1038 }
1039
1040 void MEDCouplingRemapper::BuildFieldTemplatesFrom(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method, MCAuto<MEDCouplingFieldTemplate>& src, MCAuto<MEDCouplingFieldTemplate>& target)
1041 {
1042   if(!srcMesh || !targetMesh)
1043     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::BuildFieldTemplatesFrom : presence of NULL input pointer !");
1044   std::string srcMethod,targetMethod;
1045   INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,srcMethod,targetMethod);
1046   src=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(srcMethod));
1047   src->setMesh(srcMesh);
1048   target=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(targetMethod));
1049   target->setMesh(targetMesh);
1050 }
1051
1052 void MEDCouplingRemapper::releaseData(bool matrixSuppression)
1053 {
1054   _src_ft=0;
1055   _target_ft=0;
1056   if(matrixSuppression)
1057     {
1058       _matrix.clear();
1059       _deno_multiply.clear();
1060       _deno_reverse_multiply.clear();
1061     }
1062 }
1063
1064 void MEDCouplingRemapper::restartUsing(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target)
1065 {
1066   if(!src || !target)
1067     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::restartUsingData : presence of NULL input pointer !");
1068   if(!src->getMesh() || !target->getMesh())
1069     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareEx : presence of NULL mesh pointer in given field template !");
1070   releaseData(true);
1071   _src_ft.takeRef(const_cast<MEDCouplingFieldTemplate *>(src));
1072   _target_ft.takeRef(const_cast<MEDCouplingFieldTemplate *>(target));
1073 }
1074
1075 void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue)
1076 {
1077   if(!srcField || !targetField)
1078     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferUnderground : srcField or targetField is NULL !");
1079   srcField->checkConsistencyLight();
1080   checkPrepare();
1081   if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
1082     throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
1083   if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
1084     throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
1085   if(srcField->getNature()!=targetField->getNature())
1086     throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
1087   if(srcField->getNumberOfTuplesExpected()!=_src_ft->getNumberOfTuplesExpected())
1088     {
1089       std::ostringstream oss;
1090       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();
1091       oss << " ! It appears that the source support is not the same between the prepare and the transfer !";
1092       throw INTERP_KERNEL::Exception(oss.str().c_str());
1093     }
1094   DataArrayDouble *array(targetField->getArray());
1095   std::size_t srcNbOfCompo(srcField->getNumberOfComponents());
1096   if(array)
1097     {
1098       targetField->checkConsistencyLight();
1099       if(srcNbOfCompo!=targetField->getNumberOfComponents())
1100         throw INTERP_KERNEL::Exception("Number of components mismatch !");
1101     }
1102   else
1103     {
1104       if(!isDftVal)
1105         throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !");
1106       MCAuto<DataArrayDouble> tmp(DataArrayDouble::New());
1107       tmp->alloc(targetField->getNumberOfTuples(),srcNbOfCompo);
1108       targetField->setArray(tmp);
1109     }
1110   computeDeno(srcField->getNature(),srcField,targetField);
1111   double *resPointer(targetField->getArray()->getPointer());
1112   const double *inputPointer(srcField->getArray()->getConstPointer());
1113   computeProduct(inputPointer,(int)srcNbOfCompo,isDftVal,dftValue,resPointer);
1114 }
1115
1116 void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1117 {
1118   if(nat==NoNature)
1119     return computeDenoFromScratch(nat,srcField,trgField);
1120   else if(nat!=_nature_of_deno)
1121     return computeDenoFromScratch(nat,srcField,trgField);
1122   else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis())
1123     return computeDenoFromScratch(nat,srcField,trgField);
1124 }
1125
1126 void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1127 {
1128   _nature_of_deno=nat;
1129   std::size_t _time_deno_update=getTimeOfThis();
1130   switch(_nature_of_deno)
1131   {
1132     case IntensiveMaximum:
1133       {
1134         ComputeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1135         break;
1136       }
1137     case ExtensiveMaximum:
1138       {
1139         MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1140         MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1141         const double *denoPtr=deno->getArray()->getConstPointer();
1142         const double *denoRPtr=denoR->getArray()->getConstPointer();
1143         if(trgField->getMesh()->getMeshDimension()==-1)
1144           {
1145             double *denoRPtr2=denoR->getArray()->getPointer();
1146             denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1147           }
1148         if(srcField->getMesh()->getMeshDimension()==-1)
1149           {
1150             double *denoPtr2=deno->getArray()->getPointer();
1151             denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1152           }
1153         int idx=0;
1154         for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1155           for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1156             {
1157               _deno_multiply[idx][(*iter2).first]=denoPtr[(*iter2).first];
1158               _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[idx];
1159             }
1160         deno->decrRef();
1161         denoR->decrRef();
1162         break;
1163       }
1164     case ExtensiveConservation:
1165       {
1166         ComputeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1167         break;
1168       }
1169     case IntensiveConservation:
1170       {
1171         MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1172         MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1173         const double *denoPtr=deno->getArray()->getConstPointer();
1174         const double *denoRPtr=denoR->getArray()->getConstPointer();
1175         if(trgField->getMesh()->getMeshDimension()==-1)
1176           {
1177             double *denoRPtr2=denoR->getArray()->getPointer();
1178             denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1179           }
1180         if(srcField->getMesh()->getMeshDimension()==-1)
1181           {
1182             double *denoPtr2=deno->getArray()->getPointer();
1183             denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1184           }
1185         mcIdType idx=0;
1186         for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1187           for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1188             {
1189               _deno_multiply[idx][(*iter2).first]=denoPtr[idx];
1190               _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[(*iter2).first];
1191             }
1192         deno->decrRef();
1193         denoR->decrRef();
1194         break;
1195       }
1196     case NoNature:
1197       throw INTERP_KERNEL::Exception("No nature specified ! Select one !");
1198   }
1199 }
1200
1201 void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)
1202 {
1203   int idx=0;
1204   double *tmp=new double[inputNbOfCompo];
1205   for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1206     {
1207       if((*iter1).empty())
1208         {
1209           if(isDftVal)
1210             std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1211           continue;
1212         }
1213       else
1214         std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,0.);
1215       std::map<mcIdType,double>::const_iterator iter3=_deno_multiply[idx].begin();
1216       for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter3++)
1217         {
1218           std::transform(inputPointer+(*iter2).first*inputNbOfCompo,inputPointer+((*iter2).first+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/(*iter3).second));
1219           std::transform(tmp,tmp+inputNbOfCompo,resPointer+idx*inputNbOfCompo,resPointer+idx*inputNbOfCompo,std::plus<double>());
1220         }
1221     }
1222   delete [] tmp;
1223 }
1224
1225 void MEDCouplingRemapper::computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer)
1226 {
1227   std::vector<bool> isReached(_deno_reverse_multiply.size(),false);
1228   mcIdType idx=0;
1229   double *tmp=new double[inputNbOfCompo];
1230   std::fill(resPointer,resPointer+inputNbOfCompo*_deno_reverse_multiply.size(),0.);
1231   for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1232     {
1233       for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1234         {
1235           isReached[(*iter2).first]=true;
1236           std::transform(inputPointer+idx*inputNbOfCompo,inputPointer+(idx+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/_deno_reverse_multiply[(*iter2).first][idx]));
1237           std::transform(tmp,tmp+inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,std::plus<double>());
1238         }
1239     }
1240   delete [] tmp;
1241   idx=0;
1242   for(std::vector<bool>::const_iterator iter3=isReached.begin();iter3!=isReached.end();iter3++,idx++)
1243     if(!*iter3)
1244       std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1245 }
1246
1247 void MEDCouplingRemapper::ReverseMatrix(const std::vector<std::map<mcIdType,double> >& matIn, mcIdType nbColsMatIn, std::vector<std::map<mcIdType,double> >& matOut)
1248 {
1249   matOut.resize(nbColsMatIn);
1250   mcIdType id=0;
1251   for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
1252     for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1253       matOut[(*iter2).first][id]=(*iter2).second;
1254 }
1255
1256 void MEDCouplingRemapper::ComputeRowSumAndColSum(const std::vector<std::map<mcIdType,double> >& matrixDeno,
1257                                                  std::vector<std::map<mcIdType,double> >& deno, std::vector<std::map<mcIdType,double> >& denoReverse)
1258 {
1259   std::map<mcIdType,double> values;
1260   mcIdType idx=0;
1261   for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1262     {
1263       double sum=0.;
1264       for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1265         {
1266           sum+=(*iter2).second;
1267           values[(*iter2).first]+=(*iter2).second;
1268         }
1269       for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1270         deno[idx][(*iter2).first]=sum;
1271     }
1272   idx=0;
1273   for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1274     {
1275       for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1276         denoReverse[(*iter2).first][idx]=values[(*iter2).first];
1277     }
1278 }
1279
1280 void MEDCouplingRemapper::ComputeColSumAndRowSum(const std::vector<std::map<mcIdType,double> >& matrixDeno,
1281                                                  std::vector<std::map<mcIdType,double> >& deno, std::vector<std::map<mcIdType,double> >& denoReverse)
1282 {
1283   std::map<mcIdType,double> values;
1284   mcIdType idx=0;
1285   for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1286     {
1287       double sum=0.;
1288       for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1289         {
1290           sum+=(*iter2).second;
1291           values[(*iter2).first]+=(*iter2).second;
1292         }
1293       for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1294         denoReverse[(*iter2).first][idx]=sum;
1295     }
1296   idx=0;
1297   for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1298     {
1299       for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1300         deno[idx][(*iter2).first]=values[(*iter2).first];
1301     }
1302 }
1303
1304 void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<mcIdType,double> >& m1D,
1305                                                                      const std::vector< std::map<mcIdType,double> >& m2D,
1306                                                                      const mcIdType *corrCellIdSrc, mcIdType nbOf2DCellsSrc, mcIdType nbOf1DCellsSrc,
1307                                                                      const mcIdType *corrCellIdTrg)
1308 {
1309   mcIdType nbOf2DCellsTrg=ToIdType(m2D.size());
1310   mcIdType nbOf1DCellsTrg=ToIdType(m1D.size());
1311   mcIdType nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
1312   _matrix.resize(nbOf3DCellsTrg);
1313   mcIdType id2R=0;
1314   for(std::vector< std::map<mcIdType,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
1315     {
1316       for(std::map<mcIdType,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
1317         {
1318           mcIdType id1R=0;
1319           for(std::vector< std::map<mcIdType,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
1320             {
1321               for(std::map<mcIdType,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
1322                 {
1323                   _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
1324                 }
1325             }
1326         }
1327     }
1328 }
1329
1330 void MEDCouplingRemapper::PrintMatrix(const std::vector<std::map<mcIdType,double> >& m)
1331 {
1332   mcIdType id=0;
1333   for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++)
1334     {
1335       std::cout << "Target Cell # " << id << " : ";
1336       for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1337         std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), ";
1338       std::cout << std::endl;
1339     }
1340 }
1341
1342 const std::vector<std::map<mcIdType,double> >& MEDCouplingRemapper::getCrudeMatrix() const
1343 {
1344   return _matrix;
1345 }
1346
1347 /*!
1348  * Returns the number of columns of matrix returned by MEDCouplingRemapper::getCrudeMatrix method.
1349  */
1350 mcIdType MEDCouplingRemapper::getNumberOfColsOfMatrix() const
1351 {
1352   return ToIdType(_deno_reverse_multiply.size());
1353 }
1354
1355 /*!
1356  * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1357  * If not the behaviour is unpredictable.
1358  * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than \a maxValAbs this coefficient is
1359  * 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.
1360  * 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
1361  * can occur.
1362  *
1363  * \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.
1364  * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged.
1365  * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix
1366  */
1367 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs(double maxValAbs)
1368 {
1369   int ret=0;
1370   std::vector<std::map<mcIdType,double> > matrixNew(_matrix.size());
1371   mcIdType i=0;
1372   for(std::vector<std::map<mcIdType,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++,i++)
1373     {
1374       std::map<mcIdType,double>& rowNew=matrixNew[i];
1375       for(std::map<mcIdType,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1376         {
1377           if(fabs((*it2).second)>maxValAbs)
1378             rowNew[(*it2).first]=(*it2).second;
1379           else
1380             ret++;
1381         }
1382     }
1383   if(ret>0)
1384     _matrix=matrixNew;
1385   return ret;
1386 }
1387
1388 /*!
1389  * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1390  * If not the behaviour is unpredictable.
1391  * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than delta multiplied by \a scaleFactor this coefficient is
1392  * 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.
1393  * delta is the value returned by MEDCouplingRemapper::getMaxValueInCrudeMatrix method.
1394  * 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
1395  * can occur.
1396  *
1397  * \param [in] scaleFactor is the scale factor from which coefficients lower than \a scaleFactor times range width of coefficients are set to zero.
1398  * \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
1399  *         that all coefficients are null.
1400  * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs
1401  */
1402 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix(double scaleFactor)
1403 {
1404   double maxVal=getMaxValueInCrudeMatrix();
1405   if(maxVal==0.)
1406     return -1;
1407   return nullifiedTinyCoeffInCrudeMatrixAbs(scaleFactor*maxVal);
1408 }
1409
1410 /*!
1411  * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1412  * If not the behaviour is unpredictable.
1413  * This method returns the maximum of the absolute values of coefficients into the sparse crude matrix.
1414  * The returned value is positive.
1415  */
1416 double MEDCouplingRemapper::getMaxValueInCrudeMatrix() const
1417 {
1418   double ret=0.;
1419   for(std::vector<std::map<mcIdType,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++)
1420     for(std::map<mcIdType,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1421       if(fabs((*it2).second)>ret)
1422         ret=fabs((*it2).second);
1423   return ret;
1424 }