src/MEDCoupling/MEDCouplingFieldDiscretization.cxx

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