integral in GAUSS_NE fields

diff --git a/src/INTERP_KERNEL/CellModel.cxx b/src/INTERP_KERNEL/CellModel.cxx
index 6a027a81990dc74a0211972229b7f598368e7d18..f1d517d575429bd80310c3f1f7f6a01f80293ec4 100644 (file)
--- a/src/INTERP_KERNEL/CellModel.cxx
+++ b/src/INTERP_KERNEL/CellModel.cxx
@@ -113,6 +113,7 @@ namespace INTERP_KERNEL
     _quadratic=false;
     _dyn=false;
     _extruded_type=NORM_ERROR;
+    _reverse_extruded_type=NORM_ERROR;
     _linear_type=NORM_ERROR;
     _quadratic_type=NORM_ERROR;
     switch(type)
@@ -124,7 +125,7 @@ namespace INTERP_KERNEL
         break;
       case NORM_SEG2:
         {
-          _nb_of_pts=2; _nb_of_sons=2; _dim=1; _extruded_type=NORM_QUAD4; _quadratic_type=NORM_SEG3; _is_simplex=true; _is_extruded=true;
+          _nb_of_pts=2; _nb_of_sons=2; _dim=1; _extruded_type=NORM_QUAD4; _quadratic_type=NORM_SEG3; _is_simplex=true; _is_extruded=true; _reverse_extruded_type=NORM_POINT1;
           _sons_type[0]=NORM_POINT1; _sons_type[1]=NORM_POINT1;
           _sons_con[0][0]=0; _nb_of_sons_con[0]=1;
           _sons_con[1][0]=1; _nb_of_sons_con[1]=1;
@@ -161,7 +162,7 @@ namespace INTERP_KERNEL
         break;
       case NORM_HEXA8:
         {
-          _nb_of_pts=8; _nb_of_sons=6; _dim=3; _quadratic_type=NORM_HEXA20; _is_simplex=false; _is_extruded=true;
+          _nb_of_pts=8; _nb_of_sons=6; _dim=3; _quadratic_type=NORM_HEXA20; _is_simplex=false; _is_extruded=true; _reverse_extruded_type=NORM_QUAD4;
           _sons_type[0]=NORM_QUAD4; _sons_type[1]=NORM_QUAD4; _sons_type[2]=NORM_QUAD4; _sons_type[3]=NORM_QUAD4; _sons_type[4]=NORM_QUAD4; _sons_type[5]=NORM_QUAD4;
           _sons_con[0][0]=0; _sons_con[0][1]=1; _sons_con[0][2]=2; _sons_con[0][3]=3; _nb_of_sons_con[0]=4;
           _sons_con[1][0]=4; _sons_con[1][1]=7; _sons_con[1][2]=6; _sons_con[1][3]=5; _nb_of_sons_con[1]=4;
@@ -241,7 +242,7 @@ namespace INTERP_KERNEL
         break;
       case NORM_PENTA6:
         {
-          _nb_of_pts=6; _nb_of_sons=5; _dim=3; _quadratic_type=NORM_PENTA15; _is_simplex=false; _is_extruded=true;
+          _nb_of_pts=6; _nb_of_sons=5; _dim=3; _quadratic_type=NORM_PENTA15; _is_simplex=false; _is_extruded=true; _reverse_extruded_type=NORM_TRI3;
           _sons_type[0]=NORM_TRI3; _sons_type[1]=NORM_TRI3; _sons_type[2]=NORM_QUAD4; _sons_type[3]=NORM_QUAD4; _sons_type[4]=NORM_QUAD4;
           _sons_con[0][0]=0; _sons_con[0][1]=1; _sons_con[0][2]=2; _nb_of_sons_con[0]=3;
           _sons_con[1][0]=3; _sons_con[1][1]=5; _sons_con[1][2]=4; _nb_of_sons_con[1]=3;

diff --git a/src/INTERP_KERNEL/CellModel.hxx b/src/INTERP_KERNEL/CellModel.hxx
index 27ec69f056128202b6c7583b282dce033f31b000..62c7828ee2c2f4c1540862fdd29c15d1b3ede0b5 100644 (file)
--- a/src/INTERP_KERNEL/CellModel.hxx
+++ b/src/INTERP_KERNEL/CellModel.hxx
@@ -58,6 +58,7 @@ namespace INTERP_KERNEL
     INTERPKERNEL_EXPORT unsigned getNumberOfNodesConstituentTheSon(unsigned sonId) const { return _nb_of_sons_con[sonId]; }
     INTERPKERNEL_EXPORT unsigned getNumberOfNodesConstituentTheSon2(unsigned sonId, const int *nodalConn, int lgth) const;
     INTERPKERNEL_EXPORT NormalizedCellType getExtrudedType() const { return _extruded_type; }
+    INTERPKERNEL_EXPORT NormalizedCellType getReverseExtrudedType() const { return _reverse_extruded_type; }
     INTERPKERNEL_EXPORT NormalizedCellType getLinearType() const { return _linear_type; }
     INTERPKERNEL_EXPORT NormalizedCellType getQuadraticType() const { return _quadratic_type; }
     INTERPKERNEL_EXPORT NormalizedCellType getSonType(unsigned sonId) const { return _sons_type[sonId]; }
@@ -74,6 +75,7 @@ namespace INTERP_KERNEL
     unsigned _nb_of_sons;
     NormalizedCellType _type;
     NormalizedCellType _extruded_type;
+    NormalizedCellType _reverse_extruded_type;
     NormalizedCellType _linear_type;
     NormalizedCellType _quadratic_type;
     unsigned _sons_con[MAX_NB_OF_SONS][MAX_NB_OF_NODES_PER_ELEM];

diff --git a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx
index c045279d262248bc88eade712c9501c2c29cfde6..f63772167eab5059862574aac9a13a5bd128cbc4 100644 (file)
--- a/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx
+++ b/src/INTERP_KERNEL/GaussPoints/InterpKernelGaussCoords.hxx
@@ -17,8 +17,8 @@
 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
 //
 
-#ifndef __INTERPKERNELGAUSS_HXX__
-#define __INTERPKERNELGAUSS_HXX__
+#ifndef __INTERPKERNELGAUSSCOORDS_HXX__
+#define __INTERPKERNELGAUSSCOORDS_HXX__
 
 #include "INTERPKERNELDefines.hxx"
 #include "NormalizedUnstructuredMesh.hxx"
@@ -157,4 +157,4 @@ namespace INTERP_KERNEL
     GaussInfoVector _my_gauss_info;
   };
 }
-#endif //INTERPKERNELGAUSS
+#endif //INTERPKERNELGAUSSCOORDS

diff --git a/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx b/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx
index 875889dadc19cb4d5f2d0c61752c5d8594fe2c16..e2a184545e146dae3898d684da7b477d9b878820 100644 (file)
--- a/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingExtrudedMesh.cxx
@@ -252,6 +252,41 @@ std::set<INTERP_KERNEL::NormalizedCellType> MEDCouplingExtrudedMesh::getAllGeoTy
   return ret;
 }
 
