src/MEDCoupling/MEDCouplingRemapper.cxx

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