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)
int nbOfCompo=arr->getNumberOfComponents();
std::fill(res,res+nbOfCompo,0.);
//
- MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> vol=getMeasureField(mesh,isWAbs);
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldDouble> vol=mesh->getMeasureField(isWAbs);
std::set<INTERP_KERNEL::NormalizedCellType> types=mesh->getAllGeoTypes();
MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nbOfNodesPerCell=mesh->computeNbOfNodesPerCell();
nbOfNodesPerCell->computeOffsets2();
{
if(!mesh)
throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGaussNE::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_NE);
+ ret->setMesh(mesh);
+ //
+ std::set<INTERP_KERNEL::NormalizedCellType> types=mesh->getAllGeoTypes();
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nbOfNodesPerCell=mesh->computeNbOfNodesPerCell();
+ int nbTuples=nbOfNodesPerCell->accumulate(0);
+ nbOfNodesPerCell->computeOffsets2();
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr=DataArrayDouble::New(); arr->alloc(nbTuples,1);
+ ret->setArray(arr);
+ double *arrPtr=arr->getPointer();
+ 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++)
+ for(std::size_t j=0;j<wArrSz;j++,ptIds2++)
+ arrPtr[*ptIds2]=wArr2[j]*volPtr[*ptIds];
+ }
+ ret->synchronizeTimeWithSupport();
+ return ret.retn();
}
void MEDCouplingFieldDiscretizationGaussNE::getValueOn(const DataArrayDouble *arr, const MEDCouplingMesh *mesh, const double *loc, double *res) const
#include "MEDCouplingAutoRefCountObjectPtr.hxx"
#include "MEDCouplingNatureOfField.hxx"
+#include "InterpKernelAutoPtr.hxx"
+
#include <sstream>
#include <limits>
#include <algorithm>
* \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<MEDCouplingFieldDouble> 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<DataArrayDouble> arr=getArray()->deepCpy();
+ arr->multiplyEqual(w->getArray());
+ std::transform(arr->begin(),arr->end(),arr->getPointer(),std::bind2nd(std::multiplies<double>(),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<double> res=new double[nbComps];
+ getWeightedAverageValue(res,isWAbs);
+ return res[compId];
}
/*!
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<double> res=new double[nbComps];
+ _type->normL1(_mesh,getArray(),res);
+ return res[compId];
}
/*!
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<double> res=new double[nbComps];
+ _type->normL2(_mesh,getArray(),res);
+ return res[compId];
}
/*!
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<double> res=new double[nbComps];
+ _type->integral(_mesh,getArray(),isWAbs,res);
+ return res[compId];
}
/*!
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);
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);
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);
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
%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;
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);
int sz=self->getNumberOfComponents();
INTERP_KERNEL::AutoPtr<double> 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<double> 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<double> 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<double> 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<double> 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)
{