+DataArrayInt *MEDCouplingExtrudedMesh::giveCellsWithType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception)
+{
+  const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(type);
+  INTERP_KERNEL::NormalizedCellType revExtTyp=cm.getReverseExtrudedType();
+  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
+  if(revExtTyp==INTERP_KERNEL::NORM_ERROR)
+    {
+      ret->alloc(0,1);
+      return ret.retn();
+    }
+  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> tmp=_mesh2D->giveCellsWithType(revExtTyp);
+  int nbOfLevs=_mesh1D->getNumberOfCells();
+  int nbOfCells2D=_mesh2D->getNumberOfCells();
+  int nbOfTuples=tmp->getNumberOfTuples();
+  ret->alloc(nbOfLevs*nbOfTuples,1);
+  int *pt=ret->getPointer();
+  for(int i=0;i<nbOfLevs;i++,pt+=nbOfTuples)
+    std::transform(tmp->begin(),tmp->end(),pt,std::bind2nd(std::plus<double>(),i*nbOfCells2D));
+  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret2=ret->renumberR(_mesh3D_ids->begin());
+  ret2->sort();
+  return ret2.retn();
+}
+
+DataArrayInt *MEDCouplingExtrudedMesh::computeNbOfNodesPerCell() const throw(INTERP_KERNEL::Exception)
+{
+  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret2D=_mesh2D->computeNbOfNodesPerCell();
+  int nbOfLevs=_mesh1D->getNumberOfCells();
+  int nbOfCells2D=_mesh2D->getNumberOfCells();
+  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret3D=DataArrayInt::New(); ret3D->alloc(nbOfLevs*nbOfCells2D,1);
+  int *pt=ret3D->getPointer();
+  for(int i=0;i<nbOfLevs;i++,pt+=nbOfCells2D)
+     std::copy(ret2D->begin(),ret2D->end(),pt);
+  return ret3D->renumberR(_mesh3D_ids->begin());
+}
+
 int MEDCouplingExtrudedMesh::getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const
 {
   int ret=0;

diff --git a/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx b/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx
index 8f9bc9c0918ae98377964298cd36f33ce3c8cc82..951a591a4f579fcdf6978321282d90484ce9678c 100644 (file)
--- a/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx
+++ b/src/MEDCoupling/MEDCouplingExtrudedMesh.hxx
@@ -55,6 +55,8 @@ namespace ParaMEDMEM
                                          DataArrayInt *&cellCor) const throw(INTERP_KERNEL::Exception);
     INTERP_KERNEL::NormalizedCellType getTypeOfCell(int cellId) const;
     std::set<INTERP_KERNEL::NormalizedCellType> getAllGeoTypes() const;
+    DataArrayInt *giveCellsWithType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception);
+    DataArrayInt *computeNbOfNodesPerCell() const throw(INTERP_KERNEL::Exception);
     int getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const;
     void getNodeIdsOfCell(int cellId, std::vector<int>& conn) const;
     void getCoordinatesOfNode(int nodeId, std::vector<double>& coo) const throw(INTERP_KERNEL::Exception);

diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
index 461efc9aa4161ff148172c9e0b33d99913dbbde9..b1f8912105c539ee4538b6ac25f45fb50bcea0a8 100644 (file)
--- a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
+++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
@@ -34,6 +34,7 @@
 #include <list>
 #include <limits>
 #include <sstream>
+#include <numeric>
 #include <algorithm>
 #include <functional>
 
@@ -63,6 +64,20 @@ const char MEDCouplingFieldDiscretizationKriging::REPR[]="KRIGING";
 
 const TypeOfField MEDCouplingFieldDiscretizationKriging::TYPE=ON_NODES_KR;
 
+const double MEDCouplingFieldDiscretizationGaussNE::FGP_SEG2[2]={1.,1.};
+const double MEDCouplingFieldDiscretizationGaussNE::FGP_SEG3[3]={0.5555555555555556,0.5555555555555556,0.8888888888888888};
+const double MEDCouplingFieldDiscretizationGaussNE::FGP_SEG4[4]={0.347854845137454,0.347854845137454,0.652145154862546,0.652145154862546};
+const double MEDCouplingFieldDiscretizationGaussNE::FGP_TRI3[3]={0.16666666666666666,0.16666666666666666,0.16666666666666666};
+const double MEDCouplingFieldDiscretizationGaussNE::FGP_TRI6[6]={0.0549758718227661,0.0549758718227661,0.0549758718227661,0.11169079483905,0.11169079483905,0.11169079483905};
+const double MEDCouplingFieldDiscretizationGaussNE::FGP_TRI7[7]={0.062969590272413,0.062969590272413,0.062969590272413,0.066197076394253,0.066197076394253,0.066197076394253,0.1125};
+const double MEDCouplingFieldDiscretizationGaussNE::FGP_QUAD4[4]={1.,1.,1.,1.};
+const double MEDCouplingFieldDiscretizationGaussNE::FGP_QUAD9[9]={0.30864197530864196,0.30864197530864196,0.30864197530864196,0.30864197530864196,0.49382716049382713,0.49382716049382713,0.49382716049382713,0.49382716049382713,0.7901234567901234};
+const double MEDCouplingFieldDiscretizationGaussNE::FGP_TETRA4[4]={0.041666666666666664,0.041666666666666664,0.041666666666666664,0.041666666666666664};
+const double MEDCouplingFieldDiscretizationGaussNE::FGP_PENTA6[6]={0.16666666666666666,0.16666666666666666,0.16666666666666666,0.16666666666666666,0.16666666666666666,0.16666666666666666};
+const double MEDCouplingFieldDiscretizationGaussNE::FGP_HEXA8[8]={1.,1.,1.,1.,1.,1.,1.,1.};
+const double MEDCouplingFieldDiscretizationGaussNE::FGP_HEXA27[27]={0.1714677640603567,0.1714677640603567,0.1714677640603567,0.1714677640603567,0.1714677640603567,0.1714677640603567,0.1714677640603567,0.1714677640603567,0.27434842249657065,0.27434842249657065,0.27434842249657065,0.27434842249657065,0.27434842249657065,0.27434842249657065,0.27434842249657065,0.27434842249657065,0.27434842249657065,0.27434842249657065,0.27434842249657065,0.27434842249657065,0.43895747599451296,0.43895747599451296,0.43895747599451296,0.43895747599451296,0.43895747599451296,0.43895747599451296,0.7023319615912208};
+const double MEDCouplingFieldDiscretizationGaussNE::FGP_PYRA5[5]={0.13333333333333333,0.13333333333333333,0.13333333333333333,0.13333333333333333,0.13333333333333333};
+
 MEDCouplingFieldDiscretization::MEDCouplingFieldDiscretization():_precision(DFLT_PRECISION)
 {
 }
@@ -140,7 +155,7 @@ std::size_t MEDCouplingFieldDiscretization::getHeapMemorySize() const
  */
 void MEDCouplingFieldDiscretization::normL1(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, double *res) const throw(INTERP_KERNEL::Exception)
 {
-  MEDCouplingFieldDouble *vol=getMeasureField(mesh,true);
+  MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> vol=getMeasureField(mesh,true);
   int nbOfCompo=arr->getNumberOfComponents();
   int nbOfElems=getNumberOfTuples(mesh);
   std::fill(res,res+nbOfCompo,0.);
@@ -155,7 +170,6 @@ void MEDCouplingFieldDiscretization::normL1(const MEDCouplingMesh *mesh, const D
       deno+=v;
     }
   std::transform(res,res+nbOfCompo,res,std::bind2nd(std::multiplies<double>(),1./deno));
-  vol->decrRef();
 }
 
 /*!
@@ -165,7 +179,7 @@ void MEDCouplingFieldDiscretization::normL1(const MEDCouplingMesh *mesh, const D
  */
 void MEDCouplingFieldDiscretization::normL2(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, double *res) const throw(INTERP_KERNEL::Exception)
 {
-  MEDCouplingFieldDouble *vol=getMeasureField(mesh,true);
+  MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> vol=getMeasureField(mesh,true);
   int nbOfCompo=arr->getNumberOfComponents();
   int nbOfElems=getNumberOfTuples(mesh);
   std::fill(res,res+nbOfCompo,0.);
@@ -181,7 +195,6 @@ void MEDCouplingFieldDiscretization::normL2(const MEDCouplingMesh *mesh, const D
     }
   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));
