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