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