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