refactor!: remove adm_local/ directory

[tools/medcoupling.git] / src / MEDCoupling / MEDCouplingFieldDiscretization.cxx

diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
index 52c630549cc08aaddba0a25d3e086df3ad3c910c..1e50212834c1af5cc6f15e83b60dd0cb5b25c227 100755 (executable)
--- a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
+++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2020  CEA/DEN, EDF R&D
+// Copyright (C) 2007-2024  CEA, EDF
 //
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -19,6 +19,7 @@
 // Author : Anthony Geay (EDF R&D)
 
 #include "MEDCouplingFieldDiscretization.hxx"
+#include "MEDCouplingFieldDiscretizationOnNodesFE.hxx"
 #include "MEDCouplingCMesh.hxx"
 #include "MEDCouplingUMesh.hxx"
 #include "MEDCouplingFieldDouble.hxx"
@@ -29,6 +30,7 @@
 #include "InterpKernelAutoPtr.hxx"
 #include "InterpKernelGaussCoords.hxx"
 #include "InterpKernelMatrixTools.hxx"
+#include "InterpKernelDenseMatrix.hxx"
 
 #include <set>
 #include <list>
@@ -44,26 +46,16 @@ const double MEDCouplingFieldDiscretization::DFLT_PRECISION=1.e-12;
 
 const char MEDCouplingFieldDiscretizationP0::REPR[]="P0";
 
-const TypeOfField MEDCouplingFieldDiscretizationP0::TYPE=ON_CELLS;
-
 const char MEDCouplingFieldDiscretizationP1::REPR[]="P1";
 
-const TypeOfField MEDCouplingFieldDiscretizationP1::TYPE=ON_NODES;
-
 const mcIdType MEDCouplingFieldDiscretizationPerCell::DFT_INVALID_LOCID_VALUE=-1;
 
 const char MEDCouplingFieldDiscretizationGauss::REPR[]="GAUSS";
 
-const TypeOfField MEDCouplingFieldDiscretizationGauss::TYPE=ON_GAUSS_PT;
-
 const char MEDCouplingFieldDiscretizationGaussNE::REPR[]="GSSNE";
 
-const TypeOfField MEDCouplingFieldDiscretizationGaussNE::TYPE=ON_GAUSS_NE;
-
 const char MEDCouplingFieldDiscretizationKriging::REPR[]="KRIGING";
 
-const TypeOfField MEDCouplingFieldDiscretizationKriging::TYPE=ON_NODES_KR;
-
 // doc is here http://www.code-aster.org/V2/doc/default/fr/man_r/r3/r3.01.01.pdf
 const double MEDCouplingFieldDiscretizationGaussNE::FGP_POINT1[1]={0.};
 const double MEDCouplingFieldDiscretizationGaussNE::FGP_SEG2[2]={1.,1.};
@@ -79,6 +71,7 @@ const double MEDCouplingFieldDiscretizationGaussNE::FGP_TETRA4[4]={0.04166666666
 const double MEDCouplingFieldDiscretizationGaussNE::FGP_TETRA10[10]={1.,1.,1.,1.,1.,1.,1.,1.,1.,1.};//to check
 const double MEDCouplingFieldDiscretizationGaussNE::FGP_PENTA6[6]={0.16666666666666666,0.16666666666666666,0.16666666666666666,0.16666666666666666,0.16666666666666666,0.16666666666666666};
 const double MEDCouplingFieldDiscretizationGaussNE::FGP_PENTA15[15]={1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.};//to check
+const double MEDCouplingFieldDiscretizationGaussNE::FGP_PENTA18[18]={1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.};//to check
 const double MEDCouplingFieldDiscretizationGaussNE::FGP_HEXA8[8]={1.,1.,1.,1.,1.,1.,1.,1.};
 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.};
 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.};
@@ -97,6 +90,7 @@ const double MEDCouplingFieldDiscretizationGaussNE::REF_TETRA4[12]={0.,1.,0.,0.,
 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.};
 const double MEDCouplingFieldDiscretizationGaussNE::REF_PENTA6[18]={-1.,1.,0.,-1.,0.,1.,-1.,0.,0.,1.,1.,0.,1.,0.,1.,1.,0.,0.};
 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.};
+const double MEDCouplingFieldDiscretizationGaussNE::REF_PENTA18[54]={-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.,0.,0.5,0.5,0.,0.,0.5,0.,0.5,0.};
 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.};
 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.};
 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.};
@@ -115,6 +109,7 @@ const double MEDCouplingFieldDiscretizationGaussNE::LOC_TETRA4[12]={0.1381966011
 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
 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.};
 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
