src/MEDCoupling/MEDCouplingFieldDiscretization.cxx

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