From: ageay <ageay>
Date: Thu, 21 Mar 2013 10:32:30 +0000 (+0000)
Subject: integral getMeasureField on Gauss point discretization
X-Git-Tag: V6_main_FINAL~256
X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=cfb2cbb1d6c5334973e926cf2a0746c7217bcda0;p=tools%2Fmedcoupling.git

integral getMeasureField on Gauss point discretization
---

diff --git a/src/MEDCoupling/MEDCouplingField.cxx b/src/MEDCoupling/MEDCouplingField.cxx
index 4d75b20eb..3be9efef1 100644
--- a/src/MEDCoupling/MEDCouplingField.cxx
+++ b/src/MEDCoupling/MEDCouplingField.cxx
@@ -388,6 +388,13 @@ int MEDCouplingField::getNumberOfTuplesExpected() const throw(INTERP_KERNEL::Exc
     throw INTERP_KERNEL::Exception("MEDCouplingField::getNumberOfTuplesExpected : Empty mesh !");
 }
 
+void MEDCouplingField::setDiscretization(MEDCouplingFieldDiscretization *newDisc)
+{
+  _type=newDisc;
+  if(newDisc)
+    newDisc->incrRef();
+}
+
 /*!
  * This method returns number of mesh placed expected regarding its discretization and its _mesh attribute.
  * This method expected a not null _mesh instance. If null, an exception will be thrown.
diff --git a/src/MEDCoupling/MEDCouplingField.hxx b/src/MEDCoupling/MEDCouplingField.hxx
index 22b9345ef..e2ad798c1 100644
--- a/src/MEDCoupling/MEDCouplingField.hxx
+++ b/src/MEDCoupling/MEDCouplingField.hxx
@@ -66,7 +66,7 @@ namespace ParaMEDMEM
     DataArrayInt *computeTupleIdsToSelectFromCellIds(const int *startCellIds, const int *endCellIds) const;
     const MEDCouplingFieldDiscretization *getDiscretization() const { return _type; }
     MEDCouplingFieldDiscretization *getDiscretization() { return _type; }
-    void setDiscretization(MEDCouplingFieldDiscretization *newDisc) { _type=newDisc; }
+    void setDiscretization(MEDCouplingFieldDiscretization *newDisc);
     int getNumberOfTuplesExpected() const throw(INTERP_KERNEL::Exception);
     int getNumberOfMeshPlacesExpected() const throw(INTERP_KERNEL::Exception);
     // Gauss point specific methods
diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
index b65d83cb2..16aba6e71 100644
--- a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
+++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx
@@ -232,19 +232,28 @@ 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)
 {
+  if(!mesh)
+    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretization::integral : mesh is NULL !");
+  if(!arr)
+    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretization::integral : input array is NULL !");
   MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> vol=getMeasureField(mesh,isWAbs);
   int nbOfCompo=arr->getNumberOfComponents();
   int nbOfElems=getNumberOfTuples(mesh);
+  if(nbOfElems!=arr->getNumberOfTuples())
+    {
+      std::ostringstream oss; oss << "MEDCouplingFieldDiscretization::integral : field is not correct ! number of tuples in array is " << arr->getNumberOfTuples();
+      oss << " whereas number of tuples expected is " << nbOfElems << " !";
+      throw INTERP_KERNEL::Exception(oss.str().c_str());
+    }
   std::fill(res,res+nbOfCompo,0.);
   const double *arrPtr=arr->getConstPointer();
   const double *volPtr=vol->getArray()->getConstPointer();
-  double *tmp=new double[nbOfCompo];
+  INTERP_KERNEL::AutoPtr<double> tmp=new double[nbOfCompo];
   for (int i=0;i<nbOfElems;i++)
     {
-      std::transform(arrPtr+i*nbOfCompo,arrPtr+(i+1)*nbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),volPtr[i]));
-      std::transform(tmp,tmp+nbOfCompo,res,res,std::plus<double>());
+      std::transform(arrPtr+i*nbOfCompo,arrPtr+(i+1)*nbOfCompo,(double *)tmp,std::bind2nd(std::multiplies<double>(),volPtr[i]));
+      std::transform((double *)tmp,(double *)tmp+nbOfCompo,res,res,std::plus<double>());
     }
-  delete [] tmp;
 }
 
 void MEDCouplingFieldDiscretization::getSerializationIntArray(DataArrayInt *& arr) const
@@ -923,6 +932,23 @@ void MEDCouplingFieldDiscretizationPerCell::checkNoOrphanCells() const throw(INT
     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell::checkNoOrphanCells : presence of orphan cells !");
 }
 
+/*!
+ * This method is useful when 'this' describes a field discretization with several gauss discretization on a \b same cell type.
+ * For example same NORM_TRI3 cells having 6 gauss points and others with 12 gauss points.
+ * This method returns 2 arrays with same size : the return value and 'locIds' output parameter.
+ * 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.
+ * The return vector contains a set of newly created instance to deal with.
+ * The returned vector represents a \b partition of cells ids with a gauss discretization set.
+ * 
+ * If no descretization is set in 'this' and exception will be thrown.
+ */
+std::vector<DataArrayInt *> MEDCouplingFieldDiscretizationPerCell::splitIntoSingleGaussDicrPerCellType(std::vector<int>& locIds) const throw(INTERP_KERNEL::Exception)
+{
+  if(!_discr_per_cell)
+    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell::splitIntoSingleGaussDicrPerCellType : no descretization set !");
+  return _discr_per_cell->partitionByDifferentValues(locIds);
+}
+
 const DataArrayInt *MEDCouplingFieldDiscretizationPerCell::getArrayOfDiscIds() const
 {
   return _discr_per_cell;
@@ -1264,7 +1290,50 @@ MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationGauss::getMeasureField(con
 {
   if(!mesh)
     throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getMeasureField : mesh instance specified is NULL !");
-  throw INTERP_KERNEL::Exception("Not implemented yet !");
+  MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> vol=mesh->getMeasureField(isAbs);
+  const double *volPtr=vol->getArray()->begin();
+  MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> ret=MEDCouplingFieldDouble::New(ON_GAUSS_PT);
+  ret->setMesh(mesh);
+  ret->setDiscretization(const_cast<MEDCouplingFieldDiscretizationGauss *>(this));
+  if(!_discr_per_cell)
+    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::getMeasureField : no discr per cell array not defined ! spatial localization is incorrect !");
+  _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())
+    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 !");
+  MEDCouplingAutoRefCountObjectPtr<DataArrayInt> offset=getOffsetArr(mesh);
+  MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr=DataArrayDouble::New(); arr->alloc(getNumberOfTuples(mesh),1);
+  ret->setArray(arr);
+  double *arrPtr=arr->getPointer();
+  const int *offsetPtr=offset->getConstPointer();
+  int maxGaussLoc=(int)_loc.size();
+  std::vector<int> locIds;
+  std::vector<DataArrayInt *> ids=splitIntoSingleGaussDicrPerCellType(locIds);
+  std::vector< MEDCouplingAutoRefCountObjectPtr<DataArrayInt> > ids2(ids.size()); std::copy(ids.begin(),ids.end(),ids2.begin());
+  for(std::size_t i=0;i<locIds.size();i++)
+    {
+      const DataArrayInt *curIds=ids[i];
+      int locId=locIds[i];
+      if(locId>=0 && locId<maxGaussLoc)
+        {
+          const MEDCouplingGaussLocalization& loc=_loc[locId];
+          int 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 int *cellId=curIds->begin();cellId!=curIds->end();cellId++)
+            for(int j=0;j<nbOfGaussPt;j++)
+              arrPtr[offsetPtr[*cellId]+j]=weights[j]*volPtr[*cellId];
+        }
+      else
+        {
+          std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGauss::getMeasureField : Presence of localization id " << locId << " in cell #" << curIds->getIJ(0,0) << " ! Must be in [0," << maxGaussLoc << ") !";
+          throw INTERP_KERNEL::Exception(oss.str().c_str());
+        }
+    }
+  ret->synchronizeTimeWithSupport();
+  return ret.retn();
 }
 
 void MEDCouplingFieldDiscretizationGauss::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