+const double MEDCouplingFieldDiscretizationGaussNE::LOC_PENTA18[54]={-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.,0.,0.5,0.5,0.,0.,0.5,0.,0.5,0.};//to check
 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};
 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
 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.};
@@ -139,6 +134,8 @@ MEDCouplingFieldDiscretization *MEDCouplingFieldDiscretization::New(TypeOfField
       return new MEDCouplingFieldDiscretizationGaussNE;
     case MEDCouplingFieldDiscretizationKriging::TYPE:
       return new MEDCouplingFieldDiscretizationKriging;
+    case MEDCouplingFieldDiscretizationOnNodesFE::TYPE:
+      return new MEDCouplingFieldDiscretizationOnNodesFE;
     default:
       throw INTERP_KERNEL::Exception("Chosen discretization is not implemented yet.");
   }
@@ -156,22 +153,30 @@ TypeOfField MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(const s
     return MEDCouplingFieldDiscretizationGaussNE::TYPE;
   if(repr==MEDCouplingFieldDiscretizationKriging::REPR)
     return MEDCouplingFieldDiscretizationKriging::TYPE;
+  if(repr==MEDCouplingFieldDiscretizationOnNodesFE::REPR)
+    return MEDCouplingFieldDiscretizationOnNodesFE::TYPE;
   throw INTERP_KERNEL::Exception("Representation does not match with any field discretization !");
 }
 
 std::string MEDCouplingFieldDiscretization::GetTypeOfFieldRepr(TypeOfField type)
 {
-  if(type==MEDCouplingFieldDiscretizationP0::TYPE)
-    return MEDCouplingFieldDiscretizationP0::REPR;
-  if(type==MEDCouplingFieldDiscretizationP1::TYPE)
-    return MEDCouplingFieldDiscretizationP1::REPR;
-  if(type==MEDCouplingFieldDiscretizationGauss::TYPE)
-    return MEDCouplingFieldDiscretizationGauss::REPR;
-  if(type==MEDCouplingFieldDiscretizationGaussNE::TYPE)
-    return MEDCouplingFieldDiscretizationGaussNE::REPR;
-  if(type==MEDCouplingFieldDiscretizationKriging::TYPE)
-    return MEDCouplingFieldDiscretizationKriging::REPR;
-  throw INTERP_KERNEL::Exception("GetTypeOfFieldRepr : Representation does not match with any field discretization !");
+  switch(type)
+  {
+    case MEDCouplingFieldDiscretizationP0::TYPE:
+      return MEDCouplingFieldDiscretizationP0::REPR;
+    case MEDCouplingFieldDiscretizationP1::TYPE:
+      return MEDCouplingFieldDiscretizationP1::REPR;
+    case MEDCouplingFieldDiscretizationGauss::TYPE:
+      return MEDCouplingFieldDiscretizationGauss::REPR;
+    case MEDCouplingFieldDiscretizationGaussNE::TYPE:
+      return MEDCouplingFieldDiscretizationGaussNE::REPR;
+    case MEDCouplingFieldDiscretizationKriging::TYPE:
+      return MEDCouplingFieldDiscretizationKriging::REPR;
+    case MEDCouplingFieldDiscretizationOnNodesFE::TYPE:
+      return MEDCouplingFieldDiscretizationOnNodesFE::REPR;
+    default:
+      throw INTERP_KERNEL::Exception("GetTypeOfFieldRepr : Representation does not match with any field discretization !");
+  }
 }
 
 bool MEDCouplingFieldDiscretization::isEqual(const MEDCouplingFieldDiscretization *other, double eps) const
@@ -248,7 +253,7 @@ void MEDCouplingFieldDiscretization::normL1(const MEDCouplingMesh *mesh, const D
         res[j]+=fabs(arrPtr[i*nbOfCompo+j])*v;
       deno+=v;
     }