-  vol->decrRef();
 }
 
 /*!
@@ -191,7 +204,7 @@ void MEDCouplingFieldDiscretization::normL2(const MEDCouplingMesh *mesh, const D
  */
 void MEDCouplingFieldDiscretization::integral(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, bool isWAbs, double *res) const throw(INTERP_KERNEL::Exception)
 {
-  MEDCouplingFieldDouble *vol=getMeasureField(mesh,isWAbs);
+  MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> vol=getMeasureField(mesh,isWAbs);
   int nbOfCompo=arr->getNumberOfComponents();
   int nbOfElems=getNumberOfTuples(mesh);
   std::fill(res,res+nbOfCompo,0.);
@@ -204,7 +217,6 @@ void MEDCouplingFieldDiscretization::integral(const MEDCouplingMesh *mesh, const
       std::transform(tmp,tmp+nbOfCompo,res,res,std::plus<double>());
     }
   delete [] tmp;
-  vol->decrRef();
 }
 
 void MEDCouplingFieldDiscretization::getSerializationIntArray(DataArrayInt *& arr) const
@@ -309,7 +321,7 @@ void MEDCouplingFieldDiscretization::RenumberEntitiesFromO2NArr(double eps, cons
   int oldNbOfElems=arr->getNumberOfTuples();
   int nbOfComp=arr->getNumberOfComponents();
   int newNbOfTuples=newNbOfEntity;
-  DataArrayDouble *arrCpy=arr->deepCpy();
+  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arrCpy=arr->deepCpy();
   const double *ptSrc=arrCpy->getConstPointer();
   arr->reAlloc(newNbOfTuples);
   double *ptToFill=arr->getPointer();
@@ -330,7 +342,6 @@ void MEDCouplingFieldDiscretization::RenumberEntitiesFromO2NArr(double eps, cons
               //if(!std::equal(ptSrc+i*nbOfComp,ptSrc+(i+1)*nbOfComp,ptToFill+newNb*nbOfComp))
               if(*std::max_element((double *)tmp,((double *)tmp)+nbOfComp)>eps)
                 {
-                  arrCpy->decrRef();
                   std::ostringstream oss;
                   oss << msg << " " << i << " and " << std::find(old2NewPtr,old2NewPtr+i,newNb)-old2NewPtr
                       << " have been merged and " << msg << " field on them are different !";
@@ -339,13 +350,12 @@ void MEDCouplingFieldDiscretization::RenumberEntitiesFromO2NArr(double eps, cons
             }
         }
     }
-  arrCpy->decrRef();
 }
 
 void MEDCouplingFieldDiscretization::RenumberEntitiesFromN2OArr(const int *new2OldPtr, int new2OldSz, DataArrayDouble *arr, const char *msg)
 {
   int nbOfComp=arr->getNumberOfComponents();
-  DataArrayDouble *arrCpy=arr->deepCpy();
+  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arrCpy=arr->deepCpy();
   const double *ptSrc=arrCpy->getConstPointer();
   arr->reAlloc(new2OldSz);
   double *ptToFill=arr->getPointer();
@@ -354,7 +364,6 @@ void MEDCouplingFieldDiscretization::RenumberEntitiesFromN2OArr(const int *new2O
       int oldNb=new2OldPtr[i];
       std::copy(ptSrc+oldNb*nbOfComp,ptSrc+(oldNb+1)*nbOfComp,ptToFill+i*nbOfComp);
     }
-  arrCpy->decrRef();
 }
 
 MEDCouplingFieldDiscretization::~MEDCouplingFieldDiscretization()
@@ -1593,6 +1602,91 @@ DataArrayDouble *MEDCouplingFieldDiscretizationGaussNE::getLocalizationOfDiscVal
   throw INTERP_KERNEL::Exception("Not implemented yet !");
 }
 
