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