src/MEDCoupling/MEDCouplingFieldDiscretization.cxx

   1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (CEA/DEN)
  20
  21 #include "MEDCouplingFieldDiscretization.hxx"
  22 #include "MEDCouplingCMesh.hxx"
  23 #include "MEDCouplingUMesh.hxx"
  24 #include "MEDCouplingFieldDouble.hxx"
  25 #include "MEDCouplingAutoRefCountObjectPtr.hxx"
  26
  27 #include "CellModel.hxx"
  28 #include "InterpolationUtils.hxx"
  29 #include "InterpKernelAutoPtr.hxx"
  30 #include "InterpKernelGaussCoords.hxx"
  31 #include "InterpKernelMatrixTools.hxx"
  32
  33 #include <set>
  34 #include <list>
  35 #include <limits>
  36 #include <sstream>
  37 #include <numeric>
  38 #include <algorithm>
  39 #include <functional>
  40
  41 using namespace ParaMEDMEM;
  42
  43 const double MEDCouplingFieldDiscretization::DFLT_PRECISION=1.e-12;
  44
  45 const char MEDCouplingFieldDiscretizationP0::REPR[]="P0";
  46
  47 const TypeOfField MEDCouplingFieldDiscretizationP0::TYPE=ON_CELLS;
  48
  49 const char MEDCouplingFieldDiscretizationP1::REPR[]="P1";
  50
  51 const TypeOfField MEDCouplingFieldDiscretizationP1::TYPE=ON_NODES;
  52
  53 const int MEDCouplingFieldDiscretizationPerCell::DFT_INVALID_LOCID_VALUE=-1;
  54
  55 const char MEDCouplingFieldDiscretizationGauss::REPR[]="GAUSS";
  56
  57 const TypeOfField MEDCouplingFieldDiscretizationGauss::TYPE=ON_GAUSS_PT;
  58
  59 const char MEDCouplingFieldDiscretizationGaussNE::REPR[]="GSSNE";
  60
  61 const TypeOfField MEDCouplingFieldDiscretizationGaussNE::TYPE=ON_GAUSS_NE;
  62
  63 const char MEDCouplingFieldDiscretizationKriging::REPR[]="KRIGING";
  64
  65 const TypeOfField MEDCouplingFieldDiscretizationKriging::TYPE=ON_NODES_KR;
  66
  67 // doc is here http://www.code-aster.org/V2/doc/default/fr/man_r/r3/r3.01.01.pdf
  68 const double MEDCouplingFieldDiscretizationGaussNE::FGP_SEG2[2]={1.,1.};
  69 const double MEDCouplingFieldDiscretizationGaussNE::FGP_SEG3[3]={0.5555555555555556,0.5555555555555556,0.8888888888888888};
  70 const double MEDCouplingFieldDiscretizationGaussNE::FGP_SEG4[4]={0.347854845137454,0.347854845137454,0.652145154862546,0.652145154862546};
  71 const double MEDCouplingFieldDiscretizationGaussNE::FGP_TRI3[3]={0.16666666666666666,0.16666666666666666,0.16666666666666666};
  72 const double MEDCouplingFieldDiscretizationGaussNE::FGP_TRI6[6]={0.0549758718227661,0.0549758718227661,0.0549758718227661,0.11169079483905,0.11169079483905,0.11169079483905};
  73 const double MEDCouplingFieldDiscretizationGaussNE::FGP_TRI7[7]={0.062969590272413,0.062969590272413,0.062969590272413,0.066197076394253,0.066197076394253,0.066197076394253,0.1125};
  74 const double MEDCouplingFieldDiscretizationGaussNE::FGP_QUAD4[4]={1.,1.,1.,1.};
  75 const double MEDCouplingFieldDiscretizationGaussNE::FGP_QUAD9[9]={0.30864197530864196,0.30864197530864196,0.30864197530864196,0.30864197530864196,0.49382716049382713,0.49382716049382713,0.49382716049382713,0.49382716049382713,0.7901234567901234};
  76 const double MEDCouplingFieldDiscretizationGaussNE::FGP_TETRA4[4]={0.041666666666666664,0.041666666666666664,0.041666666666666664,0.041666666666666664};
  77 const double MEDCouplingFieldDiscretizationGaussNE::FGP_PENTA6[6]={0.16666666666666666,0.16666666666666666,0.16666666666666666,0.16666666666666666,0.16666666666666666,0.16666666666666666};
  78 const double MEDCouplingFieldDiscretizationGaussNE::FGP_HEXA8[8]={1.,1.,1.,1.,1.,1.,1.,1.};
  79 const double MEDCouplingFieldDiscretizationGaussNE::FGP_HEXA27[27]={0.1714677640603567,0.1714677640603567,0.1714677640603567,0.1714677640603567,0.1714677640603567,0.1714677640603567,0.1714677640603567,0.1714677640603567,0.27434842249657065,0.27434842249657065,0.27434842249657065,0.27434842249657065,0.27434842249657065,0.27434842249657065,0.27434842249657065,0.27434842249657065,0.27434842249657065,0.27434842249657065,0.27434842249657065,0.27434842249657065,0.43895747599451296,0.43895747599451296,0.43895747599451296,0.43895747599451296,0.43895747599451296,0.43895747599451296,0.7023319615912208};
  80 const double MEDCouplingFieldDiscretizationGaussNE::FGP_PYRA5[5]={0.13333333333333333,0.13333333333333333,0.13333333333333333,0.13333333333333333,0.13333333333333333};
  81 const double MEDCouplingFieldDiscretizationGaussNE::REF_SEG2[2]={-1.,1.};
  82 const double MEDCouplingFieldDiscretizationGaussNE::REF_SEG3[3]={-1.,0.,1.};
  83 const double MEDCouplingFieldDiscretizationGaussNE::REF_SEG4[4]={-1.,1.,-0.3333333333333333,0.3333333333333333};
  84 const double MEDCouplingFieldDiscretizationGaussNE::REF_TRI3[6]={0.,0.,1.,0.,0.,1.};
  85 const double MEDCouplingFieldDiscretizationGaussNE::REF_TRI6[12]={0.,0.,1.,0.,0.,1.,0.5,0.,0.5,0.5,0.,0.5};
  86 const double MEDCouplingFieldDiscretizationGaussNE::REF_TRI7[14]={0.,0.,1.,0.,0.,1.,0.5,0.,0.5,0.5,0.,0.5,0.3333333333333333,0.3333333333333333};
  87 const double MEDCouplingFieldDiscretizationGaussNE::REF_QUAD4[8]={-1.,-1.,1.,-1.,1.,1.,-1.,1.};
  88 const double MEDCouplingFieldDiscretizationGaussNE::REF_QUAD8[16]={-1.,-1.,1.,-1.,1.,1.,-1.,1.,0.,-1.,1.,0.,0.,1.,-1.,0.};
  89 const double MEDCouplingFieldDiscretizationGaussNE::REF_QUAD9[18]={-1.,-1.,1.,-1.,1.,1.,-1.,1.,0.,-1.,1.,0.,0.,1.,-1.,0.,0.,0.};
  90 const double MEDCouplingFieldDiscretizationGaussNE::REF_TETRA4[12]={0.,1.,0.,0.,0.,1.,0.,0.,0.,1.,0.,0.};
  91 const double MEDCouplingFieldDiscretizationGaussNE::REF_TETRA10[30]={0.,1.,0.,0.,0.,1.,0.,0.,0.,1.,0.,0.,0.,0.5,0.5,0.,0.,0.5,0.,0.5,0.,0.5,0.5,0.,0.5,0.,0.5,0.5,0.,0.};
  92 const double MEDCouplingFieldDiscretizationGaussNE::REF_PENTA6[18]={-1.,1.,0.,-1.,0.,1.,-1.,0.,0.,1.,1.,0.,1.,0.,1.,1.,0.,0.};
  93 const double MEDCouplingFieldDiscretizationGaussNE::REF_PENTA15[45]={-1.,1.,0.,-1.,0.,1.,-1.,0.,0.,1.,1.,0.,1.,0.,1.,1.,0.,0.,-1.,0.5,0.5,-1.,0.,0.5,-1.,0.5,0.,0.,1.,0.,0.,0.,1.,0.,0.,0.,1.,0.5,0.5,1.,0.,0.5,1.,0.5,0.};
  94 const double MEDCouplingFieldDiscretizationGaussNE::REF_HEXA8[24]={-1.,-1.,-1.,1.,-1.,-1.,1.,1.,-1.,-1.,1.,-1.,-1.,-1.,1.,1.,-1.,1.,1.,1.,1.,-1.,1.,1.};
  95 const double MEDCouplingFieldDiscretizationGaussNE::REF_HEXA20[60]={-1.,-1.,-1.,1.,-1.,-1.,1.,1.,-1.,-1.,1.,-1.,-1.,-1.,1.,1.,-1.,1.,1.,1.,1.,-1.,1.,1.,0.,-1.,-1.,1.,0.,-1.,0.,1.,-1.,-1.,0.,-1.,-1.,-1.,0.,1.,-1.,0.,1.,1.,0.,-1.,1.,0.,0.,-1.,1.,1.,0.,1.,0.,1.,1.,-1.,0.,1.};
  96 const double MEDCouplingFieldDiscretizationGaussNE::REF_HEXA27[81]={-1.,-1.,-1.,1.,-1.,-1.,1.,1.,-1.,-1.,1.,-1.,-1.,-1.,1.,1.,-1.,1.,1.,1.,1.,-1.,1.,1.,0.,-1.,-1.,1.,0.,-1.,0.,1.,-1.,-1.,0.,-1.,-1.,-1.,0.,1.,-1.,0.,1.,1.,0.,-1.,1.,0.,0.,-1.,1.,1.,0.,1.,0.,1.,1.,-1.,0.,1.,0.,0.,-1.,0.,-1.,0.,1.,0.,0.,0.,1.,0.,-1.,0.,0.,0.,0.,1.,0.,0.,0.};
  97 const double MEDCouplingFieldDiscretizationGaussNE::REF_PYRA5[15]={1.,0.,0.,0.,1.,0.,-1.,0.,0.,0.,-1.,0.,0.,0.,1.};
  98 const double MEDCouplingFieldDiscretizationGaussNE::REF_PYRA13[39]={1.,0.,0.,0.,1.,0.,-1.,0.,0.,0.,-1.,0.,0.,0.,1.,0.5,0.5,0.,-0.5,0.5,0.,-0.5,-0.5,0.,0.5,-0.5,0.,0.5,0.,0.5,0.,0.5,0.5,-0.5,0.,0.5,0.,-0.5,0.5};
  99
 100 MEDCouplingFieldDiscretization::MEDCouplingFieldDiscretization():_precision(DFLT_PRECISION)
 101 {
 102 }
 103
 104 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretization::New(TypeOfField type)
 105 {
 106   switch(type)
 107     {
 108     case MEDCouplingFieldDiscretizationP0::TYPE:
 109       return new MEDCouplingFieldDiscretizationP0;
 110     case MEDCouplingFieldDiscretizationP1::TYPE:
 111       return new MEDCouplingFieldDiscretizationP1;
 112     case MEDCouplingFieldDiscretizationGauss::TYPE:
 113       return new MEDCouplingFieldDiscretizationGauss;
 114     case MEDCouplingFieldDiscretizationGaussNE::TYPE:
 115       return new MEDCouplingFieldDiscretizationGaussNE;
 116     case MEDCouplingFieldDiscretizationKriging::TYPE:
 117       return new MEDCouplingFieldDiscretizationKriging;
 118     default:
 119       throw INTERP_KERNEL::Exception("Choosen discretization is not implemented yet.");
 120     }
 121 }
 122
 123 TypeOfField MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(const char *repr) throw(INTERP_KERNEL::Exception)
 124 {
 125   std::string reprCpp(repr);
 126   if(reprCpp==MEDCouplingFieldDiscretizationP0::REPR)
 127     return MEDCouplingFieldDiscretizationP0::TYPE;
 128   if(reprCpp==MEDCouplingFieldDiscretizationP1::REPR)
 129     return MEDCouplingFieldDiscretizationP1::TYPE;
 130   if(reprCpp==MEDCouplingFieldDiscretizationGauss::REPR)
 131     return MEDCouplingFieldDiscretizationGauss::TYPE;
 132   if(reprCpp==MEDCouplingFieldDiscretizationGaussNE::REPR)
 133     return MEDCouplingFieldDiscretizationGaussNE::TYPE;
 134   if(reprCpp==MEDCouplingFieldDiscretizationKriging::REPR)
 135     return MEDCouplingFieldDiscretizationKriging::TYPE;
 136   throw INTERP_KERNEL::Exception("Representation does not match with any field discretization !");
 137 }
 138
 139 bool MEDCouplingFieldDiscretization::isEqual(const MEDCouplingFieldDiscretization *other, double eps) const
 140 {
 141   std::string reason;
 142   return isEqualIfNotWhy(other,eps,reason);
 143 }
 144
 145 bool MEDCouplingFieldDiscretization::isEqualWithoutConsideringStr(const MEDCouplingFieldDiscretization *other, double eps) const
 146 {
 147   return isEqual(other,eps);
 148 }
 149
 150 /*!
 151  * This method is an alias of MEDCouplingFieldDiscretization::clone. It is only here for coherency with all the remaining of MEDCoupling.
 152  * \sa MEDCouplingFieldDiscretization::clone.
 153  */
 154 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretization::deepCpy() const
 155 {
 156   return clone();
 157 }
 158
 159 /*!
 160  * For all field discretization excepted GaussPts the [ \a startCellIds, \a endCellIds ) has no impact on the cloned instance.
 161  */
 162 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretization::clonePart(const int *startCellIds, const int *endCellIds) const
 163 {
 164   return clone();
 165 }
 166
 167 /*!
 168  * For all field discretization excepted GaussPts the slice( \a beginCellId, \a endCellIds, \a stepCellId ) has no impact on the cloned instance.
 169  */
 170 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretization::clonePartRange(int beginCellIds, int endCellIds, int stepCellIds) const
 171 {
 172   return clone();
 173 }
 174
 175 /*!
 176  * Excepted for MEDCouplingFieldDiscretizationPerCell no underlying TimeLabel object : nothing to do in generally.
 177  */
 178 void MEDCouplingFieldDiscretization::updateTime() const
 179 {
 180 }
 181
 182 std::size_t MEDCouplingFieldDiscretization::getHeapMemorySizeWithoutChildren() const
 183 {
 184   return 0;
 185 }
 186
 187 std::vector<RefCountObject *> MEDCouplingFieldDiscretization::getDirectChildren() const
 188 {
 189   return std::vector<RefCountObject *>();
 190 }
 191
 192 /*!
 193  * Computes normL1 of DataArrayDouble instance arr.
 194  * @param res output parameter expected to be of size arr->getNumberOfComponents();
 195  * @throw when the field discretization fails on getMeasure fields (gauss points for example)
 196  */
 197 void MEDCouplingFieldDiscretization::normL1(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, double *res) const throw(INTERP_KERNEL::Exception)
 198 {
 199   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> vol=getMeasureField(mesh,true);
 200   int nbOfCompo=arr->getNumberOfComponents();
 201   int nbOfElems=getNumberOfTuples(mesh);
 202   std::fill(res,res+nbOfCompo,0.);
 203   const double *arrPtr=arr->getConstPointer();
 204   const double *volPtr=vol->getArray()->getConstPointer();
 205   double deno=0.;
 206   for(int i=0;i<nbOfElems;i++)
 207     {
 208       double v=fabs(volPtr[i]);
 209       for(int j=0;j<nbOfCompo;j++)
 210         res[j]+=fabs(arrPtr[i*nbOfCompo+j])*v;
 211       deno+=v;
 212     }
 213   std::transform(res,res+nbOfCompo,res,std::bind2nd(std::multiplies<double>(),1./deno));
 214 }
 215
 216 /*!
 217  * Computes normL2 of DataArrayDouble instance arr.
 218  * @param res output parameter expected to be of size arr->getNumberOfComponents();
 219  * @throw when the field discretization fails on getMeasure fields (gauss points for example)
 220  */
 221 void MEDCouplingFieldDiscretization::normL2(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, double *res) const throw(INTERP_KERNEL::Exception)
 222 {
 223   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> vol=getMeasureField(mesh,true);
 224   int nbOfCompo=arr->getNumberOfComponents();
 225   int nbOfElems=getNumberOfTuples(mesh);
 226   std::fill(res,res+nbOfCompo,0.);
 227   const double *arrPtr=arr->getConstPointer();
 228   const double *volPtr=vol->getArray()->getConstPointer();
 229   double deno=0.;
 230   for(int i=0;i<nbOfElems;i++)
 231     {
 232       double v=fabs(volPtr[i]);
 233       for(int j=0;j<nbOfCompo;j++)
 234         res[j]+=arrPtr[i*nbOfCompo+j]*arrPtr[i*nbOfCompo+j]*v;
 235       deno+=v;
 236     }
 237   std::transform(res,res+nbOfCompo,res,std::bind2nd(std::multiplies<double>(),1./deno));
 238   std::transform(res,res+nbOfCompo,res,std::ptr_fun<double,double>(std::sqrt));
 239 }
 240
 241 /*!
 242  * Computes integral of DataArrayDouble instance arr.
 243  * @param res output parameter expected to be of size arr->getNumberOfComponents();
 244  * @throw when the field discretization fails on getMeasure fields (gauss points for example)
 245  */
 246 void MEDCouplingFieldDiscretization::integral(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, bool isWAbs, double *res) const throw(INTERP_KERNEL::Exception)
 247 {
 248   if(!mesh)
 249     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretization::integral : mesh is NULL !");
 250   if(!arr)
 251     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretization::integral : input array is NULL !");
 252   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> vol=getMeasureField(mesh,isWAbs);
 253   int nbOfCompo=arr->getNumberOfComponents();
 254   int nbOfElems=getNumberOfTuples(mesh);
 255   if(nbOfElems!=arr->getNumberOfTuples())
 256     {
 257       std::ostringstream oss; oss << "MEDCouplingFieldDiscretization::integral : field is not correct ! number of tuples in array is " << arr->getNumberOfTuples();
 258       oss << " whereas number of tuples expected is " << nbOfElems << " !";
 259       throw INTERP_KERNEL::Exception(oss.str().c_str());
 260     }
 261   std::fill(res,res+nbOfCompo,0.);
 262   const double *arrPtr=arr->getConstPointer();
 263   const double *volPtr=vol->getArray()->getConstPointer();
 264   INTERP_KERNEL::AutoPtr<double> tmp=new double[nbOfCompo];
 265   for (int i=0;i<nbOfElems;i++)
 266     {
 267       std::transform(arrPtr+i*nbOfCompo,arrPtr+(i+1)*nbOfCompo,(double *)tmp,std::bind2nd(std::multiplies<double>(),volPtr[i]));
 268       std::transform((double *)tmp,(double *)tmp+nbOfCompo,res,res,std::plus<double>());
 269     }
 270 }
 271
 272 /*!
 273  * This method is strictly equivalent to MEDCouplingFieldDiscretization::buildSubMeshData except that it is optimized for input defined as a range of cell ids.
 274  *
 275  * \param [out] beginOut Valid only if \a di is NULL
 276  * \param [out] endOut Valid only if \a di is NULL
 277  * \param [out] stepOut Valid only if \a di is NULL
 278  * \param [out] di is an array returned that specifies entity ids (nodes, cells, Gauss points... ) in array if no output range is foundable.
 279  *
 280  * \sa MEDCouplingFieldDiscretization::buildSubMeshData
 281  */
 282 MEDCouplingMesh *MEDCouplingFieldDiscretization::buildSubMeshDataRange(const MEDCouplingMesh *mesh, int beginCellIds, int endCellIds, int stepCellIds, int& beginOut, int& endOut, int& stepOut, DataArrayInt *&di) const
 283 {
 284   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> da=DataArrayInt::Range(beginCellIds,endCellIds,stepCellIds);
 285   return buildSubMeshData(mesh,da->begin(),da->end(),di);
 286 }
 287
 288 void MEDCouplingFieldDiscretization::getSerializationIntArray(DataArrayInt *& arr) const
 289 {
 290   arr=0;
 291 }
 292
 293 /*!
 294  * Empty : Not a bug
 295  */
 296 void MEDCouplingFieldDiscretization::getTinySerializationIntInformation(std::vector<int>& tinyInfo) const
 297 {
 298 }
 299
 300 /*!
 301  * Empty : Not a bug
 302  */
 303 void MEDCouplingFieldDiscretization::getTinySerializationDbleInformation(std::vector<double>& tinyInfo) const
 304 {
 305 }
 306
 307 void MEDCouplingFieldDiscretization::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *& arr)
 308 {
 309   arr=0;
 310 }
 311
 312 /*!
 313  * Empty : Not a bug
 314  */
 315 void MEDCouplingFieldDiscretization::finishUnserialization(const std::vector<double>& tinyInfo)
 316 {
 317 }
 318
 319 /*!
 320  * This method is typically the first step of renumbering. The implementation is empty it is not a bug only gauss is impacted
 321  * virtualy by this method.
 322  */
 323 void MEDCouplingFieldDiscretization::renumberCells(const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
 324 {
 325 }
 326
 327 double MEDCouplingFieldDiscretization::getIJK(const MEDCouplingMesh *mesh, const DataArrayDouble *da,
 328                                               int cellId, int nodeIdInCell, int compoId) const throw(INTERP_KERNEL::Exception)
 329 {
 330   throw INTERP_KERNEL::Exception("getIJK Invalid ! only for GaussPoint and GaussNE discretizations !");
 331 }
 332
 333 void MEDCouplingFieldDiscretization::setGaussLocalizationOnType(const MEDCouplingMesh *m, INTERP_KERNEL::NormalizedCellType type, const std::vector<double>& refCoo,
 334                                                                 const std::vector<double>& gsCoo, const std::vector<double>& wg) throw(INTERP_KERNEL::Exception)
 335 {
 336   throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
 337 }
 338
 339 void MEDCouplingFieldDiscretization::setGaussLocalizationOnCells(const MEDCouplingMesh *m, const int *begin, const int *end, const std::vector<double>& refCoo,
 340                                                                  const std::vector<double>& gsCoo, const std::vector<double>& wg) throw(INTERP_KERNEL::Exception)
 341 {
 342   throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
 343 }
 344
 345 void MEDCouplingFieldDiscretization::clearGaussLocalizations() throw(INTERP_KERNEL::Exception)
 346 {
 347   throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
 348 }
 349
 350 MEDCouplingGaussLocalization& MEDCouplingFieldDiscretization::getGaussLocalization(int locId) throw(INTERP_KERNEL::Exception)
 351 {
 352   throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
 353 }
 354
 355 const MEDCouplingGaussLocalization& MEDCouplingFieldDiscretization::getGaussLocalization(int locId) const throw(INTERP_KERNEL::Exception)
 356 {
 357   throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
 358 }
 359
 360 int MEDCouplingFieldDiscretization::getNbOfGaussLocalization() const throw(INTERP_KERNEL::Exception)
 361 {
 362   throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
 363 }
 364
 365 int MEDCouplingFieldDiscretization::getGaussLocalizationIdOfOneCell(int cellId) const throw(INTERP_KERNEL::Exception)
 366 {
 367   throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
 368 }
 369
 370 int MEDCouplingFieldDiscretization::getGaussLocalizationIdOfOneType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception)
 371 {
 372   throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
 373 }
 374
 375 std::set<int> MEDCouplingFieldDiscretization::getGaussLocalizationIdsOfOneType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception)
 376 {
 377   throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
 378 }
 379
 380 void MEDCouplingFieldDiscretization::getCellIdsHavingGaussLocalization(int locId, std::vector<int>& cellIds) const throw(INTERP_KERNEL::Exception)
 381 {
 382   throw INTERP_KERNEL::Exception("Invalid method for the corresponding field discretization : available only for GaussPoint discretization !");
 383 }
 384
 385 void MEDCouplingFieldDiscretization::RenumberEntitiesFromO2NArr(double eps, const int *old2NewPtr, int newNbOfEntity, DataArrayDouble *arr, const char *msg)
 386 {
 387   if(!arr)
 388     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretization::RenumberEntitiesFromO2NArr : input array is NULL !");
 389   int oldNbOfElems=arr->getNumberOfTuples();
 390   int nbOfComp=arr->getNumberOfComponents();
 391   int newNbOfTuples=newNbOfEntity;
 392   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arrCpy=arr->deepCpy();
 393   const double *ptSrc=arrCpy->getConstPointer();
 394   arr->reAlloc(newNbOfTuples);
 395   double *ptToFill=arr->getPointer();
 396   std::fill(ptToFill,ptToFill+nbOfComp*newNbOfTuples,std::numeric_limits<double>::max());
 397   INTERP_KERNEL::AutoPtr<double> tmp=new double[nbOfComp];
 398   for(int i=0;i<oldNbOfElems;i++)
 399     {
 400       int newNb=old2NewPtr[i];
 401       if(newNb>=0)//if newNb<0 the node is considered as out.
 402         {
 403           if(std::find_if(ptToFill+newNb*nbOfComp,ptToFill+(newNb+1)*nbOfComp,std::bind2nd(std::not_equal_to<double>(),std::numeric_limits<double>::max()))
 404              ==ptToFill+(newNb+1)*nbOfComp)
 405             std::copy(ptSrc+i*nbOfComp,ptSrc+(i+1)*nbOfComp,ptToFill+newNb*nbOfComp);
 406           else
 407             {
 408               std::transform(ptSrc+i*nbOfComp,ptSrc+(i+1)*nbOfComp,ptToFill+newNb*nbOfComp,(double *)tmp,std::minus<double>());
 409               std::transform((double *)tmp,((double *)tmp)+nbOfComp,(double *)tmp,std::ptr_fun<double,double>(fabs));
 410               //if(!std::equal(ptSrc+i*nbOfComp,ptSrc+(i+1)*nbOfComp,ptToFill+newNb*nbOfComp))
 411               if(*std::max_element((double *)tmp,((double *)tmp)+nbOfComp)>eps)
 412                 {
 413                   std::ostringstream oss;
 414                   oss << msg << " " << i << " and " << std::find(old2NewPtr,old2NewPtr+i,newNb)-old2NewPtr
 415                       << " have been merged and " << msg << " field on them are different !";
 416                   throw INTERP_KERNEL::Exception(oss.str().c_str());
 417                 }
 418             }
 419         }
 420     }
 421 }
 422
 423 void MEDCouplingFieldDiscretization::RenumberEntitiesFromN2OArr(const int *new2OldPtr, int new2OldSz, DataArrayDouble *arr, const char *msg)
 424 {
 425   int nbOfComp=arr->getNumberOfComponents();
 426   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arrCpy=arr->deepCpy();
 427   const double *ptSrc=arrCpy->getConstPointer();
 428   arr->reAlloc(new2OldSz);
 429   double *ptToFill=arr->getPointer();
 430   for(int i=0;i<new2OldSz;i++)
 431     {
 432       int oldNb=new2OldPtr[i];
 433       std::copy(ptSrc+oldNb*nbOfComp,ptSrc+(oldNb+1)*nbOfComp,ptToFill+i*nbOfComp);
 434     }
 435 }
 436
 437 MEDCouplingFieldDiscretization::~MEDCouplingFieldDiscretization()
 438 {
 439 }
 440
 441 TypeOfField MEDCouplingFieldDiscretizationP0::getEnum() const
 442 {
 443   return TYPE;
 444 }
 445
 446 /*!
 447  * This method is simply called by MEDCouplingFieldDiscretization::deepCpy. It performs the deep copy of \a this.
 448  *
 449  * \sa MEDCouplingFieldDiscretization::deepCpy.
 450  */
 451 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationP0::clone() const
 452 {
 453   return new MEDCouplingFieldDiscretizationP0;
 454 }
 455
 456 std::string MEDCouplingFieldDiscretizationP0::getStringRepr() const
 457 {
 458   return std::string(REPR);
 459 }
 460
 461 const char *MEDCouplingFieldDiscretizationP0::getRepr() const
 462 {
 463   return REPR;
 464 }
 465
 466 bool MEDCouplingFieldDiscretizationP0::isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const
 467 {
 468   if(!other)
 469     {
 470       reason="other spatial discretization is NULL, and this spatial discretization (P0) is defined.";
 471       return false;
 472     }
 473   const MEDCouplingFieldDiscretizationP0 *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationP0 *>(other);
 474   bool ret=otherC!=0;
 475   if(!ret)
 476     reason="Spatial discrtization of this is ON_CELLS, which is not the case of other.";
 477   return ret;
 478 }
 479
 480 int MEDCouplingFieldDiscretizationP0::getNumberOfTuples(const MEDCouplingMesh *mesh) const throw(INTERP_KERNEL::Exception)
 481 {
 482   if(!mesh)
 483     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP0::getNumberOfTuples : NULL input mesh !");
 484   return mesh->getNumberOfCells();
 485 }
 486
 487 /*!
 488  * mesh is not used here. It is not a bug !
 489  */
 490 int MEDCouplingFieldDiscretizationP0::getNumberOfTuplesExpectedRegardingCode(const MEDCouplingMesh *mesh, const std::vector<int>& code, const std::vector<const DataArrayInt *>& idsPerType) const throw(INTERP_KERNEL::Exception)
 491 {
 492   if(code.size()%3!=0)
 493     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP0::getNumberOfTuplesExpectedRegardingCode : invalid input code !");
 494   int nbOfSplit=(int)idsPerType.size();
 495   int nbOfTypes=(int)code.size()/3;
 496   int ret=0;
 497   for(int i=0;i<nbOfTypes;i++)
 498     {
 499       int nbOfEltInChunk=code[3*i+1];
 500       if(nbOfEltInChunk<0)
 501         throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP0::getNumberOfTuplesExpectedRegardingCode : invalid input code ! presence of negative value in a type !");
 502       int pos=code[3*i+2];
 503       if(pos!=-1)
 504         {
 505           if(pos<0 || pos>=nbOfSplit)
 506             {
 507               std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationP0::getNumberOfTuplesExpectedRegardingCode : input code points to pos " << pos << " in typeid " << i << " ! Should be in [0," << nbOfSplit << ") !";
 508               throw INTERP_KERNEL::Exception(oss.str().c_str());
 509             }
 510           const DataArrayInt *ids(idsPerType[pos]);
 511           if(!ids || !ids->isAllocated() || ids->getNumberOfComponents()!=1 || ids->getNumberOfTuples()!=nbOfEltInChunk || ids->getMinValueInArray()<0)
 512             {
 513               std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationP0::getNumberOfTuplesExpectedRegardingCode : input pfl chunck at pos " << pos << " should have " << i << " tuples and one component and with ids all >=0 !";
 514               throw INTERP_KERNEL::Exception(oss.str().c_str());
 515             }
 516         }
 517       ret+=nbOfEltInChunk;
 518     }
 519   return ret;
 520 }
 521
 522 int MEDCouplingFieldDiscretizationP0::getNumberOfMeshPlaces(const MEDCouplingMesh *mesh) const
 523 {
 524   if(!mesh)
 525     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP0::getNumberOfMeshPlaces : NULL input mesh !");
 526   return mesh->getNumberOfCells();
 527 }
 528
 529 DataArrayInt *MEDCouplingFieldDiscretizationP0::getOffsetArr(const MEDCouplingMesh *mesh) const
 530 {
 531   if(!mesh)
 532     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP0::getOffsetArr : NULL input mesh !");
 533   int nbOfTuples=mesh->getNumberOfCells();
 534   DataArrayInt *ret=DataArrayInt::New();
 535   ret->alloc(nbOfTuples+1,1);
 536   ret->iota(0);
 537   return ret;
 538 }
 539
 540 void MEDCouplingFieldDiscretizationP0::renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector<DataArray *>& arrays,
 541                                                              const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
 542 {
 543   if(!mesh)
 544     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP0::renumberArraysForCell : NULL input mesh !");
 545   const int *array=old2NewBg;
 546   if(check)
 547     array=DataArrayInt::CheckAndPreparePermutation(old2NewBg,old2NewBg+mesh->getNumberOfCells());
 548   for(std::vector<DataArray *>::const_iterator it=arrays.begin();it!=arrays.end();it++)
 549     {
 550       if(*it)
 551         (*it)->renumberInPlace(array);
 552     }
 553   if(check)
 554     free(const_cast<int *>(array));
 555 }
 556
 557 DataArrayDouble *MEDCouplingFieldDiscretizationP0::getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const
 558 {
 559   if(!mesh)
 560     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP0::getLocalizationOfDiscValues : NULL input mesh !");
 561   return mesh->getBarycenterAndOwner();
 562 }
 563
 564 void MEDCouplingFieldDiscretizationP0::computeMeshRestrictionFromTupleIds(const MEDCouplingMesh *mesh, const int *tupleIdsBg, const int *tupleIdsEnd,
 565                                                                           DataArrayInt *&cellRestriction, DataArrayInt *&trueTupleRestriction) const throw(INTERP_KERNEL::Exception)
 566 {
 567   if(!mesh)
 568     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP0::computeMeshRestrictionFromTupleIds : NULL input mesh !");
 569   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp=DataArrayInt::New();
 570   tmp->alloc((int)std::distance(tupleIdsBg,tupleIdsEnd),1);
 571   std::copy(tupleIdsBg,tupleIdsEnd,tmp->getPointer());
 572   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp2(tmp->deepCpy());
 573   cellRestriction=tmp.retn();
 574   trueTupleRestriction=tmp2.retn();
 575 }
 576
 577 void MEDCouplingFieldDiscretizationP0::reprQuickOverview(std::ostream& stream) const throw(INTERP_KERNEL::Exception)
 578 {
 579   stream << "P0 spatial discretization.";
 580 }
 581
 582 void MEDCouplingFieldDiscretizationP0::checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception)
 583 {
 584 }
 585
 586 void MEDCouplingFieldDiscretizationP0::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArray *da) const throw(INTERP_KERNEL::Exception)
 587 {
 588   if(!mesh || !da)
 589     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP0::checkCoherencyBetween : NULL input mesh or DataArray !");
 590   if(mesh->getNumberOfCells()!=da->getNumberOfTuples())
 591     {
 592       std::ostringstream message;
 593       message << "Field on cells invalid because there are " << mesh->getNumberOfCells();
 594       message << " cells in mesh and " << da->getNumberOfTuples() << " tuples in field !";
 595       throw INTERP_KERNEL::Exception(message.str().c_str());
 596     }
 597 }
 598
 599 MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP0::getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const
 600 {
 601   if(!mesh)
 602     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP0::getMeasureField : mesh instance specified is NULL !");
 603   return mesh->getMeasureField(isAbs);
 604 }
 605
 606 void MEDCouplingFieldDiscretizationP0::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
 607 {
 608   if(!mesh)
 609     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP0::getValueOn : NULL input mesh !");
 610   int id=mesh->getCellContainingPoint(loc,_precision);
 611   if(id==-1)
 612     throw INTERP_KERNEL::Exception("Specified point is detected outside of mesh : unable to apply P0::getValueOn !");
 613   arr->getTuple(id,res);
 614 }
 615
 616 void MEDCouplingFieldDiscretizationP0::getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, int i, int j, int k, double *res) const
 617 {
 618   const MEDCouplingCMesh *meshC=dynamic_cast<const MEDCouplingCMesh *>(mesh);
 619   if(!meshC)
 620     throw INTERP_KERNEL::Exception("P0::getValueOnPos is only accessible for structured meshes !");
 621   int id=meshC->getCellIdFromPos(i,j,k);
 622   arr->getTuple(id,res);
 623 }
 624
 625 DataArrayDouble *MEDCouplingFieldDiscretizationP0::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, int nbOfPoints) const
 626 {
 627   if(!mesh)
 628     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP0::getValueOnMulti : NULL input mesh !");
 629   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> eltsArr,eltsIndexArr;
 630   mesh->getCellsContainingPoints(loc,nbOfPoints,_precision,eltsArr,eltsIndexArr);
 631   const int *elts(eltsArr->begin()),*eltsIndex(eltsIndexArr->begin());
 632   int spaceDim=mesh->getSpaceDimension();
 633   int nbOfComponents=arr->getNumberOfComponents();
 634   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret=DataArrayDouble::New();
 635   ret->alloc(nbOfPoints,nbOfComponents);
 636   double *ptToFill=ret->getPointer();
 637   for(int i=0;i<nbOfPoints;i++,ptToFill+=nbOfComponents)
 638     if(eltsIndex[i+1]-eltsIndex[i]>=1)
 639       arr->getTuple(elts[eltsIndex[i]],ptToFill);
 640     else
 641       {
 642         std::ostringstream oss; oss << "Point #" << i << " with coordinates : (";
 643         std::copy(loc+i*spaceDim,loc+(i+1)*spaceDim,std::ostream_iterator<double>(oss,", "));
 644         oss << ") detected outside mesh : unable to apply P0::getValueOnMulti ! ";
 645         throw INTERP_KERNEL::Exception(oss.str().c_str());
 646       }
 647   return ret.retn();
 648 }
 649
 650 /*!
 651  * Nothing to do. It's not a bug.
 652  */
 653 void MEDCouplingFieldDiscretizationP0::renumberValuesOnNodes(double , const int *, int newNbOfNodes, DataArrayDouble *) const
 654 {
 655 }
 656
 657 void MEDCouplingFieldDiscretizationP0::renumberValuesOnCells(double epsOnVals, const MEDCouplingMesh *mesh, const int *old2New, int newSz, DataArrayDouble *arr) const
 658 {
 659   RenumberEntitiesFromO2NArr(epsOnVals,old2New,newSz,arr,"Cell");
 660 }
 661
 662 void MEDCouplingFieldDiscretizationP0::renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const int *new2old, int newSz, DataArrayDouble *arr) const
 663 {
 664   RenumberEntitiesFromN2OArr(new2old,newSz,arr,"Cell");
 665 }
 666
 667 /*!
 668  * This method returns a tuple ids selection from cell ids selection [start;end).
 669  * This method is called by MEDCouplingFieldDiscretizationP0::buildSubMeshData to return parameter \b di.
 670  * Here for P0 it's very simple !
 671  *
 672  * \return a newly allocated array containing ids to select into the DataArrayDouble of the field.
 673  *
 674  */
 675 DataArrayInt *MEDCouplingFieldDiscretizationP0::computeTupleIdsToSelectFromCellIds(const MEDCouplingMesh *mesh, const int *startCellIds, const int *endCellIds) const
 676 {
 677   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
 678   ret->alloc((int)std::distance(startCellIds,endCellIds),1);
 679   std::copy(startCellIds,endCellIds,ret->getPointer());
 680   return ret.retn();
 681 }
 682
 683 /*!
 684  * This method returns a submesh of 'mesh' instance constituting cell ids contained in array defined as an interval [start;end).
 685  * @param di is an array returned that specifies entity ids (here cells ids) in mesh 'mesh' of entity in returned submesh.
 686  * Example : The first cell id of returned mesh has the (*di)[0] id in 'mesh'
 687  *
 688  * \sa MEDCouplingFieldDiscretizationP0::buildSubMeshDataRange
 689  */
 690 MEDCouplingMesh *MEDCouplingFieldDiscretizationP0::buildSubMeshData(const MEDCouplingMesh *mesh, const int *start, const int *end, DataArrayInt *&di) const
 691 {
 692   if(!mesh)
 693     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP0::buildSubMeshData : NULL input mesh !");
 694   MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> ret=mesh->buildPart(start,end);
 695   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> diSafe=DataArrayInt::New();
 696   diSafe->alloc((int)std::distance(start,end),1);
 697   std::copy(start,end,diSafe->getPointer());
 698   di=diSafe.retn();
 699   return ret.retn();
 700 }
 701
 702 /*!
 703  * This method is strictly equivalent to MEDCouplingFieldDiscretizationP0::buildSubMeshData except that it is optimized for input defined as a range of cell ids.
 704  *
 705  * \param [out] beginOut Valid only if \a di is NULL
 706  * \param [out] endOut Valid only if \a di is NULL
 707  * \param [out] stepOut Valid only if \a di is NULL
 708  * \param [out] di is an array returned that specifies entity ids (nodes, cells, Gauss points... ) in array if no output range is foundable.
 709  *
 710  * \sa MEDCouplingFieldDiscretizationP0::buildSubMeshData
 711  */
 712 MEDCouplingMesh *MEDCouplingFieldDiscretizationP0::buildSubMeshDataRange(const MEDCouplingMesh *mesh, int beginCellIds, int endCellIds, int stepCellIds, int& beginOut, int& endOut, int& stepOut, DataArrayInt *&di) const
 713 {
 714   if(!mesh)
 715     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP0::buildSubMeshDataRange : NULL input mesh !");
 716   MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> ret=mesh->buildPartRange(beginCellIds,endCellIds,stepCellIds);
 717   di=0; beginOut=beginCellIds; endOut=endCellIds; stepOut=stepCellIds;
 718   return ret.retn();
 719 }
 720
 721 int MEDCouplingFieldDiscretizationOnNodes::getNumberOfTuples(const MEDCouplingMesh *mesh) const throw(INTERP_KERNEL::Exception)
 722 {
 723   if(!mesh)
 724     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationNodes::getNumberOfTuples : NULL input mesh !");
 725   return mesh->getNumberOfNodes();
 726 }
 727
 728 /*!
 729  * mesh is not used here. It is not a bug !
 730  */
 731 int MEDCouplingFieldDiscretizationOnNodes::getNumberOfTuplesExpectedRegardingCode(const MEDCouplingMesh *mesh, const std::vector<int>& code, const std::vector<const DataArrayInt *>& idsPerType) const throw(INTERP_KERNEL::Exception)
 732 {
 733   if(code.size()%3!=0)
 734     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationOnNodes::getNumberOfTuplesExpectedRegardingCode : invalid input code !");
 735   int nbOfSplit=(int)idsPerType.size();
 736   int nbOfTypes=(int)code.size()/3;
 737   int ret=0;
 738   for(int i=0;i<nbOfTypes;i++)
 739     {
 740       int nbOfEltInChunk=code[3*i+1];
 741       if(nbOfEltInChunk<0)
 742         throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationOnNodes::getNumberOfTuplesExpectedRegardingCode : invalid input code ! presence of negative value in a type !");
 743       int pos=code[3*i+2];
 744       if(pos!=-1)
 745         {
 746           if(pos<0 || pos>=nbOfSplit)
 747             {
 748               std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationOnNodes::getNumberOfTuplesExpectedRegardingCode : input code points to pos " << pos << " in typeid " << i << " ! Should be in [0," << nbOfSplit << ") !";
 749               throw INTERP_KERNEL::Exception(oss.str().c_str());
 750             }
 751           const DataArrayInt *ids(idsPerType[pos]);
 752           if(!ids || !ids->isAllocated() || ids->getNumberOfComponents()!=1 || ids->getNumberOfTuples()!=nbOfEltInChunk || ids->getMinValueInArray()<0)
 753             {
 754               std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationOnNodes::getNumberOfTuplesExpectedRegardingCode : input pfl chunck at pos " << pos << " should have " << i << " tuples and one component and with ids all >=0 !";
 755               throw INTERP_KERNEL::Exception(oss.str().c_str());
 756             }
 757         }
 758       ret+=nbOfEltInChunk;
 759     }
 760   return ret;
 761 }
 762
 763 int MEDCouplingFieldDiscretizationOnNodes::getNumberOfMeshPlaces(const MEDCouplingMesh *mesh) const
 764 {
 765   if(!mesh)
 766     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationNodes::getNumberOfMeshPlaces : NULL input mesh !");
 767   return mesh->getNumberOfNodes();
 768 }
 769
 770 /*!
 771  * Nothing to do here.
 772  */
 773 void MEDCouplingFieldDiscretizationOnNodes::renumberArraysForCell(const MEDCouplingMesh *, const std::vector<DataArray *>& arrays,
 774                                                                   const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
 775 {
 776 }
 777
 778 DataArrayInt *MEDCouplingFieldDiscretizationOnNodes::getOffsetArr(const MEDCouplingMesh *mesh) const
 779 {
 780   if(!mesh)
 781     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationNodes::getOffsetArr : NULL input mesh !");
 782   int nbOfTuples=mesh->getNumberOfNodes();
 783   DataArrayInt *ret=DataArrayInt::New();
 784   ret->alloc(nbOfTuples+1,1);
 785   ret->iota(0);
 786   return ret;
 787 }
 788
 789 DataArrayDouble *MEDCouplingFieldDiscretizationOnNodes::getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const
 790 {
 791   if(!mesh)
 792     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationNodes::getLocalizationOfDiscValues : NULL input mesh !");
 793   return mesh->getCoordinatesAndOwner();
 794 }
 795
 796 void MEDCouplingFieldDiscretizationOnNodes::computeMeshRestrictionFromTupleIds(const MEDCouplingMesh *mesh, const int *tupleIdsBg, const int *tupleIdsEnd,
 797                                                                                DataArrayInt *&cellRestriction, DataArrayInt *&trueTupleRestriction) const throw(INTERP_KERNEL::Exception)
 798 {
 799   if(!mesh)
 800     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationOnNodes::computeMeshRestrictionFromTupleIds : NULL input mesh !");
 801   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret1=mesh->getCellIdsFullyIncludedInNodeIds(tupleIdsBg,tupleIdsEnd);
 802   const MEDCouplingUMesh *meshc=dynamic_cast<const MEDCouplingUMesh *>(mesh);
 803   if(!meshc)
 804     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationOnNodes::computeMeshRestrictionFromTupleIds : trying to subpart field on nodes by node ids ! Your mesh has to be unstructured !");
 805   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> meshPart=static_cast<MEDCouplingUMesh *>(meshc->buildPartOfMySelf(ret1->begin(),ret1->end(),true));
 806   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret2=meshPart->computeFetchedNodeIds();
 807   cellRestriction=ret1.retn();
 808   trueTupleRestriction=ret2.retn();
 809 }
 810
 811 void MEDCouplingFieldDiscretizationOnNodes::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArray *da) const throw(INTERP_KERNEL::Exception)
 812 {
 813   if(!mesh || !da)
 814     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationNodes::checkCoherencyBetween : NULL input mesh or DataArray !");
 815   if(mesh->getNumberOfNodes()!=da->getNumberOfTuples())
 816     {
 817       std::ostringstream message;
 818       message << "Field on nodes invalid because there are " << mesh->getNumberOfNodes();
 819       message << " nodes in mesh and " << da->getNumberOfTuples() << " tuples in field !";
 820       throw INTERP_KERNEL::Exception(message.str().c_str());
 821     }
 822 }
 823
 824 /*!
 825  * This method returns a submesh of 'mesh' instance constituting cell ids contained in array defined as an interval [start;end).
 826 * @param di is an array returned that specifies entity ids (here nodes ids) in mesh 'mesh' of entity in returned submesh.
 827  * Example : The first node id of returned mesh has the (*di)[0] id in 'mesh'
 828  */
 829 MEDCouplingMesh *MEDCouplingFieldDiscretizationOnNodes::buildSubMeshData(const MEDCouplingMesh *mesh, const int *start, const int *end, DataArrayInt *&di) const
 830 {
 831   if(!mesh)
 832     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationNodes::buildSubMeshData : NULL input mesh !");
 833   DataArrayInt *diTmp=0;
 834   MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> ret=mesh->buildPartAndReduceNodes(start,end,diTmp);
 835   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> diTmpSafe(diTmp);
 836   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> di2=diTmpSafe->invertArrayO2N2N2O(ret->getNumberOfNodes());
 837   di=di2.retn();
 838   return ret.retn();
 839 }
 840
 841 /*!
 842  * This method is strictly equivalent to MEDCouplingFieldDiscretizationNodes::buildSubMeshData except that it is optimized for input defined as a range of cell ids.
 843  *
 844  * \param [out] beginOut Valid only if \a di is NULL
 845  * \param [out] endOut Valid only if \a di is NULL
 846  * \param [out] stepOut Valid only if \a di is NULL
 847  * \param [out] di is an array returned that specifies entity ids (nodes, cells, Gauss points... ) in array if no output range is foundable.
 848  *
 849  * \sa MEDCouplingFieldDiscretizationNodes::buildSubMeshData
 850  */
 851 MEDCouplingMesh *MEDCouplingFieldDiscretizationOnNodes::buildSubMeshDataRange(const MEDCouplingMesh *mesh, int beginCellIds, int endCellIds, int stepCellIds, int& beginOut, int& endOut, int& stepOut, DataArrayInt *&di) const
 852 {
 853   if(!mesh)
 854     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationOnNodes::buildSubMeshDataRange : NULL input mesh !");
 855   DataArrayInt *diTmp=0;
 856   MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> ret=mesh->buildPartRangeAndReduceNodes(beginCellIds,endCellIds,stepCellIds,beginOut,endOut,stepOut,diTmp);
 857   if(diTmp)
 858     {
 859       MEDCouplingAutoRefCountObjectPtr<DataArrayInt> diTmpSafe(diTmp);
 860       MEDCouplingAutoRefCountObjectPtr<DataArrayInt> di2=diTmpSafe->invertArrayO2N2N2O(ret->getNumberOfNodes());
 861       di=di2.retn();
 862     }
 863   return ret.retn();
 864 }
 865
 866 /*!
 867  * This method returns a tuple ids selection from cell ids selection [start;end).
 868  * This method is called by MEDCouplingFieldDiscretizationOnNodes::buildSubMeshData to return parameter \b di.
 869  * Here for P1 only nodes fetched by submesh of mesh[startCellIds:endCellIds) is returned !
 870  *
 871  * \return a newly allocated array containing ids to select into the DataArrayDouble of the field.
 872  *
 873  */
 874 DataArrayInt *MEDCouplingFieldDiscretizationOnNodes::computeTupleIdsToSelectFromCellIds(const MEDCouplingMesh *mesh, const int *startCellIds, const int *endCellIds) const
 875 {
 876   if(!mesh)
 877     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP1::computeTupleIdsToSelectFromCellIds : NULL input mesh !");
 878   const MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> umesh=mesh->buildUnstructured();
 879   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> umesh2=static_cast<MEDCouplingUMesh *>(umesh->buildPartOfMySelf(startCellIds,endCellIds,true));
 880   return umesh2->computeFetchedNodeIds();
 881 }
 882
 883 void MEDCouplingFieldDiscretizationOnNodes::renumberValuesOnNodes(double epsOnVals, const int *old2NewPtr, int newNbOfNodes, DataArrayDouble *arr) const
 884 {
 885   RenumberEntitiesFromO2NArr(epsOnVals,old2NewPtr,newNbOfNodes,arr,"Node");
 886 }
 887
 888 /*!
 889  * Nothing to do it's not a bug.
 890  */
 891 void MEDCouplingFieldDiscretizationOnNodes::renumberValuesOnCells(double epsOnVals, const MEDCouplingMesh *mesh, const int *old2New, int newSz, DataArrayDouble *arr) const
 892 {
 893 }
 894
 895 /*!
 896  * Nothing to do it's not a bug.
 897  */
 898 void MEDCouplingFieldDiscretizationOnNodes::renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const int *new2old, int newSz, DataArrayDouble *arr) const
 899 {
 900 }
 901
 902 void MEDCouplingFieldDiscretizationOnNodes::getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, int i, int j, int k, double *res) const
 903 {
 904   const MEDCouplingCMesh *meshC=dynamic_cast<const MEDCouplingCMesh *>(mesh);
 905   if(!meshC)
 906     throw INTERP_KERNEL::Exception("OnNodes::getValueOnPos(i,j,k) is only accessible for structured meshes !");
 907   int id=meshC->getNodeIdFromPos(i,j,k);
 908   arr->getTuple(id,res);
 909 }
 910
 911 TypeOfField MEDCouplingFieldDiscretizationP1::getEnum() const
 912 {
 913   return TYPE;
 914 }
 915
 916 /*!
 917  * This method is simply called by MEDCouplingFieldDiscretization::deepCpy. It performs the deep copy of \a this.
 918  *
 919  * \sa MEDCouplingFieldDiscretization::deepCpy.
 920  */
 921 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationP1::clone() const
 922 {
 923   return new MEDCouplingFieldDiscretizationP1;
 924 }
 925
 926 std::string MEDCouplingFieldDiscretizationP1::getStringRepr() const
 927 {
 928   return std::string(REPR);
 929 }
 930
 931 const char *MEDCouplingFieldDiscretizationP1::getRepr() const
 932 {
 933   return REPR;
 934 }
 935
 936 bool MEDCouplingFieldDiscretizationP1::isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const
 937 {
 938   if(!other)
 939     {
 940       reason="other spatial discretization is NULL, and this spatial discretization (P1) is defined.";
 941       return false;
 942     }
 943   const MEDCouplingFieldDiscretizationP1 *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationP1 *>(other);
 944   bool ret=otherC!=0;
 945   if(!ret)
 946     reason="Spatial discrtization of this is ON_NODES, which is not the case of other.";
 947   return ret;
 948 }
 949
 950 void MEDCouplingFieldDiscretizationP1::checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception)
 951 {
 952   if(nat!=ConservativeVolumic)
 953     throw INTERP_KERNEL::Exception("Invalid nature for P1 field  : expected ConservativeVolumic !");
 954 }
 955
 956 MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationP1::getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const
 957 {
 958   if(!mesh)
 959     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP1::getMeasureField : mesh instance specified is NULL !");
 960   return mesh->getMeasureFieldOnNode(isAbs);
 961 }
 962
 963 void MEDCouplingFieldDiscretizationP1::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
 964 {
 965   if(!mesh)
 966     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP1::getValueOn : NULL input mesh !");
 967   int id=mesh->getCellContainingPoint(loc,_precision);
 968   if(id==-1)
 969     throw INTERP_KERNEL::Exception("Specified point is detected outside of mesh : unable to apply P1::getValueOn !");
 970   INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell(id);
 971   if(type!=INTERP_KERNEL::NORM_SEG2 && type!=INTERP_KERNEL::NORM_TRI3 && type!=INTERP_KERNEL::NORM_TETRA4)
 972     throw INTERP_KERNEL::Exception("P1 getValueOn is not specified for not simplex cells !");
 973   getValueInCell(mesh,id,arr,loc,res);
 974 }
 975
 976 /*!
 977  * This method localizes a point defined by 'loc' in a cell with id 'cellId' into mesh 'mesh'.
 978  * The result is put into res expected to be of size at least arr->getNumberOfComponents()
 979  */
 980 void MEDCouplingFieldDiscretizationP1::getValueInCell(const MEDCouplingMesh *mesh, int cellId, const DataArrayDouble *arr, const double *loc, double *res) const
 981 {
 982   if(!mesh)
 983     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP1::getValueInCell : NULL input mesh !");
 984   std::vector<int> conn;
 985   std::vector<double> coo;
 986   mesh->getNodeIdsOfCell(cellId,conn);
 987   for(std::vector<int>::const_iterator iter=conn.begin();iter!=conn.end();iter++)
 988     mesh->getCoordinatesOfNode(*iter,coo);
 989   int spaceDim=mesh->getSpaceDimension();
 990   std::size_t nbOfNodes=conn.size();
 991   std::vector<const double *> vec(nbOfNodes);
 992   for(std::size_t i=0;i<nbOfNodes;i++)
 993     vec[i]=&coo[i*spaceDim];
 994   INTERP_KERNEL::AutoPtr<double> tmp=new double[nbOfNodes];
 995   INTERP_KERNEL::barycentric_coords(vec,loc,tmp);
 996   int sz=arr->getNumberOfComponents();
 997   INTERP_KERNEL::AutoPtr<double> tmp2=new double[sz];
 998   std::fill(res,res+sz,0.);
 999   for(std::size_t i=0;i<nbOfNodes;i++)
1000     {
1001       arr->getTuple(conn[i],(double *)tmp2);
1002       std::transform((double *)tmp2,((double *)tmp2)+sz,(double *)tmp2,std::bind2nd(std::multiplies<double>(),tmp[i]));
1003       std::transform(res,res+sz,(double *)tmp2,res,std::plus<double>());
1004     }
1005 }
1006
1007 DataArrayDouble *MEDCouplingFieldDiscretizationP1::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, int nbOfPoints) const
1008 {
1009   if(!mesh)
1010     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationP1::getValueOnMulti : NULL input mesh !");
1011   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> eltsArr,eltsIndexArr;
1012   mesh->getCellsContainingPoints(loc,nbOfPoints,_precision,eltsArr,eltsIndexArr);
1013   const int *elts(eltsArr->begin()),*eltsIndex(eltsIndexArr->begin());
1014   int spaceDim=mesh->getSpaceDimension();
1015   int nbOfComponents=arr->getNumberOfComponents();
1016   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret=DataArrayDouble::New();
1017   ret->alloc(nbOfPoints,nbOfComponents);
1018   double *ptToFill=ret->getPointer();
1019   for(int i=0;i<nbOfPoints;i++)
1020     if(eltsIndex[i+1]-eltsIndex[i]>=1)
1021       getValueInCell(mesh,elts[eltsIndex[i]],arr,loc+i*spaceDim,ptToFill+i*nbOfComponents);
1022     else
1023       {
1024         std::ostringstream oss; oss << "Point #" << i << " with coordinates : (";
1025         std::copy(loc+i*spaceDim,loc+(i+1)*spaceDim,std::ostream_iterator<double>(oss,", "));
1026         oss << ") detected outside mesh : unable to apply P1::getValueOnMulti ! ";
1027         throw INTERP_KERNEL::Exception(oss.str().c_str());
1028       }
1029   return ret.retn();
1030 }
1031
1032 void MEDCouplingFieldDiscretizationP1::reprQuickOverview(std::ostream& stream) const throw(INTERP_KERNEL::Exception)
1033 {
1034   stream << "P1 spatial discretization.";
1035 }
1036
1037 MEDCouplingFieldDiscretizationPerCell::MEDCouplingFieldDiscretizationPerCell():_discr_per_cell(0)
1038 {
1039 }
1040
1041 MEDCouplingFieldDiscretizationPerCell::~MEDCouplingFieldDiscretizationPerCell()
1042 {
1043   if(_discr_per_cell)
1044     _discr_per_cell->decrRef();
1045 }
1046
1047 /*!
1048  * This constructor deep copies ParaMEDMEM::DataArrayInt instance from other (if any).
1049  */
1050 MEDCouplingFieldDiscretizationPerCell::MEDCouplingFieldDiscretizationPerCell(const MEDCouplingFieldDiscretizationPerCell& other, const int *startCellIds, const int *endCellIds):_discr_per_cell(0)
1051 {
1052   DataArrayInt *arr=other._discr_per_cell;
1053   if(arr)
1054     {
1055       if(startCellIds==0 && endCellIds==0)
1056         _discr_per_cell=arr->deepCpy();
1057       else
1058         _discr_per_cell=arr->selectByTupleIdSafe(startCellIds,endCellIds);
1059     }
1060 }
1061
1062 MEDCouplingFieldDiscretizationPerCell::MEDCouplingFieldDiscretizationPerCell(const MEDCouplingFieldDiscretizationPerCell& other, int beginCellIds, int endCellIds, int stepCellIds):_discr_per_cell(0)
1063 {
1064   DataArrayInt *arr=other._discr_per_cell;
1065   if(arr)
1066     {
1067       _discr_per_cell=arr->selectByTupleId2(beginCellIds,endCellIds,stepCellIds);
1068     }
1069 }
1070
1071 void MEDCouplingFieldDiscretizationPerCell::updateTime() const
1072 {
1073   if(_discr_per_cell)
1074     updateTimeWith(*_discr_per_cell);
1075 }
1076
1077 std::size_t MEDCouplingFieldDiscretizationPerCell::getHeapMemorySizeWithoutChildren() const
1078 {
1079   std::size_t ret(MEDCouplingFieldDiscretization::getHeapMemorySizeWithoutChildren());
1080   return ret;
1081 }
1082
1083 std::vector<RefCountObject *> MEDCouplingFieldDiscretizationPerCell::getDirectChildren() const
1084 {
1085   std::vector<RefCountObject *> ret(MEDCouplingFieldDiscretization::getDirectChildren());
1086   if(_discr_per_cell)
1087     ret.push_back(const_cast<DataArrayInt *>(_discr_per_cell));
1088   return ret;
1089 }
1090
1091 void MEDCouplingFieldDiscretizationPerCell::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArray *da) const throw(INTERP_KERNEL::Exception)
1092 {
1093   if(!_discr_per_cell)
1094     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell has no discretization per cell !");
1095   if(!mesh)
1096     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell::checkCoherencyBetween : NULL input mesh or DataArray !");
1097   int nbOfTuples=_discr_per_cell->getNumberOfTuples();
1098   if(nbOfTuples!=mesh->getNumberOfCells())
1099     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell has a discretization per cell but it's not matching the underlying mesh !");
1100 }
1101
1102 bool MEDCouplingFieldDiscretizationPerCell::isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const
1103 {
1104   if(!other)
1105     {
1106       reason="other spatial discretization is NULL, and this spatial discretization (PerCell) is defined.";
1107       return false;
1108     }
1109   const MEDCouplingFieldDiscretizationPerCell *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationPerCell *>(other);
1110   if(!otherC)
1111     {
1112       reason="Spatial discretization of this is ON_GAUSS, which is not the case of other.";
1113       return false;
1114     }
1115   if(_discr_per_cell==0)
1116     return otherC->_discr_per_cell==0;
1117   if(otherC->_discr_per_cell==0)
1118     return false;
1119   bool ret=_discr_per_cell->isEqualIfNotWhy(*otherC->_discr_per_cell,reason);
1120   if(!ret)
1121     reason.insert(0,"Field discretization per cell DataArrayInt given the discid per cell :");
1122   return ret;
1123 }
1124
1125 bool MEDCouplingFieldDiscretizationPerCell::isEqualWithoutConsideringStr(const MEDCouplingFieldDiscretization *other, double eps) const
1126 {
1127   const MEDCouplingFieldDiscretizationPerCell *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationPerCell *>(other);
1128   if(!otherC)
1129     return false;
1130   if(_discr_per_cell==0)
1131     return otherC->_discr_per_cell==0;
1132   if(otherC->_discr_per_cell==0)
1133     return false;
1134   return _discr_per_cell->isEqualWithoutConsideringStr(*otherC->_discr_per_cell);
1135 }
1136
1137 /*!
1138  * This method is typically the first step of renumbering. The impact on _discr_per_cell is necessary here.
1139  * virtualy by this method.
1140  */
1141 void MEDCouplingFieldDiscretizationPerCell::renumberCells(const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
1142 {
1143   int nbCells=_discr_per_cell->getNumberOfTuples();
1144   const int *array=old2NewBg;
1145   if(check)
1146     array=DataArrayInt::CheckAndPreparePermutation(old2NewBg,old2NewBg+nbCells);
1147   //
1148   DataArrayInt *dpc=_discr_per_cell->renumber(array);
1149   _discr_per_cell->decrRef();
1150   _discr_per_cell=dpc;
1151   //
1152   if(check)
1153     free(const_cast<int *>(array));
1154 }
1155
1156 void MEDCouplingFieldDiscretizationPerCell::buildDiscrPerCellIfNecessary(const MEDCouplingMesh *mesh)
1157 {
1158   if(!mesh)
1159     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell::buildDiscrPerCellIfNecessary : NULL input mesh !");
1160   if(!_discr_per_cell)
1161     {
1162       _discr_per_cell=DataArrayInt::New();
1163       int nbTuples=mesh->getNumberOfCells();
1164       _discr_per_cell->alloc(nbTuples,1);
1165       int *ptr=_discr_per_cell->getPointer();
1166       std::fill(ptr,ptr+nbTuples,DFT_INVALID_LOCID_VALUE);
1167     }
1168 }
1169
1170 void MEDCouplingFieldDiscretizationPerCell::checkNoOrphanCells() const throw(INTERP_KERNEL::Exception)
1171 {
1172   if(!_discr_per_cell)
1173     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell::checkNoOrphanCells : no discretization defined !");
1174   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> test=_discr_per_cell->getIdsEqual(DFT_INVALID_LOCID_VALUE);
1175   if(test->getNumberOfTuples()!=0)
1176     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell::checkNoOrphanCells : presence of orphan cells !");
1177 }
1178
1179 /*!
1180  * This method is useful when 'this' describes a field discretization with several gauss discretization on a \b same cell type.
1181  * For example same NORM_TRI3 cells having 6 gauss points and others with 12 gauss points.
1182  * This method returns 2 arrays with same size : the return value and 'locIds' output parameter.
1183  * 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.
1184  * The return vector contains a set of newly created instance to deal with.
1185  * The returned vector represents a \b partition of cells ids with a gauss discretization set.
1186  *
1187  * If no descretization is set in 'this' and exception will be thrown.
1188  */
1189 std::vector<DataArrayInt *> MEDCouplingFieldDiscretizationPerCell::splitIntoSingleGaussDicrPerCellType(std::vector<int>& locIds) const throw(INTERP_KERNEL::Exception)
1190 {
1191   if(!_discr_per_cell)
1192     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell::splitIntoSingleGaussDicrPerCellType : no descretization set !");
1193   return _discr_per_cell->partitionByDifferentValues(locIds);
1194 }
1195
1196 const DataArrayInt *MEDCouplingFieldDiscretizationPerCell::getArrayOfDiscIds() const
1197 {
1198   return _discr_per_cell;
1199 }
1200
1201 void MEDCouplingFieldDiscretizationPerCell::setArrayOfDiscIds(const DataArrayInt *adids) throw(INTERP_KERNEL::Exception)
1202 {
1203   if(adids!=_discr_per_cell)
1204     {
1205       if(_discr_per_cell)
1206         _discr_per_cell->decrRef();
1207       _discr_per_cell=const_cast<DataArrayInt *>(adids);
1208       if(_discr_per_cell)
1209         _discr_per_cell->incrRef();
1210       declareAsNew();
1211     }
1212 }
1213
1214 MEDCouplingFieldDiscretizationGauss::MEDCouplingFieldDiscretizationGauss()
1215 {
1216 }
1217
1218 MEDCouplingFieldDiscretizationGauss::MEDCouplingFieldDiscretizationGauss(const MEDCouplingFieldDiscretizationGauss& other, const int *startCellIds, const int *endCellIds):MEDCouplingFieldDiscretizationPerCell(other,startCellIds,endCellIds),_loc(other._loc)
1219 {
1220 }
1221
1222 MEDCouplingFieldDiscretizationGauss::MEDCouplingFieldDiscretizationGauss(const MEDCouplingFieldDiscretizationGauss& other, int beginCellIds, int endCellIds, int stepCellIds):MEDCouplingFieldDiscretizationPerCell(other,beginCellIds,endCellIds,stepCellIds),_loc(other._loc)
1223 {
1224 }
1225
1226 TypeOfField MEDCouplingFieldDiscretizationGauss::getEnum() const
1227 {
1228   return TYPE;
1229 }
1230
1231 bool MEDCouplingFieldDiscretizationGauss::isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const
1232 {
1233   if(!other)
1234     {
1235       reason="other spatial discretization is NULL, and this spatial discretization (Gauss) is defined.";
1236       return false;
1237     }
1238   const MEDCouplingFieldDiscretizationGauss *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationGauss *>(other);
1239   if(!otherC)
1240     {
1241       reason="Spatial discrtization of this is ON_GAUSS, which is not the case of other.";
1242       return false;
1243     }
1244   if(!MEDCouplingFieldDiscretizationPerCell::isEqualIfNotWhy(other,eps,reason))
1245     return false;
1246   if(_loc.size()!=otherC->_loc.size())
1247     {
1248       reason="Gauss spatial discretization : localization sizes differ";
1249       return false;
1250     }
1251   std::size_t sz=_loc.size();
1252   for(std::size_t i=0;i<sz;i++)
1253     if(!_loc[i].isEqual(otherC->_loc[i],eps))
1254       {
1255         std::ostringstream oss; oss << "Gauss spatial discretization : Localization #" << i << " differ from this to other.";
1256         reason=oss.str();
1257         return false;
1258       }
1259   return true;
1260 }
1261
1262 bool MEDCouplingFieldDiscretizationGauss::isEqualWithoutConsideringStr(const MEDCouplingFieldDiscretization *other, double eps) const
1263 {
1264   const MEDCouplingFieldDiscretizationGauss *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationGauss *>(other);
1265   if(!otherC)
1266     return false;
1267   if(!MEDCouplingFieldDiscretizationPerCell::isEqualWithoutConsideringStr(other,eps))
1268     return false;
1269   if(_loc.size()!=otherC->_loc.size())
1270     return false;
1271   std::size_t sz=_loc.size();
1272   for(std::size_t i=0;i<sz;i++)
1273     if(!_loc[i].isEqual(otherC->_loc[i],eps))
1274       return false;
1275   return true;
1276 }
1277
1278 /*!
1279  * This method is simply called by MEDCouplingFieldDiscretization::deepCpy. It performs the deep copy of \a this.
1280  *
1281  * \sa MEDCouplingFieldDiscretization::deepCpy.
1282  */
1283 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationGauss::clone() const
1284 {
1285   return new MEDCouplingFieldDiscretizationGauss(*this);
1286 }
1287
1288 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationGauss::clonePart(const int *startCellIds, const int *endCellIds) const
1289 {
1290   return new MEDCouplingFieldDiscretizationGauss(*this,startCellIds,endCellIds);
1291 }
1292
1293 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationGauss::clonePartRange(int beginCellIds, int endCellIds, int stepCellIds) const
1294 {
1295   return new MEDCouplingFieldDiscretizationGauss(*this,beginCellIds,endCellIds,stepCellIds);
1296 }
1297
1298 std::string MEDCouplingFieldDiscretizationGauss::getStringRepr() const
1299 {
1300   std::ostringstream oss; oss << REPR << "." << std::endl;
1301   if(_discr_per_cell)
1302     {
1303       if(_discr_per_cell->isAllocated())
1304         {
1305           oss << "Discretization per cell : ";
1306           std::copy(_discr_per_cell->begin(),_discr_per_cell->end(),std::ostream_iterator<int>(oss,", "));
1307           oss << std::endl;
1308         }
1309     }
1310   oss << "Presence of " << _loc.size() << " localizations." << std::endl;
1311   int i=0;
1312   for(std::vector<MEDCouplingGaussLocalization>::const_iterator it=_loc.begin();it!=_loc.end();it++,i++)
1313     {
1314       oss << "+++++ Localization #" << i << " +++++" << std::endl;
1315       oss << (*it).getStringRepr();
1316       oss << "++++++++++" << std::endl;
1317     }
1318   return oss.str();
1319 }
1320
1321 std::size_t MEDCouplingFieldDiscretizationGauss::getHeapMemorySizeWithoutChildren() const
1322 {
1323   std::size_t ret(MEDCouplingFieldDiscretizationPerCell::getHeapMemorySizeWithoutChildren());
1324   ret+=_loc.capacity()*sizeof(MEDCouplingGaussLocalization);
1325   for(std::vector<MEDCouplingGaussLocalization>::const_iterator it=_loc.begin();it!=_loc.end();it++)
1326     ret+=(*it).getHeapMemorySize();
1327   return ret;
1328 }
1329
1330 const char *MEDCouplingFieldDiscretizationGauss::getRepr() const
1331 {
1332   return REPR;
1333 }
1334
1335 /*!
1336  * mesh is not used here. It is not a bug !
1337  */
1338 int MEDCouplingFieldDiscretizationGauss::getNumberOfTuplesExpectedRegardingCode(const MEDCouplingMesh *mesh, const std::vector<int>& code, const std::vector<const DataArrayInt *>& idsPerType) const throw(INTERP_KERNEL::Exception)
1339 {
1340   if(!_discr_per_cell || !_discr_per_cell->isAllocated() || _discr_per_cell->getNumberOfComponents()!=1)
1341     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getNumberOfTuplesExpectedRegardingCode");
1342   if(code.size()%3!=0)
1343     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getNumberOfTuplesExpectedRegardingCode : invalid input code !");
1344   int nbOfSplit=(int)idsPerType.size();
1345   int nbOfTypes=(int)code.size()/3;
1346   int ret=0;
1347   for(int i=0;i<nbOfTypes;i++)
1348     {
1349       int nbOfEltInChunk=code[3*i+1];
1350       if(nbOfEltInChunk<0)
1351         throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getNumberOfTuplesExpectedRegardingCode : invalid input code ! presence of negative value in a type !");
1352       int pos=code[3*i+2];
1353       if(pos!=-1)
1354         {
1355           if(pos<0 || pos>=nbOfSplit)
1356             {
1357               std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGauss::getNumberOfTuplesExpectedRegardingCode : input code points to pos " << pos << " in typeid " << i << " ! Should be in [0," << nbOfSplit << ") !";
1358               throw INTERP_KERNEL::Exception(oss.str().c_str());
1359             }
1360           const DataArrayInt *ids(idsPerType[pos]);
1361           if(!ids || !ids->isAllocated() || ids->getNumberOfComponents()!=1 || ids->getNumberOfTuples()!=nbOfEltInChunk || ids->getMinValueInArray()<0)
1362             {
1363               std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGauss::getNumberOfTuplesExpectedRegardingCode : input pfl chunck at pos " << pos << " should have " << i << " tuples and one component and with ids all >=0 !";
1364               throw INTERP_KERNEL::Exception(oss.str().c_str());
1365             }
1366         }
1367       ret+=nbOfEltInChunk;
1368     }
1369   if(ret!=_discr_per_cell->getNumberOfTuples())
1370     {
1371       std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGauss::getNumberOfTuplesExpectedRegardingCode : input code points to " << ret << " cells whereas discretization percell array lgth is " <<  _discr_per_cell->getNumberOfTuples() << " !";
1372     }
1373   return getNumberOfTuples(0);//0 is not an error ! It is to be sure that input mesh is not used
1374 }
1375
1376 int MEDCouplingFieldDiscretizationGauss::getNumberOfTuples(const MEDCouplingMesh *) const throw(INTERP_KERNEL::Exception)
1377 {
1378   int ret=0;
1379   if (_discr_per_cell == 0)
1380     throw INTERP_KERNEL::Exception("Discretization is not initialized!");
1381   const int *dcPtr=_discr_per_cell->getConstPointer();
1382   int nbOfTuples=_discr_per_cell->getNumberOfTuples();
1383   int maxSz=(int)_loc.size();
1384   for(const int *w=dcPtr;w!=dcPtr+nbOfTuples;w++)
1385     {
1386       if(*w>=0 && *w<maxSz)
1387         ret+=_loc[*w].getNumberOfGaussPt();
1388       else
1389         {
1390           std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGauss::getNumberOfTuples : At cell #" << std::distance(dcPtr,w) << " localization id is " << *w << " should be in [0," << maxSz << ") !";
1391           throw INTERP_KERNEL::Exception(oss.str().c_str());
1392         }
1393     }
1394   return ret;
1395 }
1396
1397 int MEDCouplingFieldDiscretizationGauss::getNumberOfMeshPlaces(const MEDCouplingMesh *mesh) const
1398 {
1399   if(!mesh)
1400     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getNumberOfMeshPlaces : NULL input mesh !");
1401   return mesh->getNumberOfCells();
1402 }
1403
1404 /*!
1405  * This method is redevelopped for performance reasons, but it is equivalent to a call to MEDCouplingFieldDiscretizationGauss::buildNbOfGaussPointPerCellField
1406  * and a call to DataArrayDouble::computeOffsets2 on the returned array.
1407  */
1408 DataArrayInt *MEDCouplingFieldDiscretizationGauss::getOffsetArr(const MEDCouplingMesh *mesh) const
1409 {
1410   if(!mesh)
1411     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getOffsetArr : NULL input mesh !");
1412   int nbOfTuples=mesh->getNumberOfCells();
1413   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
1414   ret->alloc(nbOfTuples+1,1);
1415   int *retPtr=ret->getPointer();
1416   const int *start=_discr_per_cell->getConstPointer();
1417   if(_discr_per_cell->getNumberOfTuples()!=nbOfTuples)
1418     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getOffsetArr : mismatch between the mesh and the discretization ids array length !");
1419   int maxPossible=(int)_loc.size();
1420   retPtr[0]=0;
1421   for(int i=0;i<nbOfTuples;i++,start++)
1422     {
1423       if(*start>=0 && *start<maxPossible)
1424         retPtr[i+1]=retPtr[i]+_loc[*start].getNumberOfGaussPt();
1425       else
1426         {
1427           std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGauss::getOffsetArr : At position #" << i << " the locid = " << *start << " whereas it should be in [0," << maxPossible << ") !";
1428           throw INTERP_KERNEL::Exception(oss.str().c_str());
1429         }
1430     }
1431   return ret.retn();
1432 }
1433
1434 void MEDCouplingFieldDiscretizationGauss::renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector<DataArray *>& arrays,
1435                                                                 const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
1436 {
1437   if(!mesh)
1438     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::renumberArraysForCell : NULL input mesh !");
1439   const int *array=old2NewBg;
1440   if(check)
1441     array=DataArrayInt::CheckAndPreparePermutation(old2NewBg,old2NewBg+mesh->getNumberOfCells());
1442   int nbOfCells=_discr_per_cell->getNumberOfTuples();
1443   int nbOfTuples=getNumberOfTuples(0);
1444   const int *dcPtr=_discr_per_cell->getConstPointer();
1445   int *array2=new int[nbOfTuples];//stores the final conversion array old2New to give to arrays in renumberInPlace.
1446   int *array3=new int[nbOfCells];//store for each cell in present dcp array (already renumbered) the offset needed by each cell in new numbering.
1447   array3[0]=0;
1448   for(int i=1;i<nbOfCells;i++)
1449     array3[i]=array3[i-1]+_loc[dcPtr[i-1]].getNumberOfGaussPt();
1450   int j=0;
1451   for(int i=0;i<nbOfCells;i++)
1452     {
1453       int nbOfGaussPt=_loc[dcPtr[array[i]]].getNumberOfGaussPt();
1454       for(int k=0;k<nbOfGaussPt;k++,j++)
1455         array2[j]=array3[array[i]]+k;
1456     }
1457   delete [] array3;
1458   for(std::vector<DataArray *>::const_iterator it=arrays.begin();it!=arrays.end();it++)
1459     if(*it)
1460       (*it)->renumberInPlace(array2);
1461   delete [] array2;
1462   if(check)
1463     free(const_cast<int*>(array));
1464 }
1465
1466 DataArrayDouble *MEDCouplingFieldDiscretizationGauss::getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const
1467 {
1468   if(!mesh)
1469     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getLocalizationOfDiscValues : NULL input mesh !");
1470   checkNoOrphanCells();
1471   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> umesh=mesh->buildUnstructured();//in general do nothing
1472   int nbOfTuples=getNumberOfTuples(mesh);
1473   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret=DataArrayDouble::New();
1474   int spaceDim=mesh->getSpaceDimension();
1475   ret->alloc(nbOfTuples,spaceDim);
1476   std::vector< int > locIds;
1477   std::vector<DataArrayInt *> parts=splitIntoSingleGaussDicrPerCellType(locIds);
1478   std::vector< MEDCouplingAutoRefCountObjectPtr<DataArrayInt> > parts2(parts.size());
1479   std::copy(parts.begin(),parts.end(),parts2.begin());
1480   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> offsets=buildNbOfGaussPointPerCellField();
1481   offsets->computeOffsets();
1482   const int *ptrOffsets=offsets->getConstPointer();
1483   const double *coords=umesh->getCoords()->getConstPointer();
1484   const int *connI=umesh->getNodalConnectivityIndex()->getConstPointer();
1485   const int *conn=umesh->getNodalConnectivity()->getConstPointer();
1486   double *valsToFill=ret->getPointer();
1487   for(std::size_t i=0;i<parts2.size();i++)
1488     {
1489       INTERP_KERNEL::GaussCoords calculator;
1490       //
1491       const MEDCouplingGaussLocalization& cli=_loc[locIds[i]];//curLocInfo
1492       INTERP_KERNEL::NormalizedCellType typ=cli.getType();
1493       const std::vector<double>& wg=cli.getWeights();
1494       calculator.addGaussInfo(typ,INTERP_KERNEL::CellModel::GetCellModel(typ).getDimension(),
1495                                   &cli.getGaussCoords()[0],(int)wg.size(),&cli.getRefCoords()[0],
1496                                   INTERP_KERNEL::CellModel::GetCellModel(typ).getNumberOfNodes());
1497       //
1498       int nbt=parts2[i]->getNumberOfTuples();
1499       for(const int *w=parts2[i]->getConstPointer();w!=parts2[i]->getConstPointer()+nbt;w++)
1500         calculator.calculateCoords(cli.getType(),coords,spaceDim,conn+connI[*w]+1,valsToFill+spaceDim*(ptrOffsets[*w]));
1501     }
1502   ret->copyStringInfoFrom(*umesh->getCoords());
1503   return ret.retn();
1504 }
1505
1506 void MEDCouplingFieldDiscretizationGauss::computeMeshRestrictionFromTupleIds(const MEDCouplingMesh *mesh, const int *tupleIdsBg, const int *tupleIdsEnd,
1507                                                                              DataArrayInt *&cellRestriction, DataArrayInt *&trueTupleRestriction) const throw(INTERP_KERNEL::Exception)
1508 {
1509   if(!mesh)
1510     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::computeMeshRestrictionFromTupleIds : NULL input mesh !");
1511   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp=DataArrayInt::New(); tmp->alloc((int)std::distance(tupleIdsBg,tupleIdsEnd),1);
1512   std::copy(tupleIdsBg,tupleIdsEnd,tmp->getPointer());
1513   tmp->sort(true);
1514   tmp=tmp->buildUnique();
1515   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nbOfNodesPerCell=buildNbOfGaussPointPerCellField();
1516   nbOfNodesPerCell->computeOffsets2();
1517   nbOfNodesPerCell->searchRangesInListOfIds(tmp,cellRestriction,trueTupleRestriction);
1518 }
1519
1520 /*!
1521  * Empty : not a bug
1522  */
1523 void MEDCouplingFieldDiscretizationGauss::checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception)
1524 {
1525 }
1526
1527 void MEDCouplingFieldDiscretizationGauss::getTinySerializationIntInformation(std::vector<int>& tinyInfo) const
1528 {
1529   int val=-1;
1530   if(_discr_per_cell)
1531     val=_discr_per_cell->getNumberOfTuples();
1532   tinyInfo.push_back(val);
1533   tinyInfo.push_back((int)_loc.size());
1534   if(_loc.empty())
1535     tinyInfo.push_back(-1);
1536   else
1537     tinyInfo.push_back(_loc[0].getDimension());
1538   for(std::vector<MEDCouplingGaussLocalization>::const_iterator iter=_loc.begin();iter!=_loc.end();iter++)
1539     (*iter).pushTinySerializationIntInfo(tinyInfo);
1540 }
1541
1542 void MEDCouplingFieldDiscretizationGauss::getTinySerializationDbleInformation(std::vector<double>& tinyInfo) const
1543 {
1544   for(std::vector<MEDCouplingGaussLocalization>::const_iterator iter=_loc.begin();iter!=_loc.end();iter++)
1545     (*iter).pushTinySerializationDblInfo(tinyInfo);
1546 }
1547
1548 void MEDCouplingFieldDiscretizationGauss::getSerializationIntArray(DataArrayInt *& arr) const
1549 {
1550   arr=0;
1551   if(_discr_per_cell)
1552     arr=_discr_per_cell;
1553 }
1554
1555 void MEDCouplingFieldDiscretizationGauss::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *& arr)
1556 {
1557   int val=tinyInfo[0];
1558   if(val>=0)
1559     {
1560       _discr_per_cell=DataArrayInt::New();
1561       _discr_per_cell->alloc(val,1);
1562     }
1563   else
1564     _discr_per_cell=0;
1565   arr=_discr_per_cell;
1566   int nbOfLoc=tinyInfo[1];
1567   _loc.clear();
1568   int dim=tinyInfo[2];
1569   int delta=-1;
1570   if(nbOfLoc>0)
1571     delta=((int)tinyInfo.size()-3)/nbOfLoc;
1572   for(int i=0;i<nbOfLoc;i++)
1573     {
1574       std::vector<int> tmp(tinyInfo.begin()+3+i*delta,tinyInfo.begin()+3+(i+1)*delta);
1575       MEDCouplingGaussLocalization elt=MEDCouplingGaussLocalization::BuildNewInstanceFromTinyInfo(dim,tmp);
1576       _loc.push_back(elt);
1577     }
1578 }
1579
1580 void MEDCouplingFieldDiscretizationGauss::finishUnserialization(const std::vector<double>& tinyInfo)
1581 {
1582   double *tmp=new double[tinyInfo.size()];
1583   std::copy(tinyInfo.begin(),tinyInfo.end(),tmp);
1584   const double *work=tmp;
1585   for(std::vector<MEDCouplingGaussLocalization>::iterator iter=_loc.begin();iter!=_loc.end();iter++)
1586     work=(*iter).fillWithValues(work);
1587   delete [] tmp;
1588 }
1589
1590 double MEDCouplingFieldDiscretizationGauss::getIJK(const MEDCouplingMesh *mesh, const DataArrayDouble *da,
1591                                                    int cellId, int nodeIdInCell, int compoId) const throw(INTERP_KERNEL::Exception)
1592 {
1593   int offset=getOffsetOfCell(cellId);
1594   return da->getIJ(offset+nodeIdInCell,compoId);
1595 }
1596
1597 void MEDCouplingFieldDiscretizationGauss::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArray *da) const throw(INTERP_KERNEL::Exception)
1598 {
1599   if(!mesh || !da)
1600     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::checkCoherencyBetween : NULL input mesh or DataArray !");
1601   MEDCouplingFieldDiscretizationPerCell::checkCoherencyBetween(mesh,da);
1602   for(std::vector<MEDCouplingGaussLocalization>::const_iterator iter=_loc.begin();iter!=_loc.end();iter++)
1603     (*iter).checkCoherency();
1604   int nbOfDesc=(int)_loc.size();
1605   int nbOfCells=mesh->getNumberOfCells();
1606   const int *dc=_discr_per_cell->getConstPointer();
1607   for(int i=0;i<nbOfCells;i++)
1608     {
1609       if(dc[i]>=nbOfDesc)
1610         {
1611           std::ostringstream oss; oss << "Cell # " << i << " of mesh \"" << mesh->getName() << "\" has an undefined gauss location ! Should never happend !";
1612           throw INTERP_KERNEL::Exception(oss.str().c_str());
1613         }
1614       if(dc[i]<0)
1615         {
1616           std::ostringstream oss; oss << "Cell # " << i << " of mesh \"" << mesh->getName() << "\" has no gauss location !";
1617           throw INTERP_KERNEL::Exception(oss.str().c_str());
1618         }
1619       if(mesh->getTypeOfCell(i)!=_loc[dc[i]].getType())
1620         {
1621           std::ostringstream oss; oss << "Types of mesh and gauss location mismatch for cell # " << i;
1622           throw INTERP_KERNEL::Exception(oss.str().c_str());
1623         }
1624     }
1625   int nbOfTuples=getNumberOfTuples(mesh);
1626   if(nbOfTuples!=da->getNumberOfTuples())
1627     {
1628       std::ostringstream oss; oss << "Invalid number of tuples in the array : expecting " << nbOfTuples << " !";
1629       throw INTERP_KERNEL::Exception(oss.str().c_str());
1630     }
1631 }
1632
1633 MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationGauss::getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const
1634 {
1635   if(!mesh)
1636     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getMeasureField : mesh instance specified is NULL !");
1637   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> vol=mesh->getMeasureField(isAbs);
1638   const double *volPtr=vol->getArray()->begin();
1639   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(ON_GAUSS_PT);
1640   ret->setMesh(mesh);
1641   ret->setDiscretization(const_cast<MEDCouplingFieldDiscretizationGauss *>(this));
1642   if(!_discr_per_cell)
1643     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getMeasureField : no discr per cell array not defined ! spatial localization is incorrect !");
1644   _discr_per_cell->checkAllocated();
1645   if(_discr_per_cell->getNumberOfComponents()!=1)
1646     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getMeasureField : no discr per cell array defined but with nb of components different from 1 !");
1647   if(_discr_per_cell->getNumberOfTuples()!=vol->getNumberOfTuples())
1648     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getMeasureField : no discr per cell array defined but mismatch between nb of cells of mesh and size of spatial disr array !");
1649   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> offset=getOffsetArr(mesh);
1650   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr=DataArrayDouble::New(); arr->alloc(getNumberOfTuples(mesh),1);
1651   ret->setArray(arr);
1652   double *arrPtr=arr->getPointer();
1653   const int *offsetPtr=offset->getConstPointer();
1654   int maxGaussLoc=(int)_loc.size();
1655   std::vector<int> locIds;
1656   std::vector<DataArrayInt *> ids=splitIntoSingleGaussDicrPerCellType(locIds);
1657   std::vector< MEDCouplingAutoRefCountObjectPtr<DataArrayInt> > ids2(ids.size()); std::copy(ids.begin(),ids.end(),ids2.begin());
1658   for(std::size_t i=0;i<locIds.size();i++)
1659     {
1660       const DataArrayInt *curIds=ids[i];
1661       int locId=locIds[i];
1662       if(locId>=0 && locId<maxGaussLoc)
1663         {
1664           const MEDCouplingGaussLocalization& loc=_loc[locId];
1665           int nbOfGaussPt=loc.getNumberOfGaussPt();
1666           INTERP_KERNEL::AutoPtr<double> weights=new double[nbOfGaussPt];
1667           double sum=std::accumulate(loc.getWeights().begin(),loc.getWeights().end(),0.);
1668           std::transform(loc.getWeights().begin(),loc.getWeights().end(),(double *)weights,std::bind2nd(std::multiplies<double>(),1./sum));
1669           for(const int *cellId=curIds->begin();cellId!=curIds->end();cellId++)
1670             for(int j=0;j<nbOfGaussPt;j++)
1671               arrPtr[offsetPtr[*cellId]+j]=weights[j]*volPtr[*cellId];
1672         }
1673       else
1674         {
1675           std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGauss::getMeasureField : Presence of localization id " << locId << " in cell #" << curIds->getIJ(0,0) << " ! Must be in [0," << maxGaussLoc << ") !";
1676           throw INTERP_KERNEL::Exception(oss.str().c_str());
1677         }
1678     }
1679   ret->synchronizeTimeWithSupport();
1680   return ret.retn();
1681 }
1682
1683 void MEDCouplingFieldDiscretizationGauss::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
1684 {
1685   throw INTERP_KERNEL::Exception("Not implemented yet !");
1686 }
1687
1688 void MEDCouplingFieldDiscretizationGauss::getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, int i, int j, int k, double *res) const
1689 {
1690   throw INTERP_KERNEL::Exception("getValueOnPos(i,j,k) : Not applyable for Gauss points !");
1691 }
1692
1693 DataArrayDouble *MEDCouplingFieldDiscretizationGauss::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, int nbOfPoints) const
1694 {
1695   throw INTERP_KERNEL::Exception("getValueOnMulti : Not implemented yet for gauss points !");
1696 }
1697
1698 MEDCouplingMesh *MEDCouplingFieldDiscretizationGauss::buildSubMeshData(const MEDCouplingMesh *mesh, const int *start, const int *end, DataArrayInt *&di) const
1699 {
1700   if(!mesh)
1701     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::buildSubMeshData : NULL input mesh !");
1702   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> diSafe=computeTupleIdsToSelectFromCellIds(mesh,start,end);
1703   MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> ret=mesh->buildPart(start,end);
1704   di=diSafe.retn();
1705   return ret.retn();
1706 }
1707
1708 /*!
1709  * This method is strictly equivalent to MEDCouplingFieldDiscretizationGauss::buildSubMeshData except that it is optimized for input defined as a range of cell ids.
1710  *
1711  * \param [out] beginOut Valid only if \a di is NULL
1712  * \param [out] endOut Valid only if \a di is NULL
1713  * \param [out] stepOut Valid only if \a di is NULL
1714  * \param [out] di is an array returned that specifies entity ids (nodes, cells, Gauss points... ) in array if no output range is foundable.
1715  *
1716  * \sa MEDCouplingFieldDiscretizationGauss::buildSubMeshData
1717  */
1718 MEDCouplingMesh *MEDCouplingFieldDiscretizationGauss::buildSubMeshDataRange(const MEDCouplingMesh *mesh, int beginCellIds, int endCellIds, int stepCellIds, int& beginOut, int& endOut, int& stepOut, DataArrayInt *&di) const
1719 {
1720   if(stepCellIds!=1)//even for stepCellIds==-1 the output will not be a range
1721     return MEDCouplingFieldDiscretization::buildSubMeshDataRange(mesh,beginCellIds,endCellIds,stepCellIds,beginOut,endOut,stepOut,di);
1722   if(!mesh)
1723     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::buildSubMeshDataRange : NULL input mesh !");
1724   if(!_discr_per_cell)
1725     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::buildSubMeshDataRange : no discretization array set !");
1726   di=0; beginOut=0; endOut=0; stepOut=stepCellIds;
1727   const char msg[]="MEDCouplingFieldDiscretizationGauss::buildSubMeshDataRange : cell #";
1728   int nbOfTuples=_discr_per_cell->getNumberOfTuples();
1729   const int *w=_discr_per_cell->begin();
1730   int nbMaxOfLocId=(int)_loc.size();
1731   for(int i=0;i<nbOfTuples;i++,w++)
1732     {
1733       if(*w!=DFT_INVALID_LOCID_VALUE)
1734         {
1735           if(*w>=0 && *w<nbMaxOfLocId)
1736             {
1737               int delta=_loc[*w].getNumberOfGaussPt();
1738               if(i<beginCellIds)
1739                 beginOut+=delta;
1740               endOut+=delta;
1741               if(i>=endCellIds)
1742                 break;
1743             }
1744           else
1745             { std::ostringstream oss; oss << msg << i << " has invalid id (" << *w << ") ! Should be in [0," << nbMaxOfLocId << ") !"; throw INTERP_KERNEL::Exception(oss.str().c_str()); }
1746         }
1747       else
1748         { std::ostringstream oss; oss << msg << i << " is detected as orphan !"; throw INTERP_KERNEL::Exception(oss.str().c_str()); }
1749     }
1750   MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> ret=mesh->buildPartRange(beginCellIds,endCellIds,stepCellIds);
1751   return ret.retn();
1752 }
1753
1754 /*!
1755  * This method returns a tuple ids selection from cell ids selection [start;end).
1756  * This method is called by MEDCouplingFieldDiscretizationGauss::buildSubMeshData to return parameter \b di.
1757  *
1758  * \return a newly allocated array containing ids to select into the DataArrayDouble of the field.
1759  *
1760  */
1761 DataArrayInt *MEDCouplingFieldDiscretizationGauss::computeTupleIdsToSelectFromCellIds(const MEDCouplingMesh *mesh, const int *startCellIds, const int *endCellIds) const
1762 {
1763   if(!mesh)
1764     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::computeTupleIdsToSelectFromCellIds : null mesh !");
1765   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nbOfNodesPerCell=buildNbOfGaussPointPerCellField();//check of _discr_per_cell not NULL pointer
1766   int nbOfCells=mesh->getNumberOfCells();
1767   if(_discr_per_cell->getNumberOfTuples()!=nbOfCells)
1768     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::computeTupleIdsToSelectFromCellIds : mismatch of nb of tuples of cell ids array and number of cells !");
1769   nbOfNodesPerCell->computeOffsets2();
1770   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> sel=DataArrayInt::New(); sel->useArray(startCellIds,false,CPP_DEALLOC,(int)std::distance(startCellIds,endCellIds),1);
1771   return sel->buildExplicitArrByRanges(nbOfNodesPerCell);
1772 }
1773
1774 /*!
1775  * No implementation needed !
1776  */
1777 void MEDCouplingFieldDiscretizationGauss::renumberValuesOnNodes(double , const int *, int newNbOfNodes, DataArrayDouble *) const
1778 {
1779 }
1780
1781 void MEDCouplingFieldDiscretizationGauss::renumberValuesOnCells(double epsOnVals, const MEDCouplingMesh *mesh, const int *old2New, int newSz, DataArrayDouble *arr) const
1782 {
1783   throw INTERP_KERNEL::Exception("Not implemented yet !");
1784 }
1785
1786 void MEDCouplingFieldDiscretizationGauss::renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const int *new2old, int newSz, DataArrayDouble *arr) const
1787 {
1788   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 !");
1789 }
1790
1791 void MEDCouplingFieldDiscretizationGauss::setGaussLocalizationOnType(const MEDCouplingMesh *mesh, INTERP_KERNEL::NormalizedCellType type, const std::vector<double>& refCoo,
1792                                                                      const std::vector<double>& gsCoo, const std::vector<double>& wg) throw(INTERP_KERNEL::Exception)
1793 {
1794   if(!mesh)
1795     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::setGaussLocalizationOnType : NULL input mesh !");
1796   const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
1797   if((int)cm.getDimension()!=mesh->getMeshDimension())
1798     {
1799       std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGauss::setGaussLocalizationOnType : mismatch of dimensions ! MeshDim==" << mesh->getMeshDimension();
1800       oss << " whereas Type '" << cm.getRepr() << "' has dimension " << cm.getDimension() << " !";
1801       throw INTERP_KERNEL::Exception(oss.str().c_str());
1802     }
1803   buildDiscrPerCellIfNecessary(mesh);
1804   int id=(int)_loc.size();
1805   MEDCouplingGaussLocalization elt(type,refCoo,gsCoo,wg);
1806   _loc.push_back(elt);
1807   int *ptr=_discr_per_cell->getPointer();
1808   int nbCells=mesh->getNumberOfCells();
1809   for(int i=0;i<nbCells;i++)
1810     if(mesh->getTypeOfCell(i)==type)
1811       ptr[i]=id;
1812   zipGaussLocalizations();
1813 }
1814
1815 void MEDCouplingFieldDiscretizationGauss::setGaussLocalizationOnCells(const MEDCouplingMesh *mesh, const int *begin, const int *end, const std::vector<double>& refCoo,
1816                                                                       const std::vector<double>& gsCoo, const std::vector<double>& wg) throw(INTERP_KERNEL::Exception)
1817 {
1818   if(!mesh)
1819     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::setGaussLocalizationOnCells : NULL input mesh !");
1820   buildDiscrPerCellIfNecessary(mesh);
1821   if(std::distance(begin,end)<1)
1822     throw INTERP_KERNEL::Exception("Size of [begin,end) must be equal or greater than 1 !");
1823   INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell(*begin);
1824   MEDCouplingGaussLocalization elt(type,refCoo,gsCoo,wg);
1825   int id=(int)_loc.size();
1826   int *ptr=_discr_per_cell->getPointer();
1827   for(const int *w=begin+1;w!=end;w++)
1828     {
1829       if(mesh->getTypeOfCell(*w)!=type)
1830         {
1831           std::ostringstream oss; oss << "The cell with id " << *w << " has been detected to be incompatible in the [begin,end) array specified !";
1832           throw INTERP_KERNEL::Exception(oss.str().c_str());
1833         }
1834     }
1835   //
1836   for(const int *w2=begin;w2!=end;w2++)
1837     ptr[*w2]=id;
1838   //
1839   _loc.push_back(elt);
1840   zipGaussLocalizations();
1841 }
1842
1843 void MEDCouplingFieldDiscretizationGauss::clearGaussLocalizations() throw(INTERP_KERNEL::Exception)
1844 {
1845   if(_discr_per_cell)
1846     {
1847       _discr_per_cell->decrRef();
1848       _discr_per_cell=0;
1849     }
1850   _loc.clear();
1851 }
1852
1853 void MEDCouplingFieldDiscretizationGauss::setGaussLocalization(int locId, const MEDCouplingGaussLocalization& loc) throw(INTERP_KERNEL::Exception)
1854 {
1855   if(locId<0)
1856     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::setGaussLocalization : localization id has to be >=0 !");
1857   int sz=(int)_loc.size();
1858   MEDCouplingGaussLocalization gLoc(INTERP_KERNEL::NORM_ERROR);
1859   if(locId>=sz)
1860     _loc.resize(locId+1,gLoc);
1861   _loc[locId]=loc;
1862 }
1863
1864 void MEDCouplingFieldDiscretizationGauss::resizeLocalizationVector(int newSz) throw(INTERP_KERNEL::Exception)
1865 {
1866   if(newSz<0)
1867     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::resizeLocalizationVector : new size has to be >=0 !");
1868   MEDCouplingGaussLocalization gLoc(INTERP_KERNEL::NORM_ERROR);
1869   _loc.resize(newSz,gLoc);
1870 }
1871
1872 MEDCouplingGaussLocalization& MEDCouplingFieldDiscretizationGauss::getGaussLocalization(int locId) throw(INTERP_KERNEL::Exception)
1873 {
1874   checkLocalizationId(locId);
1875   return _loc[locId];
1876 }
1877
1878 int MEDCouplingFieldDiscretizationGauss::getNbOfGaussLocalization() const throw(INTERP_KERNEL::Exception)
1879 {
1880   return (int)_loc.size();
1881 }
1882
1883 int MEDCouplingFieldDiscretizationGauss::getGaussLocalizationIdOfOneCell(int cellId) const throw(INTERP_KERNEL::Exception)
1884 {
1885   if(!_discr_per_cell)
1886     throw INTERP_KERNEL::Exception("No Gauss localization still set !");
1887   int locId=_discr_per_cell->begin()[cellId];
1888   if(locId<0)
1889     throw INTERP_KERNEL::Exception("No Gauss localization set for the specified cell !");
1890   return locId;
1891 }
1892
1893 int MEDCouplingFieldDiscretizationGauss::getGaussLocalizationIdOfOneType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception)
1894 {
1895   std::set<int> ret=getGaussLocalizationIdsOfOneType(type);
1896   if(ret.empty())
1897     throw INTERP_KERNEL::Exception("No gauss discretization found for the specified type !");
1898   if(ret.size()>1)
1899     throw INTERP_KERNEL::Exception("Several gauss discretizations have been found for the specified type !");
1900   return *ret.begin();
1901 }
1902
1903 std::set<int> MEDCouplingFieldDiscretizationGauss::getGaussLocalizationIdsOfOneType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception)
1904 {
1905   if(!_discr_per_cell)
1906     throw INTERP_KERNEL::Exception("No Gauss localization still set !");
1907   std::set<int> ret;
1908   int id=0;
1909   for(std::vector<MEDCouplingGaussLocalization>::const_iterator iter=_loc.begin();iter!=_loc.end();iter++,id++)
1910     if((*iter).getType()==type)
1911       ret.insert(id);
1912   return ret;
1913 }
1914
1915 void MEDCouplingFieldDiscretizationGauss::getCellIdsHavingGaussLocalization(int locId, std::vector<int>& cellIds) const throw(INTERP_KERNEL::Exception)
1916 {
1917   if(locId<0 || locId>=(int)_loc.size())
1918     throw INTERP_KERNEL::Exception("Invalid locId given : must be in range [0:getNbOfGaussLocalization()) !");
1919   int nbOfTuples=_discr_per_cell->getNumberOfTuples();
1920   const int *ptr=_discr_per_cell->getConstPointer();
1921   for(int i=0;i<nbOfTuples;i++)
1922     if(ptr[i]==locId)
1923       cellIds.push_back(i);
1924 }
1925
1926 const MEDCouplingGaussLocalization& MEDCouplingFieldDiscretizationGauss::getGaussLocalization(int locId) const throw(INTERP_KERNEL::Exception)
1927 {
1928   checkLocalizationId(locId);
1929   return _loc[locId];
1930 }
1931
1932 void MEDCouplingFieldDiscretizationGauss::checkLocalizationId(int locId) const throw(INTERP_KERNEL::Exception)
1933 {
1934   if(locId<0 || locId>=(int)_loc.size())
1935     throw INTERP_KERNEL::Exception("Invalid locId given : must be in range [0:getNbOfGaussLocalization()) !");
1936 }
1937
1938 int MEDCouplingFieldDiscretizationGauss::getOffsetOfCell(int cellId) const throw(INTERP_KERNEL::Exception)
1939 {
1940   int ret=0;
1941   const int *start=_discr_per_cell->getConstPointer();
1942   for(const int *w=start;w!=start+cellId;w++)
1943     ret+=_loc[*w].getNumberOfGaussPt();
1944   return ret;
1945 }
1946
1947 /*!
1948  * This method do the assumption that there is no orphan cell. If there is an exception is thrown.
1949  * This method makes the assumption too that '_discr_per_cell' is defined. If not an exception is thrown.
1950  * This method returns a newly created array with number of tuples equals to '_discr_per_cell->getNumberOfTuples' and number of components equal to 1.
1951  * The i_th tuple in returned array is the number of gauss point if the corresponding cell.
1952  */
1953 DataArrayInt *MEDCouplingFieldDiscretizationGauss::buildNbOfGaussPointPerCellField() const throw(INTERP_KERNEL::Exception)
1954 {
1955   if(!_discr_per_cell)
1956     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::buildNbOfGaussPointPerCellField : no discretization array set !");
1957   int nbOfTuples=_discr_per_cell->getNumberOfTuples();
1958   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
1959   const int *w=_discr_per_cell->begin();
1960   ret->alloc(nbOfTuples,1);
1961   int *valsToFill=ret->getPointer();
1962   int nbMaxOfLocId=(int)_loc.size();
1963   for(int i=0;i<nbOfTuples;i++,w++)
1964     if(*w!=DFT_INVALID_LOCID_VALUE)
1965       {
1966         if(*w>=0 && *w<nbMaxOfLocId)
1967           valsToFill[i]=_loc[*w].getNumberOfGaussPt();
1968         else
1969           {
1970             std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGauss::buildNbOfGaussPointPerCellField : cell #" << i << " has invalid id (" << *w << ") ! Should be in [0," << nbMaxOfLocId << ") !";
1971             throw INTERP_KERNEL::Exception(oss.str().c_str());
1972           }
1973       }
1974     else
1975       {
1976         std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGauss::buildNbOfGaussPointPerCellField : cell #" << i << " is detected as orphan !";
1977         throw INTERP_KERNEL::Exception(oss.str().c_str());
1978       }
1979   return ret.retn();
1980 }
1981
1982 void MEDCouplingFieldDiscretizationGauss::reprQuickOverview(std::ostream& stream) const throw(INTERP_KERNEL::Exception)
1983 {
1984   stream << "Gauss points spatial discretization.";
1985 }
1986
1987 /*!
1988  * This method makes the assumption that _discr_per_cell is set.
1989  * This method reduces as much as possible number size of _loc.
1990  * This method is useful when several set on same cells has been done and that some Gauss Localization are no more used.
1991  */
1992 void MEDCouplingFieldDiscretizationGauss::zipGaussLocalizations()
1993 {
1994   const int *start=_discr_per_cell->begin();
1995   int nbOfTuples=_discr_per_cell->getNumberOfTuples();
1996   INTERP_KERNEL::AutoPtr<int> tmp=new int[_loc.size()];
1997   std::fill((int *)tmp,(int *)tmp+_loc.size(),-2);
1998   for(const int *w=start;w!=start+nbOfTuples;w++)
1999     if(*w>=0)
2000       tmp[*w]=1;
2001   int fid=0;
2002   for(int i=0;i<(int)_loc.size();i++)
2003     if(tmp[i]!=-2)
2004       tmp[i]=fid++;
2005   if(fid==(int)_loc.size())
2006     return;
2007   // zip needed
2008   int *start2=_discr_per_cell->getPointer();
2009   for(int *w2=start2;w2!=start2+nbOfTuples;w2++)
2010     if(*w2>=0)
2011       *w2=tmp[*w2];
2012   std::vector<MEDCouplingGaussLocalization> tmpLoc;
2013   for(int i=0;i<(int)_loc.size();i++)
2014     if(tmp[i]!=-2)
2015       tmpLoc.push_back(_loc[tmp[i]]);
2016   _loc=tmpLoc;
2017 }
2018
2019 MEDCouplingFieldDiscretizationGaussNE::MEDCouplingFieldDiscretizationGaussNE()
2020 {
2021 }
2022
2023 TypeOfField MEDCouplingFieldDiscretizationGaussNE::getEnum() const
2024 {
2025   return TYPE;
2026 }
2027
2028 /*!
2029  * This method is simply called by MEDCouplingFieldDiscretization::deepCpy. It performs the deep copy of \a this.
2030  *
2031  * \sa MEDCouplingFieldDiscretization::deepCpy.
2032  */
2033 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationGaussNE::clone() const
2034 {
2035   return new MEDCouplingFieldDiscretizationGaussNE(*this);
2036 }
2037
2038 std::string MEDCouplingFieldDiscretizationGaussNE::getStringRepr() const
2039 {
2040   return std::string(REPR);
2041 }
2042
2043 const char *MEDCouplingFieldDiscretizationGaussNE::getRepr() const
2044 {
2045   return REPR;
2046 }
2047
2048 bool MEDCouplingFieldDiscretizationGaussNE::isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const
2049 {
2050   if(!other)
2051     {
2052       reason="other spatial discretization is NULL, and this spatial discretization (GaussNE) is defined.";
2053       return false;
2054     }
2055   const MEDCouplingFieldDiscretizationGaussNE *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationGaussNE *>(other);
2056   bool ret=otherC!=0;
2057   if(!ret)
2058     reason="Spatial discrtization of this is ON_GAUSS_NE, which is not the case of other.";
2059   return ret;
2060 }
2061
2062 int MEDCouplingFieldDiscretizationGaussNE::getNumberOfTuplesExpectedRegardingCode(const MEDCouplingMesh *mesh, const std::vector<int>& code, const std::vector<const DataArrayInt *>& idsPerType) const throw(INTERP_KERNEL::Exception)
2063 {
2064   if(code.size()%3!=0)
2065     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::getNumberOfTuplesExpectedRegardingCode : invalid input code !");
2066   int nbOfSplit=(int)idsPerType.size();
2067   int nbOfTypes=(int)code.size()/3;
2068   int ret=0;
2069   for(int i=0;i<nbOfTypes;i++)
2070     {
2071       int nbOfEltInChunk=code[3*i+1];
2072       if(nbOfEltInChunk<0)
2073         throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::getNumberOfTuplesExpectedRegardingCode : invalid input code ! presence of negative value in a type !");
2074       int pos=code[3*i+2];
2075       if(pos!=-1)
2076         {
2077           if(pos<0 || pos>=nbOfSplit)
2078             {
2079               std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGaussNE::getNumberOfTuplesExpectedRegardingCode : input code points to pos " << pos << " in typeid " << i << " ! Should be in [0," << nbOfSplit << ") !";
2080               throw INTERP_KERNEL::Exception(oss.str().c_str());
2081             }
2082           const DataArrayInt *ids(idsPerType[pos]);
2083           if(!ids || !ids->isAllocated() || ids->getNumberOfComponents()!=1 || ids->getNumberOfTuples()!=nbOfEltInChunk || ids->getMinValueInArray()<0)
2084             {
2085               std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGaussNE::getNumberOfTuplesExpectedRegardingCode : input pfl chunck at pos " << pos << " should have " << i << " tuples and one component and with ids all >=0 !";
2086               throw INTERP_KERNEL::Exception(oss.str().c_str());
2087             }
2088         }
2089       ret+=nbOfEltInChunk;
2090     }
2091   if(!mesh)
2092     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::getNumberOfTuplesExpectedRegardingCode : NULL input mesh !");
2093   if(ret!=mesh->getNumberOfCells())
2094     {
2095       std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGaussNE::getNumberOfTuplesExpectedRegardingCode : input code points to " << ret << " number of cells should be " <<  mesh->getNumberOfCells() << " !";
2096     }
2097   return getNumberOfTuples(mesh);
2098 }
2099
2100 int MEDCouplingFieldDiscretizationGaussNE::getNumberOfTuples(const MEDCouplingMesh *mesh) const throw(INTERP_KERNEL::Exception)
2101 {
2102   if(!mesh)
2103     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::getNumberOfTuples : NULL input mesh !");
2104   int ret=0;
2105   int nbOfCells=mesh->getNumberOfCells();
2106   for(int i=0;i<nbOfCells;i++)
2107     {
2108       INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell(i);
2109       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
2110       if(cm.isDynamic())
2111         throw INTERP_KERNEL::Exception("Not implemented yet Gauss node on elements for polygons and polyedrons !");
2112       ret+=cm.getNumberOfNodes();
2113     }
2114   return ret;
2115 }
2116
2117 int MEDCouplingFieldDiscretizationGaussNE::getNumberOfMeshPlaces(const MEDCouplingMesh *mesh) const
2118 {
2119   if(!mesh)
2120     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::getNumberOfMeshPlaces : NULL input mesh !");
2121   return mesh->getNumberOfCells();
2122 }
2123
2124 DataArrayInt *MEDCouplingFieldDiscretizationGaussNE::getOffsetArr(const MEDCouplingMesh *mesh) const
2125 {
2126   if(!mesh)
2127     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::getOffsetArr : NULL input mesh !");
2128   int nbOfTuples=mesh->getNumberOfCells();
2129   DataArrayInt *ret=DataArrayInt::New();
2130   ret->alloc(nbOfTuples+1,1);
2131   int *retPtr=ret->getPointer();
2132   retPtr[0]=0;
2133   for(int i=0;i<nbOfTuples;i++)
2134     {
2135       INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell(i);
2136       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
2137       if(cm.isDynamic())
2138         throw INTERP_KERNEL::Exception("Not implemented yet Gauss node on elements for polygons and polyedrons !");
2139       retPtr[i+1]=retPtr[i]+cm.getNumberOfNodes();
2140     }
2141   return ret;
2142 }
2143
2144 void MEDCouplingFieldDiscretizationGaussNE::renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector<DataArray *>& arrays,
2145                                                                   const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
2146 {
2147   if(!mesh)
2148     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::renumberArraysForCell : NULL input mesh !");
2149   const int *array=old2NewBg;
2150   if(check)
2151     array=DataArrayInt::CheckAndPreparePermutation(old2NewBg,old2NewBg+mesh->getNumberOfCells());
2152   int nbOfCells=mesh->getNumberOfCells();
2153   int nbOfTuples=getNumberOfTuples(mesh);
2154   int *array2=new int[nbOfTuples];//stores the final conversion array old2New to give to arrays in renumberInPlace.
2155   int *array3=new int[nbOfCells];//store for each cell in after renumbering the offset needed by each cell in new numbering.
2156   array3[0]=0;
2157   for(int i=1;i<nbOfCells;i++)
2158     {
2159       INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell((int)std::distance(array,std::find(array,array+nbOfCells,i-1)));
2160       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
2161       array3[i]=array3[i-1]+cm.getNumberOfNodes();
2162     }
2163   int j=0;
2164   for(int i=0;i<nbOfCells;i++)
2165     {
2166       INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell(i);
2167       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
2168       for(int k=0;k<(int)cm.getNumberOfNodes();k++,j++)
2169         array2[j]=array3[array[i]]+k;
2170     }
2171   delete [] array3;
2172   for(std::vector<DataArray *>::const_iterator it=arrays.begin();it!=arrays.end();it++)
2173     if(*it)
2174       (*it)->renumberInPlace(array2);
2175   delete [] array2;
2176   if(check)
2177     free(const_cast<int *>(array));
2178 }
2179
2180 DataArrayDouble *MEDCouplingFieldDiscretizationGaussNE::getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const
2181 {
2182   if(!mesh)
2183     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::getLocalizationOfDiscValues : NULL input mesh !");
2184   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret=DataArrayDouble::New();
2185   MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> umesh=mesh->buildUnstructured();//in general do nothing
2186   int nbOfTuples=getNumberOfTuples(umesh);
2187   int spaceDim=mesh->getSpaceDimension();
2188   ret->alloc(nbOfTuples,spaceDim);
2189   const double *coords=umesh->getCoords()->begin();
2190   const int *connI=umesh->getNodalConnectivityIndex()->getConstPointer();
2191   const int *conn=umesh->getNodalConnectivity()->getConstPointer();
2192   int nbCells=umesh->getNumberOfCells();
2193   double *retPtr=ret->getPointer();
2194   for(int i=0;i<nbCells;i++,connI++)
2195     for(const int *w=conn+connI[0]+1;w!=conn+connI[1];w++)
2196       if(*w>=0)
2197         retPtr=std::copy(coords+(*w)*spaceDim,coords+((*w)+1)*spaceDim,retPtr);
2198   return ret.retn();
2199 }
2200
2201 /*!
2202  * Reimplemented from MEDCouplingFieldDiscretization::integral for performance reason. The default implementation is valid too for GAUSS_NE spatial discretization.
2203  */
2204 void MEDCouplingFieldDiscretizationGaussNE::integral(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, bool isWAbs, double *res) const throw(INTERP_KERNEL::Exception)
2205 {
2206   if(!mesh || !arr)
2207     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::integral : input mesh or array is null !");
2208   int nbOfCompo=arr->getNumberOfComponents();
2209   std::fill(res,res+nbOfCompo,0.);
2210   //
2211   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> vol=mesh->getMeasureField(isWAbs);
2212   std::set<INTERP_KERNEL::NormalizedCellType> types=mesh->getAllGeoTypes();
2213   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nbOfNodesPerCell=mesh->computeNbOfNodesPerCell();
2214   nbOfNodesPerCell->computeOffsets2();
2215   const double *arrPtr=arr->begin(),*volPtr=vol->getArray()->begin();
2216   for(std::set<INTERP_KERNEL::NormalizedCellType>::const_iterator it=types.begin();it!=types.end();it++)
2217     {
2218       std::size_t wArrSz=-1;
2219       const double *wArr=GetWeightArrayFromGeometricType(*it,wArrSz);
2220       INTERP_KERNEL::AutoPtr<double> wArr2=new double[wArrSz];
2221       double sum=std::accumulate(wArr,wArr+wArrSz,0.);
2222       std::transform(wArr,wArr+wArrSz,(double *)wArr2,std::bind2nd(std::multiplies<double>(),1./sum));
2223       MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ids=mesh->giveCellsWithType(*it);
2224       MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ids2=ids->buildExplicitArrByRanges(nbOfNodesPerCell);
2225       const int *ptIds2=ids2->begin(),*ptIds=ids->begin();
2226       int nbOfCellsWithCurGeoType=ids->getNumberOfTuples();
2227       for(int i=0;i<nbOfCellsWithCurGeoType;i++,ptIds++,ptIds2+=wArrSz)
2228         {
2229           for(int k=0;k<nbOfCompo;k++)
2230             {
2231               double tmp=0.;
2232               for(std::size_t j=0;j<wArrSz;j++)
2233                 tmp+=arrPtr[nbOfCompo*ptIds2[j]+k]*wArr2[j];
2234               res[k]+=tmp*volPtr[*ptIds];
2235             }
2236         }
2237     }
2238 }
2239
2240 const double *MEDCouplingFieldDiscretizationGaussNE::GetWeightArrayFromGeometricType(INTERP_KERNEL::NormalizedCellType geoType, std::size_t& lgth) throw(INTERP_KERNEL::Exception)
2241 {
2242   switch(geoType)
2243     {
2244     case INTERP_KERNEL::NORM_SEG2:
2245       lgth=(int)sizeof(FGP_SEG2)/sizeof(double);
2246       return FGP_SEG2;
2247     case INTERP_KERNEL::NORM_SEG3:
2248       lgth=(int)sizeof(FGP_SEG3)/sizeof(double);
2249       return FGP_SEG3;
2250     case INTERP_KERNEL::NORM_SEG4:
2251       lgth=(int)sizeof(FGP_SEG4)/sizeof(double);
2252       return FGP_SEG4;
2253     case INTERP_KERNEL::NORM_TRI3:
2254       lgth=(int)sizeof(FGP_TRI3)/sizeof(double);
2255       return FGP_TRI3;
2256     case INTERP_KERNEL::NORM_TRI6:
2257       lgth=(int)sizeof(FGP_TRI6)/sizeof(double);
2258       return FGP_TRI6;
2259     case INTERP_KERNEL::NORM_TRI7:
2260       lgth=(int)sizeof(FGP_TRI7)/sizeof(double);
2261       return FGP_TRI7;
2262     case INTERP_KERNEL::NORM_QUAD4:
2263       lgth=(int)sizeof(FGP_QUAD4)/sizeof(double);
2264       return FGP_QUAD4;
2265     case INTERP_KERNEL::NORM_QUAD9:
2266       lgth=(int)sizeof(FGP_QUAD9)/sizeof(double);
2267       return FGP_QUAD9;
2268     case INTERP_KERNEL::NORM_TETRA4:
2269       lgth=(int)sizeof(FGP_TETRA4)/sizeof(double);
2270       return FGP_TETRA4;
2271     case INTERP_KERNEL::NORM_PENTA6:
2272       lgth=(int)sizeof(FGP_PENTA6)/sizeof(double);
2273       return FGP_PENTA6;
2274     case INTERP_KERNEL::NORM_HEXA8:
2275       lgth=(int)sizeof(FGP_HEXA8)/sizeof(double);
2276       return FGP_HEXA8;
2277     case INTERP_KERNEL::NORM_HEXA27:
2278       lgth=(int)sizeof(FGP_HEXA27)/sizeof(double);
2279       return FGP_HEXA27;
2280     case INTERP_KERNEL::NORM_PYRA5:
2281       lgth=(int)sizeof(FGP_PYRA5)/sizeof(double);
2282       return FGP_PYRA5;
2283     default:
2284       throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::GetWeightArrayFromGeometricType : only SEG[2,3,4], TRI[3,6,7], QUAD[4,9], TETRA4, PENTA6, HEXA[8,27], PYRA5 supported !");
2285     }
2286 }
2287
2288 const double *MEDCouplingFieldDiscretizationGaussNE::GetRefCoordsFromGeometricType(INTERP_KERNEL::NormalizedCellType geoType, std::size_t& lgth) throw(INTERP_KERNEL::Exception)
2289 {
2290   switch(geoType)
2291     {
2292     case INTERP_KERNEL::NORM_SEG2:
2293       lgth=(int)sizeof(REF_SEG2)/sizeof(double);
2294       return REF_SEG2;
2295     case INTERP_KERNEL::NORM_SEG3:
2296       lgth=(int)sizeof(REF_SEG3)/sizeof(double);
2297       return REF_SEG3;
2298     case INTERP_KERNEL::NORM_SEG4:
2299       lgth=(int)sizeof(REF_SEG4)/sizeof(double);
2300       return REF_SEG4;
2301     case INTERP_KERNEL::NORM_TRI3:
2302       lgth=(int)sizeof(REF_TRI3)/sizeof(double);
2303       return REF_TRI3;
2304     case INTERP_KERNEL::NORM_TRI6:
2305       lgth=(int)sizeof(REF_TRI6)/sizeof(double);
2306       return REF_TRI6;
2307     case INTERP_KERNEL::NORM_TRI7:
2308       lgth=(int)sizeof(REF_TRI7)/sizeof(double);
2309       return REF_TRI7;
2310     case INTERP_KERNEL::NORM_QUAD4:
2311       lgth=(int)sizeof(REF_QUAD4)/sizeof(double);
2312       return REF_QUAD4;
2313     case INTERP_KERNEL::NORM_QUAD8:
2314       lgth=(int)sizeof(REF_QUAD8)/sizeof(double);
2315       return REF_QUAD8;
2316     case INTERP_KERNEL::NORM_QUAD9:
2317       lgth=(int)sizeof(REF_QUAD9)/sizeof(double);
2318       return REF_QUAD9;
2319     case INTERP_KERNEL::NORM_TETRA4:
2320       lgth=(int)sizeof(REF_TETRA4)/sizeof(double);
2321       return REF_TETRA4;
2322     case INTERP_KERNEL::NORM_TETRA10:
2323       lgth=(int)sizeof(REF_TETRA10)/sizeof(double);
2324       return REF_TETRA10;
2325     case INTERP_KERNEL::NORM_PENTA6:
2326       lgth=(int)sizeof(REF_PENTA6)/sizeof(double);
2327       return REF_PENTA6;
2328     case INTERP_KERNEL::NORM_PENTA15:
2329       lgth=(int)sizeof(REF_PENTA15)/sizeof(double);
2330       return REF_PENTA15;
2331     case INTERP_KERNEL::NORM_HEXA8:
2332       lgth=(int)sizeof(REF_HEXA8)/sizeof(double);
2333       return REF_HEXA8;
2334     case INTERP_KERNEL::NORM_HEXA20:
2335       lgth=(int)sizeof(REF_HEXA20)/sizeof(double);
2336       return REF_HEXA20;
2337     case INTERP_KERNEL::NORM_HEXA27:
2338       lgth=(int)sizeof(REF_HEXA27)/sizeof(double);
2339       return REF_HEXA27;
2340     case INTERP_KERNEL::NORM_PYRA5:
2341       lgth=(int)sizeof(REF_PYRA5)/sizeof(double);
2342       return REF_PYRA5;
2343     case INTERP_KERNEL::NORM_PYRA13:
2344       lgth=(int)sizeof(REF_PYRA13)/sizeof(double);
2345       return REF_PYRA13;
2346     default:
2347       throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::GetRefCoordsFromGeometricType : only SEG[2,3,4], TRI[3,6,7], QUAD[4,8,9], TETRA[4,10], PENTA[6,15], HEXA[8,20,27], PYRA[5,13] supported !");
2348     }
2349 }
2350
2351 void MEDCouplingFieldDiscretizationGaussNE::computeMeshRestrictionFromTupleIds(const MEDCouplingMesh *mesh, const int *tupleIdsBg, const int *tupleIdsEnd,
2352                                                                                DataArrayInt *&cellRestriction, DataArrayInt *&trueTupleRestriction) const throw(INTERP_KERNEL::Exception)
2353 {
2354   if(!mesh)
2355     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::computeMeshRestrictionFromTupleIds : NULL input mesh !");
2356   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp=DataArrayInt::New(); tmp->alloc((int)std::distance(tupleIdsBg,tupleIdsEnd),1);
2357   std::copy(tupleIdsBg,tupleIdsEnd,tmp->getPointer());
2358   tmp->sort(true);
2359   tmp=tmp->buildUnique();
2360   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nbOfNodesPerCell=mesh->computeNbOfNodesPerCell();
2361   nbOfNodesPerCell->computeOffsets2();
2362   nbOfNodesPerCell->searchRangesInListOfIds(tmp,cellRestriction,trueTupleRestriction);
2363 }
2364
2365 void MEDCouplingFieldDiscretizationGaussNE::checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception)
2366 {
2367 }
2368
2369 double MEDCouplingFieldDiscretizationGaussNE::getIJK(const MEDCouplingMesh *mesh, const DataArrayDouble *da,
2370                                                      int cellId, int nodeIdInCell, int compoId) const throw(INTERP_KERNEL::Exception)
2371 {
2372   if(!mesh)
2373     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::getIJK : NULL input mesh !");
2374   int offset=0;
2375   for(int i=0;i<cellId;i++)
2376     {
2377       INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell(i);
2378       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
2379       offset+=cm.getNumberOfNodes();
2380     }
2381   return da->getIJ(offset+nodeIdInCell,compoId);
2382 }
2383
2384 void MEDCouplingFieldDiscretizationGaussNE::checkCoherencyBetween(const MEDCouplingMesh *mesh, const DataArray *da) const throw(INTERP_KERNEL::Exception)
2385 {
2386   int nbOfTuples=getNumberOfTuples(mesh);
2387   if(nbOfTuples!=da->getNumberOfTuples())
2388     {
2389       std::ostringstream oss; oss << "Invalid number of tuples in the array : expecting " << nbOfTuples << " !";
2390       throw INTERP_KERNEL::Exception(oss.str().c_str());
2391     }
2392 }
2393
2394 MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationGaussNE::getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const
2395 {
2396   if(!mesh)
2397     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::getMeasureField : mesh instance specified is NULL !");
2398   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> vol=mesh->getMeasureField(isAbs);
2399   const double *volPtr=vol->getArray()->begin();
2400   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(ON_GAUSS_NE);
2401   ret->setMesh(mesh);
2402   //
2403   std::set<INTERP_KERNEL::NormalizedCellType> types=mesh->getAllGeoTypes();
2404   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nbOfNodesPerCell=mesh->computeNbOfNodesPerCell();
2405   int nbTuples=nbOfNodesPerCell->accumulate(0);
2406   nbOfNodesPerCell->computeOffsets2();
2407   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr=DataArrayDouble::New(); arr->alloc(nbTuples,1);
2408   ret->setArray(arr);
2409   double *arrPtr=arr->getPointer();
2410   for(std::set<INTERP_KERNEL::NormalizedCellType>::const_iterator it=types.begin();it!=types.end();it++)
2411     {
2412       std::size_t wArrSz=-1;
2413       const double *wArr=GetWeightArrayFromGeometricType(*it,wArrSz);
2414       INTERP_KERNEL::AutoPtr<double> wArr2=new double[wArrSz];
2415       double sum=std::accumulate(wArr,wArr+wArrSz,0.);
2416       std::transform(wArr,wArr+wArrSz,(double *)wArr2,std::bind2nd(std::multiplies<double>(),1./sum));
2417       MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ids=mesh->giveCellsWithType(*it);
2418       MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ids2=ids->buildExplicitArrByRanges(nbOfNodesPerCell);
2419       const int *ptIds2=ids2->begin(),*ptIds=ids->begin();
2420       int nbOfCellsWithCurGeoType=ids->getNumberOfTuples();
2421       for(int i=0;i<nbOfCellsWithCurGeoType;i++,ptIds++)
2422         for(std::size_t j=0;j<wArrSz;j++,ptIds2++)
2423           arrPtr[*ptIds2]=wArr2[j]*volPtr[*ptIds];
2424     }
2425   ret->synchronizeTimeWithSupport();
2426   return ret.retn();
2427 }
2428
2429 void MEDCouplingFieldDiscretizationGaussNE::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
2430 {
2431   throw INTERP_KERNEL::Exception("Not implemented yet !");
2432 }
2433
2434 void MEDCouplingFieldDiscretizationGaussNE::getValueOnPos(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, int i, int j, int k, double *res) const
2435 {
2436   throw INTERP_KERNEL::Exception("getValueOnPos(i,j,k) : Not applyable for Gauss points !");
2437 }
2438
2439 DataArrayDouble *MEDCouplingFieldDiscretizationGaussNE::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, int nbOfPoints) const
2440 {
2441   throw INTERP_KERNEL::Exception("getValueOnMulti : Not implemented for Gauss NE !");
2442 }
2443
2444 MEDCouplingMesh *MEDCouplingFieldDiscretizationGaussNE::buildSubMeshData(const MEDCouplingMesh *mesh, const int *start, const int *end, DataArrayInt *&di) const
2445 {
2446   if(!mesh)
2447     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::buildSubMeshData : NULL input mesh !");
2448   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> diSafe=computeTupleIdsToSelectFromCellIds(mesh,start,end);
2449   MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> ret=mesh->buildPart(start,end);
2450   di=diSafe.retn();
2451   return ret.retn();
2452 }
2453
2454 /*!
2455  * This method is strictly equivalent to MEDCouplingFieldDiscretizationGauss::buildSubMeshData except that it is optimized for input defined as a range of cell ids.
2456  *
2457  * \param [out] beginOut Valid only if \a di is NULL
2458  * \param [out] endOut Valid only if \a di is NULL
2459  * \param [out] stepOut Valid only if \a di is NULL
2460  * \param [out] di is an array returned that specifies entity ids (nodes, cells, Gauss points... ) in array if no output range is foundable.
2461  *
2462  * \sa MEDCouplingFieldDiscretizationGauss::buildSubMeshData
2463  */
2464 MEDCouplingMesh *MEDCouplingFieldDiscretizationGaussNE::buildSubMeshDataRange(const MEDCouplingMesh *mesh, int beginCellIds, int endCellIds, int stepCellIds, int& beginOut, int& endOut, int& stepOut, DataArrayInt *&di) const
2465 {
2466   if(stepCellIds!=1)//even for stepCellIds==-1 the output will not be a range
2467     return MEDCouplingFieldDiscretization::buildSubMeshDataRange(mesh,beginCellIds,endCellIds,stepCellIds,beginOut,endOut,stepOut,di);
2468   if(!mesh)
2469     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::buildSubMeshDataRange : NULL input mesh !");
2470   int nbOfCells=mesh->getNumberOfCells();
2471   di=0; beginOut=0; endOut=0; stepOut=stepCellIds;
2472   const char msg[]="MEDCouplingFieldDiscretizationGaussNE::buildSubMeshDataRange : cell #";
2473   for(int i=0;i<nbOfCells;i++)
2474     {
2475       INTERP_KERNEL::NormalizedCellType type=mesh->getTypeOfCell(i);
2476       const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
2477       if(cm.isDynamic())
2478         { std::ostringstream oss; oss << msg << i << " presence of dynamic cell (polygons and polyedrons) ! Not implemented !"; throw INTERP_KERNEL::Exception(oss.str().c_str()); }
2479       int delta=cm.getNumberOfNodes();
2480       if(i<beginCellIds)
2481         beginOut+=delta;
2482       endOut+=delta;
2483       if(i>=endCellIds)
2484         break;
2485     }
2486   MEDCouplingAutoRefCountObjectPtr<MEDCouplingMesh> ret=mesh->buildPartRange(beginCellIds,endCellIds,stepCellIds);
2487   return ret.retn();
2488 }
2489
2490
2491 /*!
2492  * This method returns a tuple ids selection from cell ids selection [start;end).
2493  * This method is called by MEDCouplingFieldDiscretizationGaussNE::buildSubMeshData to return parameter \b di.
2494  *
2495  * \return a newly allocated array containing ids to select into the DataArrayDouble of the field.
2496  *
2497  */
2498 DataArrayInt *MEDCouplingFieldDiscretizationGaussNE::computeTupleIdsToSelectFromCellIds(const MEDCouplingMesh *mesh, const int *startCellIds, const int *endCellIds) const
2499 {
2500   if(!mesh)
2501     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::computeTupleIdsToSelectFromCellIds : null mesh !");
2502   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nbOfNodesPerCell=mesh->computeNbOfNodesPerCell();
2503   nbOfNodesPerCell->computeOffsets2();
2504   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> sel=DataArrayInt::New(); sel->useArray(startCellIds,false,CPP_DEALLOC,(int)std::distance(startCellIds,endCellIds),1);
2505   return sel->buildExplicitArrByRanges(nbOfNodesPerCell);
2506 }
2507
2508 /*!
2509  * No implementation needed !
2510  */
2511 void MEDCouplingFieldDiscretizationGaussNE::renumberValuesOnNodes(double , const int *, int newNbOfNodes, DataArrayDouble *) const
2512 {
2513 }
2514
2515 void MEDCouplingFieldDiscretizationGaussNE::renumberValuesOnCells(double epsOnVals, const MEDCouplingMesh *mesh, const int *old2New, int newSz, DataArrayDouble *arr) const
2516 {
2517   throw INTERP_KERNEL::Exception("Not implemented yet !");
2518 }
2519
2520 void MEDCouplingFieldDiscretizationGaussNE::renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const int *new2old, int newSz, DataArrayDouble *arr) const
2521 {
2522   throw INTERP_KERNEL::Exception("Not implemented yet !");
2523 }
2524
2525 void MEDCouplingFieldDiscretizationGaussNE::reprQuickOverview(std::ostream& stream) const throw(INTERP_KERNEL::Exception)
2526 {
2527   stream << "Gauss points on nodes per element spatial discretization.";
2528 }
2529
2530 MEDCouplingFieldDiscretizationGaussNE::MEDCouplingFieldDiscretizationGaussNE(const MEDCouplingFieldDiscretizationGaussNE& other):MEDCouplingFieldDiscretization(other)
2531 {
2532 }
2533
2534 TypeOfField MEDCouplingFieldDiscretizationKriging::getEnum() const
2535 {
2536   return TYPE;
2537 }
2538
2539 const char *MEDCouplingFieldDiscretizationKriging::getRepr() const
2540 {
2541   return REPR;
2542 }
2543
2544 /*!
2545  * This method is simply called by MEDCouplingFieldDiscretization::deepCpy. It performs the deep copy of \a this.
2546  *
2547  * \sa MEDCouplingFieldDiscretization::deepCpy.
2548  */
2549 MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretizationKriging::clone() const
2550 {
2551   return new MEDCouplingFieldDiscretizationKriging;
2552 }
2553
2554 std::string MEDCouplingFieldDiscretizationKriging::getStringRepr() const
2555 {
2556   return std::string(REPR);
2557 }
2558
2559 void MEDCouplingFieldDiscretizationKriging::checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception)
2560 {
2561   if(nat!=ConservativeVolumic)
2562     throw INTERP_KERNEL::Exception("Invalid nature for Kriging field : expected ConservativeVolumic !");
2563 }
2564
2565 bool MEDCouplingFieldDiscretizationKriging::isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const
2566 {
2567   if(!other)
2568     {
2569       reason="other spatial discretization is NULL, and this spatial discretization (Kriginig) is defined.";
2570       return false;
2571     }
2572   const MEDCouplingFieldDiscretizationKriging *otherC=dynamic_cast<const MEDCouplingFieldDiscretizationKriging *>(other);
2573   bool ret=otherC!=0;
2574   if(!ret)
2575     reason="Spatial discrtization of this is ON_NODES_KR, which is not the case of other.";
2576   return ret;
2577 }
2578
2579 MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationKriging::getMeasureField(const MEDCouplingMesh *mesh, bool isAbs) const
2580 {
2581   if(!mesh)
2582     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::getMeasureField : mesh instance specified is NULL !");
2583   throw INTERP_KERNEL::Exception("getMeasureField on FieldDiscretizationKriging : not implemented yet !");
2584 }
2585
2586 void MEDCouplingFieldDiscretizationKriging::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
2587 {
2588   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> res2=MEDCouplingFieldDiscretizationKriging::getValueOnMulti(arr,mesh,loc,1);
2589   std::copy(res2->begin(),res2->end(),res);
2590 }
2591
2592 DataArrayDouble *MEDCouplingFieldDiscretizationKriging::getValueOnMulti(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, int nbOfTargetPoints) const
2593 {
2594   if(!arr || !arr->isAllocated())
2595     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::getValueOnMulti : input array is null or not allocated !");
2596   int nbOfRows(getNumberOfMeshPlaces(mesh));
2597   if(arr->getNumberOfTuples()!=nbOfRows)
2598     {
2599       std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationKriging::getValueOnMulti : input array does not have correct number of tuples ! Excepted " << nbOfRows << " having " << arr->getNumberOfTuples() << " !";
2600       throw INTERP_KERNEL::Exception(oss.str().c_str());
2601     }
2602   int nbCols(-1),nbCompo(arr->getNumberOfComponents());
2603   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> m(computeEvaluationMatrixOnGivenPts(mesh,loc,nbOfTargetPoints,nbCols));
2604   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New());
2605   ret->alloc(nbOfTargetPoints,nbCompo);
2606   INTERP_KERNEL::matrixProduct(m->begin(),nbOfTargetPoints,nbCols,arr->begin(),nbOfRows,nbCompo,ret->getPointer());
2607   return ret.retn();
2608 }
2609
2610 void MEDCouplingFieldDiscretizationKriging::reprQuickOverview(std::ostream& stream) const throw(INTERP_KERNEL::Exception)
2611 {
2612   stream << "Kriging spatial discretization.";
2613 }
2614
2615 /*!
2616  * Returns the matrix of size nbRows = \a nbOfTargetPoints and \a nbCols = \a nbCols. This matrix is useful if
2617  *
2618  * \return the new result matrix to be deallocated by the caller.
2619  */
2620 DataArrayDouble *MEDCouplingFieldDiscretizationKriging::computeEvaluationMatrixOnGivenPts(const MEDCouplingMesh *mesh, const double *loc, int nbOfTargetPoints, int& nbCols) const
2621 {
2622   int isDrift(-1),nbRows(-1);
2623   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixInv(computeInverseMatrix(mesh,isDrift,nbRows));
2624   //
2625   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coords=getLocalizationOfDiscValues(mesh);
2626   int nbOfPts(coords->getNumberOfTuples()),dimension(coords->getNumberOfComponents());
2627   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> locArr=DataArrayDouble::New();
2628   locArr->useArray(loc,false,CPP_DEALLOC,nbOfTargetPoints,dimension);
2629   nbCols=nbOfPts;
2630   //
2631   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrix2=coords->buildEuclidianDistanceDenseMatrixWith(locArr);
2632   operateOnDenseMatrix(mesh->getSpaceDimension(),nbOfTargetPoints*nbOfPts,matrix2->getPointer());
2633   //
2634   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrix3=DataArrayDouble::New();
2635   matrix3->alloc(nbOfTargetPoints*nbRows,1);
2636   double *work=matrix3->getPointer();
2637   const double *workCst(matrix2->begin()),*workCst2(loc);
2638   for(int i=0;i<nbOfTargetPoints;i++,workCst+=nbOfPts,workCst2+=isDrift-1)
2639     {
2640       for(int j=0;j<nbOfPts;j++)
2641         work[i*nbRows+j]=workCst[j];
2642       work[i*nbRows+nbOfPts]=1.0;
2643       for(int j=0;j<isDrift-1;j++)
2644         work[i*nbRows+(nbOfPts+1+j)]=workCst2[j];
2645     }
2646   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New());
2647   ret->alloc(nbOfTargetPoints,nbRows);
2648   INTERP_KERNEL::matrixProduct(matrix3->begin(),nbOfTargetPoints,nbRows,matrixInv->begin(),nbRows,nbRows,ret->getPointer());
2649   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret2(DataArrayDouble::New());
2650   ret2->alloc(nbOfTargetPoints*nbOfPts,1);
2651   workCst=ret->begin(); work=ret2->getPointer();
2652   for(int i=0;i<nbOfTargetPoints;i++,workCst+=nbRows)
2653     work=std::copy(workCst,workCst+nbOfPts,work);
2654   return ret2.retn();
2655 }
2656
2657 /*!
2658  * This method returns the square matrix of size \a matSz that is the inverse of the kriging matrix. The returned matrix can returned all the coeffs of kriging
2659  * when multiplied by the vector of values attached to each point.
2660  *
2661  * \param [out] isDrift return if drift coefficients are present in the returned vector of coefficients. If different from 0 there is presence of drift coefficients.
2662  * \param [out] matSz the size of returned square matrix
2663  * \return the new result matrix to be deallocated by the caller.
2664  */
2665 DataArrayDouble *MEDCouplingFieldDiscretizationKriging::computeInverseMatrix(const MEDCouplingMesh *mesh, int& isDrift, int& matSz) const
2666 {
2667   if(!mesh)
2668     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::computeVectorOfCoefficients : NULL input mesh !");
2669   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coords=getLocalizationOfDiscValues(mesh);
2670   int nbOfPts=coords->getNumberOfTuples();
2671   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrix=coords->buildEuclidianDistanceDenseMatrix();
2672   operateOnDenseMatrix(mesh->getSpaceDimension(),nbOfPts*nbOfPts,matrix->getPointer());
2673   // Drift
2674   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixWithDrift=performDrift(matrix,coords,isDrift);
2675   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixInv=DataArrayDouble::New();
2676   matSz=nbOfPts+isDrift;
2677   matrixInv->alloc(matSz*matSz,1);
2678   INTERP_KERNEL::inverseMatrix(matrixWithDrift->getConstPointer(),matSz,matrixInv->getPointer());
2679   return matrixInv.retn();
2680 }
2681
2682 /*!
2683  * 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
2684  * number of tuples should be equal to the number of representing points in \a mesh.
2685  *
2686  * \param [in] mesh is the sources of nodes on which kriging will be done regarding the parameters and the value of \c this->getSpaceDimension()
2687  * \param [in] arr input field DataArrayDouble whose number of tuples must be equal to the number of nodes in \a mesh
2688  * \param [out] isDrift return if drift coefficients are present in the returned vector of coefficients. If different from 0 there is presence of drift coefficients.
2689  *              Whatever the value of \a isDrift the number of tuples of returned DataArrayDouble  will be equal to \c arr->getNumberOfTuples() + \a isDrift.
2690  * \return a newly allocated array containing coefficients including or not drift coefficient at the end depending the value of \a isDrift parameter.
2691  */
2692 DataArrayDouble *MEDCouplingFieldDiscretizationKriging::computeVectorOfCoefficients(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, int& isDrift) const
2693 {
2694   int nbRows(-1);
2695   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> matrixInv(computeInverseMatrix(mesh,isDrift,nbRows));
2696   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> KnewiK=DataArrayDouble::New();
2697   KnewiK->alloc(nbRows*1,1);
2698   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr2=DataArrayDouble::New();
2699   arr2->alloc(nbRows*1,1);
2700   double *work=std::copy(arr->begin(),arr->end(),arr2->getPointer());
2701   std::fill(work,work+isDrift,0.);
2702   INTERP_KERNEL::matrixProduct(matrixInv->getConstPointer(),nbRows,nbRows,arr2->getConstPointer(),nbRows,1,KnewiK->getPointer());
2703   return KnewiK.retn();
2704 }
2705
2706 /*!
2707  * 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.
2708  *
2709  * \param [in] spaceDimension space dimension of the input mesh on which the Kriging has to be performed
2710  * \param [in] nbOfElems is the result of the product of nb of rows and the nb of columns of matrix \a matrixPtr
2711  * \param [in,out] matrixPtr is the dense matrix whose on each values the operation will be applied
2712  */
2713 void MEDCouplingFieldDiscretizationKriging::operateOnDenseMatrix(int spaceDimension, int nbOfElems, double *matrixPtr) const
2714 {
2715   switch(spaceDimension)
2716     {
2717     case 1:
2718       {
2719         for(int i=0;i<nbOfElems;i++)
2720           {
2721             double val=matrixPtr[i];
2722             matrixPtr[i]=val*val*val;
2723           }
2724         break;
2725       }
2726     case 2:
2727       {
2728         for(int i=0;i<nbOfElems;i++)
2729           {
2730             double val=matrixPtr[i];
2731             if(val!=0.)
2732               matrixPtr[i]=val*val*log(val);
2733           }
2734         break;
2735       }
2736     case 3:
2737       {
2738         //nothing here : it is not a bug g(h)=h with spaceDim 3.
2739         break;
2740       }
2741     default:
2742       throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::operateOnDenseMatrix : only dimension 1, 2 and 3 implemented !");
2743     }
2744 }
2745
2746 /*!
2747  * Starting from a square matrix \a matr, this method returns a newly allocated dense square matrix whose \a matr is included in returned matrix
2748  * in the top left corner, and in the remaining returned matrix the parameters to take into account about the kriging drift.
2749  * For the moment only linear srift is implemented.
2750  *
2751  * \param [in] arr the position of points were input mesh geometry is considered for Kriging
2752  * \param [in] matr input matrix whose drift part will be added
2753  * \param [out] delta the difference between the size of the output matrix and the input matrix \a matr.
2754  * \return a newly allocated matrix bigger than input matrix \a matr.
2755  */
2756 DataArrayDouble *MEDCouplingFieldDiscretizationKriging::performDrift(const DataArrayDouble *matr, const DataArrayDouble *arr, int& delta) const
2757 {
2758   int spaceDimension=arr->getNumberOfComponents();
2759   delta=spaceDimension+1;
2760   int szOfMatrix=arr->getNumberOfTuples();
2761   if(szOfMatrix*szOfMatrix!=matr->getNumberOfTuples())
2762     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationKriging::performDrift : invalid size");
2763   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret=DataArrayDouble::New();
2764   ret->alloc((szOfMatrix+delta)*(szOfMatrix+delta),1);
2765   const double *srcWork=matr->getConstPointer();
2766   const double *srcWork2=arr->getConstPointer();
2767   double *destWork=ret->getPointer();
2768   for(int i=0;i<szOfMatrix;i++)
2769     {
2770       destWork=std::copy(srcWork,srcWork+szOfMatrix,destWork);
2771       srcWork+=szOfMatrix;
2772       *destWork++=1.;
2773       destWork=std::copy(srcWork2,srcWork2+spaceDimension,destWork);
2774       srcWork2+=spaceDimension;
2775     }
2776   std::fill(destWork,destWork+szOfMatrix,1.); destWork+=szOfMatrix;
2777   std::fill(destWork,destWork+spaceDimension+1,0.); destWork+=spaceDimension+1;
2778   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arrNoI=arr->toNoInterlace();
2779   srcWork2=arrNoI->getConstPointer();
2780   for(int i=0;i<spaceDimension;i++)
2781     {
2782       destWork=std::copy(srcWork2,srcWork2+szOfMatrix,destWork);
2783       srcWork2+=szOfMatrix;
2784       std::fill(destWork,destWork+spaceDimension+1,0.);
2785       destWork+=spaceDimension+1;
2786     }
2787   //
2788   return ret.retn();
2789 }