+void MEDCouplingFieldDiscretizationGaussNE::integral(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, bool isWAbs, double *res) const throw(INTERP_KERNEL::Exception)
+{
+  if(!mesh || !arr)
+    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::integral : input mesh or array is null !");
+  int nbOfCompo=arr->getNumberOfComponents();
+  int nbOfTuples=arr->getNumberOfTuples();
+  std::fill(res,res+nbOfCompo,0.);
+  //
+  MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> vol=getMeasureField(mesh,isWAbs);
+  std::set<INTERP_KERNEL::NormalizedCellType> types=mesh->getAllGeoTypes();
+  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nbOfNodesPerCell=mesh->computeNbOfNodesPerCell();
+  nbOfNodesPerCell->computeOffsets2();
+  const double *arrPtr=arr->begin(),*volPtr=vol->getArray()->begin();
+  for(std::set<INTERP_KERNEL::NormalizedCellType>::const_iterator it=types.begin();it!=types.end();it++)
+    {
+      std::size_t wArrSz=-1;
+      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));      
+      MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ids=mesh->giveCellsWithType(*it);
+      MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ids2=ids->buildExplicitArrByRanges(nbOfNodesPerCell);
+      const int *ptIds2=ids2->begin(),*ptIds=ids->begin();
+      int nbOfCellsWithCurGeoType=ids->getNumberOfTuples();
+      for(int i=0;i<nbOfCellsWithCurGeoType;i++,ptIds++,ptIds2+=wArrSz)
+        {
+          for(int k=0;k<nbOfCompo;k++)
+            {
+              double tmp=0.;
+              for(std::size_t j=0;j<wArrSz;j++)
+                tmp+=arrPtr[nbOfCompo*ptIds2[j]+k]*wArr2[j];
+              res[k]+=tmp*volPtr[*ptIds];
+            }
+        }
+    }
+}
+
+const double *MEDCouplingFieldDiscretizationGaussNE::GetWeightArrayFromGeometricType(INTERP_KERNEL::NormalizedCellType geoType, std::size_t& lgth) throw(INTERP_KERNEL::Exception)
+{
+  switch(geoType)
+    {
+    case INTERP_KERNEL::NORM_SEG2:
+      lgth=(int)sizeof(FGP_SEG2)/sizeof(double);
+      return FGP_SEG2;
+    case INTERP_KERNEL::NORM_SEG3:
+      lgth=(int)sizeof(FGP_SEG3)/sizeof(double);
+      return FGP_SEG3;
+    case INTERP_KERNEL::NORM_SEG4:
+      lgth=(int)sizeof(FGP_SEG4)/sizeof(double);
+      return FGP_SEG4;
+    case INTERP_KERNEL::NORM_TRI3:
+      lgth=(int)sizeof(FGP_TRI3)/sizeof(double);
+      return FGP_TRI3;
+    case INTERP_KERNEL::NORM_TRI6:
+      lgth=(int)sizeof(FGP_TRI6)/sizeof(double);
+      return FGP_TRI6;
+    case INTERP_KERNEL::NORM_TRI7:
+      lgth=(int)sizeof(FGP_TRI7)/sizeof(double);
+      return FGP_TRI7;
+    case INTERP_KERNEL::NORM_QUAD4:
+      lgth=(int)sizeof(FGP_QUAD4)/sizeof(double);
+      return FGP_QUAD4;
+    case INTERP_KERNEL::NORM_QUAD9:
+      lgth=(int)sizeof(FGP_QUAD9)/sizeof(double);
+      return FGP_QUAD9;
+    case INTERP_KERNEL::NORM_TETRA4:
+      lgth=(int)sizeof(FGP_TETRA4)/sizeof(double);
+      return FGP_TETRA4;
+    case INTERP_KERNEL::NORM_PENTA6:
+      lgth=(int)sizeof(FGP_PENTA6)/sizeof(double);
+      return FGP_PENTA6;
+    case INTERP_KERNEL::NORM_HEXA8:
+      lgth=(int)sizeof(FGP_HEXA8)/sizeof(double);
+      return FGP_HEXA8;
+    case INTERP_KERNEL::NORM_HEXA27:
+      lgth=(int)sizeof(FGP_HEXA27)/sizeof(double);
+      return FGP_HEXA27;
+    case INTERP_KERNEL::NORM_PYRA5:
+      lgth=(int)sizeof(FGP_PYRA5)/sizeof(double);
+      return FGP_PYRA5;
+    default:
+      throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::GetWeightArrayFromGeometricType : only SEG[2,3,4], TRI[3,6,7], QUAD[4,9], TETRA4, PENTA6, HEXA[8,27], PYRA5 supported !");
+    }
+}
+
 void MEDCouplingFieldDiscretizationGaussNE::computeMeshRestrictionFromTupleIds(const MEDCouplingMesh *mesh, const int *partBg, const int *partEnd,
                                                                                DataArrayInt *&cellRest)
 {
@@ -1665,8 +1759,7 @@ DataArrayInt *MEDCouplingFieldDiscretizationGaussNE::computeTupleIdsToSelectFrom
 {
   if(!mesh)
     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::computeTupleIdsToSelectFromCellIds : null mesh !");
-  const MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> umesh=mesh->buildUnstructured();
-  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nbOfNodesPerCell=umesh->computeNbOfNodesPerCell();
+  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nbOfNodesPerCell=mesh->computeNbOfNodesPerCell();
   nbOfNodesPerCell->computeOffsets2();
   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> sel=DataArrayInt::New(); sel->useArray(startCellIds,false,CPP_DEALLOC,(int)std::distance(startCellIds,endCellIds),1);
   return sel->buildExplicitArrByRanges(nbOfNodesPerCell);

diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx
index 8e194911f71a02c4d5ce0df43a621dd2bc8c50fc..c4930dc0184e28eab707d138d60c166e8bf6c5ea 100644 (file)
--- a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx
+++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx
@@ -286,6 +286,7 @@ namespace ParaMEDMEM
     void renumberArraysForCell(const MEDCouplingMesh *mesh, const std::vector<DataArrayDouble *>& arrays,
                                const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception);
     DataArrayDouble *getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const;
+    void integral(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, bool isWAbs, double *res) const throw(INTERP_KERNEL::Exception);
     void computeMeshRestrictionFromTupleIds(const MEDCouplingMesh *mesh, const int *partBg, const int *partEnd,
                                             DataArrayInt *&cellRest);
     void checkCompatibilityWithNature(NatureOfField nat) const throw(INTERP_KERNEL::Exception);
@@ -300,11 +301,29 @@ namespace ParaMEDMEM
     void renumberValuesOnNodes(double epsOnVals, const int *old2New, int newNbOfNodes, DataArrayDouble *arr) const;
     void renumberValuesOnCells(double epsOnVals, const MEDCouplingMesh *mesh, const int *old2New, int newSz, DataArrayDouble *arr) const;
     void renumberValuesOnCellsR(const MEDCouplingMesh *mesh, const int *new2old, int newSz, DataArrayDouble *arr) const;
+    static const double *GetWeightArrayFromGeometricType(INTERP_KERNEL::NormalizedCellType geoType, std::size_t& lgth) throw(INTERP_KERNEL::Exception);
   protected:
     MEDCouplingFieldDiscretizationGaussNE(const MEDCouplingFieldDiscretizationGaussNE& other);
   public:
     static const char REPR[];
     static const TypeOfField TYPE;
+    static const double FGP_SEG2[2];
+    static const double FGP_SEG3[3];
+    static const double FGP_SEG4[4];
+    static const double FGP_TRI3[3];
+    static const double FGP_TRI6[6];
+    static const double FGP_TRI7[7];
+    static const double FGP_QUAD4[4];
+    //static const double FGP_QUAD8[8];
+    static const double FGP_QUAD9[9];
+    static const double FGP_TETRA4[4];
+    //static const double FGP_TETRA10[10];
+    static const double FGP_PENTA6[6];
+    //static const double FGP_PENTA15[15];
+    static const double FGP_HEXA8[8];
+    static const double FGP_HEXA27[27];
+    static const double FGP_PYRA5[5];
+    //static const double FGP_PYRA13[13];
   };
 
   class MEDCOUPLING_EXPORT MEDCouplingFieldDiscretizationKriging : public MEDCouplingFieldDiscretizationOnNodes

diff --git a/src/MEDCoupling/MEDCouplingMesh.hxx b/src/MEDCoupling/MEDCouplingMesh.hxx
index 9f86a7445ba6bab3ee53eaa6d65d4a7d3201c440..bf508c8ebb3b3ca104c68a7109ce81d68a44c8a2 100644 (file)
--- a/src/MEDCoupling/MEDCouplingMesh.hxx
+++ b/src/MEDCoupling/MEDCouplingMesh.hxx
@@ -84,6 +84,8 @@ namespace ParaMEDMEM
     virtual int getMeshDimension() const = 0;
     virtual DataArrayDouble *getCoordinatesAndOwner() const = 0;
     virtual DataArrayDouble *getBarycenterAndOwner() const = 0;
+    virtual DataArrayInt *giveCellsWithType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception) = 0;
+    virtual DataArrayInt *computeNbOfNodesPerCell() const throw(INTERP_KERNEL::Exception) = 0;
     virtual int getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const = 0;
     virtual INTERP_KERNEL::NormalizedCellType getTypeOfCell(int cellId) const = 0;
     virtual std::set<INTERP_KERNEL::NormalizedCellType> getAllGeoTypes() const = 0;

