src/MEDCoupling/MEDCouplingFieldDiscretization.cxx

   1 // Copyright (C) 2007-2012  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 "MEDCouplingFieldDiscretization.hxx"
  22 #include "MEDCouplingCMesh.hxx"
  23 #include "MEDCouplingUMesh.hxx"
  24 #include "MEDCouplingFieldDouble.hxx"
  25 #include "MEDCouplingAutoRefCountObjectPtr.hxx"
  26
  27 #include "CellModel.hxx"
  28 #include "InterpolationUtils.hxx"
  29 #include "InterpKernelAutoPtr.hxx"
  30 #include "InterpKernelGaussCoords.hxx"
  31 #include "InterpKernelMatrixTools.hxx"
  32
  33 #include <set>
  34 #include <list>
  35 #include <limits>
  36 #include <sstream>
  37 #include <algorithm>
  38 #include <functional>
  39
  40 using namespace ParaMEDMEM;
  41
  42 const double MEDCouplingFieldDiscretization::DFLT_PRECISION=1.e-12;
  43
  44 const char MEDCouplingFieldDiscretizationP0::REPR[]="P0";
  45
  46 const TypeOfField MEDCouplingFieldDiscretizationP0::TYPE=ON_CELLS;
  47
  48 const char MEDCouplingFieldDiscretizationP1::REPR[]="P1";
  49
  50 const TypeOfField MEDCouplingFieldDiscretizationP1::TYPE=ON_NODES;
  51
  52 const int MEDCouplingFieldDiscretizationPerCell::DFT_INVALID_LOCID_VALUE=-1;
  53
  54 const char MEDCouplingFieldDiscretizationGauss::REPR[]="GAUSS";
  55
  56 const TypeOfField MEDCouplingFieldDiscretizationGauss::TYPE=ON_GAUSS_PT;
  57
  58 const char MEDCouplingFieldDiscretizationGaussNE::REPR[]="GSSNE";
  59
  60 const TypeOfField MEDCouplingFieldDiscretizationGaussNE::TYPE=ON_GAUSS_NE;
  61
  62 const char MEDCouplingFieldDiscretizationKriging::REPR[]="KRIGING";
  63
  64 const TypeOfField MEDCouplingFieldDiscretizationKriging::TYPE=ON_NODES_KR;
  65
  66 MEDCouplingFieldDiscretization::MEDCouplingFieldDiscretization():_precision(DFLT_PRECISION)
  67 {
  68 }
  69
  70 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretization::New(TypeOfField type)
  71 {
  72   switch(type)
  73     {
  74     case MEDCouplingFieldDiscretizationP0::TYPE:
  75       return new MEDCouplingFieldDiscretizationP0;
  76     case MEDCouplingFieldDiscretizationP1::TYPE:
  77       return new MEDCouplingFieldDiscretizationP1;
  78     case MEDCouplingFieldDiscretizationGauss::TYPE:
  79       return new MEDCouplingFieldDiscretizationGauss;
  80     case MEDCouplingFieldDiscretizationGaussNE::TYPE:
  81       return new MEDCouplingFieldDiscretizationGaussNE;
  82     case MEDCouplingFieldDiscretizationKriging::TYPE:
  83       return new MEDCouplingFieldDiscretizationKriging;
  84     default:
  85       throw INTERP_KERNEL::Exception("Choosen discretization is not implemented yet.");
  86     }
  87 }
  88
  89 TypeOfField MEDCouplingFieldDiscretization::getTypeOfFieldFromStringRepr(const char *repr) throw(INTERP_KERNEL::Exception)
  90 {
  91   std::string reprCpp(repr);
  92   if(reprCpp==MEDCouplingFieldDiscretizationP0::REPR)
  93     return MEDCouplingFieldDiscretizationP0::TYPE;
  94   if(reprCpp==MEDCouplingFieldDiscretizationP1::REPR)
  95     return MEDCouplingFieldDiscretizationP1::TYPE;
  96   if(reprCpp==MEDCouplingFieldDiscretizationGauss::REPR)
  97     return MEDCouplingFieldDiscretizationGauss::TYPE;
  98   if(reprCpp==MEDCouplingFieldDiscretizationGaussNE::REPR)
  99     return MEDCouplingFieldDiscretizationGaussNE::TYPE;
 100   if(reprCpp==MEDCouplingFieldDiscretizationKriging::REPR)
 101     return MEDCouplingFieldDiscretizationKriging::TYPE;
 102   throw INTERP_KERNEL::Exception("Representation does not match with any field discretization !");
 103 }
 104
 105 bool MEDCouplingFieldDiscretization::isEqual(const MEDCouplingFieldDiscretization *other, double eps) const
 106 {
 107   std::string reason;
 108   return isEqualIfNotWhy(other,eps,reason);
 109 }
 110
 111 bool MEDCouplingFieldDiscretization::isEqualWithoutConsideringStr(const MEDCouplingFieldDiscretization *other, double eps) const
 112 {
 113   return isEqual(other,eps);
 114 }
 115
 116 /*!
 117  * Excepted for MEDCouplingFieldDiscretizationPerCell no underlying TimeLabel object : nothing to do in generally.
 118  */
 119 void MEDCouplingFieldDiscretization::updateTime() const
 120 {
 121 }
 122
 123 /*!
 124  * Computes normL1 of DataArrayDouble instance arr.
 125  * @param res output parameter expected to be of size arr->getNumberOfComponents();
 126  * @throw when the field discretization fails on getMeasure fields (gauss points for example)
 127  */
 128 void MEDCouplingFieldDiscretization::normL1(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, double *res) const throw(INTERP_KERNEL::Exception)
 129 {
 130   MEDCouplingFieldDouble *vol=getMeasureField(mesh,true);
 131   int nbOfCompo=arr->getNumberOfComponents();
 132   int nbOfElems=getNumberOfTuples(mesh);
 133   std::fill(res,res+nbOfCompo,0.);
 134   const double *arrPtr=arr->getConstPointer();
 135   const double *volPtr=vol->getArray()->getConstPointer();
 136   double deno=0.;
 137   for(int i=0;i<nbOfElems;i++)
 138     {
 139       double v=fabs(volPtr[i]);
 140       for(int j=0;j<nbOfCompo;j++)
 141         res[j]+=fabs(arrPtr[i*nbOfCompo+j])*v;
 142       deno+=v;
 143     }
 144   std::transform(res,res+nbOfCompo,res,std::bind2nd(std::multiplies<double>(),1./deno));
 145   vol->decrRef();
 146 }
 147
 148 /*!
 149  * Computes normL2 of DataArrayDouble instance arr.
 150  * @param res output parameter expected to be of size arr->getNumberOfComponents();
 151  * @throw when the field discretization fails on getMeasure fields (gauss points for example)
 152  */
 153 void MEDCouplingFieldDiscretization::normL2(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, double *res) const throw(INTERP_KERNEL::Exception)
 154 {
 155   MEDCouplingFieldDouble *vol=getMeasureField(mesh,true);
 156   int nbOfCompo=arr->getNumberOfComponents();
 157   int nbOfElems=getNumberOfTuples(mesh);
 158   std::fill(res,res+nbOfCompo,0.);
 159   const double *arrPtr=arr->getConstPointer();
 160   const double *volPtr=vol->getArray()->getConstPointer();
 161   double deno=0.;
 162   for(int i=0;i<nbOfElems;i++)
 163     {
 164       double v=fabs(volPtr[i]);
 165       for(int j=0;j<nbOfCompo;j++)
 166         res[j]+=arrPtr[i*nbOfCompo+j]*arrPtr[i*nbOfCompo+j]*v;
 167       deno+=v;
 168     }
 169   std::transform(res,res+nbOfCompo,res,std::bind2nd(std::multiplies<double>(),1./deno));
 170   std::transform(res,res+nbOfCompo,res,std::ptr_fun<double,double>(std::sqrt));
 171   vol->decrRef();
 172 }
 173
 174 /*!
 175  * Computes integral of DataArrayDouble instance arr.
 176  * @param res output parameter expected to be of size arr->getNumberOfComponents();
 177  * @throw when the field discretization fails on getMeasure fields (gauss points for example)
 178  */
 179 void MEDCouplingFieldDiscretization::integral(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, bool isWAbs, double *res) const throw(INTERP_KERNEL::Exception)
 180 {
 181   MEDCouplingFieldDouble *vol=getMeasureField(mesh,isWAbs);
 182   int nbOfCompo=arr->getNumberOfComponents();
 183   int nbOfElems=getNumberOfTuples(mesh);
 184   std::fill(res,res+nbOfCompo,0.);
 185   const double *arrPtr=arr->getConstPointer();
 186   const double *volPtr=vol->getArray()->getConstPointer();
 187   double *tmp=new double[nbOfCompo];
 188   for (int i=0;i<nbOfElems;i++)
 189     {
 190       std::transform(arrPtr+i*nbOfCompo,arrPtr+(i+1)*nbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),volPtr[i]));
 191       std::transform(tmp,tmp+nbOfCompo,res,res,std::plus<double>());
 192     }
 193   delete [] tmp;
 194   vol->decrRef();
 195 }
 196
 197 void MEDCouplingFieldDiscretization::getSerializationIntArray(DataArrayInt *& arr) const
 198 {
 199   arr=0;
 200 }
 201
 202 /*!
 203  * Empty : Not a bug
 204  */
 205 void MEDCouplingFieldDiscretization::getTinySerializationIntInformation(std::vector<int>& tinyInfo) const
 206 {
 207 }
 208
 209 /*!
 210  * Empty : Not a bug
 211  */
 212 void MEDCouplingFieldDiscretization::getTinySerializationDbleInformation(std::vector<double>& tinyInfo) const
 213 {
 214 }
 215
 216 void MEDCouplingFieldDiscretization::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *& arr)
 217 {
 218   arr=0;
 219 }
 220
 221 /*!
 222  * Empty : Not a bug
 223  */
 224 void MEDCouplingFieldDiscretization::finishUnserialization(const std::vector<double>& tinyInfo)
 225 {
 226 }
 227
 228 /*!
 229  * This method is typically the first step of renumbering. The implementation is empty it is not a bug only gauss is impacted
 230  * virtualy by this method.
 231  */
 232 void MEDCouplingFieldDiscretization::renumberCells(const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
 233 {
 234 }
 235
 236 double MEDCouplingFieldDiscretization::getIJK(const MEDCouplingMesh *mesh, const DataArrayDouble *da,
 237                                               int cellId, int nodeIdInCell, int compoId) const throw(INTERP_KERNEL::Exception)
 238 {
 239   throw INTERP_KERNEL::Exception("getIJK Invalid ! only for GaussPoint and GaussNE discretizations !");
 240 }
 241
 242 void MEDCouplingFieldDiscretization::setGaussLocalizationOnType(const MEDCouplingMesh *m, INTERP_KERNEL::NormalizedCellType type, const std::vector<double>& refCoo,
 243                                                                 const std::vector<double>& gsCoo, const std::vector<double>& wg) throw(INTERP_KERNEL::Exception)
 244 {
 245   throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
 246 }
 247
 248 void MEDCouplingFieldDiscretization::setGaussLocalizationOnCells(const MEDCouplingMesh *m, const int *begin, const int *end, const std::vector<double>& refCoo,
 249                                                                  const std::vector<double>& gsCoo, const std::vector<double>& wg) throw(INTERP_KERNEL::Exception)
 250 {
 251   throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
 252 }
 253
 254 void MEDCouplingFieldDiscretization::clearGaussLocalizations() throw(INTERP_KERNEL::Exception)
 255 {
 256   throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
 257 }
 258
 259 MEDCouplingGaussLocalization& MEDCouplingFieldDiscretization::getGaussLocalization(int locId) throw(INTERP_KERNEL::Exception)
 260 {
 261   throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
 262 }
 263
 264 const MEDCouplingGaussLocalization& MEDCouplingFieldDiscretization::getGaussLocalization(int locId) const throw(INTERP_KERNEL::Exception)
 265 {
 266   throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
 267 }
 268
 269 int MEDCouplingFieldDiscretization::getNbOfGaussLocalization() const throw(INTERP_KERNEL::Exception)
 270 {
 271   throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
 272 }
 273
 274 int MEDCouplingFieldDiscretization::getGaussLocalizationIdOfOneCell(int cellId) const throw(INTERP_KERNEL::Exception)
 275 {
 276   throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
 277 }
 278
 279 int MEDCouplingFieldDiscretization::getGaussLocalizationIdOfOneType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception)
 280 {
 281   throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
 282 }
 283
 284 void MEDCouplingFieldDiscretization::getCellIdsHavingGaussLocalization(int locId, std::vector<int>& cellIds) const throw(INTERP_KERNEL::Exception)
 285 {
 286   throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
 287 }
 288
 289 void MEDCouplingFieldDiscretization::renumberEntitiesFromO2NArr(double eps, const int *old2NewPtr, DataArrayDouble *arr, const char *msg)
 290 {
 291   int oldNbOfElems=arr->getNumberOfTuples();
 292   int nbOfComp=arr->getNumberOfComponents();
 293   int newNbOfTuples=(*std::max_element(old2NewPtr,old2NewPtr+oldNbOfElems))+1;
 294   DataArrayDouble *arrCpy=arr->deepCpy();
 295   const double *ptSrc=arrCpy->getConstPointer();
 296   arr->reAlloc(newNbOfTuples);
 297   double *ptToFill=arr->getPointer();
 298   std::fill(ptToFill,ptToFill+nbOfComp*newNbOfTuples,std::numeric_limits<double>::max());
 299   INTERP_KERNEL::AutoPtr<double> tmp=new double[nbOfComp];
 300   for(int i=0;i<oldNbOfElems;i++)
 301     {
 302       int newNb=old2NewPtr[i];
 303       if(newNb>=0)//if newNb<0 the node is considered as out.
 304         {
 305           if(std::find_if(ptToFill+newNb*nbOfComp,ptToFill+(newNb+1)*nbOfComp,std::bind2nd(std::not_equal_to<double>(),std::numeric_limits<double>::max()))
 306              ==ptToFill+(newNb+1)*nbOfComp)
 307             std::copy(ptSrc+i*nbOfComp,ptSrc+(i+1)*nbOfComp,ptToFill+newNb*nbOfComp);
 308           else
 309             {
 310               std::transform(ptSrc+i*nbOfComp,ptSrc+(i+1)*nbOfComp,ptToFill+newNb*nbOfComp,(double *)tmp,std::minus<double>());
 311               std::transform((double *)tmp,((double *)tmp)+nbOfComp,(double *)tmp,std::ptr_fun<double,double>(fabs));
 312               //if(!std::equal(ptSrc+i*nbOfComp,ptSrc+(i+1)*nbOfComp,ptToFill+newNb*nbOfComp))
 313               if(*std::max_element((double *)tmp,((double *)tmp)+nbOfComp)>eps)
 314                 {
 315                   arrCpy->decrRef();
 316                   std::ostringstream oss;
 317                   oss << msg << " " << i << " and " << std::find(old2NewPtr,old2NewPtr+i,newNb)-old2NewPtr
 318                       << " have been merged and " << msg << " field on them are different !";
 319                   throw INTERP_KERNEL::Exception(oss.str().c_str());
 320                 }
 321             }
 322         }
 323     }
 324   arrCpy->decrRef();
 325 }
 326
 327 void MEDCouplingFieldDiscretization::renumberEntitiesFromN2OArr(const int *new2OldPtr, int new2OldSz, DataArrayDouble *arr, const char *msg)
 328 {
 329   int nbOfComp=arr->getNumberOfComponents();
 330   DataArrayDouble *arrCpy=arr->deepCpy();
 331   const double *ptSrc=arrCpy->getConstPointer();
 332   arr->reAlloc(new2OldSz);
 333   double *ptToFill=arr->getPointer();
 334   for(int i=0;i<new2OldSz;i++)
 335     {
 336       int oldNb=new2OldPtr[i];
 337       std::copy(ptSrc+oldNb*nbOfComp,ptSrc+(oldNb+1)*nbOfComp,ptToFill+i*nbOfComp);
 338     }
 339   arrCpy->decrRef();
 340 }
 341
 342 MEDCouplingFieldDiscretization::~MEDCouplingFieldDiscretization()
 343 {
 344 }
 345
 346 TypeOfField MEDCouplingFieldDiscretizationP0::getEnum() const
 347 {
 348   return TYPE;
 349 }
 350
 351 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationP0::clone() const
 352 {
 353   return new MEDCouplingFieldDiscretizationP0;
 354 }
 355
 356 std::string MEDCouplingFieldDiscretizationP0::getStringRepr() const
 357 {
 358   return std::string(REPR);
 359 }
 360
 361 const char *MEDCouplingFieldDiscretizationP0::getRepr() const
 362 {
 363   return REPR;
 364 }
 365
 366 bool MEDCouplingFieldDiscretizationP0::isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const
 367 {
 368   const MEDCouplingFieldDiscretizationP0 *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationP0 *>(other);
 369   bool ret=otherC!=0;
 370   if(!ret)
 371     reason="Spatial discrtization of this is ON_CELLS, which is not the case of other.";
 372   return ret;
 373 }
 374
 375 int MEDCouplingFieldDiscretizationP0::getNumberOfTuples(const MEDCouplingMesh *mesh) const
 376 {
 377   return mesh->getNumberOfCells();
 378 }
 379
 380 int MEDCouplingFieldDiscretizationP0::getNumberOfMeshPlaces(const MEDCouplingMesh *mesh) const
 381 {
 382   return mesh->getNumberOfCells();
 383 }
 384
 385 DataArrayInt *MEDCouplingFieldDiscretizationP0::getOffsetArr(const MEDCouplingMesh *mesh) const
 386 {
 387   int nbOfTuples=mesh->getNumberOfCells();
 388   DataArrayInt *ret=DataArrayInt::New();
 389   ret->alloc(nbOfTuples+1,1);
 390   ret->iota(0);
 391   return ret;
 392 }
 393
 394 void MEDCouplingFieldDiscretizationP0::renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector<DataArrayDouble *>& arrays,
 395                                                              const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
 396 {
 397   const int *array=old2NewBg;
 398   if(check)
 399     array=DataArrayInt::CheckAndPreparePermutation(old2NewBg,old2NewBg+mesh->getNumberOfCells());
 400   for(std::vector<DataArrayDouble *>::const_iterator it=arrays.begin();it!=arrays.end();it++)
 401     {
 402       if(*it)
 403         (*it)->renumberInPlace(array);
 404     }
 405   if(check)
 406     delete [] array;
 407 }
 408
 409 DataArrayDouble *MEDCouplingFieldDiscretizationP0::getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const
 410 {
 411   return mesh->getBarycenterAndOwner();
 412 }
 413
 414 void MEDCouplingFieldDiscretizationP0::computeMeshRestrictionFromTupleIds(const MEDCouplingMesh *mesh, const int *partBg, const int *partEnd,
 415                                                                           DataArrayInt *&cellRest)
 416 {
 417   cellRest=DataArrayInt::New();
 418   cellRest->alloc((int)std::distance(partBg,partEnd),1);
 419   std::copy(partBg,partEnd,cellRest->getPointer());
 420 }
 421
 422 void MEDCouplingFieldDiscretizationP0::checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception)
 423 {
 424 }
 425
 426 void MEDCouplingFieldDiscretizationP0::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception)
 427 {
 428   if(mesh->getNumberOfCells()!=da->getNumberOfTuples())
 429     {
 430       std::ostringstream message;
 431       message << "Field on cells invalid because there are " << mesh->getNumberOfCells();
 432       message << " cells in mesh and " << da->getNumberOfTuples() << " tuples in field !";
 433       throw INTERP_KERNEL::Exception(message.str().c_str());
 434     }
 435 }
 436
 437 MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP0::getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const
 438 {
 439   return mesh->getMeasureField(isAbs);
 440 }
 441
 442 void MEDCouplingFieldDiscretizationP0::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
 443 {
 444   int id=mesh->getCellContainingPoint(loc,_precision);
 445   if(id==-1)
 446     throw INTERP_KERNEL::Exception("Specified point is detected outside of mesh : unable to apply P0::getValueOn !");
 447   arr->getTuple(id,res);
 448 }
 449
 450 void MEDCouplingFieldDiscretizationP0::getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, int i, int j, int k, double *res) const
 451 {
 452   const MEDCouplingCMesh *meshC=dynamic_cast<const MEDCouplingCMesh *>(mesh);
 453   if(!meshC)
 454     throw INTERP_KERNEL::Exception("P0::getValueOnPos is only accessible for structured meshes !");
 455   int id=meshC->getCellIdFromPos(i,j,k);
 456   arr->getTuple(id,res);
 457 }
 458
 459 DataArrayDouble *MEDCouplingFieldDiscretizationP0::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, int nbOfPoints) const
 460 {
 461   std::vector<int> elts,eltsIndex;
 462   mesh->getCellsContainingPoints(loc,nbOfPoints,_precision,elts,eltsIndex);
 463   int spaceDim=mesh->getSpaceDimension();
 464   int nbOfComponents=arr->getNumberOfComponents();
 465   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret=DataArrayDouble::New();
 466   ret->alloc(nbOfPoints,nbOfComponents);
 467   double *ptToFill=ret->getPointer();
 468   for(int i=0;i<nbOfPoints;i++,ptToFill+=nbOfComponents)
 469     if(eltsIndex[i+1]-eltsIndex[i]>=1)
 470       arr->getTuple(elts[eltsIndex[i]],ptToFill);
 471     else
 472       {
 473         std::ostringstream oss; oss << "Point #" << i << " with coordinates : (";
 474         std::copy(loc+i*spaceDim,loc+(i+1)*spaceDim,std::ostream_iterator<double>(oss,", "));
 475         oss << ") detected outside mesh : unable to apply P0::getValueOnMulti ! ";
 476         throw INTERP_KERNEL::Exception(oss.str().c_str());
 477       }
 478   ret->incrRef();
 479   return ret;
 480 }
 481
 482 /*!
 483  * Nothing to do. It's not a bug.
 484  */
 485 void MEDCouplingFieldDiscretizationP0::renumberValuesOnNodes(double , const int *, DataArrayDouble *) const
 486 {
 487 }
 488
 489 void MEDCouplingFieldDiscretizationP0::renumberValuesOnCells(double epsOnVals, const MEDCouplingMesh *mesh, const int *old2New, DataArrayDouble *arr) const
 490 {
 491   renumberEntitiesFromO2NArr(epsOnVals,old2New,arr,"Cell");
 492 }
 493
 494 void MEDCouplingFieldDiscretizationP0::renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const int *new2old, int newSz, DataArrayDouble *arr) const
 495 {
 496   renumberEntitiesFromN2OArr(new2old,newSz,arr,"Cell");
 497 }
 498
 499 /*!
 500  * This method returns a tuple ids selection from cell ids selection [start;end).
 501  * This method is called by MEDCouplingFieldDiscretizationP0::buildSubMeshData to return parameter \b di.
 502  * Here for P0 it's very simple !
 503  *
 504  * \return a newly allocated array containing ids to select into the DataArrayDouble of the field.
 505  *
 506  */
 507 DataArrayInt *MEDCouplingFieldDiscretizationP0::computeTupleIdsToSelectFromCellIds(const MEDCouplingMesh *mesh, const int *startCellIds, const int *endCellIds) const
 508 {
 509   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
 510   ret->alloc((int)std::distance(startCellIds,endCellIds),1);
 511   std::copy(startCellIds,endCellIds,ret->getPointer());
 512   ret->incrRef(); return ret;
 513 }
 514
 515 /*!
 516  * This method returns a submesh of 'mesh' instance constituting cell ids contained in array defined as an interval [start;end).
 517  * @param di is an array returned that specifies entity ids (here cells ids) in mesh 'mesh' of entity in returned submesh.
 518  * Example : The first cell id of returned mesh has the (*di)[0] id in 'mesh'
 519  */
 520 MEDCouplingMesh *MEDCouplingFieldDiscretizationP0::buildSubMeshData(const MEDCouplingMesh *mesh, const int *start, const int *end, DataArrayInt *&di) const
 521 {
 522   MEDCouplingMesh *ret=mesh->buildPart(start,end);
 523   di=DataArrayInt::New();
 524   di->alloc((int)std::distance(start,end),1);
 525   int *pt=di->getPointer();
 526   std::copy(start,end,pt);
 527   return ret;
 528 }
 529
 530 int MEDCouplingFieldDiscretizationOnNodes::getNumberOfTuples(const MEDCouplingMesh *mesh) const
 531 {
 532   return mesh->getNumberOfNodes();
 533 }
 534
 535 int MEDCouplingFieldDiscretizationOnNodes::getNumberOfMeshPlaces(const MEDCouplingMesh *mesh) const
 536 {
 537   return mesh->getNumberOfNodes();
 538 }
 539
 540 /*!
 541  * Nothing to do here.
 542  */
 543 void MEDCouplingFieldDiscretizationOnNodes::renumberArraysForCell(const MEDCouplingMesh *, const std::vector<DataArrayDouble *>& arrays,
 544                                                                   const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
 545 {
 546 }
 547
 548 DataArrayInt *MEDCouplingFieldDiscretizationOnNodes::getOffsetArr(const MEDCouplingMesh *mesh) const
 549 {
 550   int nbOfTuples=mesh->getNumberOfNodes();
 551   DataArrayInt *ret=DataArrayInt::New();
 552   ret->alloc(nbOfTuples+1,1);
 553   ret->iota(0);
 554   return ret;
 555 }
 556
 557 DataArrayDouble *MEDCouplingFieldDiscretizationOnNodes::getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const
 558 {
 559   return mesh->getCoordinatesAndOwner();
 560 }
 561
 562 void MEDCouplingFieldDiscretizationOnNodes::computeMeshRestrictionFromTupleIds(const MEDCouplingMesh *mesh, const int *partBg, const int *partEnd,
 563                                                                                DataArrayInt *&cellRest)
 564 {
 565   cellRest=mesh->getCellIdsFullyIncludedInNodeIds(partBg,partEnd);
 566 }
 567
 568 void MEDCouplingFieldDiscretizationOnNodes::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception)
 569 {
 570   if(mesh->getNumberOfNodes()!=da->getNumberOfTuples())
 571     {
 572       std::ostringstream message;
 573       message << "Field on nodes invalid because there are " << mesh->getNumberOfNodes();
 574       message << " nodes in mesh and " << da->getNumberOfTuples() << " tuples in field !";
 575       throw INTERP_KERNEL::Exception(message.str().c_str());
 576     }
 577 }
 578
 579 /*!
 580  * This method returns a submesh of 'mesh' instance constituting cell ids contained in array defined as an interval [start;end).
 581 * @param di is an array returned that specifies entity ids (here nodes ids) in mesh 'mesh' of entity in returned submesh.
 582  * Example : The first node id of returned mesh has the (*di)[0] id in 'mesh'
 583  */
 584 MEDCouplingMesh *MEDCouplingFieldDiscretizationOnNodes::buildSubMeshData(const MEDCouplingMesh *mesh, const int *start, const int *end, DataArrayInt *&di) const
 585 {
 586   MEDCouplingMesh *ret=mesh->buildPartAndReduceNodes(start,end,di);
 587   DataArrayInt *di2=di->invertArrayO2N2N2O(ret->getNumberOfNodes());
 588   di->decrRef();
 589   di=di2;
 590   return ret;
 591 }
 592
 593 /*!
 594  * This method returns a tuple ids selection from cell ids selection [start;end).
 595  * This method is called by MEDCouplingFieldDiscretizationP0::buildSubMeshData to return parameter \b di.
 596  * Here for P1 only nodes fetched by submesh of mesh[startCellIds:endCellIds) is returned !
 597  *
 598  * \return a newly allocated array containing ids to select into the DataArrayDouble of the field.
 599  *
 600  */
 601 DataArrayInt *MEDCouplingFieldDiscretizationOnNodes::computeTupleIdsToSelectFromCellIds(const MEDCouplingMesh *mesh, const int *startCellIds, const int *endCellIds) const
 602 {
 603   if(!mesh)
 604     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP1::computeTupleIdsToSelectFromCellIds : null mesh !");
 605   const MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> umesh=mesh->buildUnstructured();
 606   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> umesh2=static_cast<MEDCouplingUMesh *>(umesh->buildPartOfMySelf(startCellIds,endCellIds,true));
 607   return umesh2->computeFetchedNodeIds();
 608 }
 609
 610 void MEDCouplingFieldDiscretizationOnNodes::renumberValuesOnNodes(double epsOnVals, const int *old2NewPtr, DataArrayDouble *arr) const
 611 {
 612   renumberEntitiesFromO2NArr(epsOnVals,old2NewPtr,arr,"Node");
 613 }
 614
 615 /*!
 616  * Nothing to do it's not a bug.
 617  */
 618 void MEDCouplingFieldDiscretizationOnNodes::renumberValuesOnCells(double epsOnVals, const MEDCouplingMesh *mesh, const int *old2New, DataArrayDouble *arr) const
 619 {
 620 }
 621
 622 /*!
 623  * Nothing to do it's not a bug.
 624  */
 625 void MEDCouplingFieldDiscretizationOnNodes::renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const int *new2old, int newSz, DataArrayDouble *arr) const
 626 {
 627 }
 628
 629 void MEDCouplingFieldDiscretizationOnNodes::getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, int i, int j, int k, double *res) const
 630 {
 631   const MEDCouplingCMesh *meshC=dynamic_cast<const MEDCouplingCMesh *>(mesh);
 632   if(!meshC)
 633     throw INTERP_KERNEL::Exception("OnNodes::getValueOnPos(i,j,k) is only accessible for structured meshes !");
 634   int id=meshC->getNodeIdFromPos(i,j,k);
 635   arr->getTuple(id,res);
 636 }
 637
 638 TypeOfField MEDCouplingFieldDiscretizationP1::getEnum() const
 639 {
 640   return TYPE;
 641 }
 642
 643 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationP1::clone() const
 644 {
 645   return new MEDCouplingFieldDiscretizationP1;
 646 }
 647
 648 std::string MEDCouplingFieldDiscretizationP1::getStringRepr() const
 649 {
 650   return std::string(REPR);
 651 }
 652
 653 const char *MEDCouplingFieldDiscretizationP1::getRepr() const
 654 {
 655   return REPR;
 656 }
 657
 658 bool MEDCouplingFieldDiscretizationP1::isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const
 659 {
 660   const MEDCouplingFieldDiscretizationP1 *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationP1 *>(other);
 661   bool ret=otherC!=0;
 662   if(!ret)
 663     reason="Spatial discrtization of this is ON_NODES, which is not the case of other.";
 664   return ret;
 665 }
 666
 667 void MEDCouplingFieldDiscretizationP1::checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception)
 668 {
 669   if(nat!=ConservativeVolumic)
 670     throw INTERP_KERNEL::Exception("Invalid nature for P1 field  : expected ConservativeVolumic !");
 671 }
 672
 673 MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP1::getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const
 674 {
 675   return mesh->getMeasureFieldOnNode(isAbs);
 676 }
 677
 678 void MEDCouplingFieldDiscretizationP1::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
 679 {
 680   int id=mesh->getCellContainingPoint(loc,_precision);
 681   if(id==-1)
 682     throw INTERP_KERNEL::Exception("Specified point is detected outside of mesh : unable to apply P1::getValueOn !");
 683   INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell(id);
 684   if(type!=INTERP_KERNEL::NORM_SEG2 && type!=INTERP_KERNEL::NORM_TRI3 && type!=INTERP_KERNEL::NORM_TETRA4)
 685     throw INTERP_KERNEL::Exception("P1 getValueOn is not specified for not simplex cells !");
 686   getValueInCell(mesh,id,arr,loc,res);
 687 }
 688
 689 /*!
 690  * This method localizes a point defined by 'loc' in a cell with id 'cellId' into mesh 'mesh'.
 691  * The result is put into res expected to be of size at least arr->getNumberOfComponents()
 692  */
 693 void MEDCouplingFieldDiscretizationP1::getValueInCell(const MEDCouplingMesh *mesh, int cellId, const DataArrayDouble *arr, const double *loc, double *res) const
 694 {
 695   std::vector<int> conn;
 696   std::vector<double> coo;
 697   mesh->getNodeIdsOfCell(cellId,conn);
 698   for(std::vector<int>::const_iterator iter=conn.begin();iter!=conn.end();iter++)
 699     mesh->getCoordinatesOfNode(*iter,coo);
 700   int spaceDim=mesh->getSpaceDimension();
 701   std::size_t nbOfNodes=conn.size();
 702   std::vector<const double *> vec(nbOfNodes);
 703   for(std::size_t i=0;i<nbOfNodes;i++)
 704     vec[i]=&coo[i*spaceDim];
 705   INTERP_KERNEL::AutoPtr<double> tmp=new double[nbOfNodes];
 706   INTERP_KERNEL::barycentric_coords(vec,loc,tmp);
 707   int sz=arr->getNumberOfComponents();
 708   INTERP_KERNEL::AutoPtr<double> tmp2=new double[sz];
 709   std::fill(res,res+sz,0.);
 710   for(std::size_t i=0;i<nbOfNodes;i++)
 711     {
 712       arr->getTuple(conn[i],(double *)tmp2);
 713       std::transform((double *)tmp2,((double *)tmp2)+sz,(double *)tmp2,std::bind2nd(std::multiplies<double>(),tmp[i]));
 714       std::transform(res,res+sz,(double *)tmp2,res,std::plus<double>());
 715     }
 716 }
 717
 718 DataArrayDouble *MEDCouplingFieldDiscretizationP1::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, int nbOfPoints) const
 719 {
 720   std::vector<int> elts,eltsIndex;
 721   mesh->getCellsContainingPoints(loc,nbOfPoints,_precision,elts,eltsIndex);
 722   int spaceDim=mesh->getSpaceDimension();
 723   int nbOfComponents=arr->getNumberOfComponents();
 724   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret=DataArrayDouble::New();
 725   ret->alloc(nbOfPoints,nbOfComponents);
 726   double *ptToFill=ret->getPointer();
 727   for(int i=0;i<nbOfPoints;i++)
 728     if(eltsIndex[i+1]-eltsIndex[i]>=1)
 729       getValueInCell(mesh,elts[eltsIndex[i]],arr,loc+i*spaceDim,ptToFill+i*nbOfComponents);
 730     else
 731       {
 732         std::ostringstream oss; oss << "Point #" << i << " with coordinates : (";
 733         std::copy(loc+i*spaceDim,loc+(i+1)*spaceDim,std::ostream_iterator<double>(oss,", "));
 734         oss << ") detected outside mesh : unable to apply P1::getValueOnMulti ! ";
 735         throw INTERP_KERNEL::Exception(oss.str().c_str());
 736       }
 737   ret->incrRef();
 738   return ret;
 739 }
 740
 741 MEDCouplingFieldDiscretizationPerCell::MEDCouplingFieldDiscretizationPerCell():_discr_per_cell(0)
 742 {
 743 }
 744
 745 MEDCouplingFieldDiscretizationPerCell::~MEDCouplingFieldDiscretizationPerCell()
 746 {
 747   if(_discr_per_cell)
 748     _discr_per_cell->decrRef();
 749 }
 750
 751 MEDCouplingFieldDiscretizationPerCell::MEDCouplingFieldDiscretizationPerCell(const MEDCouplingFieldDiscretizationPerCell& other):_discr_per_cell(0)
 752 {
 753   DataArrayInt *arr=other._discr_per_cell;
 754   if(arr)
 755     _discr_per_cell=arr->deepCpy();
 756 }
 757
 758 void MEDCouplingFieldDiscretizationPerCell::updateTime() const
 759 {
 760   if(_discr_per_cell)
 761     updateTimeWith(*_discr_per_cell);
 762 }
 763
 764 void MEDCouplingFieldDiscretizationPerCell::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception)
 765 {
 766   if(!_discr_per_cell)
 767     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell has no discretization per cell !");
 768   int nbOfTuples=_discr_per_cell->getNumberOfTuples();
 769   if(nbOfTuples!=mesh->getNumberOfCells())
 770     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell has a discretization per cell but it's not matching the underlying mesh !");
 771 }
 772
 773 bool MEDCouplingFieldDiscretizationPerCell::isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const
 774 {
 775   const MEDCouplingFieldDiscretizationPerCell *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationPerCell *>(other);
 776   if(!otherC)
 777     {
 778       reason="Spatial discrtization of this is ON_GAUSS, which is not the case of other.";
 779       return false;
 780     }
 781   if(_discr_per_cell==0)
 782     return otherC->_discr_per_cell==0;
 783   if(otherC->_discr_per_cell==0)
 784     return false;
 785   bool ret=_discr_per_cell->isEqualIfNotWhy(*otherC->_discr_per_cell,reason);
 786   if(!ret)
 787     reason.insert(0,"Field discretization per cell DataArrayInt given the discid per cell :");
 788   return ret;
 789 }
 790
 791 bool MEDCouplingFieldDiscretizationPerCell::isEqualWithoutConsideringStr(const MEDCouplingFieldDiscretization *other, double eps) const
 792 {
 793   const MEDCouplingFieldDiscretizationPerCell *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationPerCell *>(other);
 794   if(!otherC)
 795     return false;
 796   if(_discr_per_cell==0)
 797     return otherC->_discr_per_cell==0;
 798   if(otherC->_discr_per_cell==0)
 799     return false;
 800   return _discr_per_cell->isEqualWithoutConsideringStr(*otherC->_discr_per_cell);
 801 }
 802
 803 /*!
 804  * This method is typically the first step of renumbering. The impact on _discr_per_cell is necessary here.
 805  * virtualy by this method.
 806  */
 807 void MEDCouplingFieldDiscretizationPerCell::renumberCells(const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
 808 {
 809   int nbCells=_discr_per_cell->getNumberOfTuples();
 810   const int *array=old2NewBg;
 811   if(check)
 812     array=DataArrayInt::CheckAndPreparePermutation(old2NewBg,old2NewBg+nbCells);
 813   //
 814   DataArrayInt *dpc=_discr_per_cell->renumber(array);
 815   _discr_per_cell->decrRef();
 816   _discr_per_cell=dpc;
 817   //
 818   if(check)
 819     delete [] const_cast<int *>(array);
 820 }
 821
 822 void MEDCouplingFieldDiscretizationPerCell::buildDiscrPerCellIfNecessary(const MEDCouplingMesh *m)
 823 {
 824   if(!_discr_per_cell)
 825     {
 826       _discr_per_cell=DataArrayInt::New();
 827       int nbTuples=m->getNumberOfCells();
 828       _discr_per_cell->alloc(nbTuples,1);
 829       int *ptr=_discr_per_cell->getPointer();
 830       std::fill(ptr,ptr+nbTuples,DFT_INVALID_LOCID_VALUE);
 831     }
 832 }
 833
 834 void MEDCouplingFieldDiscretizationPerCell::checkNoOrphanCells() const throw(INTERP_KERNEL::Exception)
 835 {
 836   if(!_discr_per_cell)
 837     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell::checkNoOrphanCells : no discretization defined !");
 838   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> test=_discr_per_cell->getIdsEqual(DFT_INVALID_LOCID_VALUE);
 839   if(test->getNumberOfTuples()!=0)
 840     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell::checkNoOrphanCells : presence of orphan cells !");
 841 }
 842
 843 const DataArrayInt *MEDCouplingFieldDiscretizationPerCell::getArrayOfDiscIds() const
 844 {
 845   return _discr_per_cell;
 846 }
 847
 848 MEDCouplingFieldDiscretizationGauss::MEDCouplingFieldDiscretizationGauss()
 849 {
 850 }
 851
 852 MEDCouplingFieldDiscretizationGauss::MEDCouplingFieldDiscretizationGauss(const MEDCouplingFieldDiscretizationGauss& other):MEDCouplingFieldDiscretizationPerCell(other),_loc(other._loc)
 853 {
 854 }
 855
 856 TypeOfField MEDCouplingFieldDiscretizationGauss::getEnum() const
 857 {
 858   return TYPE;
 859 }
 860
 861 bool MEDCouplingFieldDiscretizationGauss::isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const
 862 {
 863   const MEDCouplingFieldDiscretizationGauss *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationGauss *>(other);
 864   if(!otherC)
 865     {
 866       reason="Spatial discrtization of this is ON_GAUSS, which is not the case of other.";
 867       return false;
 868     }
 869   if(!MEDCouplingFieldDiscretizationPerCell::isEqualIfNotWhy(other,eps,reason))
 870     return false;
 871   if(_loc.size()!=otherC->_loc.size())
 872     {
 873       reason="Gauss spatial discretization : localization sizes differ";
 874       return false;
 875     }
 876   std::size_t sz=_loc.size();
 877   for(std::size_t i=0;i<sz;i++)
 878     if(!_loc[i].isEqual(otherC->_loc[i],eps))
 879       {
 880         std::ostringstream oss; oss << "Gauss spatial discretization : Localization #" << i << " differ from this to other.";
 881         reason=oss.str();
 882         return false;
 883       }
 884   return true;
 885 }
 886
 887 bool MEDCouplingFieldDiscretizationGauss::isEqualWithoutConsideringStr(const MEDCouplingFieldDiscretization *other, double eps) const
 888 {
 889   const MEDCouplingFieldDiscretizationGauss *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationGauss *>(other);
 890   if(!otherC)
 891     return false;
 892   if(!MEDCouplingFieldDiscretizationPerCell::isEqualWithoutConsideringStr(other,eps))
 893     return false;
 894   if(_loc.size()!=otherC->_loc.size())
 895     return false;
 896   std::size_t sz=_loc.size();
 897   for(std::size_t i=0;i<sz;i++)
 898     if(!_loc[i].isEqual(otherC->_loc[i],eps))
 899       return false;
 900   return true;
 901 }
 902
 903 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationGauss::clone() const
 904 {
 905   return new MEDCouplingFieldDiscretizationGauss(*this);
 906 }
 907
 908 std::string MEDCouplingFieldDiscretizationGauss::getStringRepr() const
 909 {
 910   std::ostringstream oss; oss << REPR << "." << std::endl;
 911   if(_discr_per_cell)
 912     {
 913       if(_discr_per_cell->isAllocated())
 914         {
 915           oss << "Discretization per cell : ";
 916           std::copy(_discr_per_cell->begin(),_discr_per_cell->end(),std::ostream_iterator<int>(oss,", "));
 917           oss << std::endl;
 918         }
 919     }
 920   oss << "Presence of " << _loc.size() << " localizations." << std::endl;
 921   int i=0;
 922   for(std::vector<MEDCouplingGaussLocalization>::const_iterator it=_loc.begin();it!=_loc.end();it++,i++)
 923     {
 924       oss << "+++++ Localization #" << i << " +++++" << std::endl;
 925       oss << (*it).getStringRepr();
 926       oss << "++++++++++" << std::endl;
 927     }
 928   return oss.str();
 929 }
 930
 931 const char *MEDCouplingFieldDiscretizationGauss::getRepr() const
 932 {
 933   return REPR;
 934 }
 935
 936 int MEDCouplingFieldDiscretizationGauss::getNumberOfTuples(const MEDCouplingMesh *) const
 937 {
 938   int ret=0;
 939   const int *dcPtr=_discr_per_cell->getConstPointer();
 940   int nbOfTuples=_discr_per_cell->getNumberOfTuples();
 941   for(const int *w=dcPtr;w!=dcPtr+nbOfTuples;w++)
 942     ret+=_loc[*w].getNumberOfGaussPt();
 943   return ret;
 944 }
 945
 946 int MEDCouplingFieldDiscretizationGauss::getNumberOfMeshPlaces(const MEDCouplingMesh *mesh) const
 947 {
 948   return mesh->getNumberOfCells();
 949 }
 950
 951 DataArrayInt *MEDCouplingFieldDiscretizationGauss::getOffsetArr(const MEDCouplingMesh *mesh) const
 952 {
 953   int nbOfTuples=mesh->getNumberOfCells();
 954   DataArrayInt *ret=DataArrayInt::New();
 955   ret->alloc(nbOfTuples+1,1);
 956   int *retPtr=ret->getPointer();
 957   const int *start=_discr_per_cell->getConstPointer();
 958   retPtr[0]=0;
 959   for(int i=0;i<nbOfTuples;i++,start++)
 960     retPtr[i+1]=retPtr[i]+_loc[*start].getNumberOfGaussPt();
 961   return ret;
 962 }
 963
 964 void MEDCouplingFieldDiscretizationGauss::renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector<DataArrayDouble *>& arrays,
 965                                                                 const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
 966 {
 967   const int *array=old2NewBg;
 968   if(check)
 969     array=DataArrayInt::CheckAndPreparePermutation(old2NewBg,old2NewBg+mesh->getNumberOfCells());
 970   int nbOfCells=_discr_per_cell->getNumberOfTuples();
 971   int nbOfTuples=getNumberOfTuples(0);
 972   const int *dcPtr=_discr_per_cell->getConstPointer();
 973   int *array2=new int[nbOfTuples];//stores the final conversion array old2New to give to arrays in renumberInPlace.
 974   int *array3=new int[nbOfCells];//store for each cell in present dcp array (already renumbered) the offset needed by each cell in new numbering.
 975   array3[0]=0;
 976   for(int i=1;i<nbOfCells;i++)
 977     array3[i]=array3[i-1]+_loc[dcPtr[i-1]].getNumberOfGaussPt();
 978   int j=0;
 979   for(int i=0;i<nbOfCells;i++)
 980     {
 981       int nbOfGaussPt=_loc[dcPtr[array[i]]].getNumberOfGaussPt();
 982       for(int k=0;k<nbOfGaussPt;k++,j++)
 983         array2[j]=array3[array[i]]+k;
 984     }
 985   delete [] array3;
 986   for(std::vector<DataArrayDouble *>::const_iterator it=arrays.begin();it!=arrays.end();it++)
 987     if(*it)
 988       (*it)->renumberInPlace(array2);
 989   delete [] array2;
 990   if(check)
 991     delete [] const_cast<int*>(array);
 992 }
 993
 994 DataArrayDouble *MEDCouplingFieldDiscretizationGauss::getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const
 995 {
 996   checkNoOrphanCells();
 997   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> umesh=mesh->buildUnstructured();//in general do nothing
 998   int nbOfTuples=getNumberOfTuples(mesh);
 999   DataArrayDouble *ret=DataArrayDouble::New();
1000   int spaceDim=mesh->getSpaceDimension();
1001   ret->alloc(nbOfTuples,spaceDim);
1002   std::vector< std::vector<int> > locIds;
1003   std::vector<DataArrayInt *> parts=splitIntoSingleGaussDicrPerCellType(locIds);
1004   std::vector< MEDCouplingAutoRefCountObjectPtr<DataArrayInt> > parts2(parts.size());
1005   std::copy(parts.begin(),parts.end(),parts2.begin());
1006   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> offsets=buildNbOfGaussPointPerCellField();
1007   offsets->computeOffsets();
1008   const int *ptrOffsets=offsets->getConstPointer();
1009   const double *coords=umesh->getCoords()->getConstPointer();
1010   const int *connI=umesh->getNodalConnectivityIndex()->getConstPointer();
1011   const int *conn=umesh->getNodalConnectivity()->getConstPointer();
1012   double *valsToFill=ret->getPointer();
1013   for(std::size_t i=0;i<parts2.size();i++)
1014     {
1015       INTERP_KERNEL::GaussCoords calculator;
1016       for(std::vector<int>::const_iterator it=locIds[i].begin();it!=locIds[i].end();it++)
1017         {
1018           const MEDCouplingGaussLocalization& cli=_loc[*it];//curLocInfo
1019           INTERP_KERNEL::NormalizedCellType typ=cli.getType();
1020           const std::vector<double>& wg=cli.getWeights();
1021           calculator.addGaussInfo(typ,INTERP_KERNEL::CellModel::GetCellModel(typ).getDimension(),
1022                                   &cli.getGaussCoords()[0],(int)wg.size(),&cli.getRefCoords()[0],
1023                                   INTERP_KERNEL::CellModel::GetCellModel(typ).getNumberOfNodes());
1024         }
1025       int nbt=parts2[i]->getNumberOfTuples();
1026       for(const int *w=parts2[i]->getConstPointer();w!=parts2[i]->getConstPointer()+nbt;w++)
1027         {
1028           const MEDCouplingGaussLocalization& cli=_loc[*w];
1029           calculator.calculateCoords(cli.getType(),coords,spaceDim,conn+connI[*w]+1,valsToFill+spaceDim*(ptrOffsets[*w]));
1030         }
1031     }
1032   ret->copyStringInfoFrom(*umesh->getCoords());
1033   return ret;
1034 }
1035
1036 void MEDCouplingFieldDiscretizationGauss::computeMeshRestrictionFromTupleIds(const MEDCouplingMesh *mesh, const int *partBg, const int *partEnd,
1037                                                                              DataArrayInt *&cellRest)
1038 {
1039   throw INTERP_KERNEL::Exception("Not implemented yet !");
1040 }
1041
1042 /*!
1043  * Empty : not a bug
1044  */
1045 void MEDCouplingFieldDiscretizationGauss::checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception)
1046 {
1047 }
1048
1049 void MEDCouplingFieldDiscretizationGauss::getTinySerializationIntInformation(std::vector<int>& tinyInfo) const
1050 {
1051   int val=-1;
1052   if(_discr_per_cell)
1053     val=_discr_per_cell->getNumberOfTuples();
1054   tinyInfo.push_back(val);
1055   tinyInfo.push_back((int)_loc.size());
1056   if(_loc.empty())
1057     tinyInfo.push_back(-1);
1058   else
1059     tinyInfo.push_back(_loc[0].getDimension());
1060   for(std::vector<MEDCouplingGaussLocalization>::const_iterator iter=_loc.begin();iter!=_loc.end();iter++)
1061     (*iter).pushTinySerializationIntInfo(tinyInfo);
1062 }
1063
1064 void MEDCouplingFieldDiscretizationGauss::getTinySerializationDbleInformation(std::vector<double>& tinyInfo) const
1065 {
1066   for(std::vector<MEDCouplingGaussLocalization>::const_iterator iter=_loc.begin();iter!=_loc.end();iter++)
1067     (*iter).pushTinySerializationDblInfo(tinyInfo);
1068 }
1069
1070 void MEDCouplingFieldDiscretizationGauss::getSerializationIntArray(DataArrayInt *& arr) const
1071 {
1072   arr=0;
1073   if(_discr_per_cell)
1074     arr=_discr_per_cell;
1075 }
1076
1077 void MEDCouplingFieldDiscretizationGauss::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *& arr)
1078 {
1079   int val=tinyInfo[0];
1080   if(val>=0)
1081     {
1082       _discr_per_cell=DataArrayInt::New();
1083       _discr_per_cell->alloc(val,1);
1084     }
1085   else
1086     _discr_per_cell=0;
1087   arr=_discr_per_cell;
1088   int nbOfLoc=tinyInfo[1];
1089   _loc.clear();
1090   int dim=tinyInfo[2];
1091   int delta=-1;
1092   if(nbOfLoc>0)
1093     delta=((int)tinyInfo.size()-3)/nbOfLoc;
1094   for(int i=0;i<nbOfLoc;i++)
1095     {
1096       std::vector<int> tmp(tinyInfo.begin()+3+i*delta,tinyInfo.begin()+3+(i+1)*delta);
1097       MEDCouplingGaussLocalization elt=MEDCouplingGaussLocalization::BuildNewInstanceFromTinyInfo(dim,tmp);
1098       _loc.push_back(elt);
1099     }
1100 }
1101
1102 void MEDCouplingFieldDiscretizationGauss::finishUnserialization(const std::vector<double>& tinyInfo)
1103 {
1104   double *tmp=new double[tinyInfo.size()];
1105   std::copy(tinyInfo.begin(),tinyInfo.end(),tmp);
1106   const double *work=tmp;
1107   for(std::vector<MEDCouplingGaussLocalization>::iterator iter=_loc.begin();iter!=_loc.end();iter++)
1108     work=(*iter).fillWithValues(work);
1109   delete [] tmp;
1110 }
1111
1112 double MEDCouplingFieldDiscretizationGauss::getIJK(const MEDCouplingMesh *mesh, const DataArrayDouble *da,
1113                                                    int cellId, int nodeIdInCell, int compoId) const throw(INTERP_KERNEL::Exception)
1114 {
1115   int offset=getOffsetOfCell(cellId);
1116   return da->getIJ(offset+nodeIdInCell,compoId);
1117 }
1118
1119 void MEDCouplingFieldDiscretizationGauss::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception)
1120 {
1121   MEDCouplingFieldDiscretizationPerCell::checkCoherencyBetween(mesh,da);
1122   for(std::vector<MEDCouplingGaussLocalization>::const_iterator iter=_loc.begin();iter!=_loc.end();iter++)
1123     (*iter).checkCoherency();
1124   int nbOfDesc=(int)_loc.size();
1125   int nbOfCells=mesh->getNumberOfCells();
1126   const int *dc=_discr_per_cell->getConstPointer();
1127   for(int i=0;i<nbOfCells;i++)
1128     {
1129       if(dc[i]>=nbOfDesc)
1130         {
1131           std::ostringstream oss; oss << "Cell # " << i << " of mesh \"" << mesh->getName() << "\" has an undefined gauss location ! Should never happend !";
1132           throw INTERP_KERNEL::Exception(oss.str().c_str());
1133         }
1134       if(dc[i]<0)
1135         {
1136           std::ostringstream oss; oss << "Cell # " << i << " of mesh \"" << mesh->getName() << "\" has no gauss location !";
1137           throw INTERP_KERNEL::Exception(oss.str().c_str());
1138         }
1139       if(mesh->getTypeOfCell(i)!=_loc[dc[i]].getType())
1140         {
1141           std::ostringstream oss; oss << "Types of mesh and gauss location mismatch for cell # " << i;
1142           throw INTERP_KERNEL::Exception(oss.str().c_str());
1143         }
1144     }
1145   int nbOfTuples=getNumberOfTuples(mesh);
1146   if(nbOfTuples!=da->getNumberOfTuples())
1147     {
1148       std::ostringstream oss; oss << "Invalid number of tuples in the array : expecting " << nbOfTuples << " !";
1149       throw INTERP_KERNEL::Exception(oss.str().c_str());
1150     }
1151 }
1152
1153 MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationGauss::getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const
1154 {
1155   throw INTERP_KERNEL::Exception("Not implemented yet !");
1156 }
1157
1158 void MEDCouplingFieldDiscretizationGauss::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
1159 {
1160   throw INTERP_KERNEL::Exception("Not implemented yet !");
1161 }
1162
1163 void MEDCouplingFieldDiscretizationGauss::getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, int i, int j, int k, double *res) const
1164 {
1165   throw INTERP_KERNEL::Exception("getValueOnPos(i,j,k) : Not applyable for Gauss points !");
1166 }
1167
1168 DataArrayDouble *MEDCouplingFieldDiscretizationGauss::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, int nbOfPoints) const
1169 {
1170   throw INTERP_KERNEL::Exception("getValueOnMulti : Not implemented yet for gauss points !");
1171 }
1172
1173 MEDCouplingMesh *MEDCouplingFieldDiscretizationGauss::buildSubMeshData(const MEDCouplingMesh *mesh, const int *start, const int *end, DataArrayInt *&di) const
1174 {
1175   di=computeTupleIdsToSelectFromCellIds(mesh,start,end);
1176   return mesh->buildPart(start,end);
1177 }
1178
1179 /*!
1180  * This method returns a tuple ids selection from cell ids selection [start;end).
1181  * This method is called by MEDCouplingFieldDiscretizationGauss::buildSubMeshData to return parameter \b di.
1182  *
1183  * \return a newly allocated array containing ids to select into the DataArrayDouble of the field.
1184  *
1185  */
1186 DataArrayInt *MEDCouplingFieldDiscretizationGauss::computeTupleIdsToSelectFromCellIds(const MEDCouplingMesh *mesh, const int *startCellIds, const int *endCellIds) const
1187 {
1188   if(!mesh)
1189     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::computeTupleIdsToSelectFromCellIds : null mesh !");
1190   if(!_discr_per_cell)
1191     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::computeTupleIdsToSelectFromCellIds : null discretization ids !");
1192   int nbOfCells=mesh->getNumberOfCells();
1193   if(_discr_per_cell->getNumberOfTuples()!=nbOfCells)
1194     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::computeTupleIdsToSelectFromCellIds : mismatch of nb of tuples of cell ids array and number of cells !");
1195   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nbOfNodesPerCell=DataArrayInt::New(); nbOfNodesPerCell->alloc(nbOfCells,1);
1196   int *retPtr=nbOfNodesPerCell->getPointer();
1197   const int *pt=_discr_per_cell->getConstPointer();
1198   int nbMaxOfLocId=(int)_loc.size();
1199   for(int i=0;i<nbOfCells;i++,retPtr++)
1200     {
1201       if(*pt>=0 && *pt<nbMaxOfLocId)
1202         *retPtr=_loc[*pt].getNumberOfGaussPt();
1203     }
1204   nbOfNodesPerCell->computeOffsets2();
1205   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> sel=DataArrayInt::New(); sel->useArray(startCellIds,false,CPP_DEALLOC,(int)std::distance(startCellIds,endCellIds),1);
1206   return sel->buildExplicitArrByRanges(nbOfNodesPerCell);
1207 }
1208
1209 /*!
1210  * No implementation needed !
1211  */
1212 void MEDCouplingFieldDiscretizationGauss::renumberValuesOnNodes(double , const int *, DataArrayDouble *) const
1213 {
1214 }
1215
1216 void MEDCouplingFieldDiscretizationGauss::renumberValuesOnCells(double epsOnVals, const MEDCouplingMesh *mesh, const int *old2New, DataArrayDouble *arr) const
1217 {
1218   throw INTERP_KERNEL::Exception("Not implemented yet !");
1219 }
1220
1221 void MEDCouplingFieldDiscretizationGauss::renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const int *new2old, int newSz, DataArrayDouble *arr) const
1222 {
1223   throw INTERP_KERNEL::Exception("Number of cells has changed and becomes higher with some cells that have been split ! Unable to conserve the Gauss field !");
1224 }
1225
1226 void MEDCouplingFieldDiscretizationGauss::setGaussLocalizationOnType(const MEDCouplingMesh *m, INTERP_KERNEL::NormalizedCellType type, const std::vector<double>& refCoo,
1227                                                                      const std::vector<double>& gsCoo, const std::vector<double>& wg) throw(INTERP_KERNEL::Exception)
1228 {
1229   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
1230   if((int)cm.getDimension()!=m->getMeshDimension())
1231     {
1232       std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGauss::setGaussLocalizationOnType : mismatch of dimensions ! MeshDim==" << m->getMeshDimension();
1233       oss << " whereas Type '" << cm.getRepr() << "' has dimension " << cm.getDimension() << " !";
1234       throw INTERP_KERNEL::Exception(oss.str().c_str());
1235     }
1236   buildDiscrPerCellIfNecessary(m);
1237   int id=(int)_loc.size();
1238   MEDCouplingGaussLocalization elt(type,refCoo,gsCoo,wg);
1239   _loc.push_back(elt);
1240   int *ptr=_discr_per_cell->getPointer();
1241   int nbCells=m->getNumberOfCells();
1242   for(int i=0;i<nbCells;i++)
1243     if(m->getTypeOfCell(i)==type)
1244       ptr[i]=id;
1245   zipGaussLocalizations();
1246 }
1247
1248 void MEDCouplingFieldDiscretizationGauss::setGaussLocalizationOnCells(const MEDCouplingMesh *m, const int *begin, const int *end, const std::vector<double>& refCoo,
1249                                                                       const std::vector<double>& gsCoo, const std::vector<double>& wg) throw(INTERP_KERNEL::Exception)
1250 {
1251   buildDiscrPerCellIfNecessary(m);
1252   if(std::distance(begin,end)<1)
1253     throw INTERP_KERNEL::Exception("Size of [begin,end) must be equal or greater than 1 !");
1254   INTERP_KERNEL::NormalizedCellType type=m->getTypeOfCell(*begin);
1255   MEDCouplingGaussLocalization elt(type,refCoo,gsCoo,wg);
1256   int id=(int)_loc.size();
1257   int *ptr=_discr_per_cell->getPointer();
1258   for(const int *w=begin+1;w!=end;w++)
1259     {
1260       if(m->getTypeOfCell(*w)!=type)
1261         {
1262           std::ostringstream oss; oss << "The cell with id " << *w << " has been detected to be incompatible in the [begin,end) array specified !";
1263           throw INTERP_KERNEL::Exception(oss.str().c_str());
1264         }
1265     }
1266   //
1267   for(const int *w2=begin;w2!=end;w2++)
1268     ptr[*w2]=id;
1269   //
1270   _loc.push_back(elt);
1271   zipGaussLocalizations();
1272 }
1273
1274 void MEDCouplingFieldDiscretizationGauss::clearGaussLocalizations() throw(INTERP_KERNEL::Exception)
1275 {
1276   if(_discr_per_cell)
1277     {
1278       _discr_per_cell->decrRef();
1279       _discr_per_cell=0;
1280     }
1281   _loc.clear();
1282 }
1283
1284 MEDCouplingGaussLocalization& MEDCouplingFieldDiscretizationGauss::getGaussLocalization(int locId) throw(INTERP_KERNEL::Exception)
1285 {
1286   checkLocalizationId(locId);
1287   return _loc[locId];
1288 }
1289
1290 int MEDCouplingFieldDiscretizationGauss::getNbOfGaussLocalization() const throw(INTERP_KERNEL::Exception)
1291 {
1292   return (int)_loc.size();
1293 }
1294
1295 int MEDCouplingFieldDiscretizationGauss::getGaussLocalizationIdOfOneCell(int cellId) const throw(INTERP_KERNEL::Exception)
1296 {
1297   if(!_discr_per_cell)
1298     throw INTERP_KERNEL::Exception("No Gauss localization still set !");
1299   int locId=_discr_per_cell->getConstPointer()[cellId];
1300   if(locId<0)
1301     throw INTERP_KERNEL::Exception("No Gauss localization set for the specified cell !");
1302   return locId;
1303 }
1304
1305 int MEDCouplingFieldDiscretizationGauss::getGaussLocalizationIdOfOneType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception)
1306 {
1307   if(!_discr_per_cell)
1308     throw INTERP_KERNEL::Exception("No Gauss localization still set !");
1309   std::set<int> ret;
1310   int id=0;
1311   for(std::vector<MEDCouplingGaussLocalization>::const_iterator iter=_loc.begin();iter!=_loc.end();iter++,id++)
1312     if((*iter).getType()==type)
1313       ret.insert(id);
1314   if(ret.empty())
1315     throw INTERP_KERNEL::Exception("No gauss discretization found for the specified type !");
1316   if(ret.size()>1)
1317     throw INTERP_KERNEL::Exception("Several gauss discretizations have been found for the specified type !");
1318   return *ret.begin();
1319 }
1320
1321 void MEDCouplingFieldDiscretizationGauss::getCellIdsHavingGaussLocalization(int locId, std::vector<int>& cellIds) const throw(INTERP_KERNEL::Exception)
1322 {
1323   if(locId<0 || locId>=(int)_loc.size())
1324     throw INTERP_KERNEL::Exception("Invalid locId given : must be in range [0:getNbOfGaussLocalization()) !");
1325   int nbOfTuples=_discr_per_cell->getNumberOfTuples();
1326   const int *ptr=_discr_per_cell->getConstPointer();
1327   for(int i=0;i<nbOfTuples;i++)
1328     if(ptr[i]==locId)
1329       cellIds.push_back(i);
1330 }
1331
1332 const MEDCouplingGaussLocalization& MEDCouplingFieldDiscretizationGauss::getGaussLocalization(int locId) const throw(INTERP_KERNEL::Exception)
1333 {
1334   checkLocalizationId(locId);
1335   return _loc[locId];
1336 }
1337
1338 void MEDCouplingFieldDiscretizationGauss::checkLocalizationId(int locId) const throw(INTERP_KERNEL::Exception)
1339 {
1340   if(locId<0 || locId>=(int)_loc.size())
1341     throw INTERP_KERNEL::Exception("Invalid locId given : must be in range [0:getNbOfGaussLocalization()) !");
1342 }
1343
1344 int MEDCouplingFieldDiscretizationGauss::getOffsetOfCell(int cellId) const throw(INTERP_KERNEL::Exception)
1345 {
1346   int ret=0;
1347   const int *start=_discr_per_cell->getConstPointer();
1348   for(const int *w=start;w!=start+cellId;w++)
1349     ret+=_loc[*w].getNumberOfGaussPt();
1350   return ret;
1351 }
1352
1353 /*!
1354  * This method do the assumption that there is no orphan cell. If there is an exception is thrown.
1355  * This method makes the assumption too that '_discr_per_cell' is defined. If not an exception is thrown.
1356  * This method returns a newly created array with number of tuples equals to '_discr_per_cell->getNumberOfTuples' and number of components equal to 1.
1357  * The i_th tuple in returned array is the number of gauss point if the corresponding cell.
1358  */
1359 DataArrayInt *MEDCouplingFieldDiscretizationGauss::buildNbOfGaussPointPerCellField() const throw(INTERP_KERNEL::Exception)
1360 {
1361   if(!_discr_per_cell)
1362     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::buildNbOfGaussPointPerCellField : no discretization array set !");
1363   int nbOfTuples=_discr_per_cell->getNumberOfTuples();
1364   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
1365   const int *w=_discr_per_cell->getConstPointer();
1366   ret->alloc(nbOfTuples,1);
1367   int *valsToFill=ret->getPointer();
1368   for(int i=0;i<nbOfTuples;i++,w++)
1369     if(*w!=DFT_INVALID_LOCID_VALUE)
1370       valsToFill[i]=_loc[*w].getNumberOfGaussPt();
1371     else
1372       throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::buildNbOfGaussPointPerCellField : orphan cell detected !");
1373   ret->incrRef();
1374   return ret;
1375 }
1376
1377 /*!
1378  * This method makes the assumption that _discr_per_cell is set.
1379  * This method reduces as much as possible number size of _loc.
1380  * This method is usefull when several set on same cells has been done and that some Gauss Localization are no more used.
1381  */
1382 void MEDCouplingFieldDiscretizationGauss::zipGaussLocalizations()
1383 {
1384   const int *start=_discr_per_cell->getConstPointer();
1385   int nbOfTuples=_discr_per_cell->getNumberOfTuples();
1386   int *tmp=new int[_loc.size()];
1387   std::fill(tmp,tmp+_loc.size(),-2);
1388   for(const int *w=start;w!=start+nbOfTuples;w++)
1389     if(*w>=0)
1390       tmp[*w]=1;
1391   int fid=0;
1392   for(int i=0;i<(int)_loc.size();i++)
1393     if(tmp[i]!=-2)
1394       tmp[i]=fid++;
1395   if(fid==(int)_loc.size())
1396     {//no zip needed
1397       delete [] tmp;
1398       return;
1399     }
1400   // zip needed
1401   int *start2=_discr_per_cell->getPointer();
1402   for(int *w2=start2;w2!=start2+nbOfTuples;w2++)
1403     *w2=tmp[*w2];
1404   std::vector<MEDCouplingGaussLocalization> tmpLoc;
1405   for(int i=0;i<(int)_loc.size();i++)
1406     if(tmp[i]!=-2)
1407       tmpLoc.push_back(_loc[tmp[i]]);
1408   delete [] tmp;
1409   _loc=tmpLoc;
1410 }
1411
1412 /*!
1413  * This method is usefull when 'this' describes a field discretization with several gauss discretization on a \b same cell type.
1414  * For example same NORM_TRI3 cells having 6 gauss points and others with 12 gauss points.
1415  * This method returns 2 arrays with same size : the return value and 'locIds' output parameter.
1416  * For a given i into [0,locIds.size) ret[i] represents the set of cell ids of i_th set an locIds[i] represents the set of discretisation of the set.
1417  * The return vector contains a set of newly created instance to deal with.
1418  * The returned vector represents a \b partition of cells ids with a gauss discretization set.
1419  *
1420  * If no descretization is set in 'this' and exception will be thrown.
1421  */
1422 std::vector<DataArrayInt *> MEDCouplingFieldDiscretizationGauss::splitIntoSingleGaussDicrPerCellType(std::vector< std::vector<int> >& locIds) const throw(INTERP_KERNEL::Exception)
1423 {
1424   if(!_discr_per_cell)
1425     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::splitIntoSingleGaussDicrPerCellType : no descretization set !");
1426   locIds.clear();
1427   std::vector<DataArrayInt *> ret;
1428   const int *discrPerCell=_discr_per_cell->getConstPointer();
1429   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret2=_discr_per_cell->getIdsNotEqual(-1);
1430   int nbOfTuplesSet=ret2->getNumberOfTuples();
1431   std::list<int> idsRemaining(ret2->getConstPointer(),ret2->getConstPointer()+nbOfTuplesSet);
1432   std::list<int>::iterator it=idsRemaining.begin();
1433   while(it!=idsRemaining.end())
1434     {
1435       std::vector<int> ids;
1436       std::set<int> curLocIds;
1437       std::set<INTERP_KERNEL::NormalizedCellType> curCellTypes;
1438       while(it!=idsRemaining.end())
1439         {
1440           int curDiscrId=discrPerCell[*it];
1441           INTERP_KERNEL::NormalizedCellType typ=_loc[curDiscrId].getType();
1442           if(curCellTypes.find(typ)!=curCellTypes.end())
1443             {
1444               if(curLocIds.find(curDiscrId)!=curLocIds.end())
1445                 {
1446                   curLocIds.insert(curDiscrId);
1447                   curCellTypes.insert(typ);
1448                   ids.push_back(*it);
1449                   it=idsRemaining.erase(it);
1450                 }
1451               else
1452                 it++;
1453             }
1454           else
1455             {
1456               curLocIds.insert(curDiscrId);
1457               curCellTypes.insert(typ);
1458               ids.push_back(*it);
1459               it=idsRemaining.erase(it);
1460             }
1461         }
1462       it=idsRemaining.begin();
1463       ret.resize(ret.size()+1);
1464       DataArrayInt *part=DataArrayInt::New();
1465       part->alloc((int)ids.size(),1);
1466       std::copy(ids.begin(),ids.end(),part->getPointer());
1467       ret.back()=part;
1468       locIds.resize(locIds.size()+1);
1469       locIds.back().insert(locIds.back().end(),curLocIds.begin(),curLocIds.end());
1470     }
1471   return ret;
1472 }
1473
1474 MEDCouplingFieldDiscretizationGaussNE::MEDCouplingFieldDiscretizationGaussNE()
1475 {
1476 }
1477
1478 TypeOfField MEDCouplingFieldDiscretizationGaussNE::getEnum() const
1479 {
1480   return TYPE;
1481 }
1482
1483 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationGaussNE::clone() const
1484 {
1485   return new MEDCouplingFieldDiscretizationGaussNE(*this);
1486 }
1487
1488 std::string MEDCouplingFieldDiscretizationGaussNE::getStringRepr() const
1489 {
1490   return std::string(REPR);
1491 }
1492
1493 const char *MEDCouplingFieldDiscretizationGaussNE::getRepr() const
1494 {
1495   return REPR;
1496 }
1497
1498 bool MEDCouplingFieldDiscretizationGaussNE::isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const
1499 {
1500   const MEDCouplingFieldDiscretizationGaussNE *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationGaussNE *>(other);
1501   bool ret=otherC!=0;
1502   if(!ret)
1503     reason="Spatial discrtization of this is ON_GAUSS_NE, which is not the case of other.";
1504   return ret;
1505 }
1506
1507 int MEDCouplingFieldDiscretizationGaussNE::getNumberOfTuples(const MEDCouplingMesh *mesh) const
1508 {
1509   int ret=0;
1510   int nbOfCells=mesh->getNumberOfCells();
1511   for(int i=0;i<nbOfCells;i++)
1512     {
1513       INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell(i);
1514       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
1515       if(cm.isDynamic())
1516         throw INTERP_KERNEL::Exception("Not implemented yet Gauss node on elements for polygons and polyedrons !");
1517       ret+=cm.getNumberOfNodes();
1518     }
1519   return ret;
1520 }
1521
1522 int MEDCouplingFieldDiscretizationGaussNE::getNumberOfMeshPlaces(const MEDCouplingMesh *mesh) const
1523 {
1524   return mesh->getNumberOfCells();
1525 }
1526
1527 DataArrayInt *MEDCouplingFieldDiscretizationGaussNE::getOffsetArr(const MEDCouplingMesh *mesh) const
1528 {
1529   int nbOfTuples=mesh->getNumberOfCells();
1530   DataArrayInt *ret=DataArrayInt::New();
1531   ret->alloc(nbOfTuples+1,1);
1532   int *retPtr=ret->getPointer();
1533   retPtr[0]=0;
1534   for(int i=0;i<nbOfTuples;i++)
1535     {
1536       INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell(i);
1537       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
1538       if(cm.isDynamic())
1539         throw INTERP_KERNEL::Exception("Not implemented yet Gauss node on elements for polygons and polyedrons !");
1540       retPtr[i+1]=retPtr[i]+cm.getNumberOfNodes();
1541     }
1542   return ret;
1543 }
1544
1545 void MEDCouplingFieldDiscretizationGaussNE::renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector<DataArrayDouble *>& arrays,
1546                                                                   const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
1547 {
1548   const int *array=old2NewBg;
1549   if(check)
1550     array=DataArrayInt::CheckAndPreparePermutation(old2NewBg,old2NewBg+mesh->getNumberOfCells());
1551   int nbOfCells=mesh->getNumberOfCells();
1552   int nbOfTuples=getNumberOfTuples(mesh);
1553   int *array2=new int[nbOfTuples];//stores the final conversion array old2New to give to arrays in renumberInPlace.
1554   int *array3=new int[nbOfCells];//store for each cell in after renumbering the offset needed by each cell in new numbering.
1555   array3[0]=0;
1556   for(int i=1;i<nbOfCells;i++)
1557     {
1558       INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell((int)std::distance(array,std::find(array,array+nbOfCells,i-1)));
1559       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
1560       array3[i]=array3[i-1]+cm.getNumberOfNodes();
1561     }
1562   int j=0;
1563   for(int i=0;i<nbOfCells;i++)
1564     {
1565       INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell(i);
1566       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
1567       for(int k=0;k<(int)cm.getNumberOfNodes();k++,j++)
1568         array2[j]=array3[array[i]]+k;
1569     }
1570   delete [] array3;
1571   for(std::vector<DataArrayDouble *>::const_iterator it=arrays.begin();it!=arrays.end();it++)
1572     if(*it)
1573       (*it)->renumberInPlace(array2);
1574   delete [] array2;
1575   if(check)
1576     delete [] const_cast<int *>(array);
1577 }
1578
1579 DataArrayDouble *MEDCouplingFieldDiscretizationGaussNE::getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const
1580 {
1581   throw INTERP_KERNEL::Exception("Not implemented yet !");
1582 }
1583
1584 void MEDCouplingFieldDiscretizationGaussNE::computeMeshRestrictionFromTupleIds(const MEDCouplingMesh *mesh, const int *partBg, const int *partEnd,
1585                                                                                DataArrayInt *&cellRest)
1586 {
1587   throw INTERP_KERNEL::Exception("Not implemented yet !");
1588 }
1589
1590 void MEDCouplingFieldDiscretizationGaussNE::checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception)
1591 {
1592 }
1593
1594 double MEDCouplingFieldDiscretizationGaussNE::getIJK(const MEDCouplingMesh *mesh, const DataArrayDouble *da,
1595                                                      int cellId, int nodeIdInCell, int compoId) const throw(INTERP_KERNEL::Exception)
1596 {
1597   int offset=0;
1598   for(int i=0;i<cellId;i++)
1599     {
1600       INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell(i);
1601       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
1602       offset+=cm.getNumberOfNodes();
1603     }
1604   return da->getIJ(offset+nodeIdInCell,compoId);
1605 }
1606
1607 void MEDCouplingFieldDiscretizationGaussNE::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception)
1608 {
1609   int nbOfTuples=getNumberOfTuples(mesh);
1610   if(nbOfTuples!=da->getNumberOfTuples())
1611     {
1612       std::ostringstream oss; oss << "Invalid number of tuples in the array : expecting " << nbOfTuples << " !";
1613       throw INTERP_KERNEL::Exception(oss.str().c_str());
1614     }
1615 }
1616
1617 MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationGaussNE::getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const
1618 {
1619   throw INTERP_KERNEL::Exception("Not implemented yet !");
1620 }
1621
1622 void MEDCouplingFieldDiscretizationGaussNE::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
1623 {
1624   throw INTERP_KERNEL::Exception("Not implemented yet !");
1625 }
1626
1627 void MEDCouplingFieldDiscretizationGaussNE::getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, int i, int j, int k, double *res) const
1628 {
1629   throw INTERP_KERNEL::Exception("getValueOnPos(i,j,k) : Not applyable for Gauss points !");
1630 }
1631
1632 DataArrayDouble *MEDCouplingFieldDiscretizationGaussNE::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, int nbOfPoints) const
1633 {
1634   throw INTERP_KERNEL::Exception("getValueOnMulti : Not implemented for Gauss NE !");
1635 }
1636
1637 MEDCouplingMesh *MEDCouplingFieldDiscretizationGaussNE::buildSubMeshData(const MEDCouplingMesh *mesh, const int *start, const int *end, DataArrayInt *&di) const
1638 {
1639   di=computeTupleIdsToSelectFromCellIds(mesh,start,end);
1640   return mesh->buildPart(start,end);
1641 }
1642
1643 /*!
1644  * This method returns a tuple ids selection from cell ids selection [start;end).
1645  * This method is called by MEDCouplingFieldDiscretizationGaussNE::buildSubMeshData to return parameter \b di.
1646  *
1647  * \return a newly allocated array containing ids to select into the DataArrayDouble of the field.
1648  *
1649  */
1650 DataArrayInt *MEDCouplingFieldDiscretizationGaussNE::computeTupleIdsToSelectFromCellIds(const MEDCouplingMesh *mesh, const int *startCellIds, const int *endCellIds) const
1651 {
1652   if(!mesh)
1653     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::computeTupleIdsToSelectFromCellIds : null mesh !");
1654   const MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> umesh=mesh->buildUnstructured();
1655   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nbOfNodesPerCell=umesh->computeNbOfNodesPerCell();
1656   nbOfNodesPerCell->computeOffsets2();
1657   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> sel=DataArrayInt::New(); sel->useArray(startCellIds,false,CPP_DEALLOC,(int)std::distance(startCellIds,endCellIds),1);
1658   return sel->buildExplicitArrByRanges(nbOfNodesPerCell);
1659 }
1660
1661 /*!
1662  * No implementation needed !
1663  */
1664 void MEDCouplingFieldDiscretizationGaussNE::renumberValuesOnNodes(double , const int *, DataArrayDouble *) const
1665 {
1666 }
1667
1668 void MEDCouplingFieldDiscretizationGaussNE::renumberValuesOnCells(double epsOnVals, const MEDCouplingMesh *mesh, const int *old2New, DataArrayDouble *arr) const
1669 {
1670   throw INTERP_KERNEL::Exception("Not implemented yet !");
1671 }
1672
1673 void MEDCouplingFieldDiscretizationGaussNE::renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const int *new2old, int newSz, DataArrayDouble *arr) const
1674 {
1675   throw INTERP_KERNEL::Exception("Not implemented yet !");
1676 }
1677
1678 MEDCouplingFieldDiscretizationGaussNE::MEDCouplingFieldDiscretizationGaussNE(const MEDCouplingFieldDiscretizationGaussNE& other):MEDCouplingFieldDiscretization(other)
1679 {
1680 }
1681
1682 TypeOfField MEDCouplingFieldDiscretizationKriging::getEnum() const
1683 {
1684   return TYPE;
1685 }
1686
1687 const char *MEDCouplingFieldDiscretizationKriging::getRepr() const
1688 {
1689   return REPR;
1690 }
1691
1692 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationKriging::clone() const
1693 {
1694   return new MEDCouplingFieldDiscretizationKriging;
1695 }
1696
1697 std::string MEDCouplingFieldDiscretizationKriging::getStringRepr() const
1698 {
1699   return std::string(REPR);
1700 }
1701
1702 void MEDCouplingFieldDiscretizationKriging::checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception)
1703 {
1704   if(nat!=ConservativeVolumic)
1705     throw INTERP_KERNEL::Exception("Invalid nature for Kriging field : expected ConservativeVolumic !");
1706 }
1707
1708 bool MEDCouplingFieldDiscretizationKriging::isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const
1709 {
1710   const MEDCouplingFieldDiscretizationKriging *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationKriging *>(other);
1711   bool ret=otherC!=0;
1712   if(!ret)
1713     reason="Spatial discrtization of this is ON_NODES_KR, which is not the case of other.";
1714   return ret;
1715 }
1716
1717 MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationKriging::getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const
1718 {
1719   throw INTERP_KERNEL::Exception("getMeasureField on FieldDiscretizationKriging : not implemented yet !");
1720 }
1721
1722 void MEDCouplingFieldDiscretizationKriging::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
1723 {
1724   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> res2=MEDCouplingFieldDiscretizationKriging::getValueOnMulti(arr,mesh,loc,1);
1725   std::copy(res2->begin(),res2->end(),res);
1726 }
1727
1728 DataArrayDouble *MEDCouplingFieldDiscretizationKriging::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, int nbOfTargetPoints) const
1729 {
1730   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coords=getLocalizationOfDiscValues(mesh);
1731   int nbOfPts=coords->getNumberOfTuples();
1732   int dimension=coords->getNumberOfComponents();
1733   //
1734   int delta=0;
1735   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> KnewiK=computeVectorOfCoefficients(mesh,arr,delta);
1736   //
1737   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> locArr=DataArrayDouble::New();
1738   locArr->useArray(loc,false,CPP_DEALLOC,nbOfTargetPoints,dimension);
1739   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrix2=coords->buildEuclidianDistanceDenseMatrixWith(locArr);
1740   operateOnDenseMatrix(mesh->getSpaceDimension(),nbOfPts*nbOfTargetPoints,matrix2->getPointer());
1741   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrix3=DataArrayDouble::New();
1742   matrix3->alloc((nbOfPts+delta)*nbOfTargetPoints,1);
1743   double *work=matrix3->getPointer();
1744   const double *workCst=matrix2->getConstPointer();
1745   const double *workCst2=loc;
1746   for(int i=0;i<nbOfTargetPoints;i++,workCst+=nbOfPts,workCst2+=delta-1)
1747     {
1748       for(int j=0;j<nbOfPts;j++)
1749         work[j*nbOfTargetPoints+i]=workCst[j];
1750       work[nbOfPts*nbOfTargetPoints+i]=1.0;
1751       for(int j=0;j<delta-1;j++)
1752         work[(nbOfPts+1+j)*nbOfTargetPoints+i]=workCst2[j];
1753     }
1754   //
1755   int nbOfCompo=arr->getNumberOfComponents();
1756   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret=DataArrayDouble::New();
1757   ret->alloc(nbOfTargetPoints,nbOfCompo);
1758   INTERP_KERNEL::matrixProduct(KnewiK->getConstPointer(),1,nbOfPts+delta,matrix3->getConstPointer(),nbOfPts+delta,nbOfTargetPoints*nbOfCompo,ret->getPointer());
1759   ret->incrRef();
1760   return ret;
1761 }
1762
1763 /*!
1764  * This method computes coefficients to apply to each representing points of \a mesh, that is to say the nodes of \a mesh given a field array \a arr whose
1765  * number of tuples should be equal to the number of representing points in \a mesh.
1766  *
1767  * \param [in] mesh is the sources of nodes on which kriging will be done regarding the parameters and the value of \c this->getSpaceDimension()
1768  * \param [in] arr input field DataArrayDouble whose number of tuples must be equal to the number of nodes in \a mesh
1769  * \param [out] isDrift return if drift coefficients are present in the returned vector of coefficients, and if. If different from 0 there is presence of drift coefficients.
1770  *              Whatever the value of \a isDrift the number of tuples of returned DataArrayDouble  will be equal to \c arr->getNumberOfTuples() + \a isDrift.
1771  * \return a newly allocated array containing coefficients including or not drift coefficient at the end depending the value of \a isDrift parameter.
1772  */
1773 DataArrayDouble *MEDCouplingFieldDiscretizationKriging::computeVectorOfCoefficients(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, int& isDrift) const
1774 {
1775   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coords=getLocalizationOfDiscValues(mesh);
1776   int nbOfPts=coords->getNumberOfTuples();
1777   int dimension=coords->getNumberOfComponents();
1778   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrix=coords->buildEuclidianDistanceDenseMatrix();
1779   operateOnDenseMatrix(mesh->getSpaceDimension(),nbOfPts*nbOfPts,matrix->getPointer());
1780   // Drift
1781   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixWithDrift=performDrift(matrix,coords,isDrift);
1782   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixInv=DataArrayDouble::New();
1783   matrixInv->alloc((nbOfPts+isDrift)*(nbOfPts+isDrift),1);
1784   INTERP_KERNEL::inverseMatrix(matrixWithDrift->getConstPointer(),nbOfPts+isDrift,matrixInv->getPointer());
1785   //
1786   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> KnewiK=DataArrayDouble::New();
1787   KnewiK->alloc((nbOfPts+isDrift)*1,1);
1788   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr2=DataArrayDouble::New();
1789   arr2->alloc((nbOfPts+isDrift)*1,1);
1790   double *work=std::copy(arr->begin(),arr->end(),arr2->getPointer());
1791   std::fill(work,work+isDrift,0.);
1792   INTERP_KERNEL::matrixProduct(matrixInv->getConstPointer(),nbOfPts+isDrift,nbOfPts+isDrift,arr2->getConstPointer(),nbOfPts+isDrift,1,KnewiK->getPointer());
1793   KnewiK->incrRef();
1794   return KnewiK;
1795 }
1796
1797 /*!
1798  * Apply \f f(x) on each element x in \a matrixPtr. \a matrixPtr is expected to be a dense matrix represented by a chunck of memory of size at least equal to \a nbOfElems.
1799  *
1800  * \param [in] spaceDimension space dimension of the input mesh on which the Kriging has to be performed
1801  * \param [in] nbOfElems is the result of the product of nb of rows and the nb of columns of matrix \a matrixPtr
1802  * \param [in,out] matrixPtr is the dense matrix whose on each values the operation will be applied
1803  */
1804 void MEDCouplingFieldDiscretizationKriging::operateOnDenseMatrix(int spaceDimension, int nbOfElems, double *matrixPtr) const
1805 {
1806   switch(spaceDimension)
1807     {
1808     case 1:
1809       {
1810         for(int i=0;i<nbOfElems;i++)
1811           {
1812             double val=matrixPtr[i];
1813             matrixPtr[i]=val*val*val;
1814           }
1815         break;
1816       }
1817     default:
1818       throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::operateOnDenseMatrix : only dimension 1 implemented !");
1819     }
1820 }
1821
1822 /*!
1823  * Starting from a square matrix \a matr, this method returns a newly allocated dense square matrix whose \a matr is included in returned matrix
1824  * in the top left corner, and in the remaining returned matrix the parameters to take into account about the kriging drift.
1825  * For the moment only linear srift is implemented.
1826  *
1827  * \param [in] arr the position of points were input mesh geometry is considered for Kriging
1828  * \param [in] matr input matrix whose drift part will be added
1829  * \param [out] delta the difference between the size of the output matrix and the input matrix \a matr.
1830  * \return a newly allocated matrix bigger than input matrix \a matr.
1831  */
1832 DataArrayDouble *MEDCouplingFieldDiscretizationKriging::performDrift(const DataArrayDouble *matr, const DataArrayDouble *arr, int& delta) const
1833 {
1834   int spaceDimension=arr->getNumberOfComponents();
1835   delta=spaceDimension+1;
1836   int szOfMatrix=arr->getNumberOfTuples();
1837   if(szOfMatrix*szOfMatrix!=matr->getNumberOfTuples())
1838     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::performDrift : invalid size");
1839   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret=DataArrayDouble::New();
1840   ret->alloc((szOfMatrix+delta)*(szOfMatrix+delta),1);
1841   const double *srcWork=matr->getConstPointer();
1842   const double *srcWork2=arr->getConstPointer();
1843   double *destWork=ret->getPointer();
1844   for(int i=0;i<szOfMatrix;i++)
1845     {
1846       destWork=std::copy(srcWork,srcWork+szOfMatrix,destWork);
1847       srcWork+=szOfMatrix;
1848       *destWork++=1.;
1849       destWork=std::copy(srcWork2,srcWork2+spaceDimension,destWork);
1850       srcWork2+=spaceDimension;
1851     }
1852   std::fill(destWork,destWork+szOfMatrix,1.); destWork+=szOfMatrix;
1853   std::fill(destWork,destWork+spaceDimension+1,0.); destWork+=spaceDimension+1;
1854   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arrNoI=arr->toNoInterlace();
1855   srcWork2=arrNoI->getConstPointer();
1856   for(int i=0;i<spaceDimension;i++)
1857     {
1858       destWork=std::copy(srcWork2,srcWork2+szOfMatrix,destWork);
1859       srcWork2+=szOfMatrix;
1860       std::fill(destWork,destWork+spaceDimension+1,0.);
1861       destWork+=spaceDimension;
1862     }
1863   //
1864   ret->incrRef();
1865   return ret;
1866 }
1867