src/MEDCoupling/MEDCouplingRemapper.cxx

   1 // Copyright (C) 2007-2013  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.
   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       nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
 373       ReverseMatrix(matrixTmp,nbCols,_matrix);
 374       nbCols=matrixTmp.size();
 375     }
 376   else if(srcMeshDim==2 && trgMeshDim==1 && srcSpaceDim==2)
 377     {
 378       if(getIntersectionType()==INTERP_KERNEL::PointLocator)
 379         {
 380           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 381           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
 382           INTERP_KERNEL::Interpolation2D interpolation(*this);
 383           nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 384         }
 385       else
 386         {
 387           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 388           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
 389           INTERP_KERNEL::Interpolation2D1D interpolation(*this);
 390           std::vector<std::map<int,double> > matrixTmp;
 391           nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
 392           ReverseMatrix(matrixTmp,nbCols,_matrix);
 393           nbCols=matrixTmp.size();
 394           INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
 395           if(!duplicateFaces.empty())
 396             {
 397               std::ostringstream oss; oss << "An unexpected situation happend ! For the following 1D Cells are part of edges shared by 2D cells :\n";
 398               for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
 399                 {
 400                   oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
 401                   std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
 402                   oss << std::endl;
 403                 }
 404             }
 405         }
 406     }
 407   else if(srcMeshDim==1 && trgMeshDim==2 && srcSpaceDim==2)
 408     {
 409       if(getIntersectionType()==INTERP_KERNEL::PointLocator)
 410         {
 411           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 412           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
 413           INTERP_KERNEL::Interpolation2D interpolation(*this);
 414           std::vector<std::map<int,double> > matrixTmp;
 415           nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
 416           ReverseMatrix(matrixTmp,nbCols,_matrix);
 417           nbCols=matrixTmp.size();
 418         }
 419       else
 420         {
 421           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 422           MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
 423           INTERP_KERNEL::Interpolation2D1D interpolation(*this);
 424           nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 425           INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
 426           if(!duplicateFaces.empty())
 427             {
 428               std::ostringstream oss; oss << "An unexpected situation happend ! For the following 1D Cells are part of edges shared by 2D cells :\n";
 429               for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
 430                 {
 431                   oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
 432                   std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
 433                   oss << std::endl;
 434                 }
 435             }
 436         }
 437     }
 438   else if(srcMeshDim==2 && trgMeshDim==3 && srcSpaceDim==3)
 439     {
 440       MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 441       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
 442       INTERP_KERNEL::Interpolation3D2D interpolation(*this);
 443       nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
 444       INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
 445       if(!duplicateFaces.empty())
 446         {
 447           std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
 448           for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
 449             {
 450               oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
 451               std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
 452               oss << std::endl;
 453             }
 454         }
 455     }
 456   else if(srcMeshDim==3 && trgMeshDim==2 && srcSpaceDim==3)
 457     {
 458       MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 459       MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
 460       INTERP_KERNEL::Interpolation3D2D interpolation(*this);
 461       std::vector<std::map<int,double> > matrixTmp;
 462       nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
 463       ReverseMatrix(matrixTmp,nbCols,_matrix);
 464       nbCols=matrixTmp.size();
 465       INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
 466       if(!duplicateFaces.empty())
 467         {
 468           std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
 469           for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
 470             {
 471               oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
 472               std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
 473               oss << std::endl;
 474             }
 475         }
 476     }
 477   else if(trgMeshDim==-1)
 478     {
 479       if(srcMeshDim==2 && srcSpaceDim==2)
 480         {
 481           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
 482           INTERP_KERNEL::Interpolation2D interpolation(*this);
 483           nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
 484         }
 485       else if(srcMeshDim==3 && srcSpaceDim==3)
 486         {
 487           MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
 488           INTERP_KERNEL::Interpolation3D interpolation(*this);
 489           nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
 490         }
 491       else if(srcMeshDim==2 && srcSpaceDim==3)
 492         {
 493           MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
 494           INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
 495           nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
 496         }
 497       else
 498         throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
 499     }
 500   else if(srcMeshDim==-1)
 501     {
 502       if(trgMeshDim==2 && trgSpaceDim==2)
 503         {
 504           MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(target_mesh);
 505           INTERP_KERNEL::Interpolation2D interpolation(*this);
 506           nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
 507         }
 508       else if(trgMeshDim==3 && trgSpaceDim==3)
 509         {
 510           MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(target_mesh);
 511           INTERP_KERNEL::Interpolation3D interpolation(*this);
 512           nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
 513         }
 514       else if(trgMeshDim==2 && trgSpaceDim==3)
 515         {
 516           MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(target_mesh);
 517           INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
 518           nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
 519         }
 520       else
 521         throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh from -1D sourceMesh");
 522     }
 523   else
 524     throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension");
 525   _deno_multiply.clear();
 526   _deno_multiply.resize(_matrix.size());
 527   _deno_reverse_multiply.clear();
 528   _deno_reverse_multiply.resize(nbCols);
 529   declareAsNew();
 530   return 1;
 531 }
 532
 533 int MEDCouplingRemapper::prepareInterpKernelOnlyEE()
 534 {
 535   std::string srcMeth,trgMeth;
 536   std::string methC=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
 537   const MEDCouplingExtrudedMesh *src_mesh=static_cast<const MEDCouplingExtrudedMesh *>(_src_ft->getMesh());
 538   const MEDCouplingExtrudedMesh *target_mesh=static_cast<const MEDCouplingExtrudedMesh *>(_target_ft->getMesh());
 539   if(methC!="P0P0")
 540     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyEE : Only P0P0 method implemented for Extruded/Extruded meshes !");
 541   MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh->getMesh2D());
 542   MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh->getMesh2D());
 543   INTERP_KERNEL::Interpolation3DSurf interpolation2D(*this);
 544   std::vector<std::map<int,double> > matrix2D;
 545   int nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,methC);
 546   MEDCouplingUMesh *s1D,*t1D;
 547   double v[3];
 548   MEDCouplingExtrudedMesh::Project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v);
 549   MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D);
 550   MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
 551   std::vector<std::map<int,double> > matrix1D;
 552   INTERP_KERNEL::Interpolation1D interpolation1D(*this);
 553   int nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,methC);
 554   s1D->decrRef();
 555   t1D->decrRef();
 556   buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
 557                                              target_mesh->getMesh3DIds()->getConstPointer());
 558   //
 559   _deno_multiply.clear();
 560   _deno_multiply.resize(_matrix.size());
 561   _deno_reverse_multiply.clear();
 562   _deno_reverse_multiply.resize(nbCols2D*nbCols1D);
 563   declareAsNew();
 564   return 1;
 565 }
 566
 567 int MEDCouplingRemapper::prepareInterpKernelOnlyUC()
 568 {
 569   std::string srcMeth,trgMeth;
 570   std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
 571   if(methodCpp!="P0P0")
 572     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only P0P0 interpolation supported for the moment !");
 573   const MEDCouplingUMesh *src_mesh=static_cast<const MEDCouplingUMesh *>(_src_ft->getMesh());
 574   const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
 575   const int srcMeshDim=src_mesh->getMeshDimension();
 576   const int srcSpceDim=src_mesh->getSpaceDimension();
 577   const int trgMeshDim=target_mesh->getMeshDimension();
 578   if(srcMeshDim!=srcSpceDim || srcMeshDim!=trgMeshDim)
 579     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 !");
 580   std::vector<std::map<int,double> > res;
 581   switch(srcMeshDim)
 582     {
 583     case 1:
 584       {
 585         MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
 586         MEDCouplingNormalizedUnstructuredMesh<1,1> sourceWrapper(src_mesh);
 587         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 588         myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
 589         break;
 590       }
 591     case 2:
 592       {
 593         MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
 594         MEDCouplingNormalizedUnstructuredMesh<2,2> sourceWrapper(src_mesh);
 595         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 596         myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
 597         break;
 598       }
 599     case 3:
 600       {
 601         MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
 602         MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(src_mesh);
 603         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 604         myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
 605         break;
 606       }
 607     default:
 608       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only dimension 1 2 or 3 supported !");
 609     }
 610   ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix);
 611   nullifiedTinyCoeffInCrudeMatrixAbs(0.);
 612   //
 613   _deno_multiply.clear();
 614   _deno_multiply.resize(_matrix.size());
 615   _deno_reverse_multiply.clear();
 616   _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
 617   declareAsNew();
 618   return 1;
 619 }
 620
 621 int MEDCouplingRemapper::prepareInterpKernelOnlyCU()
 622 {
 623   std::string srcMeth,trgMeth;
 624   std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
 625   if(methodCpp!="P0P0")
 626     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only P0P0 interpolation supported for the moment !");
 627   const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
 628   const MEDCouplingUMesh *target_mesh=static_cast<const MEDCouplingUMesh *>(_target_ft->getMesh());
 629   const int srcMeshDim=src_mesh->getMeshDimension();
 630   const int trgMeshDim=target_mesh->getMeshDimension();
 631   const int trgSpceDim=target_mesh->getSpaceDimension();
 632   if(trgMeshDim!=trgSpceDim || trgMeshDim!=srcMeshDim)
 633     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 !");
 634   switch(srcMeshDim)
 635     {
 636     case 1:
 637       {
 638         MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
 639         MEDCouplingNormalizedUnstructuredMesh<1,1> targetWrapper(target_mesh);
 640         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 641         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 642         break;
 643       }
 644     case 2:
 645       {
 646         MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
 647         MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(target_mesh);
 648         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 649         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 650         break;
 651       }
 652     case 3:
 653       {
 654         MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
 655         MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(target_mesh);
 656         INTERP_KERNEL::InterpolationCU myInterpolator(*this);
 657         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 658         break;
 659       }
 660     default:
 661       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only dimension 1 2 or 3 supported !");
 662     }
 663   nullifiedTinyCoeffInCrudeMatrixAbs(0.);
 664   //
 665   _deno_multiply.clear();
 666   _deno_multiply.resize(_matrix.size());
 667   _deno_reverse_multiply.clear();
 668   _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
 669   declareAsNew();
 670   return 1;
 671 }
 672
 673 int MEDCouplingRemapper::prepareInterpKernelOnlyCC()
 674 {
 675   std::string srcMeth,trgMeth;
 676   std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
 677   if(methodCpp!="P0P0")
 678     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only P0P0 interpolation supported for the moment !");
 679   const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
 680   const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
 681   const int srcMeshDim=src_mesh->getMeshDimension();
 682   const int trgMeshDim=target_mesh->getMeshDimension();
 683   if(trgMeshDim!=srcMeshDim)
 684     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : dim of target cartesian should be equal to dim of source cartesian dimension !");
 685   switch(srcMeshDim)
 686     {
 687     case 1:
 688       {
 689         MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
 690         MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
 691         INTERP_KERNEL::InterpolationCC myInterpolator(*this);
 692         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 693         break;
 694       }
 695     case 2:
 696       {
 697         MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
 698         MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
 699         INTERP_KERNEL::InterpolationCC myInterpolator(*this);
 700         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 701         break;
 702       }
 703     case 3:
 704       {
 705         MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
 706         MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
 707         INTERP_KERNEL::InterpolationCC myInterpolator(*this);
 708         myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
 709         break;
 710       }
 711     default:
 712       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only dimension 1 2 or 3 supported !");
 713     }
 714   nullifiedTinyCoeffInCrudeMatrixAbs(0.);
 715   //
 716   _deno_multiply.clear();
 717   _deno_multiply.resize(_matrix.size());
 718   _deno_reverse_multiply.clear();
 719   _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
 720   declareAsNew();
 721   return 1;
 722 }
 723
 724 int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss()
 725 {
 726   if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
 727       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 !");
 728   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> trgLoc=_target_ft->getLocalizationOfDiscr();
 729   const double *trgLocPtr=trgLoc->begin();
 730   int trgSpaceDim=trgLoc->getNumberOfComponents();
 731   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> srcOffsetArr=_src_ft->getDiscretization()->getOffsetArr(_src_ft->getMesh());
 732   if(trgSpaceDim!=_src_ft->getMesh()->getSpaceDimension())
 733     {
 734       std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss : space dimensions mismatch between source and target !";
 735       oss << " Target discretization localization has dimension " << trgSpaceDim << ", whereas the space dimension of source is equal to ";
 736       oss << _src_ft->getMesh()->getSpaceDimension() << " !";
 737       throw INTERP_KERNEL::Exception(oss.str().c_str());
 738     }
 739   const int *srcOffsetArrPtr=srcOffsetArr->begin();
 740   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> srcLoc=_src_ft->getLocalizationOfDiscr();
 741   const double *srcLocPtr=srcLoc->begin();
 742   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> eltsArr,eltsIndexArr;
 743   int trgNbOfGaussPts=trgLoc->getNumberOfTuples();
 744   _matrix.resize(trgNbOfGaussPts);
 745   _src_ft->getMesh()->getCellsContainingPoints(trgLoc->begin(),trgNbOfGaussPts,getPrecision(),eltsArr,eltsIndexArr);
 746   const int *elts(eltsArr->begin()),*eltsIndex(eltsIndexArr->begin());
 747   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nbOfSrcCellsShTrgPts(eltsIndexArr->deltaShiftIndex());
 748   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ids0=nbOfSrcCellsShTrgPts->getIdsNotEqual(0);
 749   for(const int *trgId=ids0->begin();trgId!=ids0->end();trgId++)
 750     {
 751       const double *ptTrg=trgLocPtr+trgSpaceDim*(*trgId);
 752       int srcCellId=elts[eltsIndex[*trgId]];
 753       double dist=std::numeric_limits<double>::max();
 754       int srcEntry=-1;
 755       for(int srcId=srcOffsetArrPtr[srcCellId];srcId<srcOffsetArrPtr[srcCellId+1];srcId++)
 756         {
 757           const double *ptSrc=srcLocPtr+trgSpaceDim*srcId;
 758           double tmp=0.;
 759           for(int i=0;i<trgSpaceDim;i++)
 760             tmp+=(ptTrg[i]-ptSrc[i])*(ptTrg[i]-ptSrc[i]);
 761           if(tmp<dist)
 762             { dist=tmp; srcEntry=srcId; }
 763         }
 764       _matrix[*trgId][srcEntry]=1.;
 765     }
 766   if(ids0->getNumberOfTuples()!=trgNbOfGaussPts)
 767     {
 768       MEDCouplingAutoRefCountObjectPtr<DataArrayInt> orphanTrgIds=nbOfSrcCellsShTrgPts->getIdsEqual(0);
 769       MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> orphanTrg=trgLoc->selectByTupleId(orphanTrgIds->begin(),orphanTrgIds->end());
 770       MEDCouplingAutoRefCountObjectPtr<DataArrayInt> srcIdPerTrg=srcLoc->findClosestTupleId(orphanTrg);
 771       const int *srcIdPerTrgPtr=srcIdPerTrg->begin();
 772       for(const int *orphanTrgId=orphanTrgIds->begin();orphanTrgId!=orphanTrgIds->end();orphanTrgId++,srcIdPerTrgPtr++)
 773         _matrix[*orphanTrgId][*srcIdPerTrgPtr]=2.;
 774     }
 775   _deno_multiply.clear();
 776   _deno_multiply.resize(_matrix.size());
 777   _deno_reverse_multiply.clear();
 778   _deno_reverse_multiply.resize(srcLoc->getNumberOfTuples());
 779   declareAsNew();
 780   return 1;
 781 }
 782
 783 /*!
 784  * This method checks that the input interpolation \a method is managed by not INTERP_KERNEL only methods.
 785  * If no an INTERP_KERNEL::Exception will be thrown. If yes, a magic number will be returned to switch in the MEDCouplingRemapper::prepareNotInterpKernelOnly method.
 786  */
 787 int MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel(const std::string& method)
 788 {
 789   if(method=="GAUSSGAUSS")
 790     return 0;
 791   std::ostringstream oss; oss << "MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel : ";
 792   oss << "The method \"" << method << "\" is not manageable by not INTERP_KERNEL only method.";
 793   oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS !";
 794   throw INTERP_KERNEL::Exception(oss.str().c_str());
 795 }
 796
 797 /*!
 798  * This method determines regarding \c _interp_matrix_pol attribute ( set by MEDCouplingRemapper::setInterpolationMatrixPolicy and by default equal
 799  * 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
 800  * only INTERP_KERNEL method should be applied.
 801  */
 802 bool MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly() const
 803 {
 804   std::string srcm,trgm,method;
 805   method=checkAndGiveInterpolationMethodStr(srcm,trgm);
 806   switch(_interp_matrix_pol)
 807     {
 808     case IK_ONLY_PREFERED:
 809       {
 810         try
 811           {
 812             std::string tmp1,tmp2;
 813             INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,tmp1,tmp2);
 814             return true;
 815           }
 816         catch(INTERP_KERNEL::Exception& /*e*/)
 817           {
 818             return false;
 819           }
 820       }
 821     case NOT_IK_ONLY_PREFERED:
 822       {
 823         try
 824           {
 825             CheckInterpolationMethodManageableByNotOnlyInterpKernel(method);
 826             return false;
 827           }
 828         catch(INTERP_KERNEL::Exception& /*e*/)
 829           {
 830             return true;
 831           }
 832       }
 833     case IK_ONLY_FORCED:
 834       return true;
 835     case NOT_IK_ONLY_FORCED:
 836       return false;
 837     default:
 838       throw INTERP_KERNEL::Exception("MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly : internal error ! The interpolation matrix policy is not managed ! Try to change it using MEDCouplingRemapper::setInterpolationMatrixPolicy !");
 839     }
 840 }
 841
 842 void MEDCouplingRemapper::updateTime() const
 843 {
 844 }
 845
 846 void MEDCouplingRemapper::checkPrepare() const
 847 {
 848   const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
 849   if(!s || !t)
 850     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that MEDCouplingRemapper::prepare(Ex) has not been called !");
 851   if(!s->getMesh() || !t->getMesh())
 852     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that no all field templates have their mesh set !");
 853 }
 854
 855 /*!
 856  * This method builds a code considering already set field discretization int \a this : \a _src_ft and \a _target_ft.
 857  * This method returns 3 informations (2 in ouput parameters and 1 in return).
 858  *
 859  * \param [out] srcMeth the string code of the discretization of source field template
 860  * \param [out] trgMeth the string code of the discretization of target field template
 861  * \return the standardized string code (compatible with INTERP_KERNEL) for matrix of numerators (in \a _matrix)
 862  */
 863 std::string MEDCouplingRemapper::checkAndGiveInterpolationMethodStr(std::string& srcMeth, std::string& trgMeth) const
 864 {
 865   const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
 866   if(!s || !t)
 867     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have been set !");
 868   if(!s->getMesh() || !t->getMesh())
 869     throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have their mesh set !");
 870   srcMeth=_src_ft->getDiscretization()->getRepr();
 871   trgMeth=_target_ft->getDiscretization()->getRepr();
 872   std::string method(srcMeth); method+=trgMeth;
 873   return method;
 874 }
 875
 876 void MEDCouplingRemapper::releaseData(bool matrixSuppression)
 877 {
 878   _src_ft=0;
 879   _target_ft=0;
 880   if(matrixSuppression)
 881     {
 882       _matrix.clear();
 883       _deno_multiply.clear();
 884       _deno_reverse_multiply.clear();
 885     }
 886 }
 887
 888 void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue)
 889 {
 890   checkPrepare();
 891   if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
 892     throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
 893   if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
 894     throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
 895   if(srcField->getNature()!=targetField->getNature())
 896     throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
 897   DataArrayDouble *array=targetField->getArray();
 898   int srcNbOfCompo=srcField->getNumberOfComponents();
 899   if(array)
 900     {
 901       if(srcNbOfCompo!=targetField->getNumberOfComponents())
 902         throw INTERP_KERNEL::Exception("Number of components mismatch !");
 903     }
 904   else
 905     {
 906       if(!isDftVal)
 907         throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !");
 908       array=DataArrayDouble::New();
 909       array->alloc(targetField->getNumberOfTuples(),srcNbOfCompo);
 910       targetField->setArray(array);
 911       array->decrRef();
 912     }
 913   computeDeno(srcField->getNature(),srcField,targetField);
 914   double *resPointer=array->getPointer();
 915   const double *inputPointer=srcField->getArray()->getConstPointer();
 916   computeProduct(inputPointer,srcNbOfCompo,isDftVal,dftValue,resPointer);
 917 }
 918
 919 void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
 920 {
 921   if(nat==NoNature)
 922     return computeDenoFromScratch(nat,srcField,trgField);
 923   else if(nat!=_nature_of_deno)
 924    return computeDenoFromScratch(nat,srcField,trgField);
 925   else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis())
 926     return computeDenoFromScratch(nat,srcField,trgField);
 927 }
 928
 929 void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
 930 {
 931   _nature_of_deno=nat;
 932   _time_deno_update=getTimeOfThis();
 933   switch(_nature_of_deno)
 934     {
 935     case ConservativeVolumic:
 936       {
 937         ComputeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
 938         break;
 939       }
 940     case Integral:
 941       {
 942         MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
 943         MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
 944         const double *denoPtr=deno->getArray()->getConstPointer();
 945         const double *denoRPtr=denoR->getArray()->getConstPointer();
 946         if(trgField->getMesh()->getMeshDimension()==-1)
 947           {
 948             double *denoRPtr2=denoR->getArray()->getPointer();
 949             denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
 950           }
 951         if(srcField->getMesh()->getMeshDimension()==-1)
 952           {
 953             double *denoPtr2=deno->getArray()->getPointer();
 954             denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
 955           }
 956         int idx=0;
 957         for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
 958           for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 959             {
 960               _deno_multiply[idx][(*iter2).first]=denoPtr[(*iter2).first];
 961               _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[idx];
 962             }
 963         deno->decrRef();
 964         denoR->decrRef();
 965         break;
 966       }
 967     case IntegralGlobConstraint:
 968       {
 969         ComputeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
 970         break;
 971       }
 972     case RevIntegral:
 973       {
 974         MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
 975         MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
 976         const double *denoPtr=deno->getArray()->getConstPointer();
 977         const double *denoRPtr=denoR->getArray()->getConstPointer();
 978         if(trgField->getMesh()->getMeshDimension()==-1)
 979           {
 980             double *denoRPtr2=denoR->getArray()->getPointer();
 981             denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
 982           }
 983         if(srcField->getMesh()->getMeshDimension()==-1)
 984           {
 985             double *denoPtr2=deno->getArray()->getPointer();
 986             denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
 987           }
 988         int idx=0;
 989         for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
 990           for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
 991             {
 992               _deno_multiply[idx][(*iter2).first]=denoPtr[idx];
 993               _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[(*iter2).first];
 994             }
 995         deno->decrRef();
 996         denoR->decrRef();
 997         break;
 998       }
 999     case NoNature:
1000       throw INTERP_KERNEL::Exception("No nature specified ! Select one !");
1001     }
1002 }
1003
1004 void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)
1005 {
1006   int idx=0;
1007   double *tmp=new double[inputNbOfCompo];
1008   for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1009     {
1010       if((*iter1).empty())
1011         {
1012           if(isDftVal)
1013             std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1014           continue;
1015         }
1016       else
1017         std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,0.);
1018       std::map<int,double>::const_iterator iter3=_deno_multiply[idx].begin();
1019       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter3++)
1020         {
1021           std::transform(inputPointer+(*iter2).first*inputNbOfCompo,inputPointer+((*iter2).first+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/(*iter3).second));
1022           std::transform(tmp,tmp+inputNbOfCompo,resPointer+idx*inputNbOfCompo,resPointer+idx*inputNbOfCompo,std::plus<double>());
1023         }
1024     }
1025   delete [] tmp;
1026 }
1027
1028 void MEDCouplingRemapper::computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer)
1029 {
1030   std::vector<bool> isReached(_deno_reverse_multiply.size(),false);
1031   int idx=0;
1032   double *tmp=new double[inputNbOfCompo];
1033   std::fill(resPointer,resPointer+inputNbOfCompo*_deno_reverse_multiply.size(),0.);
1034   for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1035     {
1036       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1037         {
1038           isReached[(*iter2).first]=true;
1039           std::transform(inputPointer+idx*inputNbOfCompo,inputPointer+(idx+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/_deno_reverse_multiply[(*iter2).first][idx]));
1040           std::transform(tmp,tmp+inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,std::plus<double>());
1041         }
1042     }
1043   delete [] tmp;
1044   idx=0;
1045   for(std::vector<bool>::const_iterator iter3=isReached.begin();iter3!=isReached.end();iter3++,idx++)
1046     if(!*iter3)
1047       std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1048 }
1049
1050 void MEDCouplingRemapper::ReverseMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
1051 {
1052   matOut.resize(nbColsMatIn);
1053   int id=0;
1054   for(std::vector<std::map<int,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
1055     for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1056       matOut[(*iter2).first][id]=(*iter2).second;
1057 }
1058
1059 void MEDCouplingRemapper::ComputeRowSumAndColSum(const std::vector<std::map<int,double> >& matrixDeno,
1060                                                  std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1061 {
1062   std::map<int,double> values;
1063   int idx=0;
1064   for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1065     {
1066       double sum=0.;
1067       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1068         {
1069           sum+=(*iter2).second;
1070           values[(*iter2).first]+=(*iter2).second;
1071         }
1072       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1073         deno[idx][(*iter2).first]=sum;
1074     }
1075   idx=0;
1076   for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1077     {
1078       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1079         denoReverse[(*iter2).first][idx]=values[(*iter2).first];
1080     }
1081 }
1082
1083 void MEDCouplingRemapper::ComputeColSumAndRowSum(const std::vector<std::map<int,double> >& matrixDeno,
1084                                                  std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1085 {
1086   std::map<int,double> values;
1087   int idx=0;
1088   for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1089     {
1090       double sum=0.;
1091       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1092         {
1093           sum+=(*iter2).second;
1094           values[(*iter2).first]+=(*iter2).second;
1095         }
1096       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1097         denoReverse[(*iter2).first][idx]=sum;
1098     }
1099   idx=0;
1100   for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1101     {
1102       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1103         deno[idx][(*iter2).first]=values[(*iter2).first];
1104     }
1105 }
1106
1107 void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<int,double> >& m1D,
1108                                                                      const std::vector< std::map<int,double> >& m2D,
1109                                                                      const int *corrCellIdSrc, int nbOf2DCellsSrc, int nbOf1DCellsSrc,
1110                                                                      const int *corrCellIdTrg)
1111 {
1112   int nbOf2DCellsTrg=m2D.size();
1113   int nbOf1DCellsTrg=m1D.size();
1114   int nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
1115   _matrix.resize(nbOf3DCellsTrg);
1116   int id2R=0;
1117   for(std::vector< std::map<int,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
1118     {
1119       for(std::map<int,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
1120         {
1121           int id1R=0;
1122           for(std::vector< std::map<int,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
1123             {
1124               for(std::map<int,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
1125                 {
1126                   _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
1127                 }
1128             }
1129         }
1130     }
1131 }
1132
1133 void MEDCouplingRemapper::PrintMatrix(const std::vector<std::map<int,double> >& m)
1134 {
1135   int id=0;
1136   for(std::vector<std::map<int,double> >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++)
1137     {
1138       std::cout << "Target Cell # " << id << " : ";
1139       for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1140         std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), ";
1141       std::cout << std::endl;
1142     }
1143 }
1144
1145 const std::vector<std::map<int,double> >& MEDCouplingRemapper::getCrudeMatrix() const
1146 {
1147   return _matrix;
1148 }
1149
1150 /*!
1151  * Returns the number of columns of matrix returned by MEDCouplingRemapper::getCrudeMatrix method.
1152  */
1153 int MEDCouplingRemapper::getNumberOfColsOfMatrix() const
1154 {
1155   return (int)_deno_reverse_multiply.size();
1156 }
1157
1158 /*!
1159  * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1160  * If not the behaviour is unpredictable.
1161  * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than \a maxValAbs this coefficient is
1162  * 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.
1163  * 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
1164  * can occur.
1165  *
1166  * \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.
1167  * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged.
1168  * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix
1169  */
1170 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs(double maxValAbs)
1171 {
1172   int ret=0;
1173   std::vector<std::map<int,double> > matrixNew(_matrix.size());
1174   int i=0;
1175   for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++,i++)
1176     {
1177       std::map<int,double>& rowNew=matrixNew[i];
1178       for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1179         {
1180           if(fabs((*it2).second)>maxValAbs)
1181             rowNew[(*it2).first]=(*it2).second;
1182           else
1183             ret++;
1184         }
1185     }
1186   if(ret>0)
1187     _matrix=matrixNew;
1188   return ret;
1189 }
1190
1191 /*!
1192  * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1193  * If not the behaviour is unpredictable.
1194  * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than delta multiplied by \a scaleFactor this coefficient is
1195  * 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.
1196  * delta is the value returned by MEDCouplingRemapper::getMaxValueInCrudeMatrix method.
1197  * 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
1198  * can occur.
1199  *
1200  * \param [in] scaleFactor is the scale factor from which coefficients lower than \a scaleFactor times range width of coefficients are set to zero.
1201  * \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
1202  *         that all coefficients are null.
1203  * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs
1204  */
1205 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix(double scaleFactor)
1206 {
1207   double maxVal=getMaxValueInCrudeMatrix();
1208   if(maxVal==0.)
1209     return -1;
1210   return nullifiedTinyCoeffInCrudeMatrixAbs(scaleFactor*maxVal);
1211 }
1212
1213 /*!
1214  * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1215  * If not the behaviour is unpredictable.
1216  * This method returns the maximum of the absolute values of coefficients into the sparse crude matrix.
1217  * The returned value is positive.
1218  */
1219 double MEDCouplingRemapper::getMaxValueInCrudeMatrix() const
1220 {
1221   double ret=0.;
1222   for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++)
1223     for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1224       if(fabs((*it2).second)>ret)
1225         ret=fabs((*it2).second);
1226   return ret;
1227 }