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