src/MEDCoupling/MEDCouplingRemapper.cxx

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