@@ -1523,23 +1592,6 @@ void MEDCouplingFieldDiscretizationGauss::zipGaussLocalizations()
   _loc=tmpLoc;
 }
 
-/*!
- * This method is useful when 'this' describes a field discretization with several gauss discretization on a \b same cell type.
- * For example same NORM_TRI3 cells having 6 gauss points and others with 12 gauss points.
- * This method returns 2 arrays with same size : the return value and 'locIds' output parameter.
- * 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.
- * The return vector contains a set of newly created instance to deal with.
- * The returned vector represents a \b partition of cells ids with a gauss discretization set.
- * 
- * If no descretization is set in 'this' and exception will be thrown.
- */
-std::vector<DataArrayInt *> MEDCouplingFieldDiscretizationGauss::splitIntoSingleGaussDicrPerCellType(std::vector<int>& locIds) const throw(INTERP_KERNEL::Exception)
-{
-  if(!_discr_per_cell)
-    throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::splitIntoSingleGaussDicrPerCellType : no descretization set !");
-  return _discr_per_cell->partitionByDifferentValues(locIds);
-}
-
 MEDCouplingFieldDiscretizationGaussNE::MEDCouplingFieldDiscretizationGaussNE()
 {
 }
diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx
index 20781780e..7609ead28 100644
--- a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx
+++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx
@@ -188,6 +188,8 @@ namespace ParaMEDMEM
   {
   public:
     const DataArrayInt *getArrayOfDiscIds() const;
+    void checkNoOrphanCells() const throw(INTERP_KERNEL::Exception);
+    std::vector<DataArrayInt *> splitIntoSingleGaussDicrPerCellType(std::vector< int >& locIds) const throw(INTERP_KERNEL::Exception);
   protected:
     MEDCouplingFieldDiscretizationPerCell();
     MEDCouplingFieldDiscretizationPerCell(const MEDCouplingFieldDiscretizationPerCell& other, const int *startCellIds, const int *endCellIds);
@@ -198,7 +200,6 @@ namespace ParaMEDMEM
     bool isEqualIfNotWhy(const MEDCouplingFieldDiscretization *other, double eps, std::string& reason) const;
     bool isEqualWithoutConsideringStr(const MEDCouplingFieldDiscretization *other, double eps) const;
     void renumberCells(const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception);
-    void checkNoOrphanCells() const throw(INTERP_KERNEL::Exception);
   protected:
     void buildDiscrPerCellIfNecessary(const MEDCouplingMesh *m);
   protected:
@@ -255,7 +256,6 @@ namespace ParaMEDMEM
     std::set<int> getGaussLocalizationIdsOfOneType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception);
     void getCellIdsHavingGaussLocalization(int locId, std::vector<int>& cellIds) const throw(INTERP_KERNEL::Exception);
     const MEDCouplingGaussLocalization& getGaussLocalization(int locId) const throw(INTERP_KERNEL::Exception);
