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