diff --git a/src/MEDCoupling/MEDCouplingStructuredMesh.cxx b/src/MEDCoupling/MEDCouplingStructuredMesh.cxx
index 4d02d732b9c8bb04d9be523b481080742e434786..f473cf7bd0731ff79081933e6a255e7df4cdd427 100644 (file)
--- a/src/MEDCoupling/MEDCouplingStructuredMesh.cxx
+++ b/src/MEDCoupling/MEDCouplingStructuredMesh.cxx
@@ -71,45 +71,42 @@ INTERP_KERNEL::NormalizedCellType MEDCouplingStructuredMesh::getTypeOfCell(int c
 
 std::set<INTERP_KERNEL::NormalizedCellType> MEDCouplingStructuredMesh::getAllGeoTypes() const
 {
-  INTERP_KERNEL::NormalizedCellType ret;
-  switch(getMeshDimension())
-    {
-    case 3:
-      ret=INTERP_KERNEL::NORM_HEXA8;
-      break;
-    case 2:
-      ret=INTERP_KERNEL::NORM_QUAD4;
-      break;
-    case 1:
-      ret=INTERP_KERNEL::NORM_SEG2;
-      break;
-    default:
-      throw INTERP_KERNEL::Exception("Unexpected dimension for MEDCouplingStructuredMesh::getAllGeoTypes !");
-    }
   std::set<INTERP_KERNEL::NormalizedCellType> ret2;
-  ret2.insert(ret);
+  ret2.insert(getTypeOfCell(0));
   return ret2;
 }
 
 int MEDCouplingStructuredMesh::getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const
 {
   int ret=getNumberOfCells();
-  int dim=getMeshDimension();
-  switch(type)
+  if(type==getTypeOfCell(0))
+    return ret;
+  const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(getTypeOfCell(0));
+  std::ostringstream oss; oss << "MEDCouplingStructuredMesh::getNumberOfCellsWithType : no specified type ! Type available is " << cm.getRepr() << " !";
+  throw INTERP_KERNEL::Exception(oss.str().c_str());
+}
+
+DataArrayInt *MEDCouplingStructuredMesh::giveCellsWithType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception)
+{
+  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
+  if(getTypeOfCell(0)==type)
     {
-    case INTERP_KERNEL::NORM_HEXA8:
-      if(dim==3)
-        return ret;
-    case INTERP_KERNEL::NORM_QUAD4:
-      if(dim==2)
-        return ret;
-    case INTERP_KERNEL::NORM_SEG2:
-      if(dim==1)
-        return ret;
-    default:
-      throw INTERP_KERNEL::Exception("Unexpected dimension for MEDCouplingStructuredMesh::getTypeOfCell !");
+      ret->alloc(getNumberOfCells(),1);
+      ret->iota(0);
     }
-  return 0;
+  else
+    ret->alloc(0,1);
+  return ret.retn();
+}
+
+DataArrayInt *MEDCouplingStructuredMesh::computeNbOfNodesPerCell() const throw(INTERP_KERNEL::Exception)
+{
+  int nbCells=getNumberOfCells();
+  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ret=DataArrayInt::New();
+  ret->alloc(nbCells,1);
+  const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel(getTypeOfCell(0));
+  ret->fillWithValue((int)cm.getNumberOfNodes());
+  return ret.retn();
 }
 
 void MEDCouplingStructuredMesh::getNodeIdsOfCell(int cellId, std::vector<int>& conn) const

