From 7edaf7aa3d424acb45a52788c88b2bfb30e5e88f Mon Sep 17 00:00:00 2001 From: ageay Date: Fri, 8 Mar 2013 11:16:26 +0000 Subject: [PATCH] Implementation of integral and getWeightedAverageValue for fields with spatial discretization set to GAUSS_NE Modification of API of MEDCouplingFieldDouble::getWeightedAverageValue --- .../MEDCouplingFieldDiscretization.cxx | 35 ++++++- src/MEDCoupling/MEDCouplingFieldDouble.cxx | 98 ++++++++++--------- src/MEDCoupling/MEDCouplingFieldDouble.hxx | 3 +- .../Test/MEDCouplingBasicsTest2.cxx | 2 +- src/MEDCoupling_Swig/MEDCouplingBasicsTest.py | 74 +++++++++++++- src/MEDCoupling_Swig/MEDCouplingCommon.i | 28 ++++-- 6 files changed, 181 insertions(+), 59 deletions(-) diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx index e9bd76c9a..65b9f7e4e 100644 --- a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx +++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx @@ -1604,6 +1604,9 @@ DataArrayDouble *MEDCouplingFieldDiscretizationGaussNE::getLocalizationOfDiscVal throw INTERP_KERNEL::Exception("Not implemented yet !"); } +/*! + * Reimplemented from MEDCouplingFieldDiscretization::integral for performance reason. The default implementation is valid too for GAUSS_NE spatial discretization. + */ void MEDCouplingFieldDiscretizationGaussNE::integral(const MEDCouplingMesh *mesh, const DataArrayDouble *arr, bool isWAbs, double *res) const throw(INTERP_KERNEL::Exception) { if(!mesh || !arr) @@ -1611,7 +1614,7 @@ void MEDCouplingFieldDiscretizationGaussNE::integral(const MEDCouplingMesh *mesh int nbOfCompo=arr->getNumberOfComponents(); std::fill(res,res+nbOfCompo,0.); // - MEDCouplingAutoRefCountObjectPtr vol=getMeasureField(mesh,isWAbs); + MEDCouplingAutoRefCountObjectPtr vol=mesh->getMeasureField(isWAbs); std::set types=mesh->getAllGeoTypes(); MEDCouplingAutoRefCountObjectPtr nbOfNodesPerCell=mesh->computeNbOfNodesPerCell(); nbOfNodesPerCell->computeOffsets2(); @@ -1725,7 +1728,35 @@ MEDCouplingFieldDouble *MEDCouplingFieldDiscretizationGaussNE::getMeasureField(c { if(!mesh) throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::getMeasureField : mesh instance specified is NULL !"); - throw INTERP_KERNEL::Exception("Not implemented yet !"); + MEDCouplingAutoRefCountObjectPtr vol=mesh->getMeasureField(isAbs); + const double *volPtr=vol->getArray()->begin(); + MEDCouplingAutoRefCountObjectPtr ret=MEDCouplingFieldDouble::New(ON_GAUSS_NE); + ret->setMesh(mesh); + // + std::set types=mesh->getAllGeoTypes(); + MEDCouplingAutoRefCountObjectPtr nbOfNodesPerCell=mesh->computeNbOfNodesPerCell(); + int nbTuples=nbOfNodesPerCell->accumulate(0); + nbOfNodesPerCell->computeOffsets2(); + MEDCouplingAutoRefCountObjectPtr arr=DataArrayDouble::New(); arr->alloc(nbTuples,1); + ret->setArray(arr); + double *arrPtr=arr->getPointer(); + for(std::set::const_iterator it=types.begin();it!=types.end();it++) + { + std::size_t wArrSz=-1; + const double *wArr=GetWeightArrayFromGeometricType(*it,wArrSz); + INTERP_KERNEL::AutoPtr wArr2=new double[wArrSz]; + double sum=std::accumulate(wArr,wArr+wArrSz,0.); + std::transform(wArr,wArr+wArrSz,(double *)wArr2,std::bind2nd(std::multiplies(),1./sum)); + MEDCouplingAutoRefCountObjectPtr ids=mesh->giveCellsWithType(*it); + MEDCouplingAutoRefCountObjectPtr ids2=ids->buildExplicitArrByRanges(nbOfNodesPerCell); + const int *ptIds2=ids2->begin(),*ptIds=ids->begin(); + int nbOfCellsWithCurGeoType=ids->getNumberOfTuples(); + for(int i=0;isynchronizeTimeWithSupport(); + return ret.retn(); } void MEDCouplingFieldDiscretizationGaussNE::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const diff --git a/src/MEDCoupling/MEDCouplingFieldDouble.cxx b/src/MEDCoupling/MEDCouplingFieldDouble.cxx index 03d44c8d0..f2875bf54 100644 --- a/src/MEDCoupling/MEDCouplingFieldDouble.cxx +++ b/src/MEDCoupling/MEDCouplingFieldDouble.cxx @@ -26,6 +26,8 @@ #include "MEDCouplingAutoRefCountObjectPtr.hxx" #include "MEDCouplingNatureOfField.hxx" +#include "InterpKernelAutoPtr.hxx" + #include #include #include @@ -751,17 +753,44 @@ double MEDCouplingFieldDouble::normMax() const throw(INTERP_KERNEL::Exception) * \a this is expected to be a field with exactly \b one component. If not an exception will be thrown. * To getAverageValue on vector field applyFunc is needed before. This method looks only \b default array \b and \b only \b default. * If default array does not exist, an exception will be thrown. + * + * \param [out] res the location where the result will be stored. \a res is expected to be a location with \c this->getNumberOfComponents() places available. + * \param [in] isWAbs specifies if abs is applied on measure on underlying mesh before performing computation. For a user already sure that all cells of its underlying mesh + * are all well oriented this parameter can be set to false to be 'faster'. By default this parameter is true. */ -double MEDCouplingFieldDouble::getWeightedAverageValue() const throw(INTERP_KERNEL::Exception) +void MEDCouplingFieldDouble::getWeightedAverageValue(double *res, bool isWAbs) const throw(INTERP_KERNEL::Exception) { if(getArray()==0) throw INTERP_KERNEL::Exception("MEDCouplingFieldDouble::getWeightedAverageValue : no default array defined !"); - MEDCouplingFieldDouble *w=buildMeasureField(true); + MEDCouplingAutoRefCountObjectPtr w=buildMeasureField(isWAbs); double deno=w->getArray()->accumulate(0); - w->getArray()->multiplyEqual(getArray()); - double res=w->getArray()->accumulate(0); - w->decrRef(); - return res/deno; + MEDCouplingAutoRefCountObjectPtr arr=getArray()->deepCpy(); + arr->multiplyEqual(w->getArray()); + std::transform(arr->begin(),arr->end(),arr->getPointer(),std::bind2nd(std::multiplies(),1./deno)); + arr->accumulate(res); +} + +/*! + * This method returns the average value in \a this weighted by ParaMEDMEM::MEDCouplingField::buildMeasureField. + * \a this is expected to be a field with exactly \b one component. If not an exception will be thrown. + * To getAverageValue on vector field applyFunc is needed before. This method looks only \b default array \b and \b only \b default. + * If default array does not exist, an exception will be thrown. + * + * \param [in] compId The component id that should be in [0, \c this->getNumberOfComponents() ). If not an INTERP_KERNEL::Exception will be thrown. + * \param [in] isWAbs specifies if abs is applied on measure on underlying mesh before performing computation. For a user already sure that all cells of its underlying mesh + * are all well oriented this parameter can be set to false to be 'faster'. By default this parameter is true in C++ not in python (overloading confusion). + */ +double MEDCouplingFieldDouble::getWeightedAverageValue(int compId, bool isWAbs) const throw(INTERP_KERNEL::Exception) +{ + int nbComps=getArray()->getNumberOfComponents(); + if(compId<0 || compId>=nbComps) + { + std::ostringstream oss; oss << "MEDCouplingFieldDouble::getWeightedAverageValue : Invalid compId specified : No such nb of components ! Should be in [0," << nbComps << ") !"; + throw INTERP_KERNEL::Exception(oss.str().c_str()); + } + INTERP_KERNEL::AutoPtr res=new double[nbComps]; + getWeightedAverageValue(res,isWAbs); + return res[compId]; } /*! @@ -776,21 +805,14 @@ double MEDCouplingFieldDouble::normL1(int compId) const throw(INTERP_KERNEL::Exc if(!_mesh) throw INTERP_KERNEL::Exception("No mesh underlying this field to perform normL1"); int nbComps=getArray()->getNumberOfComponents(); - if(compId>=nbComps) - throw INTERP_KERNEL::Exception("Invalid compId specified : No such nb of components !"); - double *res=new double[nbComps]; - try - { - _type->normL1(_mesh,getArray(),res); - } - catch(INTERP_KERNEL::Exception& e) + if(compId<0 || compId>=nbComps) { - delete [] res; - throw e; + std::ostringstream oss; oss << "MEDCouplingFieldDouble::normL1 : Invalid compId specified : No such nb of components ! Should be in [0," << nbComps << ") !"; + throw INTERP_KERNEL::Exception(oss.str().c_str()); } - double ret=res[compId]; - delete [] res; - return ret; + INTERP_KERNEL::AutoPtr res=new double[nbComps]; + _type->normL1(_mesh,getArray(),res); + return res[compId]; } /*! @@ -819,21 +841,14 @@ double MEDCouplingFieldDouble::normL2(int compId) const throw(INTERP_KERNEL::Exc if(!_mesh) throw INTERP_KERNEL::Exception("No mesh underlying this field to perform normL2"); int nbComps=getArray()->getNumberOfComponents(); - if(compId>=nbComps) - throw INTERP_KERNEL::Exception("Invalid compId specified : No such nb of components !"); - double *res=new double[nbComps]; - try - { - _type->normL2(_mesh,getArray(),res); - } - catch(INTERP_KERNEL::Exception& e) + if(compId<0 || compId>=nbComps) { - delete [] res; - throw e; + std::ostringstream oss; oss << "MEDCouplingFieldDouble::normL2 : Invalid compId specified : No such nb of components ! Should be in [0," << nbComps << ") !"; + throw INTERP_KERNEL::Exception(oss.str().c_str()); } - double ret=res[compId]; - delete [] res; - return ret; + INTERP_KERNEL::AutoPtr res=new double[nbComps]; + _type->normL2(_mesh,getArray(),res); + return res[compId]; } /*! @@ -860,21 +875,14 @@ double MEDCouplingFieldDouble::integral(int compId, bool isWAbs) const throw(INT if(!_mesh) throw INTERP_KERNEL::Exception("No mesh underlying this field to perform integral"); int nbComps=getArray()->getNumberOfComponents(); - if(compId>=nbComps) - throw INTERP_KERNEL::Exception("Invalid compId specified : No such nb of components !"); - double *res=new double[nbComps]; - try - { - _type->integral(_mesh,getArray(),isWAbs,res); - } - catch(INTERP_KERNEL::Exception& e) + if(compId<0 || compId>=nbComps) { - delete [] res; - throw e; + std::ostringstream oss; oss << "MEDCouplingFieldDouble::integral : Invalid compId specified : No such nb of components ! Should be in [0," << nbComps << ") !"; + throw INTERP_KERNEL::Exception(oss.str().c_str()); } - double ret=res[compId]; - delete [] res; - return ret; + INTERP_KERNEL::AutoPtr res=new double[nbComps]; + _type->integral(_mesh,getArray(),isWAbs,res); + return res[compId]; } /*! diff --git a/src/MEDCoupling/MEDCouplingFieldDouble.hxx b/src/MEDCoupling/MEDCouplingFieldDouble.hxx index da5e01ba0..493bb4cca 100644 --- a/src/MEDCoupling/MEDCouplingFieldDouble.hxx +++ b/src/MEDCoupling/MEDCouplingFieldDouble.hxx @@ -99,7 +99,8 @@ namespace ParaMEDMEM double getAverageValue() const throw(INTERP_KERNEL::Exception); double norm2() const throw(INTERP_KERNEL::Exception); double normMax() const throw(INTERP_KERNEL::Exception); - double getWeightedAverageValue() const throw(INTERP_KERNEL::Exception); + void getWeightedAverageValue(double *res, bool isWAbs=true) const throw(INTERP_KERNEL::Exception); + double getWeightedAverageValue(int compId, bool isWAbs=true) const throw(INTERP_KERNEL::Exception); double normL1(int compId) const throw(INTERP_KERNEL::Exception); void normL1(double *res) const throw(INTERP_KERNEL::Exception); double normL2(int compId) const throw(INTERP_KERNEL::Exception); diff --git a/src/MEDCoupling/Test/MEDCouplingBasicsTest2.cxx b/src/MEDCoupling/Test/MEDCouplingBasicsTest2.cxx index 05bf33476..f441cbffe 100644 --- a/src/MEDCoupling/Test/MEDCouplingBasicsTest2.cxx +++ b/src/MEDCoupling/Test/MEDCouplingBasicsTest2.cxx @@ -1094,7 +1094,7 @@ void MEDCouplingBasicsTest2::testGetMaxValue1() CPPUNIT_ASSERT_DOUBLES_EQUAL(8.,f->getMaxValue(),1e-14); CPPUNIT_ASSERT_DOUBLES_EQUAL(0.,f->getMinValue(),1e-14); CPPUNIT_ASSERT_DOUBLES_EQUAL(5.,f->getAverageValue(),1e-14); - CPPUNIT_ASSERT_DOUBLES_EQUAL(5.125,f->getWeightedAverageValue(),1e-14); + CPPUNIT_ASSERT_DOUBLES_EQUAL(5.125,f->getWeightedAverageValue(0),1e-14); a1->setIJ(0,2,9.5); CPPUNIT_ASSERT_DOUBLES_EQUAL(9.5,f->getMaxValue(),1e-14); CPPUNIT_ASSERT_DOUBLES_EQUAL(0.,f->getMinValue(),1e-14); diff --git a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py index c583167a7..06111e164 100644 --- a/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py +++ b/src/MEDCoupling_Swig/MEDCouplingBasicsTest.py @@ -2898,7 +2898,7 @@ class MEDCouplingBasicsTest(unittest.TestCase): self.assertAlmostEqual(8.,f.getMaxValue(),14); self.assertAlmostEqual(0.,f.getMinValue(),14); self.assertAlmostEqual(5.,f.getAverageValue(),14); - self.assertAlmostEqual(5.125,f.getWeightedAverageValue(),14); + self.assertAlmostEqual(5.125,f.getWeightedAverageValue(0,True),14); a1.setIJ(0,2,9.5); self.assertAlmostEqual(9.5,f.getMaxValue(),14); self.assertAlmostEqual(0.,f.getMinValue(),14); @@ -11373,6 +11373,78 @@ class MEDCouplingBasicsTest(unittest.TestCase): self.assertTrue(m3D.getNodalConnectivityIndex().isEqual(DataArrayInt([0,28,56,84,112]))) pass + def testSwig2GaussNEIntegral1(self): + m2D=MEDCouplingDataForTest.build2DTargetMesh_1() + m0=m2D[0] ; m0.zipCoords() + m1=m2D[[1,2]] ; m1.zipCoords() + m2=m2D[[3,4]] ; m2.zipCoords() + m0.convertLinearCellsToQuadratic(1) + m1.convertLinearCellsToQuadratic(0) + m2.convertLinearCellsToQuadratic(1) + m=MEDCouplingUMesh.MergeUMeshes([m0,m1,m2]) + m.mergeNodes(1e-12) + f=MEDCouplingFieldDouble(ON_GAUSS_NE) + f.setMesh(m) + arr=DataArrayDouble([1.1,2.2,3.3,4.4,5.5,6.6,7.7,8.8,9.9, + 11.1,12.2,13.3,14.4,15.5,16.6, + 21.1,22.2,23.3,24.4,25.5,26.6, + 31.1,32.2,33.3,34.4,35.5,36.6,37.7,38.8,39.9, + 41.1,42.2,43.3,44.4,45.5,46.6,47.7,48.8,49.9]) + arr2=DataArrayDouble(len(arr),2) + arr2[:,0]=arr ; arr2[:,1]=arr+100 + f.setArray(arr2) + f.checkCoherency() + res=f.integral(False) + # a=25./81 ; b=40./81 ; c=64./81 + # p1=0.11169079483905 ; p2=0.0549758718227661 + # 1st compo + # c0=(a*(1.1+2.2+3.3+4.4)+b*(5.5+6.6+7.7+8.8)+c*9.9)*0.25/3.9999999999999978 ; c0=1.5837962962962973 + # c1=(p2*(11.1+12.2+13.3)+p1*(14.4+15.5+16.6))*0.125/0.4999999999854482 ; c1=1.8014347172346943 + # c2=(p2*(21.1+22.2+23.3)+p1*(24.4+25.5+26.6))*0.125/0.4999999999854482 ; c2=3.0514347172346943 + # c3=(a*(31.1+32.2+33.3+34.4)+b*(35.5+36.6+37.7+38.8)+c*39.9)*0.25/3.9999999999999978 ; c3=9.0837962962963 + # c4=(a*(41.1+42.2+43.3+44.4)+b*(45.5+46.6+47.7+48.8)+c*49.9)*0.25/3.9999999999999978 ; c4=11.583796296296303 + # c0+c1+c2+c3+c4=27.104258323358287 + integExp0=27.104258323358287 + self.assertAlmostEqual(res[0],integExp0,13) + # 2nd compo + # c0=(a*(101.1+102.2+103.3+104.4)+b*(105.5+106.6+107.7+108.8)+c*109.9)*0.25/3.9999999999999978 ; c0=26.58379629629631 + # c1=(p2*(111.1+112.2+113.3)+p1*(114.4+115.5+116.6))*0.125/0.4999999999854482 ; c1=14.301434717234699 + # c2=(p2*(121.1+122.2+123.3)+p1*(124.4+125.5+126.6))*0.125/0.4999999999854482 ; c2=15.5514347172347 + # c3=(a*(131.1+132.2+133.3+134.4)+b*(135.5+136.6+137.7+138.8)+c*139.9)*0.25/3.9999999999999978 ; c3=34.08379629629631 + # c4=(a*(141.1+142.2+143.3+144.4)+b*(145.5+146.6+147.7+148.8)+c*149.9)*0.25/3.9999999999999978 ; c4=36.58379629629632 + # c0+c1+c2+c3+c4=127.10425832335835 + integExp1=127.10425832335835 + self.assertAlmostEqual(res[1],integExp1,12) + meas=f.getDiscretization().getMeasureField(f.getMesh(),False) + intPerTuple=meas*f + res2=intPerTuple.accumulate() + self.assertAlmostEqual(res2[0],integExp0,13) + self.assertAlmostEqual(res2[1],integExp1,12) + # + meas2=f.buildMeasureField(False) + intPerTuple=meas2*f + res3=intPerTuple.accumulate() + self.assertAlmostEqual(res3[0],integExp0,13) + self.assertAlmostEqual(res3[1],integExp1,12) + # + res4=f.getWeightedAverageValue(False) # res4==res2 because sum of area of mesh is equal to 1 + self.assertAlmostEqual(res4[0],integExp0,13) + self.assertAlmostEqual(res4[1],integExp1,12) + # + m.scale([0,0],2.) + # + res5=f.getWeightedAverageValue() # res4==res4 because weighted average is not sensitive to the scaling + self.assertAlmostEqual(res5[0],integExp0,13) + self.assertAlmostEqual(res5[1],integExp1,12) + meas3=f.buildMeasureField(False) + delta=4*meas2.getArray()-meas3.getArray() + delta.abs() + self.assertTrue(delta.isUniform(0.,1e-16)) + res6=f.integral(False) + self.assertAlmostEqual(res6[0],4.*integExp0,12) + self.assertAlmostEqual(res6[1],4.*integExp1,11) + pass + def setUp(self): pass pass diff --git a/src/MEDCoupling_Swig/MEDCouplingCommon.i b/src/MEDCoupling_Swig/MEDCouplingCommon.i index a006ed889..ffe999e44 100644 --- a/src/MEDCoupling_Swig/MEDCouplingCommon.i +++ b/src/MEDCoupling_Swig/MEDCouplingCommon.i @@ -86,6 +86,12 @@ using namespace INTERP_KERNEL; %newobject ParaMEDMEM::MEDCouplingFieldDiscretization::getOffsetArr; %newobject ParaMEDMEM::MEDCouplingFieldDiscretization::clone; %newobject ParaMEDMEM::MEDCouplingFieldDiscretization::clonePart; +%newobject ParaMEDMEM::MEDCouplingFieldDiscretization::getMeasureField; +%newobject ParaMEDMEM::MEDCouplingFieldDiscretization::getOffsetArr; +%newobject ParaMEDMEM::MEDCouplingFieldDiscretization::getLocalizationOfDiscValues; +%newobject ParaMEDMEM::MEDCouplingFieldDiscretization::getValueOnMulti; +%newobject ParaMEDMEM::MEDCouplingFieldDiscretization::computeTupleIdsToSelectFromCellIds; +%newobject ParaMEDMEM::MEDCouplingFieldDiscretization::buildSubMeshData; %newobject ParaMEDMEM::MEDCouplingField::buildMeasureField; %newobject ParaMEDMEM::MEDCouplingField::getLocalizationOfDiscr; %newobject ParaMEDMEM::MEDCouplingField::computeTupleIdsToSelectFromCellIds; @@ -6576,7 +6582,8 @@ namespace ParaMEDMEM double getAverageValue() const throw(INTERP_KERNEL::Exception); double norm2() const throw(INTERP_KERNEL::Exception); double normMax() const throw(INTERP_KERNEL::Exception); - double getWeightedAverageValue() const throw(INTERP_KERNEL::Exception); + //do not put a default value to isWAbs because confusion in python with overloaded getWeightedAverageValue method + double getWeightedAverageValue(int compId, bool isWAbs) const throw(INTERP_KERNEL::Exception); double integral(int compId, bool isWAbs) const throw(INTERP_KERNEL::Exception); double normL1(int compId) const throw(INTERP_KERNEL::Exception); double normL2(int compId) const throw(INTERP_KERNEL::Exception); @@ -6798,32 +6805,35 @@ namespace ParaMEDMEM int sz=self->getNumberOfComponents(); INTERP_KERNEL::AutoPtr tmp=new double[sz]; self->accumulate(tmp); - PyObject *ret=convertDblArrToPyList(tmp,sz); - return ret; + return convertDblArrToPyList(tmp,sz); } PyObject *integral(bool isWAbs) const throw(INTERP_KERNEL::Exception) { int sz=self->getNumberOfComponents(); INTERP_KERNEL::AutoPtr tmp=new double[sz]; self->integral(isWAbs,tmp); - PyObject *ret=convertDblArrToPyList(tmp,sz); - return ret; + return convertDblArrToPyList(tmp,sz); + } + PyObject *getWeightedAverageValue(bool isWAbs=true) const throw(INTERP_KERNEL::Exception) + { + int sz=self->getNumberOfComponents(); + INTERP_KERNEL::AutoPtr tmp=new double[sz]; + self->getWeightedAverageValue(tmp,isWAbs); + return convertDblArrToPyList(tmp,sz); } PyObject *normL1() const throw(INTERP_KERNEL::Exception) { int sz=self->getNumberOfComponents(); INTERP_KERNEL::AutoPtr tmp=new double[sz]; self->normL1(tmp); - PyObject *ret=convertDblArrToPyList(tmp,sz); - return ret; + return convertDblArrToPyList(tmp,sz); } PyObject *normL2() const throw(INTERP_KERNEL::Exception) { int sz=self->getNumberOfComponents(); INTERP_KERNEL::AutoPtr tmp=new double[sz]; self->normL2(tmp); - PyObject *ret=convertDblArrToPyList(tmp,sz); - return ret; + return convertDblArrToPyList(tmp,sz); } void renumberCells(PyObject *li, bool check=true) throw(INTERP_KERNEL::Exception) { -- 2.39.2