-    std::vector<DataArrayInt *> splitIntoSingleGaussDicrPerCellType(std::vector< int >& locIds) const throw(INTERP_KERNEL::Exception);
     DataArrayInt *buildNbOfGaussPointPerCellField() const throw(INTERP_KERNEL::Exception);
   protected:
     MEDCouplingFieldDiscretizationGauss(const MEDCouplingFieldDiscretizationGauss& other, const int *startCellIds=0, const int *endCellIds=0);
diff --git a/src/MEDCoupling/MEDCouplingFieldDouble.cxx b/src/MEDCoupling/MEDCouplingFieldDouble.cxx
index a0107f665..be8873b06 100644
--- a/src/MEDCoupling/MEDCouplingFieldDouble.cxx
+++ b/src/MEDCoupling/MEDCouplingFieldDouble.cxx
@@ -516,7 +516,7 @@ MEDCouplingFieldDouble *MEDCouplingFieldDouble::buildSubPart(const int *partBg,
   MEDCouplingFieldDouble *ret=clone(false);//quick shallow copy.
   const MEDCouplingFieldDiscretization *disc=getDiscretization();
   if(disc)
-    ret->setDiscretization(disc->clonePart(partBg,partEnd));
+    ret->setDiscretization(MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDiscretization>(disc->clonePart(partBg,partEnd)));
   ret->setMesh(m);
   std::vector<DataArrayDouble *> arrays;
   _time_discr->getArrays(arrays);
diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
index 42d87603d..a7228bd95 100644
--- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
+++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py
@@ -11626,6 +11626,24 @@ class MEDCouplingBasicsTest(unittest.TestCase):
         self.assertTrue(loc.isEqual(DataArrayDouble([0.,0.,7.,1.,0.,7.,0.5,1.,7.,0.,0.,7.,1.,0.,7.,1.,1.,7.,0.,1.,7.,0.,0.,7.,1.,0.,7.,0.5,1.,7.,0.5,0.,7.,0.75,0.5,7.,0.25,0.5,7.,0.,0.,7.,1.,0.,7.,1.,1.,7.,0.,1.,7.,0.5,0.,7.,1.,0.5,7.,0.5,1.,7.,0.,0.5,7.,0.,0.,7.,0.5,1.,7.,1.,0.,7.,0.,0.,7.,0.,1.,7.,1.,1.,7.,1.,0.,7.,0.,0.,7.,0.5,1.,7.,1.,0.,7.,0.25,0.5,7.,0.75,0.5,7.,0.5,0.,7.,0.,0.,7.,0.,1.,7.,1.,1.,7.,1.,0.,7.,0.,0.5,7.,0.5,1.,7.,1.,0.5,7.,0.5,0.,7.],42,3),1e-13))
         pass
 
+    def testSwig2GaussMeasureAndIntegral(self):
+        ft=MEDCouplingDataForTest.buildFieldOnGauss_1()
+        mea=ft.buildMeasureField(False)
+        mea.checkCoherency()
+        self.assertTrue(mea.getArray().isEqual(DataArrayDouble([-0.08504076274779823,-0.06378057206084897,-0.08504076274779869,-0.10630095343474463,-0.12756114412169625,-0.10630095343474734,-0.0637805720608491,-0.0850407627477968,-0.1063009534347449,-0.0850407627477994,-0.10630095343474809,-0.1275611441216954,-0.037205333702161475,-0.037205333702161475,-0.037205333702161475,-0.037205333702161475,-0.047835429045636084,-0.047835429045636084,-0.047835429045636084,-0.047835429045636084,-0.05846552438911087,-0.05846552438911087,-0.05846552438911087,-0.05846552438911087,-0.037205333702161725,-0.037205333702161725,-0.037205333702161725,-0.037205333702161725,-0.047835429045635834,-0.047835429045635834,-0.047835429045635834,-0.047835429045635834,-0.05846552438911058,-0.05846552438911058,-0.05846552438911058,-0.05846552438911058,-0.03879154890291829,-0.03879154890291829,-0.03879154890291829,-0.04120270848015563,-0.04120270848015563,-0.04120270848015563,-0.03393028948486933,-0.03393028948486933,-0.03393028948486933,-0.03151955746491709,-0.03151955746491709,-0.03151955746491709,-0.02424752187358276,-0.02424752187358276,-0.02424752187358276,-0.026657914642918758,-0.026657914642918758,-0.026657914642918758,-0.04120270848015456,-0.04120270848015456,-0.04120270848015456,-0.03879154890291757,-0.03879154890291757,-0.03879154890291757,-0.031519557464916595,-0.031519557464916595,-0.031519557464916595,-0.03393028948487046,-0.03393028948487046,-0.03393028948487046,-0.0266579146429191,-0.0266579146429191,-0.0266579146429191,-0.024247521873582645,-0.024247521873582645,-0.024247521873582645,-0.01851718920904466,-0.01851718920904466,-0.01851718920904466,-0.01851718920904466,-0.029627502734471456,-0.029627502734471456,-0.029627502734471456,-0.029627502734471456,-0.04740400437515433,-0.015150427534672922,-0.015150427534672922,-0.015150427534672922,-0.015150427534672922,-0.024240684055476674,-0.024240684055476674,-0.024240684055476674,-0.024240684055476674,-0.038785094488762675,-0.011783665860301345,-0.011783665860301345,-0.011783665860301345,-0.011783665860301345,-0.018853865376482152,-0.018853865376482152,-0.018853865376482152,-0.018853865376482152,-0.030166184602371443,-0.018517189209044892,-0.018517189209044892,-0.018517189209044892,-0.018517189209044892,-0.029627502734471827,-0.029627502734471827,-0.029627502734471827,-0.029627502734471827,-0.04740400437515492,-0.015150427534672776,-0.015150427534672776,-0.015150427534672776,-0.015150427534672776,-0.02424068405547644,-0.02424068405547644,-0.02424068405547644,-0.02424068405547644,-0.03878509448876231,-0.011783665860301277,-0.011783665860301277,-0.011783665860301277,-0.011783665860301277,-0.01885386537648204,-0.01885386537648204,-0.01885386537648204,-0.01885386537648204,-0.030166184602371266]),1e-14))
+        f=MEDCouplingFieldDouble(ft)
+        arr=DataArrayDouble(126,2)
+        arr[:,0]=range(126)
+        arr[:,1]=range(126)
+        arr[:,1]+=1000
+        f.setArray(arr)
+        f.checkCoherency()
+        self.assertTrue(DataArrayDouble(f.integral(False)).isEqual(DataArrayDouble([-211.66121638700983,-4863.9563007698835]),1e-12))
+        self.assertTrue(DataArrayDouble(f.getWeightedAverageValue()).isEqual(DataArrayDouble([45.496085813113595, 1045.496085813114]),1e-12))
+        self.assertTrue(DataArrayDouble(f.normL1()).isEqual(DataArrayDouble([45.49608581311362,1045.496085813114]),1e-12))
+        self.assertTrue(DataArrayDouble(f.normL2()).isEqual(DataArrayDouble([58.16846378340898,1046.1241521947334]),1e-12))
+        pass
+
     def setUp(self):
         pass
     pass
diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i
index f00797c75..30cc0beb3 100644
--- a/src/MEDCoupling_Swig/MEDCouplingCommon.i
+++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i
@@ -2007,6 +2007,24 @@ namespace ParaMEDMEM
       ret->incrRef();
     return SWIG_NewPointerObj(SWIG_as_voidptr(ret),SWIGTYPE_p_ParaMEDMEM__DataArrayInt, SWIG_POINTER_OWN | 0 );
   }
+
+  PyObject *splitIntoSingleGaussDicrPerCellType() const throw(INTERP_KERNEL::Exception)
+  {
+    std::vector<int> ret1;
+    std::vector<DataArrayInt *> ret0=self->splitIntoSingleGaussDicrPerCellType(ret1);
+    std::size_t sz=ret0.size();
+    PyObject *pyRet=PyTuple_New(2);
+    PyObject *pyRet0=PyList_New((int)sz);
+    PyObject *pyRet1=PyList_New((int)sz);
+    for(std::size_t i=0;i<sz;i++)
+      {
+        PyList_SetItem(pyRet0,i,SWIG_NewPointerObj(SWIG_as_voidptr(ret0[i]),SWIGTYPE_p_ParaMEDMEM__DataArrayInt, SWIG_POINTER_OWN | 0 ));
+        PyList_SetItem(pyRet1,i,PyInt_FromLong(ret1[i]));
+      }
+    PyTuple_SetItem(pyRet,0,pyRet0);
+    PyTuple_SetItem(pyRet,1,pyRet1);
+    return pyRet;
+  }
 }
 
 %ignore ParaMEDMEM::DataArray::getInfoOnComponents;