medtool/src/MEDCoupling/MEDCouplingRemapper.cxx

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