src/MEDCoupling/MEDCouplingFieldDiscretization.cxx

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