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