src/MEDCoupling/MEDCouplingRemapper.cxx

   1 // Copyright (C) 2007-2022  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 #include "MEDCouplingFieldDiscretizationOnNodesFE.hxx"
  31
  32 #include "Interpolation1D.txx"
  33 #include "Interpolation2DCurve.hxx"
  34 #include "Interpolation2D.txx"
  35 #include "Interpolation3D.txx"
  36 #include "Interpolation3DSurf.hxx"
  37 #include "Interpolation2D1D.txx"
  38 #include "Interpolation2D3D.txx"
  39 #include "Interpolation3D1D.txx"
  40 #include "Interpolation1D0D.txx"
  41 #include "InterpolationCU.txx"
  42 #include "InterpolationCC.txx"
  43
  44 using namespace MEDCoupling;
  45
  46 MEDCouplingRemapper::MEDCouplingRemapper():_src_ft(0),_target_ft(0),_interp_matrix_pol(IK_ONLY_PREFERED),_nature_of_deno(NoNature),_time_deno_update(0)
  47 {
  48 }
  49
  50 MEDCouplingRemapper::~MEDCouplingRemapper()
  51 {
  52   releaseData(false);
  53 }
  54
  55 /*!
  56  * This method is the second step of the remapping process. The remapping process works in three phases :
  57  *
  58  * - Set remapping options appropriately
  59  * - The computation of remapping matrix
  60  * - Apply the matrix vector multiply to obtain the result of the remapping
  61  *
  62  * 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.
  63  * 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).
  64  * So this method is less precise but a more user friendly way to compute a remapping matrix.
  65  *
  66  * \param [in] srcMesh the source mesh
  67  * \param [in] targetMesh the target mesh
  68  * \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.
  69  *             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).
  70  *
  71  * \sa MEDCouplingRemapper::prepareEx
  72  */
  73 int MEDCouplingRemapper::prepare(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method)
  74 {
  75   MCAuto<MEDCouplingFieldTemplate> src,target;
  76   BuildFieldTemplatesFrom(srcMesh,targetMesh,method,src,target);
  77   return prepareEx(src,target);
  78 }
  79
  80 /*!
  81  * This method is the generalization of MEDCouplingRemapper::prepare. Indeed, MEDCouplingFieldTemplate instances gives all required information to compute the remapping matrix.
  82  * This method must be used instead of MEDCouplingRemapper::prepare if Gauss point to Gauss point must be applied.
  83  *
  84  * \param [in] src is the field template source side.
  85  * \param [in] target is the field template target side.
  86  *
  87  * \sa MEDCouplingRemapper::prepare
  88  */
  89 int MEDCouplingRemapper::prepareEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target)
  90 {
  91   restartUsing(src,target);
  92   if(isInterpKernelOnlyOrNotOnly())
  93     return prepareInterpKernelOnly();
  94   else
  95     return prepareNotInterpKernelOnly();
  96 }
  97
  98 void MEDCouplingRemapper::setCrudeMatrix(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method, const std::vector<std::map<mcIdType,double> >& m)
  99 {
 100   MCAuto<MEDCouplingFieldTemplate> src,target;
 101   BuildFieldTemplatesFrom(srcMesh,targetMesh,method,src,target);
 102   setCrudeMatrixEx(src,target,m);
 103 }
 104
 105 void MEDCouplingRemapper::setCrudeMatrixEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target, const std::vector<std::map<mcIdType,double> >& m)
 106 {
 107   restartUsing(src,target);
 108   if(ToIdType(m.size())!=target->getNumberOfTuplesExpected())
 109     {
 110       std::ostringstream oss; oss << "MEDCouplingRemapper::setMatrixEx : input matrix has " << m.size() << " rows whereas there are " << target->getNumberOfTuplesExpected() << " expected !";
 111       throw INTERP_KERNEL::Exception(oss.str());
 112     }
 113   auto srcNbElem(src->getNumberOfTuplesExpected());
 114   for(auto it: m)
 115     {
 116       for(auto it2: it)
 117         {
 118           auto idToTest(it2.first);
 119           if(idToTest<0 || idToTest>=srcNbElem)
 120             {
 121               std::ostringstream oss; oss << "MEDCouplingRemapper::setMatrixEx : presence of elt #" << idToTest << " ! not in [0," << srcNbElem << ") !";
 122               throw INTERP_KERNEL::Exception(oss.str());
 123             }
 124         }
 125     }
 126   _matrix=m;
 127   synchronizeSizeOfSideMatricesAfterMatrixComputation(srcNbElem);
 128 }
 129
 130 int MEDCouplingRemapper::prepareInterpKernelOnly()
 131 {
 132   int meshInterpType=((int)_src_ft->getMesh()->getType()*16)+(int)_target_ft->getMesh()->getType();
 133   // *** Remember:
 134 //  typedef enum
 135 //    {
 136 //      UNSTRUCTURED = 5,
 137 //      CARTESIAN = 7,
 138 //      EXTRUDED = 8,
 139 //      CURVE_LINEAR = 9,
 140 //      SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED = 10,
 141 //      SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED = 11,
 142 //      IMAGE_GRID = 12
 143 //    } MEDCouplingMeshType;
 144
 145   switch(meshInterpType)
 146   {
 147     case 90:   // UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
 148     case 91:   // UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
 149     case 165:  // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - UNSTRUCTURED
 150     case 181:  // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - UNSTRUCTURED
 151     case 170:  // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
 152     case 171:  // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
 153     case 186:  // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
 154     case 187:  // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
 155     case 85:   // UNSTRUCTURED - UNSTRUCTURED
 156       return prepareInterpKernelOnlyUU();
 157     case 167:  // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - CARTESIAN
 158     case 183:  // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - CARTESIAN
 159     case 87:   // UNSTRUCTURED - CARTESIAN
 160       return prepareInterpKernelOnlyUC();
 161     case 122:  // CARTESIAN - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
 162     case 123:  // CARTESIAN - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
 163     case 117:  // CARTESIAN - UNSTRUCTURED
 164       return prepareInterpKernelOnlyCU();
 165     case 119:  // CARTESIAN - CARTESIAN
 166       return prepareInterpKernelOnlyCC();
 167     case 136:  // EXTRUDED - EXTRUDED
 168       return prepareInterpKernelOnlyEE();
 169     default:
 170       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnly : Not managed type of meshes ! Dealt meshes type are : Unstructured<->Unstructured, Unstructured<->Cartesian, Cartesian<->Cartesian, Extruded<->Extruded !");
 171   }
 172 }
 173
 174 int MEDCouplingRemapper::prepareNotInterpKernelOnly()
 175 {
 176   std::string srcm,trgm,method;
 177   method=checkAndGiveInterpolationMethodStr(srcm,trgm);
 178   switch(CheckInterpolationMethodManageableByNotOnlyInterpKernel(method))
 179   {
 180     case 0:
 181       return prepareNotInterpKernelOnlyGaussGauss();
 182     case 1:
 183       return prepareNotInterpKernelOnlyFEFE();
 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(trgNbOfCompo != srcField->getNumberOfComponents())
 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   synchronizeSizeOfSideMatricesAfterMatrixComputation(nbCols);
 693   return 1;
 694 }
 695
 696 int MEDCouplingRemapper::prepareInterpKernelOnlyEE()
 697 {
 698   std::string srcMeth,trgMeth;
 699   std::string methC=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
 700   const MEDCouplingMappedExtrudedMesh *src_mesh=static_cast<const MEDCouplingMappedExtrudedMesh *>(_src_ft->getMesh());
 701   const MEDCouplingMappedExtrudedMesh *target_mesh=static_cast<const MEDCouplingMappedExtrudedMesh *>(_target_ft->getMesh());
 702   if(methC!="P0P0")
 703     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyEE : Only P0P0 method implemented for Extruded/Extruded meshes !");
 704   MCAuto<MEDCouplingUMesh> src2D(src_mesh->getMesh2D()->clone(false)); src2D->changeSpaceDimension(2,0.);
 705   MCAuto<MEDCouplingUMesh> trg2D(target_mesh->getMesh2D()->clone(false)); trg2D->changeSpaceDimension(2,0.);
 706   MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src2D);
 707   MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(trg2D);
 708   INTERP_KERNEL::Interpolation2D interpolation2D(*this);
 709   std::vector<std::map<mcIdType,double> > matrix2D;
 710   mcIdType nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,methC);
 711   MEDCouplingUMesh *s1D,*t1D;
 712   double v[3];
 713   MEDCouplingMappedExtrudedMesh::Project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v);
 714   MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D);
 715   MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
 716   std::vector<std::map<mcIdType,double> > matrix1D;
 717   INTERP_KERNEL::Interpolation1D interpolation1D(*this);
 718   if(interpolation1D.getIntersectionType()==INTERP_KERNEL::Geometric2D)// For intersection type of 1D, Geometric2D do not deal with it ! -> make interpolation1D not inherite from this
 719     interpolation1D.setIntersectionType(INTERP_KERNEL::Triangulation);//
 720     mcIdType nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,methC);
 721   s1D->decrRef();
 722   t1D->decrRef();
 723   buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
 724                                              target_mesh->getMesh3DIds()->getConstPointer());
 725   //
 726   synchronizeSizeOfSideMatricesAfterMatrixComputation(nbCols2D*nbCols1D);
 727   return 1;
 728 }
 729
 730 int MEDCouplingRemapper::prepareInterpKernelOnlyUC()
 731 {
 732   std::string srcMeth,trgMeth;
 733   std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
 734   if(methodCpp!="P0P0")
 735     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC: only P0P0 interpolation supported for the moment !");
 736   if(InterpolationOptions::getIntersectionType()!=INTERP_KERNEL::Triangulation)
 737       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC: only 'Triangulation' intersection type supported!");
 738   const MEDCouplingUMesh *src_mesh=static_cast<const MEDCouplingUMesh *>(_src_ft->getMesh());
 739   const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
 740   const int srcMeshDim=src_mesh->getMeshDimension();
 741   const int srcSpceDim=src_mesh->getSpaceDimension();
 742   const int trgMeshDim=target_mesh->getMeshDimension();
 743   if(srcMeshDim!=srcSpceDim || srcMeshDim!=trgMeshDim)
 744     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!");
 745   std::vector<std::map<mcIdType,double> > res;
 746   switch(srcMeshDim)
 747   {
 748     case 1:
 749       {
 750         MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
 751         MEDCouplingNormalizedUnstructuredMesh<1,1> sourceWrapper(src_mesh);
 752         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 753         myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
 754         break;
 755       }
 756     case 2:
 757       {
 758         MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
 759         MEDCouplingNormalizedUnstructuredMesh<2,2> sourceWrapper(src_mesh);
 760         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 761         myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
 762         break;
 763       }
 764     case 3:
 765       {
 766         MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
 767         MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(src_mesh);
 768         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 769         myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
 770         break;
 771       }
 772     default:
 773       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only dimension 1 2 or 3 supported !");
 774   }
 775   ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix);
 776   nullifiedTinyCoeffInCrudeMatrixAbs(0.);
 777   //
 778   synchronizeSizeOfSideMatricesAfterMatrixComputation(src_mesh->getNumberOfCells());
 779   return 1;
 780 }
 781
 782 int MEDCouplingRemapper::prepareInterpKernelOnlyCU()
 783 {
 784   std::string srcMeth,trgMeth;
 785   std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
 786   if(methodCpp!="P0P0")
 787     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only P0P0 interpolation supported for the moment !");
 788   if(InterpolationOptions::getIntersectionType()!=INTERP_KERNEL::Triangulation)
 789     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU: only 'Triangulation' intersection type supported!");
 790   const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
 791   const MEDCouplingUMesh *target_mesh=static_cast<const MEDCouplingUMesh *>(_target_ft->getMesh());
 792   const int srcMeshDim=src_mesh->getMeshDimension();
 793   const int trgMeshDim=target_mesh->getMeshDimension();
 794   const int trgSpceDim=target_mesh->getSpaceDimension();
 795   if(trgMeshDim!=trgSpceDim || trgMeshDim!=srcMeshDim)
 796     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!");
 797   switch(srcMeshDim)
 798   {
 799     case 1:
 800       {
 801         MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
 802         MEDCouplingNormalizedUnstructuredMesh<1,1> targetWrapper(target_mesh);
 803         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 804         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 805         break;
 806       }
 807     case 2:
 808       {
 809         MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
 810         MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(target_mesh);
 811         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 812         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 813         break;
 814       }
 815     case 3:
 816       {
 817         MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
 818         MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(target_mesh);
 819         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 820         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 821         break;
 822       }
 823     default:
 824       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only dimension 1 2 or 3 supported !");
 825   }
 826   nullifiedTinyCoeffInCrudeMatrixAbs(0.);
 827   //
 828   synchronizeSizeOfSideMatricesAfterMatrixComputation(src_mesh->getNumberOfCells());
 829   return 1;
 830 }
 831
 832 int MEDCouplingRemapper::prepareInterpKernelOnlyCC()
 833 {
 834   std::string srcMeth,trgMeth;
 835   std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
 836   if(methodCpp!="P0P0")
 837     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only P0P0 interpolation supported for the moment !");
 838   if(InterpolationOptions::getIntersectionType()!=INTERP_KERNEL::Triangulation)
 839     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC: only 'Triangulation' intersection type supported!");
 840   const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
 841   const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
 842   const int srcMeshDim=src_mesh->getMeshDimension();
 843   const int trgMeshDim=target_mesh->getMeshDimension();
 844   if(trgMeshDim!=srcMeshDim)
 845     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : dimension of target cartesian mesh should be equal to dimension of source cartesian mesh !");
 846   switch(srcMeshDim)
 847   {
 848     case 1:
 849       {
 850         MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
 851         MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
 852         INTERP_KERNEL::InterpolationCC myInterpolator(*this);
 853         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 854         break;
 855       }
 856     case 2:
 857       {
 858         MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
 859         MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
 860         INTERP_KERNEL::InterpolationCC myInterpolator(*this);
 861         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 862         break;
 863       }
 864     case 3:
 865       {
 866         MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
 867         MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
 868         INTERP_KERNEL::InterpolationCC myInterpolator(*this);
 869         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 870         break;
 871       }
 872     default:
 873       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only dimension 1 2 or 3 supported !");
 874   }
 875   nullifiedTinyCoeffInCrudeMatrixAbs(0.);
 876   //
 877   synchronizeSizeOfSideMatricesAfterMatrixComputation(src_mesh->getNumberOfCells());
 878   return 1;
 879 }
 880
 881 int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss()
 882 {
 883   if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
 884     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 !");
 885   MCAuto<DataArrayDouble> trgLoc=_target_ft->getLocalizationOfDiscr();
 886   const double *trgLocPtr=trgLoc->begin();
 887   mcIdType trgSpaceDim=ToIdType(trgLoc->getNumberOfComponents());
 888   MCAuto<DataArrayIdType> srcOffsetArr=_src_ft->getDiscretization()->getOffsetArr(_src_ft->getMesh());
 889   if(trgSpaceDim!=_src_ft->getMesh()->getSpaceDimension())
 890     {
 891       std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss : space dimensions mismatch between source and target !";
 892       oss << " Target discretization localization has dimension " << trgSpaceDim << ", whereas the space dimension of source is equal to ";
 893       oss << _src_ft->getMesh()->getSpaceDimension() << " !";
 894       throw INTERP_KERNEL::Exception(oss.str().c_str());
 895     }
 896   const mcIdType *srcOffsetArrPtr=srcOffsetArr->begin();
 897   MCAuto<DataArrayDouble> srcLoc=_src_ft->getLocalizationOfDiscr();
 898   const double *srcLocPtr=srcLoc->begin();
 899   MCAuto<DataArrayIdType> eltsArr,eltsIndexArr;
 900   mcIdType trgNbOfGaussPts=trgLoc->getNumberOfTuples();
 901   _matrix.resize(trgNbOfGaussPts);
 902   _src_ft->getMesh()->getCellsContainingPointsLinearPartOnlyOnNonDynType(trgLoc->begin(),trgNbOfGaussPts,getPrecision(),eltsArr,eltsIndexArr);
 903   const mcIdType *elts(eltsArr->begin()),*eltsIndex(eltsIndexArr->begin());
 904   MCAuto<DataArrayIdType> nbOfSrcCellsShTrgPts(eltsIndexArr->deltaShiftIndex());
 905   MCAuto<DataArrayIdType> ids0=nbOfSrcCellsShTrgPts->findIdsNotEqual(0);
 906   for(const mcIdType *trgId=ids0->begin();trgId!=ids0->end();trgId++)
 907     {
 908       const double *ptTrg=trgLocPtr+trgSpaceDim*(*trgId);
 909       mcIdType srcCellId=elts[eltsIndex[*trgId]];
 910       double dist=std::numeric_limits<double>::max();
 911       mcIdType srcEntry=-1;
 912       for(mcIdType srcId=srcOffsetArrPtr[srcCellId];srcId<srcOffsetArrPtr[srcCellId+1];srcId++)
 913         {
 914           const double *ptSrc=srcLocPtr+trgSpaceDim*srcId;
 915           double tmp=0.;
 916           for(mcIdType i=0;i<trgSpaceDim;i++)
 917             tmp+=(ptTrg[i]-ptSrc[i])*(ptTrg[i]-ptSrc[i]);
 918           if(tmp<dist)
 919             { dist=tmp; srcEntry=srcId; }
 920         }
 921       _matrix[*trgId][srcEntry]=1.;
 922     }
 923   if(ids0->getNumberOfTuples()!=trgNbOfGaussPts)
 924     {
 925       MCAuto<DataArrayIdType> orphanTrgIds=nbOfSrcCellsShTrgPts->findIdsEqual(0);
 926       MCAuto<DataArrayDouble> orphanTrg=trgLoc->selectByTupleId(orphanTrgIds->begin(),orphanTrgIds->end());
 927       MCAuto<DataArrayIdType> srcIdPerTrg=srcLoc->findClosestTupleId(orphanTrg);
 928       const mcIdType *srcIdPerTrgPtr=srcIdPerTrg->begin();
 929       for(const mcIdType *orphanTrgId=orphanTrgIds->begin();orphanTrgId!=orphanTrgIds->end();orphanTrgId++,srcIdPerTrgPtr++)
 930         _matrix[*orphanTrgId][*srcIdPerTrgPtr]=2.;
 931     }
 932   synchronizeSizeOfSideMatricesAfterMatrixComputation( srcLoc->getNumberOfTuples() );
 933   return 1;
 934 }
 935
 936 int MEDCouplingRemapper::prepareNotInterpKernelOnlyFEFE()
 937 {
 938   if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
 939     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareNotInterpKernelOnlyFEFE : The intersection type is not supported ! Only PointLocator is supported for FE->FE interpolation ! Please invoke setIntersectionType(PointLocator) on the MEDCouplingRemapper instance !");
 940   MCAuto<DataArrayDouble> trgLoc=_target_ft->getLocalizationOfDiscr();
 941   mcIdType trgSpaceDim=ToIdType(trgLoc->getNumberOfComponents());
 942   if(trgSpaceDim!=3)
 943     THROW_IK_EXCEPTION("prepareNotInterpKernelOnlyFEFE : only spacedim 3 supported for target !")
 944   if(_src_ft->getMesh()->getSpaceDimension() != 3)
 945     THROW_IK_EXCEPTION("prepareNotInterpKernelOnlyFEFE : only spacedim 3 supported for source !")
 946   const MEDCouplingUMesh *srcUMesh( dynamic_cast<const MEDCouplingUMesh *>(_src_ft->getMesh()) );
 947   const MEDCouplingPointSet *trgMesh( dynamic_cast<const MEDCouplingPointSet *>(_target_ft->getMesh()) );
 948   if( !srcUMesh )
 949     THROW_IK_EXCEPTION("prepareNotInterpKernelOnlyFEFE : only 3D UMesh supported as source !");
 950   if( !trgMesh )
 951     THROW_IK_EXCEPTION("prepareNotInterpKernelOnlyFEFE : only 3D PointSet mesh supported as target !");
 952
 953   _matrix.clear();
 954   _matrix.resize(trgMesh->getNumberOfNodes());
 955   mcIdType rowId(0);
 956
 957   auto matrixFeeder = [this,&rowId](const MEDCouplingGaussLocalization& gl, const std::vector<mcIdType>& conn)
 958   {
 959     auto& row = this->_matrix[rowId++];
 960     MCAuto<DataArrayDouble> resVector( gl.getShapeFunctionValues() );
 961     for(int iPt = 0 ; iPt < gl.getNumberOfPtsInRefCell(); ++iPt)
 962     {
 963       row[ conn[iPt] ] = resVector->getIJ(0,iPt);
 964     }
 965   };
 966
 967   MEDCouplingFieldDiscretizationOnNodesFE::GetRefCoordOfListOf3DPtsIn3D(srcUMesh,trgMesh->getCoords()->begin(),trgMesh->getNumberOfNodes(),matrixFeeder);
 968   synchronizeSizeOfSideMatricesAfterMatrixComputation( srcUMesh->getNumberOfNodes() );
 969   return 1;
 970 }
 971
 972 /*!
 973  * This method checks that the input interpolation \a method is managed by not INTERP_KERNEL only methods.
 974  * If no an INTERP_KERNEL::Exception will be thrown. If yes, a magic number will be returned to switch in the MEDCouplingRemapper::prepareNotInterpKernelOnly method.
 975  */
 976 int MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel(const std::string& method)
 977 {
 978   if(method=="GAUSSGAUSS")
 979     return 0;
 980   if(method=="FEFE")
 981     return 1;
 982   std::ostringstream oss; oss << "MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel : ";
 983   oss << "The method \"" << method << "\" is not manageable by not INTERP_KERNEL only method.";
 984   oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS FEFE !";
 985   throw INTERP_KERNEL::Exception(oss.str().c_str());
 986 }
 987
 988 /*!
 989  * This method determines regarding \c _interp_matrix_pol attribute ( set by MEDCouplingRemapper::setInterpolationMatrixPolicy and by default equal
 990  * 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
 991  * only INTERP_KERNEL method should be applied.
 992  */
 993 bool MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly() const
 994 {
 995   std::string srcm,trgm,method;
 996   method=checkAndGiveInterpolationMethodStr(srcm,trgm);
 997   switch(_interp_matrix_pol)
 998   {
 999     case IK_ONLY_PREFERED:
1000       {
1001         try
1002         {
1003             std::string tmp1,tmp2;
1004             INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,tmp1,tmp2);
1005             return true;
1006         }
1007         catch(INTERP_KERNEL::Exception& /*e*/)
1008         {
1009             return false;
1010         }
1011       }
1012     case NOT_IK_ONLY_PREFERED:
1013       {
1014         try
1015         {
1016             CheckInterpolationMethodManageableByNotOnlyInterpKernel(method);
1017             return false;
1018         }
1019         catch(INTERP_KERNEL::Exception& /*e*/)
1020         {
1021             return true;
1022         }
1023       }
1024     case IK_ONLY_FORCED:
1025       return true;
1026     case NOT_IK_ONLY_FORCED:
1027       return false;
1028     default:
1029       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly : internal error ! The interpolation matrix policy is not managed ! Try to change it using MEDCouplingRemapper::setInterpolationMatrixPolicy !");
1030   }
1031 }
1032
1033 void MEDCouplingRemapper::updateTime() const
1034 {
1035 }
1036
1037 void MEDCouplingRemapper::checkPrepare() const
1038 {
1039   const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
1040   if(!s || !t)
1041     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that MEDCouplingRemapper::prepare(Ex) has not been called !");
1042   if(!s->getMesh() || !t->getMesh())
1043     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that no all field templates have their mesh set !");
1044 }
1045
1046 void MEDCouplingRemapper::synchronizeSizeOfSideMatricesAfterMatrixComputation(mcIdType nbOfColsInMatrix)
1047 {
1048   _deno_multiply.clear();
1049   _deno_multiply.resize(_matrix.size());
1050   _deno_reverse_multiply.clear();
1051   _deno_reverse_multiply.resize(nbOfColsInMatrix);
1052   declareAsNew();
1053 }
1054
1055 /*!
1056  * This method builds a code considering already set field discretization int \a this : \a _src_ft and \a _target_ft.
1057  * This method returns 3 information (2 in output parameters and 1 in return).
1058  *
1059  * \param [out] srcMeth the string code of the discretization of source field template
1060  * \param [out] trgMeth the string code of the discretization of target field template
1061  * \return the standardized string code (compatible with INTERP_KERNEL) for matrix of numerators (in \a _matrix)
1062  */
1063 std::string MEDCouplingRemapper::checkAndGiveInterpolationMethodStr(std::string& srcMeth, std::string& trgMeth) const
1064 {
1065   const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
1066   if(!s || !t)
1067     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have been set !");
1068   if(!s->getMesh() || !t->getMesh())
1069     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have their mesh set !");
1070   srcMeth=_src_ft->getDiscretization()->getRepr();
1071   trgMeth=_target_ft->getDiscretization()->getRepr();
1072   return BuildMethodFrom(srcMeth,trgMeth);
1073 }
1074
1075 std::string MEDCouplingRemapper::BuildMethodFrom(const std::string& meth1, const std::string& meth2)
1076 {
1077   std::string method(meth1); method+=meth2;
1078   return method;
1079 }
1080
1081 void MEDCouplingRemapper::BuildFieldTemplatesFrom(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method, MCAuto<MEDCouplingFieldTemplate>& src, MCAuto<MEDCouplingFieldTemplate>& target)
1082 {
1083   if(!srcMesh || !targetMesh)
1084     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::BuildFieldTemplatesFrom : presence of NULL input pointer !");
1085   std::string srcMethod,targetMethod;
1086   INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,srcMethod,targetMethod);
1087   src=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(srcMethod));
1088   src->setMesh(srcMesh);
1089   target=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(targetMethod));
1090   target->setMesh(targetMesh);
1091 }
1092
1093 void MEDCouplingRemapper::releaseData(bool matrixSuppression)
1094 {
1095   _src_ft=0;
1096   _target_ft=0;
1097   if(matrixSuppression)
1098     {
1099       _matrix.clear();
1100       _deno_multiply.clear();
1101       _deno_reverse_multiply.clear();
1102     }
1103 }
1104
1105 void MEDCouplingRemapper::restartUsing(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target)
1106 {
1107   if(!src || !target)
1108     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::restartUsingData : presence of NULL input pointer !");
1109   if(!src->getMesh() || !target->getMesh())
1110     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareEx : presence of NULL mesh pointer in given field template !");
1111   releaseData(true);
1112   _src_ft.takeRef(const_cast<MEDCouplingFieldTemplate *>(src));
1113   _target_ft.takeRef(const_cast<MEDCouplingFieldTemplate *>(target));
1114 }
1115
1116 void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue)
1117 {
1118   if(!srcField || !targetField)
1119     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferUnderground : srcField or targetField is NULL !");
1120   srcField->checkConsistencyLight();
1121   checkPrepare();
1122   if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
1123     throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
1124   if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
1125     throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
1126   if(srcField->getNature()!=targetField->getNature())
1127     throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
1128   if(srcField->getNumberOfTuplesExpected()!=_src_ft->getNumberOfTuplesExpected())
1129     {
1130       std::ostringstream oss;
1131       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();
1132       oss << " ! It appears that the source support is not the same between the prepare and the transfer !";
1133       throw INTERP_KERNEL::Exception(oss.str().c_str());
1134     }
1135   DataArrayDouble *array(targetField->getArray());
1136   std::size_t srcNbOfCompo(srcField->getNumberOfComponents());
1137   if(array)
1138     {
1139       targetField->checkConsistencyLight();
1140       if(srcNbOfCompo!=targetField->getNumberOfComponents())
1141         throw INTERP_KERNEL::Exception("Number of components mismatch !");
1142     }
1143   else
1144     {
1145       if(!isDftVal)
1146         throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !");
1147       MCAuto<DataArrayDouble> tmp(DataArrayDouble::New());
1148       tmp->alloc(targetField->getNumberOfTuples(),srcNbOfCompo);
1149       targetField->setArray(tmp);
1150     }
1151   computeDeno(srcField->getNature(),srcField,targetField);
1152   double *resPointer(targetField->getArray()->getPointer());
1153   const double *inputPointer(srcField->getArray()->getConstPointer());
1154   computeProduct(inputPointer,(int)srcNbOfCompo,isDftVal,dftValue,resPointer);
1155 }
1156
1157 void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1158 {
1159   if(nat==NoNature)
1160     return computeDenoFromScratch(nat,srcField,trgField);
1161   else if(nat!=_nature_of_deno)
1162     return computeDenoFromScratch(nat,srcField,trgField);
1163   else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis())
1164     return computeDenoFromScratch(nat,srcField,trgField);
1165 }
1166
1167 void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1168 {
1169   _nature_of_deno=nat;
1170   switch(_nature_of_deno)
1171   {
1172     case IntensiveMaximum:
1173       {
1174         ComputeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1175         break;
1176       }
1177     case ExtensiveMaximum:
1178       {
1179         MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1180         MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1181         const double *denoPtr=deno->getArray()->getConstPointer();
1182         const double *denoRPtr=denoR->getArray()->getConstPointer();
1183         if(trgField->getMesh()->getMeshDimension()==-1)
1184           {
1185             double *denoRPtr2=denoR->getArray()->getPointer();
1186             denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1187           }
1188         if(srcField->getMesh()->getMeshDimension()==-1)
1189           {
1190             double *denoPtr2=deno->getArray()->getPointer();
1191             denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1192           }
1193         int idx=0;
1194         for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1195           for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1196             {
1197               _deno_multiply[idx][(*iter2).first]=denoPtr[(*iter2).first];
1198               _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[idx];
1199             }
1200         deno->decrRef();
1201         denoR->decrRef();
1202         break;
1203       }
1204     case ExtensiveConservation:
1205       {
1206         ComputeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1207         break;
1208       }
1209     case IntensiveConservation:
1210       {
1211         MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1212         MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1213         const double *denoPtr=deno->getArray()->getConstPointer();
1214         const double *denoRPtr=denoR->getArray()->getConstPointer();
1215         if(trgField->getMesh()->getMeshDimension()==-1)
1216           {
1217             double *denoRPtr2=denoR->getArray()->getPointer();
1218             denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1219           }
1220         if(srcField->getMesh()->getMeshDimension()==-1)
1221           {
1222             double *denoPtr2=deno->getArray()->getPointer();
1223             denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1224           }
1225         mcIdType idx=0;
1226         for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1227           for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1228             {
1229               _deno_multiply[idx][(*iter2).first]=denoPtr[idx];
1230               _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[(*iter2).first];
1231             }
1232         deno->decrRef();
1233         denoR->decrRef();
1234         break;
1235       }
1236     case NoNature:
1237       throw INTERP_KERNEL::Exception("No nature specified ! Select one !");
1238   }
1239 }
1240
1241 void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)
1242 {
1243   int idx=0;
1244   double *tmp=new double[inputNbOfCompo];
1245   for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1246     {
1247       if((*iter1).empty())
1248         {
1249           if(isDftVal)
1250             std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1251           continue;
1252         }
1253       else
1254         std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,0.);
1255       std::map<mcIdType,double>::const_iterator iter3=_deno_multiply[idx].begin();
1256       for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter3++)
1257         {
1258           std::transform(inputPointer+(*iter2).first*inputNbOfCompo,inputPointer+((*iter2).first+1)*inputNbOfCompo,tmp,std::bind(std::multiplies<double>(),std::placeholders::_1,(*iter2).second/(*iter3).second));
1259           std::transform(tmp,tmp+inputNbOfCompo,resPointer+idx*inputNbOfCompo,resPointer+idx*inputNbOfCompo,std::plus<double>());
1260         }
1261     }
1262   delete [] tmp;
1263 }
1264
1265 void MEDCouplingRemapper::computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer)
1266 {
1267   std::vector<bool> isReached(_deno_reverse_multiply.size(),false);
1268   mcIdType idx=0;
1269   double *tmp=new double[inputNbOfCompo];
1270   std::fill(resPointer,resPointer+inputNbOfCompo*_deno_reverse_multiply.size(),0.);
1271   for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1272     {
1273       for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1274         {
1275           isReached[(*iter2).first]=true;
1276           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]));
1277           std::transform(tmp,tmp+inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,std::plus<double>());
1278         }
1279     }
1280   delete [] tmp;
1281   idx=0;
1282   for(std::vector<bool>::const_iterator iter3=isReached.begin();iter3!=isReached.end();iter3++,idx++)
1283     if(!*iter3)
1284       std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1285 }
1286
1287 void MEDCouplingRemapper::ReverseMatrix(const std::vector<std::map<mcIdType,double> >& matIn, mcIdType nbColsMatIn, std::vector<std::map<mcIdType,double> >& matOut)
1288 {
1289   matOut.resize(nbColsMatIn);
1290   mcIdType id=0;
1291   for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
1292     for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1293       matOut[(*iter2).first][id]=(*iter2).second;
1294 }
1295
1296 void MEDCouplingRemapper::ComputeRowSumAndColSum(const std::vector<std::map<mcIdType,double> >& matrixDeno,
1297                                                  std::vector<std::map<mcIdType,double> >& deno, std::vector<std::map<mcIdType,double> >& denoReverse)
1298 {
1299   std::map<mcIdType,double> values;
1300   mcIdType idx=0;
1301   for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1302     {
1303       double sum=0.;
1304       for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1305         {
1306           sum+=(*iter2).second;
1307           values[(*iter2).first]+=(*iter2).second;
1308         }
1309       for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1310         deno[idx][(*iter2).first]=sum;
1311     }
1312   idx=0;
1313   for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1314     {
1315       for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1316         denoReverse[(*iter2).first][idx]=values[(*iter2).first];
1317     }
1318 }
1319
1320 void MEDCouplingRemapper::ComputeColSumAndRowSum(const std::vector<std::map<mcIdType,double> >& matrixDeno,
1321                                                  std::vector<std::map<mcIdType,double> >& deno, std::vector<std::map<mcIdType,double> >& denoReverse)
1322 {
1323   std::map<mcIdType,double> values;
1324   mcIdType idx=0;
1325   for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1326     {
1327       double sum=0.;
1328       for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1329         {
1330           sum+=(*iter2).second;
1331           values[(*iter2).first]+=(*iter2).second;
1332         }
1333       for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1334         denoReverse[(*iter2).first][idx]=sum;
1335     }
1336   idx=0;
1337   for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1338     {
1339       for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1340         deno[idx][(*iter2).first]=values[(*iter2).first];
1341     }
1342 }
1343
1344 void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<mcIdType,double> >& m1D,
1345                                                                      const std::vector< std::map<mcIdType,double> >& m2D,
1346                                                                      const mcIdType *corrCellIdSrc, mcIdType nbOf2DCellsSrc, mcIdType nbOf1DCellsSrc,
1347                                                                      const mcIdType *corrCellIdTrg)
1348 {
1349   mcIdType nbOf2DCellsTrg=ToIdType(m2D.size());
1350   mcIdType nbOf1DCellsTrg=ToIdType(m1D.size());
1351   mcIdType nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
1352   _matrix.resize(nbOf3DCellsTrg);
1353   mcIdType id2R=0;
1354   for(std::vector< std::map<mcIdType,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
1355     {
1356       for(std::map<mcIdType,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
1357         {
1358           mcIdType id1R=0;
1359           for(std::vector< std::map<mcIdType,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
1360             {
1361               for(std::map<mcIdType,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
1362                 {
1363                   _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
1364                 }
1365             }
1366         }
1367     }
1368 }
1369
1370 void MEDCouplingRemapper::PrintMatrix(const std::vector<std::map<mcIdType,double> >& m)
1371 {
1372   mcIdType id=0;
1373   for(std::vector<std::map<mcIdType,double> >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++)
1374     {
1375       std::cout << "Target Cell # " << id << " : ";
1376       for(std::map<mcIdType,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1377         std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), ";
1378       std::cout << std::endl;
1379     }
1380 }
1381
1382 const std::vector<std::map<mcIdType,double> >& MEDCouplingRemapper::getCrudeMatrix() const
1383 {
1384   return _matrix;
1385 }
1386
1387 /*!
1388  * Returns the number of columns of matrix returned by MEDCouplingRemapper::getCrudeMatrix method.
1389  */
1390 mcIdType MEDCouplingRemapper::getNumberOfColsOfMatrix() const
1391 {
1392   return ToIdType(_deno_reverse_multiply.size());
1393 }
1394
1395 /*!
1396  * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1397  * If not the behaviour is unpredictable.
1398  * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than \a maxValAbs this coefficient is
1399  * 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.
1400  * 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
1401  * can occur.
1402  *
1403  * \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.
1404  * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged.
1405  * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix
1406  */
1407 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs(double maxValAbs)
1408 {
1409   int ret=0;
1410   std::vector<std::map<mcIdType,double> > matrixNew(_matrix.size());
1411   mcIdType i=0;
1412   for(std::vector<std::map<mcIdType,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++,i++)
1413     {
1414       std::map<mcIdType,double>& rowNew=matrixNew[i];
1415       for(std::map<mcIdType,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1416         {
1417           if(fabs((*it2).second)>maxValAbs)
1418             rowNew[(*it2).first]=(*it2).second;
1419           else
1420             ret++;
1421         }
1422     }
1423   if(ret>0)
1424     _matrix=matrixNew;
1425   return ret;
1426 }
1427
1428 /*!
1429  * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1430  * If not the behaviour is unpredictable.
1431  * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than delta multiplied by \a scaleFactor this coefficient is
1432  * 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.
1433  * delta is the value returned by MEDCouplingRemapper::getMaxValueInCrudeMatrix method.
1434  * 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
1435  * can occur.
1436  *
1437  * \param [in] scaleFactor is the scale factor from which coefficients lower than \a scaleFactor times range width of coefficients are set to zero.
1438  * \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
1439  *         that all coefficients are null.
1440  * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs
1441  */
1442 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix(double scaleFactor)
1443 {
1444   double maxVal=getMaxValueInCrudeMatrix();
1445   if(maxVal==0.)
1446     return -1;
1447   return nullifiedTinyCoeffInCrudeMatrixAbs(scaleFactor*maxVal);
1448 }
1449
1450 /*!
1451  * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1452  * If not the behaviour is unpredictable.
1453  * This method returns the maximum of the absolute values of coefficients into the sparse crude matrix.
1454  * The returned value is positive.
1455  */
1456 double MEDCouplingRemapper::getMaxValueInCrudeMatrix() const
1457 {
1458   double ret=0.;
1459   for(std::vector<std::map<mcIdType,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++)
1460     for(std::map<mcIdType,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1461       if(fabs((*it2).second)>ret)
1462         ret=fabs((*it2).second);
1463   return ret;
1464 }