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