src/MEDCoupling/MEDCouplingRemapper.cxx

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