src/MEDCoupling/MEDCouplingRemapper.cxx

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