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