src/MEDCoupling/MEDCouplingRemapper.cxx

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