diff --git a/src/MEDCoupling/MEDCouplingStructuredMesh.hxx b/src/MEDCoupling/MEDCouplingStructuredMesh.hxx
index 80f18ec9fdfb130cd408063bcbb165104f82059b..786d4d2d58c09a3a8a921c1fe67e046be9ca1a64 100644 (file)
--- a/src/MEDCoupling/MEDCouplingStructuredMesh.hxx
+++ b/src/MEDCoupling/MEDCouplingStructuredMesh.hxx
@@ -32,6 +32,8 @@ namespace ParaMEDMEM
     INTERP_KERNEL::NormalizedCellType getTypeOfCell(int cellId) const;
     std::set<INTERP_KERNEL::NormalizedCellType> getAllGeoTypes() const;
     int getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const;
+    DataArrayInt *giveCellsWithType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception);
+    DataArrayInt *computeNbOfNodesPerCell() const throw(INTERP_KERNEL::Exception);
     static void GetPosFromId(int nodeId, int meshDim, const int *split, int *res);
     void getNodeIdsOfCell(int cellId, std::vector<int>& conn) const;
     std::size_t getHeapMemorySize() const;

diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx
index 998a2a0434b766f80cac2332a772f58dc9c1c41c..f492da9f1d8bdb7cdb327097b906ca1ac6b90b0f 100644 (file)
--- a/src/MEDCoupling/MEDCouplingUMeshDesc.cxx
+++ b/src/MEDCoupling/MEDCouplingUMeshDesc.cxx
@@ -158,6 +158,16 @@ std::set<INTERP_KERNEL::NormalizedCellType> MEDCouplingUMeshDesc::getAllGeoTypes
   return _types;
 }
 
+DataArrayInt *MEDCouplingUMeshDesc::giveCellsWithType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception)
+{
+  throw INTERP_KERNEL::Exception("MEDCouplingUMeshDesc::giveCellsWithType : not implemented yet !");
+}
+
+DataArrayInt *MEDCouplingUMeshDesc::computeNbOfNodesPerCell() const throw(INTERP_KERNEL::Exception)
+{
+  throw INTERP_KERNEL::Exception("MEDCouplingUMeshDesc::computeNbOfNodesPerCell : not implemented yet !");
+}
+
 int MEDCouplingUMeshDesc::getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const
 {
   const int *desc_connec=_desc_connec->getConstPointer();

diff --git a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx
index 18d3eb74a0e0d116916117cc5044a75be5c5f2b8..86ef28f73d18a9ca91f8306ee67c1ca736c596a3 100644 (file)
--- a/src/MEDCoupling/MEDCouplingUMeshDesc.hxx
+++ b/src/MEDCoupling/MEDCouplingUMeshDesc.hxx
@@ -51,6 +51,8 @@ namespace ParaMEDMEM
     MEDCOUPLING_EXPORT int getMeshDimension() const { return _mesh_dim; }
     MEDCOUPLING_EXPORT INTERP_KERNEL::NormalizedCellType getTypeOfCell(int cellId) const;
     MEDCOUPLING_EXPORT std::set<INTERP_KERNEL::NormalizedCellType> getAllGeoTypes() const;
+    MEDCOUPLING_EXPORT DataArrayInt *giveCellsWithType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception);
+    MEDCOUPLING_EXPORT DataArrayInt *computeNbOfNodesPerCell() const throw(INTERP_KERNEL::Exception);
     MEDCOUPLING_EXPORT int getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const;
     MEDCOUPLING_EXPORT void getNodeIdsOfCell(int cellId, std::vector<int>& conn) const;
     MEDCOUPLING_EXPORT std::string simpleRepr() const;

diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i
index 47ae7c96d35528ab09035e26f0a673d9fe7f1d9a..5cd1591ef7bb74d80843c9f11bed6063d01d89f0 100644 (file)
--- a/src/MEDCoupling_Swig/MEDCouplingCommon.i
+++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i
@@ -251,6 +251,8 @@ using namespace INTERP_KERNEL;
 %newobject ParaMEDMEM::DataArrayDoubleTuple::buildDADouble;
 %newobject ParaMEDMEM::MEDCouplingMesh::deepCpy;
 %newobject ParaMEDMEM::MEDCouplingMesh::checkTypeConsistencyAndContig;
+%newobject ParaMEDMEM::MEDCouplingMesh::computeNbOfNodesPerCell;
+%newobject ParaMEDMEM::MEDCouplingMesh::giveCellsWithType;
 %newobject ParaMEDMEM::MEDCouplingMesh::getCoordinatesAndOwner;
 %newobject ParaMEDMEM::MEDCouplingMesh::getBarycenterAndOwner;
 %newobject ParaMEDMEM::MEDCouplingMesh::buildOrthogonalField;
@@ -277,13 +279,11 @@ using namespace INTERP_KERNEL;
 %newobject ParaMEDMEM::MEDCouplingUMesh::__iter__;
 %newobject ParaMEDMEM::MEDCouplingUMesh::__getitem__;
 %newobject ParaMEDMEM::MEDCouplingUMesh::cellsByType;
-%newobject ParaMEDMEM::MEDCouplingUMesh::giveCellsWithType;
 %newobject ParaMEDMEM::MEDCouplingUMesh::zipConnectivityTraducer;
 %newobject ParaMEDMEM::MEDCouplingUMesh::buildDescendingConnectivity;
 %newobject ParaMEDMEM::MEDCouplingUMesh::buildDescendingConnectivity2;
 %newobject ParaMEDMEM::MEDCouplingUMesh::buildExtrudedMesh;
 %newobject ParaMEDMEM::MEDCouplingUMesh::buildSpreadZonesWithPoly;
-%newobject ParaMEDMEM::MEDCouplingUMesh::computeNbOfNodesPerCell;
 %newobject ParaMEDMEM::MEDCouplingUMesh::MergeUMeshes;
 %newobject ParaMEDMEM::MEDCouplingUMesh::MergeUMeshesOnSameCoords;
 %newobject ParaMEDMEM::MEDCouplingUMesh::ComputeSpreadZoneGradually;
@@ -491,6 +491,8 @@ namespace ParaMEDMEM
     virtual int getMeshDimension() const throw(INTERP_KERNEL::Exception);
     virtual DataArrayDouble *getCoordinatesAndOwner() const throw(INTERP_KERNEL::Exception);
     virtual DataArrayDouble *getBarycenterAndOwner() const throw(INTERP_KERNEL::Exception);
+    virtual DataArrayInt *giveCellsWithType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception);
+    virtual DataArrayInt *computeNbOfNodesPerCell() const throw(INTERP_KERNEL::Exception);
     virtual int getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception);
     virtual INTERP_KERNEL::NormalizedCellType getTypeOfCell(int cellId) const throw(INTERP_KERNEL::Exception);
     virtual std::string simpleRepr() const;
@@ -1375,7 +1377,6 @@ namespace ParaMEDMEM
     int getNumberOfNodesInCell(int cellId) const throw(INTERP_KERNEL::Exception);
     int getMeshLength() const throw(INTERP_KERNEL::Exception);
     void computeTypes() throw(INTERP_KERNEL::Exception);
-    DataArrayInt *giveCellsWithType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception);
     std::string reprConnectivityOfThis() const throw(INTERP_KERNEL::Exception);
     MEDCouplingUMesh *buildSetInstanceFromThis(int spaceDim) const throw(INTERP_KERNEL::Exception);
     //tools
@@ -1390,7 +1391,6 @@ namespace ParaMEDMEM
     DataArrayInt *convertCellArrayPerGeoType(const DataArrayInt *da) const throw(INTERP_KERNEL::Exception);
     DataArrayInt *computeFetchedNodeIds() const throw(INTERP_KERNEL::Exception);
     DataArrayInt *zipConnectivityTraducer(int compType, int startCellId=0) throw(INTERP_KERNEL::Exception);
-    DataArrayInt *computeNbOfNodesPerCell() const throw(INTERP_KERNEL::Exception);
     MEDCouplingUMesh *buildDescendingConnectivity(DataArrayInt *desc, DataArrayInt *descIndx, DataArrayInt *revDesc, DataArrayInt *revDescIndx) const throw(INTERP_KERNEL::Exception);
     MEDCouplingUMesh *buildDescendingConnectivity2(DataArrayInt *desc, DataArrayInt *descIndx, DataArrayInt *revDesc, DataArrayInt *revDescIndx) const throw(INTERP_KERNEL::Exception);
     void orientCorrectlyPolyhedrons() throw(INTERP_KERNEL::Exception);