-  std::transform(res,res+nbOfCompo,res,std::bind2nd(std::multiplies<double>(),1./deno));
+  std::transform(res,res+nbOfCompo,res,std::bind(std::multiplies<double>(),std::placeholders::_1,1./deno));
 }
 
 /*!
@@ -272,8 +277,8 @@ void MEDCouplingFieldDiscretization::normL2(const MEDCouplingMesh *mesh, const D
         res[j]+=arrPtr[i*nbOfCompo+j]*arrPtr[i*nbOfCompo+j]*v;
       deno+=v;
     }
-  std::transform(res,res+nbOfCompo,res,std::bind2nd(std::multiplies<double>(),1./deno));
-  std::transform(res,res+nbOfCompo,res,std::ptr_fun<double,double>(std::sqrt));
+  std::transform(res,res+nbOfCompo,res,std::bind(std::multiplies<double>(),std::placeholders::_1,1./deno));
+  std::transform(res,res+nbOfCompo,res,[](double c){return sqrt(c);});
 }
 
 /*!
@@ -301,7 +306,7 @@ void MEDCouplingFieldDiscretization::integral(const MEDCouplingMesh *mesh, const
   INTERP_KERNEL::AutoPtr<double> tmp=new double[nbOfCompo];
   for(mcIdType i=0;i<nbOfElems;i++)
     {
-      std::transform(arrPtr+i*nbOfCompo,arrPtr+(i+1)*nbOfCompo,(double *)tmp,std::bind2nd(std::multiplies<double>(),volPtr[i]));
+      std::transform(arrPtr+i*nbOfCompo,arrPtr+(i+1)*nbOfCompo,(double *)tmp,std::bind(std::multiplies<double>(),std::placeholders::_1,volPtr[i]));
       std::transform((double *)tmp,(double *)tmp+nbOfCompo,res,res,std::plus<double>());
     }
 }
@@ -443,13 +448,13 @@ void MEDCouplingFieldDiscretization::RenumberEntitiesFromO2NArr(double eps, cons
       mcIdType newNb=old2NewPtr[i];
       if(newNb>=0)//if newNb<0 the node is considered as out.
         {
-          if(std::find_if(ptToFill+newNb*nbOfComp,ptToFill+(newNb+1)*nbOfComp,std::bind2nd(std::not_equal_to<double>(),std::numeric_limits<double>::max()))
+          if(std::find_if(ptToFill+newNb*nbOfComp,ptToFill+(newNb+1)*nbOfComp,std::bind(std::not_equal_to<double>(),std::placeholders::_1,std::numeric_limits<double>::max()))
           ==ptToFill+(newNb+1)*nbOfComp)
             std::copy(ptSrc+i*nbOfComp,ptSrc+(i+1)*nbOfComp,ptToFill+newNb*nbOfComp);
           else
             {
               std::transform(ptSrc+i*nbOfComp,ptSrc+(i+1)*nbOfComp,ptToFill+newNb*nbOfComp,(double *)tmp,std::minus<double>());
-              std::transform((double *)tmp,((double *)tmp)+nbOfComp,(double *)tmp,std::ptr_fun<double,double>(fabs));
+              std::transform((double *)tmp,((double *)tmp)+nbOfComp,(double *)tmp,[](double c){return fabs(c);});
               //if(!std::equal(ptSrc+i*nbOfComp,ptSrc+(i+1)*nbOfComp,ptToFill+newNb*nbOfComp))
               if(*std::max_element((double *)tmp,((double *)tmp)+nbOfComp)>eps)
                 {
@@ -477,20 +482,6 @@ void MEDCouplingFieldDiscretization::RenumberEntitiesFromN2OArr(const mcIdType *
     }
 }
 
-template<class FIELD_DISC>
-MCAuto<MEDCouplingFieldDiscretization> MEDCouplingFieldDiscretization::EasyAggregate(std::vector<const MEDCouplingFieldDiscretization *>& fds)
-{
-  if(fds.empty())
-    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretization::aggregate : input array is empty");
-  for(const MEDCouplingFieldDiscretization * it : fds)
-    {
-      const FIELD_DISC *itc(dynamic_cast<const FIELD_DISC *>(it));
-      if(!itc)
-        throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretization::aggregate : same field discretization expected for all input discretizations !");
-    }
-  return fds[0]->clone();
-}
-
 MEDCouplingFieldDiscretization::~MEDCouplingFieldDiscretization()
 {
 }
@@ -1068,7 +1059,7 @@ void MEDCouplingFieldDiscretizationP1::getValueInCell(const MEDCouplingMesh *mes
   for(std::size_t i=0;i<nbOfNodes;i++)
     {
       arr->getTuple(conn[i],(double *)tmp2);
-      std::transform((double *)tmp2,((double *)tmp2)+sz,(double *)tmp2,std::bind2nd(std::multiplies<double>(),tmp[i]));
+      std::transform((double *)tmp2,((double *)tmp2)+sz,(double *)tmp2,std::bind(std::multiplies<double>(),std::placeholders::_1,tmp[i]));
       std::transform(res,res+sz,(double *)tmp2,res,std::plus<double>());
     }
 }
@@ -1721,8 +1712,10 @@ MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationGauss::getMeasureField(con
 {
   if(!mesh)
     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getMeasureField : mesh instance specified is NULL !");
-  MCAuto<MEDCouplingFieldDouble> vol=mesh->getMeasureField(isAbs);
-  const double *volPtr=vol->getArray()->begin();
+  MCAuto<MEDCouplingUMesh> umesh(mesh->buildUnstructured());
+  const double *coordsOfMesh( umesh->getCoords()->begin() );
+  auto spaceDim(mesh->getSpaceDimension());
+  auto meshDim(mesh->getMeshDimension());
   MCAuto<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(ON_GAUSS_PT);
   ret->setMesh(mesh);
   ret->setDiscretization(const_cast<MEDCouplingFieldDiscretizationGauss *>(this));
@@ -1731,7 +1724,7 @@ MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationGauss::getMeasureField(con
   _discr_per_cell->checkAllocated();
   if(_discr_per_cell->getNumberOfComponents()!=1)
     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getMeasureField : no discr per cell array defined but with nb of components different from 1 !");
-  if(_discr_per_cell->getNumberOfTuples()!=vol->getNumberOfTuples())
+  if(_discr_per_cell->getNumberOfTuples()!=mesh->getNumberOfCells())
     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 !");
   MCAuto<DataArrayIdType> offset=getOffsetArr(mesh);
   MCAuto<DataArrayDouble> arr=DataArrayDouble::New(); arr->alloc(getNumberOfTuples(mesh),1);
@@ -1751,11 +1744,28 @@ MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationGauss::getMeasureField(con
           const MEDCouplingGaussLocalization& loc=_loc[locId];
           mcIdType nbOfGaussPt=loc.getNumberOfGaussPt();
           INTERP_KERNEL::AutoPtr<double> weights=new double[nbOfGaussPt];
-          double sum=std::accumulate(loc.getWeights().begin(),loc.getWeights().end(),0.);
-          std::transform(loc.getWeights().begin(),loc.getWeights().end(),(double *)weights,std::bind2nd(std::multiplies<double>(),1./sum));
           for(const mcIdType *cellId=curIds->begin();cellId!=curIds->end();cellId++)
-            for(mcIdType j=0;j<nbOfGaussPt;j++)
-              arrPtr[offsetPtr[*cellId]+j]=weights[j]*volPtr[*cellId];
+          {
+            std::vector<mcIdType> conn;
+            umesh->getNodeIdsOfCell(*cellId,conn);
+            std::vector<double> ptsInCell; ptsInCell.reserve(conn.size()*loc.getDimension());
+            std::for_each( conn.cbegin(), conn.cend(), [spaceDim,coordsOfMesh,&ptsInCell](mcIdType c) { ptsInCell.insert(ptsInCell.end(),coordsOfMesh+c*spaceDim,coordsOfMesh+(c+1)*spaceDim); } );
+            std::size_t nbPtsInCell(ptsInCell.size()/spaceDim);
+            INTERP_KERNEL::DenseMatrix jacobian(spaceDim,meshDim);
+            MCAuto<DataArrayDouble> shapeFunc = loc.getDerivativeOfShapeFunctionValues();
+            for(mcIdType iGPt = 0 ; iGPt < nbOfGaussPt ; ++iGPt)
+            {
+              for(auto i = 0 ; i < spaceDim ; ++i)
+                for(auto j = 0 ; j < meshDim ; ++j)
+                {
+                  double res = 0.0;
+                  for( std::size_t k = 0 ; k < nbPtsInCell ; ++k )
+                    res += ptsInCell[spaceDim*k+i] * shapeFunc->getIJ(iGPt,meshDim*k+j);
+                  jacobian[ i ][ j ] = res;
+                }
+              arrPtr[offsetPtr[*cellId]+iGPt]=std::abs( jacobian.toJacobian() )*loc.getWeight(FromIdType<int>(iGPt));
+            }
+          }
         }
       else
         {
@@ -2372,7 +2382,7 @@ void MEDCouplingFieldDiscretizationGaussNE::integral(const MEDCouplingMesh *mesh
       const double *wArr=GetWeightArrayFromGeometricType(*it,wArrSz);
       INTERP_KERNEL::AutoPtr<double> wArr2=new double[wArrSz];
       double sum=std::accumulate(wArr,wArr+wArrSz,0.);
-      std::transform(wArr,wArr+wArrSz,(double *)wArr2,std::bind2nd(std::multiplies<double>(),1./sum));      
+      std::transform(wArr,wArr+wArrSz,(double *)wArr2,std::bind(std::multiplies<double>(),std::placeholders::_1,1./sum));        
       MCAuto<DataArrayIdType> ids=mesh->giveCellsWithType(*it);
       MCAuto<DataArrayIdType> ids2=ids->buildExplicitArrByRanges(nbOfNodesPerCell);
       const mcIdType *ptIds2=ids2->begin(),*ptIds=ids->begin();
@@ -2436,6 +2446,9 @@ const double *MEDCouplingFieldDiscretizationGaussNE::GetWeightArrayFromGeometric
     case INTERP_KERNEL::NORM_PENTA15:
       lgth=sizeof(FGP_PENTA15)/sizeof(double);
       return FGP_PENTA15;
+    case INTERP_KERNEL::NORM_PENTA18:
+      lgth=sizeof(FGP_PENTA18)/sizeof(double);
+      return FGP_PENTA18;
     case INTERP_KERNEL::NORM_HEXA8:
       lgth=sizeof(FGP_HEXA8)/sizeof(double);
       return FGP_HEXA8;
@@ -2452,7 +2465,7 @@ const double *MEDCouplingFieldDiscretizationGaussNE::GetWeightArrayFromGeometric
       lgth=sizeof(FGP_PYRA13)/sizeof(double);
       return FGP_PYRA13;
     default:
-      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 !");
+      throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::GetWeightArrayFromGeometricType : only SEG[2,3,4], TRI[3,6,7], QUAD[4,9], TETRA[4,10], PENTA[6,15,18], HEXA[8,20,27], PYRA[5,13] supported !");
   }
 }
 
@@ -2502,6 +2515,9 @@ const double *MEDCouplingFieldDiscretizationGaussNE::GetRefCoordsFromGeometricTy
     case INTERP_KERNEL::NORM_PENTA15:
       lgth=sizeof(REF_PENTA15)/sizeof(double);
       return REF_PENTA15;
+    case INTERP_KERNEL::NORM_PENTA18:
+      lgth=sizeof(REF_PENTA18)/sizeof(double);
+      return REF_PENTA18;
     case INTERP_KERNEL::NORM_HEXA8:
       lgth=sizeof(REF_HEXA8)/sizeof(double);
       return REF_HEXA8;
@@ -2518,7 +2534,7 @@ const double *MEDCouplingFieldDiscretizationGaussNE::GetRefCoordsFromGeometricTy
       lgth=sizeof(REF_PYRA13)/sizeof(double);
       return REF_PYRA13;
     default:
-      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 !");
+      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,18], HEXA[8,20,27], PYRA[5,13] supported !");
   }
 }
 
@@ -2596,6 +2612,11 @@ const double *MEDCouplingFieldDiscretizationGaussNE::GetLocsFromGeometricType(IN
         lgth=sizeof(LOC_PENTA15)/sizeof(double);
         return LOC_PENTA15;
       }
+    case INTERP_KERNEL::NORM_PENTA18:
+      {
+        lgth=sizeof(LOC_PENTA18)/sizeof(double);
+        return LOC_PENTA18;
+      }
     case INTERP_KERNEL::NORM_HEXA8:
       {
         lgth=sizeof(LOC_HEXA8)/sizeof(double);
@@ -2622,7 +2643,7 @@ const double *MEDCouplingFieldDiscretizationGaussNE::GetLocsFromGeometricType(IN
         return LOC_PYRA13;
       }
     default:
-      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 !");
+      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,18], HEXA[8,20,27], PYRA[5,13] supported !");
   }
 }
 
@@ -2690,7 +2711,7 @@ MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationGaussNE::getMeasureField(c
       const double *wArr=GetWeightArrayFromGeometricType(*it,wArrSz);
       INTERP_KERNEL::AutoPtr<double> wArr2=new double[wArrSz];
       double sum=std::accumulate(wArr,wArr+wArrSz,0.);
-      std::transform(wArr,wArr+wArrSz,(double *)wArr2,std::bind2nd(std::multiplies<double>(),1./sum));      
+      std::transform(wArr,wArr+wArrSz,(double *)wArr2,std::bind(std::multiplies<double>(),std::placeholders::_1,1./sum));     
       MCAuto<DataArrayIdType> ids=mesh->giveCellsWithType(*it);
       MCAuto<DataArrayIdType> ids2=ids->buildExplicitArrByRanges(nbOfNodesPerCell);
       const mcIdType *ptIds2=ids2->begin(),*ptIds=ids->begin();