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
 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
 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
 992 {
 993   int ret=0;
 994   const int *dcPtr=_discr_per_cell->getConstPointer();
 995   int nbOfTuples=_discr_per_cell->getNumberOfTuples();
 996   for(const int *w=dcPtr;w!=dcPtr+nbOfTuples;w++)
 997     ret+=_loc[*w].getNumberOfGaussPt();
 998   return ret;
 999 }
1000
1001 int MEDCouplingFieldDiscretizationGauss::getNumberOfMeshPlaces(const MEDCouplingMesh *mesh) const
1002 {
1003   return mesh->getNumberOfCells();
1004 }
1005
1006 DataArrayInt *MEDCouplingFieldDiscretizationGauss::getOffsetArr(const MEDCouplingMesh *mesh) const
1007 {
1008   int nbOfTuples=mesh->getNumberOfCells();
1009   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
1010   ret->alloc(nbOfTuples+1,1);
1011   int *retPtr=ret->getPointer();
1012   const int *start=_discr_per_cell->getConstPointer();
1013   if(_discr_per_cell->getNumberOfTuples()!=nbOfTuples)
1014     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getOffsetArr : mismatch between the mesh and the discretization ids array length !");
1015   int maxPossible=(int)_loc.size();
1016   retPtr[0]=0;
1017   for(int i=0;i<nbOfTuples;i++,start++)
1018     {
1019       if(*start>=0 && *start<maxPossible)
1020         retPtr[i+1]=retPtr[i]+_loc[*start].getNumberOfGaussPt();
1021       else
1022         {
1023           std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGauss::getOffsetArr : At position #" << i << " the locid = " << *start << " whereas it should be in [0," << maxPossible << ") !";
1024           throw INTERP_KERNEL::Exception(oss.str().c_str());
1025         }
1026     }
1027   return ret.retn();
1028 }
1029
1030 void MEDCouplingFieldDiscretizationGauss::renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector<DataArrayDouble *>& arrays,
1031                                                                 const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
1032 {
1033   const int *array=old2NewBg;
1034   if(check)
1035     array=DataArrayInt::CheckAndPreparePermutation(old2NewBg,old2NewBg+mesh->getNumberOfCells());
1036   int nbOfCells=_discr_per_cell->getNumberOfTuples();
1037   int nbOfTuples=getNumberOfTuples(0);
1038   const int *dcPtr=_discr_per_cell->getConstPointer();
1039   int *array2=new int[nbOfTuples];//stores the final conversion array old2New to give to arrays in renumberInPlace.
1040   int *array3=new int[nbOfCells];//store for each cell in present dcp array (already renumbered) the offset needed by each cell in new numbering.
1041   array3[0]=0;
1042   for(int i=1;i<nbOfCells;i++)
1043     array3[i]=array3[i-1]+_loc[dcPtr[i-1]].getNumberOfGaussPt();
1044   int j=0;
1045   for(int i=0;i<nbOfCells;i++)
1046     {
1047       int nbOfGaussPt=_loc[dcPtr[array[i]]].getNumberOfGaussPt();
1048       for(int k=0;k<nbOfGaussPt;k++,j++)
1049         array2[j]=array3[array[i]]+k;
1050     }
1051   delete [] array3;
1052   for(std::vector<DataArrayDouble *>::const_iterator it=arrays.begin();it!=arrays.end();it++)
1053     if(*it)
1054       (*it)->renumberInPlace(array2);
1055   delete [] array2;
1056   if(check)
1057     delete [] const_cast<int*>(array);
1058 }
1059
1060 DataArrayDouble *MEDCouplingFieldDiscretizationGauss::getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const
1061 {
1062   checkNoOrphanCells();
1063   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> umesh=mesh->buildUnstructured();//in general do nothing
1064   int nbOfTuples=getNumberOfTuples(mesh);
1065   DataArrayDouble *ret=DataArrayDouble::New();
1066   int spaceDim=mesh->getSpaceDimension();
1067   ret->alloc(nbOfTuples,spaceDim);
1068   std::vector< int > locIds;
1069   std::vector<DataArrayInt *> parts=splitIntoSingleGaussDicrPerCellType(locIds);
1070   std::vector< MEDCouplingAutoRefCountObjectPtr<DataArrayInt> > parts2(parts.size());
1071   std::copy(parts.begin(),parts.end(),parts2.begin());
1072   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> offsets=buildNbOfGaussPointPerCellField();
1073   offsets->computeOffsets();
1074   const int *ptrOffsets=offsets->getConstPointer();
1075   const double *coords=umesh->getCoords()->getConstPointer();
1076   const int *connI=umesh->getNodalConnectivityIndex()->getConstPointer();
1077   const int *conn=umesh->getNodalConnectivity()->getConstPointer();
1078   double *valsToFill=ret->getPointer();
1079   for(std::size_t i=0;i<parts2.size();i++)
1080     {
1081       INTERP_KERNEL::GaussCoords calculator;
1082       //
1083       const MEDCouplingGaussLocalization& cli=_loc[locIds[i]];//curLocInfo
1084       INTERP_KERNEL::NormalizedCellType typ=cli.getType();
1085       const std::vector<double>& wg=cli.getWeights();
1086       calculator.addGaussInfo(typ,INTERP_KERNEL::CellModel::GetCellModel(typ).getDimension(),
1087                                   &cli.getGaussCoords()[0],(int)wg.size(),&cli.getRefCoords()[0],
1088                                   INTERP_KERNEL::CellModel::GetCellModel(typ).getNumberOfNodes());
1089       //
1090       int nbt=parts2[i]->getNumberOfTuples();
1091       for(const int *w=parts2[i]->getConstPointer();w!=parts2[i]->getConstPointer()+nbt;w++)
1092         calculator.calculateCoords(cli.getType(),coords,spaceDim,conn+connI[*w]+1,valsToFill+spaceDim*(ptrOffsets[*w]));
1093     }
1094   ret->copyStringInfoFrom(*umesh->getCoords());
1095   return ret;
1096 }
1097
1098 void MEDCouplingFieldDiscretizationGauss::computeMeshRestrictionFromTupleIds(const MEDCouplingMesh *mesh, const int *partBg, const int *partEnd,
1099                                                                              DataArrayInt *&cellRest)
1100 {
1101   throw INTERP_KERNEL::Exception("Not implemented yet !");
1102 }
1103
1104 /*!
1105  * Empty : not a bug
1106  */
1107 void MEDCouplingFieldDiscretizationGauss::checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception)
1108 {
1109 }
1110
1111 void MEDCouplingFieldDiscretizationGauss::getTinySerializationIntInformation(std::vector<int>& tinyInfo) const
1112 {
1113   int val=-1;
1114   if(_discr_per_cell)
1115     val=_discr_per_cell->getNumberOfTuples();
1116   tinyInfo.push_back(val);
1117   tinyInfo.push_back((int)_loc.size());
1118   if(_loc.empty())
1119     tinyInfo.push_back(-1);
1120   else
1121     tinyInfo.push_back(_loc[0].getDimension());
1122   for(std::vector<MEDCouplingGaussLocalization>::const_iterator iter=_loc.begin();iter!=_loc.end();iter++)
1123     (*iter).pushTinySerializationIntInfo(tinyInfo);
1124 }
1125
1126 void MEDCouplingFieldDiscretizationGauss::getTinySerializationDbleInformation(std::vector<double>& tinyInfo) const
1127 {
1128   for(std::vector<MEDCouplingGaussLocalization>::const_iterator iter=_loc.begin();iter!=_loc.end();iter++)
1129     (*iter).pushTinySerializationDblInfo(tinyInfo);
1130 }
1131
1132 void MEDCouplingFieldDiscretizationGauss::getSerializationIntArray(DataArrayInt *& arr) const
1133 {
1134   arr=0;
1135   if(_discr_per_cell)
1136     arr=_discr_per_cell;
1137 }
1138
1139 void MEDCouplingFieldDiscretizationGauss::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *& arr)
1140 {
1141   int val=tinyInfo[0];
1142   if(val>=0)
1143     {
1144       _discr_per_cell=DataArrayInt::New();
1145       _discr_per_cell->alloc(val,1);
1146     }
1147   else
1148     _discr_per_cell=0;
1149   arr=_discr_per_cell;
1150   int nbOfLoc=tinyInfo[1];
1151   _loc.clear();
1152   int dim=tinyInfo[2];
1153   int delta=-1;
1154   if(nbOfLoc>0)
1155     delta=((int)tinyInfo.size()-3)/nbOfLoc;
1156   for(int i=0;i<nbOfLoc;i++)
1157     {
1158       std::vector<int> tmp(tinyInfo.begin()+3+i*delta,tinyInfo.begin()+3+(i+1)*delta);
1159       MEDCouplingGaussLocalization elt=MEDCouplingGaussLocalization::BuildNewInstanceFromTinyInfo(dim,tmp);
1160       _loc.push_back(elt);
1161     }
1162 }
1163
1164 void MEDCouplingFieldDiscretizationGauss::finishUnserialization(const std::vector<double>& tinyInfo)
1165 {
1166   double *tmp=new double[tinyInfo.size()];
1167   std::copy(tinyInfo.begin(),tinyInfo.end(),tmp);
1168   const double *work=tmp;
1169   for(std::vector<MEDCouplingGaussLocalization>::iterator iter=_loc.begin();iter!=_loc.end();iter++)
1170     work=(*iter).fillWithValues(work);
1171   delete [] tmp;
1172 }
1173
1174 double MEDCouplingFieldDiscretizationGauss::getIJK(const MEDCouplingMesh *mesh, const DataArrayDouble *da,
1175                                                    int cellId, int nodeIdInCell, int compoId) const throw(INTERP_KERNEL::Exception)
1176 {
1177   int offset=getOffsetOfCell(cellId);
1178   return da->getIJ(offset+nodeIdInCell,compoId);
1179 }
1180
1181 void MEDCouplingFieldDiscretizationGauss::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception)
1182 {
1183   MEDCouplingFieldDiscretizationPerCell::checkCoherencyBetween(mesh,da);
1184   for(std::vector<MEDCouplingGaussLocalization>::const_iterator iter=_loc.begin();iter!=_loc.end();iter++)
1185     (*iter).checkCoherency();
1186   int nbOfDesc=(int)_loc.size();
1187   int nbOfCells=mesh->getNumberOfCells();
1188   const int *dc=_discr_per_cell->getConstPointer();
1189   for(int i=0;i<nbOfCells;i++)
1190     {
1191       if(dc[i]>=nbOfDesc)
1192         {
1193           std::ostringstream oss; oss << "Cell # " << i << " of mesh \"" << mesh->getName() << "\" has an undefined gauss location ! Should never happend !";
1194           throw INTERP_KERNEL::Exception(oss.str().c_str());
1195         }
1196       if(dc[i]<0)
1197         {
1198           std::ostringstream oss; oss << "Cell # " << i << " of mesh \"" << mesh->getName() << "\" has no gauss location !";
1199           throw INTERP_KERNEL::Exception(oss.str().c_str());
1200         }
1201       if(mesh->getTypeOfCell(i)!=_loc[dc[i]].getType())
1202         {
1203           std::ostringstream oss; oss << "Types of mesh and gauss location mismatch for cell # " << i;
1204           throw INTERP_KERNEL::Exception(oss.str().c_str());
1205         }
1206     }
1207   int nbOfTuples=getNumberOfTuples(mesh);
1208   if(nbOfTuples!=da->getNumberOfTuples())
1209     {
1210       std::ostringstream oss; oss << "Invalid number of tuples in the array : expecting " << nbOfTuples << " !";
1211       throw INTERP_KERNEL::Exception(oss.str().c_str());
1212     }
1213 }
1214
1215 MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationGauss::getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const
1216 {
1217   if(!mesh)
1218     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getMeasureField : mesh instance specified is NULL !");
1219   throw INTERP_KERNEL::Exception("Not implemented yet !");
1220 }
1221
1222 void MEDCouplingFieldDiscretizationGauss::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
1223 {
1224   throw INTERP_KERNEL::Exception("Not implemented yet !");
1225 }
1226
1227 void MEDCouplingFieldDiscretizationGauss::getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, int i, int j, int k, double *res) const
1228 {
1229   throw INTERP_KERNEL::Exception("getValueOnPos(i,j,k) : Not applyable for Gauss points !");
1230 }
1231
1232 DataArrayDouble *MEDCouplingFieldDiscretizationGauss::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, int nbOfPoints) const
1233 {
1234   throw INTERP_KERNEL::Exception("getValueOnMulti : Not implemented yet for gauss points !");
1235 }
1236
1237 MEDCouplingMesh *MEDCouplingFieldDiscretizationGauss::buildSubMeshData(const MEDCouplingMesh *mesh, const int *start, const int *end, DataArrayInt *&di) const
1238 {
1239   di=computeTupleIdsToSelectFromCellIds(mesh,start,end);
1240   return mesh->buildPart(start,end);
1241 }
1242
1243 /*!
1244  * This method returns a tuple ids selection from cell ids selection [start;end).
1245  * This method is called by MEDCouplingFieldDiscretizationGauss::buildSubMeshData to return parameter \b di.
1246  *
1247  * \return a newly allocated array containing ids to select into the DataArrayDouble of the field.
1248  *
1249  */
1250 DataArrayInt *MEDCouplingFieldDiscretizationGauss::computeTupleIdsToSelectFromCellIds(const MEDCouplingMesh *mesh, const int *startCellIds, const int *endCellIds) const
1251 {
1252   if(!mesh)
1253     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::computeTupleIdsToSelectFromCellIds : null mesh !");
1254   if(!_discr_per_cell)
1255     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::computeTupleIdsToSelectFromCellIds : null discretization ids !");
1256   int nbOfCells=mesh->getNumberOfCells();
1257   if(_discr_per_cell->getNumberOfTuples()!=nbOfCells)
1258     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::computeTupleIdsToSelectFromCellIds : mismatch of nb of tuples of cell ids array and number of cells !");
1259   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nbOfNodesPerCell=DataArrayInt::New(); nbOfNodesPerCell->alloc(nbOfCells,1);
1260   int *retPtr=nbOfNodesPerCell->getPointer();
1261   const int *pt=_discr_per_cell->getConstPointer();
1262   int nbMaxOfLocId=(int)_loc.size();
1263   for(int i=0;i<nbOfCells;i++,retPtr++,pt++)
1264     {
1265       if(*pt>=0 && *pt<nbMaxOfLocId)
1266         *retPtr=_loc[*pt].getNumberOfGaussPt();
1267     }
1268   nbOfNodesPerCell->computeOffsets2();
1269   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> sel=DataArrayInt::New(); sel->useArray(startCellIds,false,CPP_DEALLOC,(int)std::distance(startCellIds,endCellIds),1);
1270   return sel->buildExplicitArrByRanges(nbOfNodesPerCell);
1271 }
1272
1273 /*!
1274  * No implementation needed !
1275  */
1276 void MEDCouplingFieldDiscretizationGauss::renumberValuesOnNodes(double , const int *, int newNbOfNodes, DataArrayDouble *) const
1277 {
1278 }
1279
1280 void MEDCouplingFieldDiscretizationGauss::renumberValuesOnCells(double epsOnVals, const MEDCouplingMesh *mesh, const int *old2New, int newSz, DataArrayDouble *arr) const
1281 {
1282   throw INTERP_KERNEL::Exception("Not implemented yet !");
1283 }
1284
1285 void MEDCouplingFieldDiscretizationGauss::renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const int *new2old, int newSz, DataArrayDouble *arr) const
1286 {
1287   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 !");
1288 }
1289
1290 void MEDCouplingFieldDiscretizationGauss::setGaussLocalizationOnType(const MEDCouplingMesh *m, INTERP_KERNEL::NormalizedCellType type, const std::vector<double>& refCoo,
1291                                                                      const std::vector<double>& gsCoo, const std::vector<double>& wg) throw(INTERP_KERNEL::Exception)
1292 {
1293   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
1294   if((int)cm.getDimension()!=m->getMeshDimension())
1295     {
1296       std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGauss::setGaussLocalizationOnType : mismatch of dimensions ! MeshDim==" << m->getMeshDimension();
1297       oss << " whereas Type '" << cm.getRepr() << "' has dimension " << cm.getDimension() << " !";
1298       throw INTERP_KERNEL::Exception(oss.str().c_str());
1299     }
1300   buildDiscrPerCellIfNecessary(m);
1301   int id=(int)_loc.size();
1302   MEDCouplingGaussLocalization elt(type,refCoo,gsCoo,wg);
1303   _loc.push_back(elt);
1304   int *ptr=_discr_per_cell->getPointer();
1305   int nbCells=m->getNumberOfCells();
1306   for(int i=0;i<nbCells;i++)
1307     if(m->getTypeOfCell(i)==type)
1308       ptr[i]=id;
1309   zipGaussLocalizations();
1310 }
1311
1312 void MEDCouplingFieldDiscretizationGauss::setGaussLocalizationOnCells(const MEDCouplingMesh *m, const int *begin, const int *end, const std::vector<double>& refCoo,
1313                                                                       const std::vector<double>& gsCoo, const std::vector<double>& wg) throw(INTERP_KERNEL::Exception)
1314 {
1315   buildDiscrPerCellIfNecessary(m);
1316   if(std::distance(begin,end)<1)
1317     throw INTERP_KERNEL::Exception("Size of [begin,end) must be equal or greater than 1 !");
1318   INTERP_KERNEL::NormalizedCellType type=m->getTypeOfCell(*begin);
1319   MEDCouplingGaussLocalization elt(type,refCoo,gsCoo,wg);
1320   int id=(int)_loc.size();
1321   int *ptr=_discr_per_cell->getPointer();
1322   for(const int *w=begin+1;w!=end;w++)
1323     {
1324       if(m->getTypeOfCell(*w)!=type)
1325         {
1326           std::ostringstream oss; oss << "The cell with id " << *w << " has been detected to be incompatible in the [begin,end) array specified !";
1327           throw INTERP_KERNEL::Exception(oss.str().c_str());
1328         }
1329     }
1330   //
1331   for(const int *w2=begin;w2!=end;w2++)
1332     ptr[*w2]=id;
1333   //
1334   _loc.push_back(elt);
1335   zipGaussLocalizations();
1336 }
1337
1338 void MEDCouplingFieldDiscretizationGauss::clearGaussLocalizations() throw(INTERP_KERNEL::Exception)
1339 {
1340   if(_discr_per_cell)
1341     {
1342       _discr_per_cell->decrRef();
1343       _discr_per_cell=0;
1344     }
1345   _loc.clear();
1346 }
1347
1348 MEDCouplingGaussLocalization& MEDCouplingFieldDiscretizationGauss::getGaussLocalization(int locId) throw(INTERP_KERNEL::Exception)
1349 {
1350   checkLocalizationId(locId);
1351   return _loc[locId];
1352 }
1353
1354 int MEDCouplingFieldDiscretizationGauss::getNbOfGaussLocalization() const throw(INTERP_KERNEL::Exception)
1355 {
1356   return (int)_loc.size();
1357 }
1358
1359 int MEDCouplingFieldDiscretizationGauss::getGaussLocalizationIdOfOneCell(int cellId) const throw(INTERP_KERNEL::Exception)
1360 {
1361   if(!_discr_per_cell)
1362     throw INTERP_KERNEL::Exception("No Gauss localization still set !");
1363   int locId=_discr_per_cell->getConstPointer()[cellId];
1364   if(locId<0)
1365     throw INTERP_KERNEL::Exception("No Gauss localization set for the specified cell !");
1366   return locId;
1367 }
1368
1369 int MEDCouplingFieldDiscretizationGauss::getGaussLocalizationIdOfOneType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception)
1370 {
1371   std::set<int> ret=getGaussLocalizationIdsOfOneType(type);
1372   if(ret.empty())
1373     throw INTERP_KERNEL::Exception("No gauss discretization found for the specified type !");
1374   if(ret.size()>1)
1375     throw INTERP_KERNEL::Exception("Several gauss discretizations have been found for the specified type !");
1376   return *ret.begin();
1377 }
1378
1379 std::set<int> MEDCouplingFieldDiscretizationGauss::getGaussLocalizationIdsOfOneType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception)
1380 {
1381   if(!_discr_per_cell)
1382     throw INTERP_KERNEL::Exception("No Gauss localization still set !");
1383   std::set<int> ret;
1384   int id=0;
1385   for(std::vector<MEDCouplingGaussLocalization>::const_iterator iter=_loc.begin();iter!=_loc.end();iter++,id++)
1386     if((*iter).getType()==type)
1387       ret.insert(id);
1388   return ret;
1389 }
1390
1391 void MEDCouplingFieldDiscretizationGauss::getCellIdsHavingGaussLocalization(int locId, std::vector<int>& cellIds) const throw(INTERP_KERNEL::Exception)
1392 {
1393   if(locId<0 || locId>=(int)_loc.size())
1394     throw INTERP_KERNEL::Exception("Invalid locId given : must be in range [0:getNbOfGaussLocalization()) !");
1395   int nbOfTuples=_discr_per_cell->getNumberOfTuples();
1396   const int *ptr=_discr_per_cell->getConstPointer();
1397   for(int i=0;i<nbOfTuples;i++)
1398     if(ptr[i]==locId)
1399       cellIds.push_back(i);
1400 }
1401
1402 const MEDCouplingGaussLocalization& MEDCouplingFieldDiscretizationGauss::getGaussLocalization(int locId) const throw(INTERP_KERNEL::Exception)
1403 {
1404   checkLocalizationId(locId);
1405   return _loc[locId];
1406 }
1407
1408 void MEDCouplingFieldDiscretizationGauss::checkLocalizationId(int locId) const throw(INTERP_KERNEL::Exception)
1409 {
1410   if(locId<0 || locId>=(int)_loc.size())
1411     throw INTERP_KERNEL::Exception("Invalid locId given : must be in range [0:getNbOfGaussLocalization()) !");
1412 }
1413
1414 int MEDCouplingFieldDiscretizationGauss::getOffsetOfCell(int cellId) const throw(INTERP_KERNEL::Exception)
1415 {
1416   int ret=0;
1417   const int *start=_discr_per_cell->getConstPointer();
1418   for(const int *w=start;w!=start+cellId;w++)
1419     ret+=_loc[*w].getNumberOfGaussPt();
1420   return ret;
1421 }
1422
1423 /*!
1424  * This method do the assumption that there is no orphan cell. If there is an exception is thrown.
1425  * This method makes the assumption too that '_discr_per_cell' is defined. If not an exception is thrown.
1426  * This method returns a newly created array with number of tuples equals to '_discr_per_cell->getNumberOfTuples' and number of components equal to 1.
1427  * The i_th tuple in returned array is the number of gauss point if the corresponding cell.
1428  */
1429 DataArrayInt *MEDCouplingFieldDiscretizationGauss::buildNbOfGaussPointPerCellField() const throw(INTERP_KERNEL::Exception)
1430 {
1431   if(!_discr_per_cell)
1432     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::buildNbOfGaussPointPerCellField : no discretization array set !");
1433   int nbOfTuples=_discr_per_cell->getNumberOfTuples();
1434   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
1435   const int *w=_discr_per_cell->getConstPointer();
1436   ret->alloc(nbOfTuples,1);
1437   int *valsToFill=ret->getPointer();
1438   for(int i=0;i<nbOfTuples;i++,w++)
1439     if(*w!=DFT_INVALID_LOCID_VALUE)
1440       valsToFill[i]=_loc[*w].getNumberOfGaussPt();
1441     else
1442       throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::buildNbOfGaussPointPerCellField : orphan cell detected !");
1443   return ret.retn();
1444 }
1445
1446 /*!
1447  * This method makes the assumption that _discr_per_cell is set.
1448  * This method reduces as much as possible number size of _loc.
1449  * This method is usefull when several set on same cells has been done and that some Gauss Localization are no more used.
1450  */
1451 void MEDCouplingFieldDiscretizationGauss::zipGaussLocalizations()
1452 {
1453   const int *start=_discr_per_cell->getConstPointer();
1454   int nbOfTuples=_discr_per_cell->getNumberOfTuples();
1455   INTERP_KERNEL::AutoPtr<int> tmp=new int[_loc.size()];
1456   std::fill((int *)tmp,(int *)tmp+_loc.size(),-2);
1457   for(const int *w=start;w!=start+nbOfTuples;w++)
1458     if(*w>=0)
1459       tmp[*w]=1;
1460   int fid=0;
1461   for(int i=0;i<(int)_loc.size();i++)
1462     if(tmp[i]!=-2)
1463       tmp[i]=fid++;
1464   if(fid==(int)_loc.size())
1465     return;
1466   // zip needed
1467   int *start2=_discr_per_cell->getPointer();
1468   for(int *w2=start2;w2!=start2+nbOfTuples;w2++)
1469     if(*w2>=0)
1470       *w2=tmp[*w2];
1471   std::vector<MEDCouplingGaussLocalization> tmpLoc;
1472   for(int i=0;i<(int)_loc.size();i++)
1473     if(tmp[i]!=-2)
1474       tmpLoc.push_back(_loc[tmp[i]]);
1475   _loc=tmpLoc;
1476 }
1477
1478 /*!
1479  * This method is usefull when 'this' describes a field discretization with several gauss discretization on a \b same cell type.
1480  * For example same NORM_TRI3 cells having 6 gauss points and others with 12 gauss points.
1481  * This method returns 2 arrays with same size : the return value and 'locIds' output parameter.
1482  * 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.
1483  * The return vector contains a set of newly created instance to deal with.
1484  * The returned vector represents a \b partition of cells ids with a gauss discretization set.
1485  *
1486  * If no descretization is set in 'this' and exception will be thrown.
1487  */
1488 std::vector<DataArrayInt *> MEDCouplingFieldDiscretizationGauss::splitIntoSingleGaussDicrPerCellType(std::vector<int>& locIds) const throw(INTERP_KERNEL::Exception)
1489 {
1490   if(!_discr_per_cell)
1491     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::splitIntoSingleGaussDicrPerCellType : no descretization set !");
1492   return _discr_per_cell->partitionByDifferentValues(locIds);
1493 }
1494
1495 MEDCouplingFieldDiscretizationGaussNE::MEDCouplingFieldDiscretizationGaussNE()
1496 {
1497 }
1498
1499 TypeOfField MEDCouplingFieldDiscretizationGaussNE::getEnum() const
1500 {
1501   return TYPE;
1502 }
1503
1504 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationGaussNE::clone() const
1505 {
1506   return new MEDCouplingFieldDiscretizationGaussNE(*this);
1507 }
1508
1509 std::string MEDCouplingFieldDiscretizationGaussNE::getStringRepr() const
1510 {
1511   return std::string(REPR);
1512 }
1513
1514 const char *MEDCouplingFieldDiscretizationGaussNE::getRepr() const
1515 {
1516   return REPR;
1517 }
1518
1519 bool MEDCouplingFieldDiscretizationGaussNE::isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const
1520 {
1521   const MEDCouplingFieldDiscretizationGaussNE *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationGaussNE *>(other);
1522   bool ret=otherC!=0;
1523   if(!ret)
1524     reason="Spatial discrtization of this is ON_GAUSS_NE, which is not the case of other.";
1525   return ret;
1526 }
1527
1528 int MEDCouplingFieldDiscretizationGaussNE::getNumberOfTuples(const MEDCouplingMesh *mesh) const
1529 {
1530   int ret=0;
1531   int nbOfCells=mesh->getNumberOfCells();
1532   for(int i=0;i<nbOfCells;i++)
1533     {
1534       INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell(i);
1535       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
1536       if(cm.isDynamic())
1537         throw INTERP_KERNEL::Exception("Not implemented yet Gauss node on elements for polygons and polyedrons !");
1538       ret+=cm.getNumberOfNodes();
1539     }
1540   return ret;
1541 }
1542
1543 int MEDCouplingFieldDiscretizationGaussNE::getNumberOfMeshPlaces(const MEDCouplingMesh *mesh) const
1544 {
1545   return mesh->getNumberOfCells();
1546 }
1547
1548 DataArrayInt *MEDCouplingFieldDiscretizationGaussNE::getOffsetArr(const MEDCouplingMesh *mesh) const
1549 {
1550   int nbOfTuples=mesh->getNumberOfCells();
1551   DataArrayInt *ret=DataArrayInt::New();
1552   ret->alloc(nbOfTuples+1,1);
1553   int *retPtr=ret->getPointer();
1554   retPtr[0]=0;
1555   for(int i=0;i<nbOfTuples;i++)
1556     {
1557       INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell(i);
1558       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
1559       if(cm.isDynamic())
1560         throw INTERP_KERNEL::Exception("Not implemented yet Gauss node on elements for polygons and polyedrons !");
1561       retPtr[i+1]=retPtr[i]+cm.getNumberOfNodes();
1562     }
1563   return ret;
1564 }
1565
1566 void MEDCouplingFieldDiscretizationGaussNE::renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector<DataArrayDouble *>& arrays,
1567                                                                   const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
1568 {
1569   const int *array=old2NewBg;
1570   if(check)
1571     array=DataArrayInt::CheckAndPreparePermutation(old2NewBg,old2NewBg+mesh->getNumberOfCells());
1572   int nbOfCells=mesh->getNumberOfCells();
1573   int nbOfTuples=getNumberOfTuples(mesh);
1574   int *array2=new int[nbOfTuples];//stores the final conversion array old2New to give to arrays in renumberInPlace.
1575   int *array3=new int[nbOfCells];//store for each cell in after renumbering the offset needed by each cell in new numbering.
1576   array3[0]=0;
1577   for(int i=1;i<nbOfCells;i++)
1578     {
1579       INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell((int)std::distance(array,std::find(array,array+nbOfCells,i-1)));
1580       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
1581       array3[i]=array3[i-1]+cm.getNumberOfNodes();
1582     }
1583   int j=0;
1584   for(int i=0;i<nbOfCells;i++)
1585     {
1586       INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell(i);
1587       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
1588       for(int k=0;k<(int)cm.getNumberOfNodes();k++,j++)
1589         array2[j]=array3[array[i]]+k;
1590     }
1591   delete [] array3;
1592   for(std::vector<DataArrayDouble *>::const_iterator it=arrays.begin();it!=arrays.end();it++)
1593     if(*it)
1594       (*it)->renumberInPlace(array2);
1595   delete [] array2;
1596   if(check)
1597     delete [] const_cast<int *>(array);
1598 }
1599
1600 DataArrayDouble *MEDCouplingFieldDiscretizationGaussNE::getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const
1601 {
1602   throw INTERP_KERNEL::Exception("Not implemented yet !");
1603 }
1604
1605 void MEDCouplingFieldDiscretizationGaussNE::integral(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, bool isWAbs, double *res) const throw(INTERP_KERNEL::Exception)
1606 {
1607   if(!mesh || !arr)
1608     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::integral : input mesh or array is null !");
1609   int nbOfCompo=arr->getNumberOfComponents();
1610   std::fill(res,res+nbOfCompo,0.);
1611   //
1612   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> vol=getMeasureField(mesh,isWAbs);
1613   std::set<INTERP_KERNEL::NormalizedCellType> types=mesh->getAllGeoTypes();
1614   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nbOfNodesPerCell=mesh->computeNbOfNodesPerCell();
1615   nbOfNodesPerCell->computeOffsets2();
1616   const double *arrPtr=arr->begin(),*volPtr=vol->getArray()->begin();
1617   for(std::set<INTERP_KERNEL::NormalizedCellType>::const_iterator it=types.begin();it!=types.end();it++)
1618     {
1619       std::size_t wArrSz=-1;
1620       const double *wArr=GetWeightArrayFromGeometricType(*it,wArrSz);
1621       INTERP_KERNEL::AutoPtr<double> wArr2=new double[wArrSz];
1622       double sum=std::accumulate(wArr,wArr+wArrSz,0.);
1623       std::transform(wArr,wArr+wArrSz,(double *)wArr2,std::bind2nd(std::multiplies<double>(),1./sum));
1624       MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ids=mesh->giveCellsWithType(*it);
1625       MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ids2=ids->buildExplicitArrByRanges(nbOfNodesPerCell);
1626       const int *ptIds2=ids2->begin(),*ptIds=ids->begin();
1627       int nbOfCellsWithCurGeoType=ids->getNumberOfTuples();
1628       for(int i=0;i<nbOfCellsWithCurGeoType;i++,ptIds++,ptIds2+=wArrSz)
1629         {
1630           for(int k=0;k<nbOfCompo;k++)
1631             {
1632               double tmp=0.;
1633               for(std::size_t j=0;j<wArrSz;j++)
1634                 tmp+=arrPtr[nbOfCompo*ptIds2[j]+k]*wArr2[j];
1635               res[k]+=tmp*volPtr[*ptIds];
1636             }
1637         }
1638     }
1639 }
1640
1641 const double *MEDCouplingFieldDiscretizationGaussNE::GetWeightArrayFromGeometricType(INTERP_KERNEL::NormalizedCellType geoType, std::size_t& lgth) throw(INTERP_KERNEL::Exception)
1642 {
1643   switch(geoType)
1644     {
1645     case INTERP_KERNEL::NORM_SEG2:
1646       lgth=(int)sizeof(FGP_SEG2)/sizeof(double);
1647       return FGP_SEG2;
1648     case INTERP_KERNEL::NORM_SEG3:
1649       lgth=(int)sizeof(FGP_SEG3)/sizeof(double);
1650       return FGP_SEG3;
1651     case INTERP_KERNEL::NORM_SEG4:
1652       lgth=(int)sizeof(FGP_SEG4)/sizeof(double);
1653       return FGP_SEG4;
1654     case INTERP_KERNEL::NORM_TRI3:
1655       lgth=(int)sizeof(FGP_TRI3)/sizeof(double);
1656       return FGP_TRI3;
1657     case INTERP_KERNEL::NORM_TRI6:
1658       lgth=(int)sizeof(FGP_TRI6)/sizeof(double);
1659       return FGP_TRI6;
1660     case INTERP_KERNEL::NORM_TRI7:
1661       lgth=(int)sizeof(FGP_TRI7)/sizeof(double);
1662       return FGP_TRI7;
1663     case INTERP_KERNEL::NORM_QUAD4:
1664       lgth=(int)sizeof(FGP_QUAD4)/sizeof(double);
1665       return FGP_QUAD4;
1666     case INTERP_KERNEL::NORM_QUAD9:
1667       lgth=(int)sizeof(FGP_QUAD9)/sizeof(double);
1668       return FGP_QUAD9;
1669     case INTERP_KERNEL::NORM_TETRA4:
1670       lgth=(int)sizeof(FGP_TETRA4)/sizeof(double);
1671       return FGP_TETRA4;
1672     case INTERP_KERNEL::NORM_PENTA6:
1673       lgth=(int)sizeof(FGP_PENTA6)/sizeof(double);
1674       return FGP_PENTA6;
1675     case INTERP_KERNEL::NORM_HEXA8:
1676       lgth=(int)sizeof(FGP_HEXA8)/sizeof(double);
1677       return FGP_HEXA8;
1678     case INTERP_KERNEL::NORM_HEXA27:
1679       lgth=(int)sizeof(FGP_HEXA27)/sizeof(double);
1680       return FGP_HEXA27;
1681     case INTERP_KERNEL::NORM_PYRA5:
1682       lgth=(int)sizeof(FGP_PYRA5)/sizeof(double);
1683       return FGP_PYRA5;
1684     default:
1685       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 !");
1686     }
1687 }
1688
1689 void MEDCouplingFieldDiscretizationGaussNE::computeMeshRestrictionFromTupleIds(const MEDCouplingMesh *mesh, const int *partBg, const int *partEnd,
1690                                                                                DataArrayInt *&cellRest)
1691 {
1692   throw INTERP_KERNEL::Exception("Not implemented yet !");
1693 }
1694
1695 void MEDCouplingFieldDiscretizationGaussNE::checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception)
1696 {
1697 }
1698
1699 double MEDCouplingFieldDiscretizationGaussNE::getIJK(const MEDCouplingMesh *mesh, const DataArrayDouble *da,
1700                                                      int cellId, int nodeIdInCell, int compoId) const throw(INTERP_KERNEL::Exception)
1701 {
1702   int offset=0;
1703   for(int i=0;i<cellId;i++)
1704     {
1705       INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell(i);
1706       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
1707       offset+=cm.getNumberOfNodes();
1708     }
1709   return da->getIJ(offset+nodeIdInCell,compoId);
1710 }
1711
1712 void MEDCouplingFieldDiscretizationGaussNE::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArrayDouble *da) const throw(INTERP_KERNEL::Exception)
1713 {
1714   int nbOfTuples=getNumberOfTuples(mesh);
1715   if(nbOfTuples!=da->getNumberOfTuples())
1716     {
1717       std::ostringstream oss; oss << "Invalid number of tuples in the array : expecting " << nbOfTuples << " !";
1718       throw INTERP_KERNEL::Exception(oss.str().c_str());
1719     }
1720 }
1721
1722 MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationGaussNE::getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const
1723 {
1724   if(!mesh)
1725     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::getMeasureField : mesh instance specified is NULL !");
1726   throw INTERP_KERNEL::Exception("Not implemented yet !");
1727 }
1728
1729 void MEDCouplingFieldDiscretizationGaussNE::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
1730 {
1731   throw INTERP_KERNEL::Exception("Not implemented yet !");
1732 }
1733
1734 void MEDCouplingFieldDiscretizationGaussNE::getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, int i, int j, int k, double *res) const
1735 {
1736   throw INTERP_KERNEL::Exception("getValueOnPos(i,j,k) : Not applyable for Gauss points !");
1737 }
1738
1739 DataArrayDouble *MEDCouplingFieldDiscretizationGaussNE::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, int nbOfPoints) const
1740 {
1741   throw INTERP_KERNEL::Exception("getValueOnMulti : Not implemented for Gauss NE !");
1742 }
1743
1744 MEDCouplingMesh *MEDCouplingFieldDiscretizationGaussNE::buildSubMeshData(const MEDCouplingMesh *mesh, const int *start, const int *end, DataArrayInt *&di) const
1745 {
1746   di=computeTupleIdsToSelectFromCellIds(mesh,start,end);
1747   return mesh->buildPart(start,end);
1748 }
1749
1750 /*!
1751  * This method returns a tuple ids selection from cell ids selection [start;end).
1752  * This method is called by MEDCouplingFieldDiscretizationGaussNE::buildSubMeshData to return parameter \b di.
1753  *
1754  * \return a newly allocated array containing ids to select into the DataArrayDouble of the field.
1755  *
1756  */
1757 DataArrayInt *MEDCouplingFieldDiscretizationGaussNE::computeTupleIdsToSelectFromCellIds(const MEDCouplingMesh *mesh, const int *startCellIds, const int *endCellIds) const
1758 {
1759   if(!mesh)
1760     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::computeTupleIdsToSelectFromCellIds : null mesh !");
1761   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nbOfNodesPerCell=mesh->computeNbOfNodesPerCell();
1762   nbOfNodesPerCell->computeOffsets2();
1763   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> sel=DataArrayInt::New(); sel->useArray(startCellIds,false,CPP_DEALLOC,(int)std::distance(startCellIds,endCellIds),1);
1764   return sel->buildExplicitArrByRanges(nbOfNodesPerCell);
1765 }
1766
1767 /*!
1768  * No implementation needed !
1769  */
1770 void MEDCouplingFieldDiscretizationGaussNE::renumberValuesOnNodes(double , const int *, int newNbOfNodes, DataArrayDouble *) const
1771 {
1772 }
1773
1774 void MEDCouplingFieldDiscretizationGaussNE::renumberValuesOnCells(double epsOnVals, const MEDCouplingMesh *mesh, const int *old2New, int newSz, DataArrayDouble *arr) const
1775 {
1776   throw INTERP_KERNEL::Exception("Not implemented yet !");
1777 }
1778
1779 void MEDCouplingFieldDiscretizationGaussNE::renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const int *new2old, int newSz, DataArrayDouble *arr) const
1780 {
1781   throw INTERP_KERNEL::Exception("Not implemented yet !");
1782 }
1783
1784 MEDCouplingFieldDiscretizationGaussNE::MEDCouplingFieldDiscretizationGaussNE(const MEDCouplingFieldDiscretizationGaussNE& other):MEDCouplingFieldDiscretization(other)
1785 {
1786 }
1787
1788 TypeOfField MEDCouplingFieldDiscretizationKriging::getEnum() const
1789 {
1790   return TYPE;
1791 }
1792
1793 const char *MEDCouplingFieldDiscretizationKriging::getRepr() const
1794 {
1795   return REPR;
1796 }
1797
1798 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationKriging::clone() const
1799 {
1800   return new MEDCouplingFieldDiscretizationKriging;
1801 }
1802
1803 std::string MEDCouplingFieldDiscretizationKriging::getStringRepr() const
1804 {
1805   return std::string(REPR);
1806 }
1807
1808 void MEDCouplingFieldDiscretizationKriging::checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception)
1809 {
1810   if(nat!=ConservativeVolumic)
1811     throw INTERP_KERNEL::Exception("Invalid nature for Kriging field : expected ConservativeVolumic !");
1812 }
1813
1814 bool MEDCouplingFieldDiscretizationKriging::isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const
1815 {
1816   const MEDCouplingFieldDiscretizationKriging *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationKriging *>(other);
1817   bool ret=otherC!=0;
1818   if(!ret)
1819     reason="Spatial discrtization of this is ON_NODES_KR, which is not the case of other.";
1820   return ret;
1821 }
1822
1823 MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationKriging::getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const
1824 {
1825   if(!mesh)
1826     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::getMeasureField : mesh instance specified is NULL !");
1827   throw INTERP_KERNEL::Exception("getMeasureField on FieldDiscretizationKriging : not implemented yet !");
1828 }
1829
1830 void MEDCouplingFieldDiscretizationKriging::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
1831 {
1832   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> res2=MEDCouplingFieldDiscretizationKriging::getValueOnMulti(arr,mesh,loc,1);
1833   std::copy(res2->begin(),res2->end(),res);
1834 }
1835
1836 DataArrayDouble *MEDCouplingFieldDiscretizationKriging::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, int nbOfTargetPoints) const
1837 {
1838   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coords=getLocalizationOfDiscValues(mesh);
1839   int nbOfPts=coords->getNumberOfTuples();
1840   int dimension=coords->getNumberOfComponents();
1841   //
1842   int delta=0;
1843   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> KnewiK=computeVectorOfCoefficients(mesh,arr,delta);
1844   //
1845   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> locArr=DataArrayDouble::New();
1846   locArr->useArray(loc,false,CPP_DEALLOC,nbOfTargetPoints,dimension);
1847   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrix2=coords->buildEuclidianDistanceDenseMatrixWith(locArr);
1848   operateOnDenseMatrix(mesh->getSpaceDimension(),nbOfPts*nbOfTargetPoints,matrix2->getPointer());
1849   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrix3=DataArrayDouble::New();
1850   matrix3->alloc((nbOfPts+delta)*nbOfTargetPoints,1);
1851   double *work=matrix3->getPointer();
1852   const double *workCst=matrix2->getConstPointer();
1853   const double *workCst2=loc;
1854   for(int i=0;i<nbOfTargetPoints;i++,workCst+=nbOfPts,workCst2+=delta-1)
1855     {
1856       for(int j=0;j<nbOfPts;j++)
1857         work[j*nbOfTargetPoints+i]=workCst[j];
1858       work[nbOfPts*nbOfTargetPoints+i]=1.0;
1859       for(int j=0;j<delta-1;j++)
1860         work[(nbOfPts+1+j)*nbOfTargetPoints+i]=workCst2[j];
1861     }
1862   //
1863   int nbOfCompo=arr->getNumberOfComponents();
1864   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret=DataArrayDouble::New();
1865   ret->alloc(nbOfTargetPoints,nbOfCompo);
1866   INTERP_KERNEL::matrixProduct(KnewiK->getConstPointer(),1,nbOfPts+delta,matrix3->getConstPointer(),nbOfPts+delta,nbOfTargetPoints*nbOfCompo,ret->getPointer());
1867   return ret.retn();
1868 }
1869
1870 /*!
1871  * 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
1872  * number of tuples should be equal to the number of representing points in \a mesh.
1873  *
1874  * \param [in] mesh is the sources of nodes on which kriging will be done regarding the parameters and the value of \c this->getSpaceDimension()
1875  * \param [in] arr input field DataArrayDouble whose number of tuples must be equal to the number of nodes in \a mesh
1876  * \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.
1877  *              Whatever the value of \a isDrift the number of tuples of returned DataArrayDouble  will be equal to \c arr->getNumberOfTuples() + \a isDrift.
1878  * \return a newly allocated array containing coefficients including or not drift coefficient at the end depending the value of \a isDrift parameter.
1879  */
1880 DataArrayDouble *MEDCouplingFieldDiscretizationKriging::computeVectorOfCoefficients(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, int& isDrift) const
1881 {
1882   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coords=getLocalizationOfDiscValues(mesh);
1883   int nbOfPts=coords->getNumberOfTuples();
1884   //int dimension=coords->getNumberOfComponents();
1885   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrix=coords->buildEuclidianDistanceDenseMatrix();
1886   operateOnDenseMatrix(mesh->getSpaceDimension(),nbOfPts*nbOfPts,matrix->getPointer());
1887   // Drift
1888   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixWithDrift=performDrift(matrix,coords,isDrift);
1889   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixInv=DataArrayDouble::New();
1890   matrixInv->alloc((nbOfPts+isDrift)*(nbOfPts+isDrift),1);
1891   INTERP_KERNEL::inverseMatrix(matrixWithDrift->getConstPointer(),nbOfPts+isDrift,matrixInv->getPointer());
1892   //
1893   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> KnewiK=DataArrayDouble::New();
1894   KnewiK->alloc((nbOfPts+isDrift)*1,1);
1895   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr2=DataArrayDouble::New();
1896   arr2->alloc((nbOfPts+isDrift)*1,1);
1897   double *work=std::copy(arr->begin(),arr->end(),arr2->getPointer());
1898   std::fill(work,work+isDrift,0.);
1899   INTERP_KERNEL::matrixProduct(matrixInv->getConstPointer(),nbOfPts+isDrift,nbOfPts+isDrift,arr2->getConstPointer(),nbOfPts+isDrift,1,KnewiK->getPointer());
1900   return KnewiK.retn();
1901 }
1902
1903 /*!
1904  * 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.
1905  *
1906  * \param [in] spaceDimension space dimension of the input mesh on which the Kriging has to be performed
1907  * \param [in] nbOfElems is the result of the product of nb of rows and the nb of columns of matrix \a matrixPtr
1908  * \param [in,out] matrixPtr is the dense matrix whose on each values the operation will be applied
1909  */
1910 void MEDCouplingFieldDiscretizationKriging::operateOnDenseMatrix(int spaceDimension, int nbOfElems, double *matrixPtr) const
1911 {
1912   switch(spaceDimension)
1913     {
1914     case 1:
1915       {
1916         for(int i=0;i<nbOfElems;i++)
1917           {
1918             double val=matrixPtr[i];
1919             matrixPtr[i]=val*val*val;
1920           }
1921         break;
1922       }
1923     default:
1924       throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::operateOnDenseMatrix : only dimension 1 implemented !");
1925     }
1926 }
1927
1928 /*!
1929  * Starting from a square matrix \a matr, this method returns a newly allocated dense square matrix whose \a matr is included in returned matrix
1930  * in the top left corner, and in the remaining returned matrix the parameters to take into account about the kriging drift.
1931  * For the moment only linear srift is implemented.
1932  *
1933  * \param [in] arr the position of points were input mesh geometry is considered for Kriging
1934  * \param [in] matr input matrix whose drift part will be added
1935  * \param [out] delta the difference between the size of the output matrix and the input matrix \a matr.
1936  * \return a newly allocated matrix bigger than input matrix \a matr.
1937  */
1938 DataArrayDouble *MEDCouplingFieldDiscretizationKriging::performDrift(const DataArrayDouble *matr, const DataArrayDouble *arr, int& delta) const
1939 {
1940   int spaceDimension=arr->getNumberOfComponents();
1941   delta=spaceDimension+1;
1942   int szOfMatrix=arr->getNumberOfTuples();
1943   if(szOfMatrix*szOfMatrix!=matr->getNumberOfTuples())
1944     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::performDrift : invalid size");
1945   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret=DataArrayDouble::New();
1946   ret->alloc((szOfMatrix+delta)*(szOfMatrix+delta),1);
1947   const double *srcWork=matr->getConstPointer();
1948   const double *srcWork2=arr->getConstPointer();
1949   double *destWork=ret->getPointer();
1950   for(int i=0;i<szOfMatrix;i++)
1951     {
1952       destWork=std::copy(srcWork,srcWork+szOfMatrix,destWork);
1953       srcWork+=szOfMatrix;
1954       *destWork++=1.;
1955       destWork=std::copy(srcWork2,srcWork2+spaceDimension,destWork);
1956       srcWork2+=spaceDimension;
1957     }
1958   std::fill(destWork,destWork+szOfMatrix,1.); destWork+=szOfMatrix;
1959   std::fill(destWork,destWork+spaceDimension+1,0.); destWork+=spaceDimension+1;
1960   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arrNoI=arr->toNoInterlace();
1961   srcWork2=arrNoI->getConstPointer();
1962   for(int i=0;i<spaceDimension;i++)
1963     {
1964       destWork=std::copy(srcWork2,srcWork2+szOfMatrix,destWork);
1965       srcWork2+=szOfMatrix;
1966       std::fill(destWork,destWork+spaceDimension+1,0.);
1967       destWork+=spaceDimension;
1968     }
1969   //
1970   return ret.retn();
1971 }
1972