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