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