From cfb2cbb1d6c5334973e926cf2a0746c7217bcda0 Mon Sep 17 00:00:00 2001 From: ageay Date: Thu, 21 Mar 2013 10:32:30 +0000 Subject: [PATCH] integral getMeasureField on Gauss point discretization --- src/MEDCoupling/MEDCouplingField.cxx | 7 ++ src/MEDCoupling/MEDCouplingField.hxx | 2 +- .../MEDCouplingFieldDiscretization.cxx | 96 ++++++++++++++----- .../MEDCouplingFieldDiscretization.hxx | 4 +- src/MEDCoupling/MEDCouplingFieldDouble.cxx | 2 +- src/MEDCoupling_Swig/MEDCouplingBasicsTest.py | 18 ++++ src/MEDCoupling_Swig/MEDCouplingCommon.i | 18 ++++ 7 files changed, 121 insertions(+), 26 deletions(-) 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 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 tmp=new double[nbOfCompo]; for (int i=0;i(),volPtr[i])); - std::transform(tmp,tmp+nbOfCompo,res,res,std::plus()); + std::transform(arrPtr+i*nbOfCompo,arrPtr+(i+1)*nbOfCompo,(double *)tmp,std::bind2nd(std::multiplies(),volPtr[i])); + std::transform((double *)tmp,(double *)tmp+nbOfCompo,res,res,std::plus()); } - 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 MEDCouplingFieldDiscretizationPerCell::splitIntoSingleGaussDicrPerCellType(std::vector& 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 vol=mesh->getMeasureField(isAbs); + const double *volPtr=vol->getArray()->begin(); + MEDCouplingAutoRefCountObjectPtr ret=MEDCouplingFieldDouble::New(ON_GAUSS_PT); + ret->setMesh(mesh); + ret->setDiscretization(const_cast(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 offset=getOffsetArr(mesh); + MEDCouplingAutoRefCountObjectPtr 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 locIds; + std::vector ids=splitIntoSingleGaussDicrPerCellType(locIds); + std::vector< MEDCouplingAutoRefCountObjectPtr > ids2(ids.size()); std::copy(ids.begin(),ids.end(),ids2.begin()); + for(std::size_t i=0;i=0 && locId 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(),1./sum)); + for(const int *cellId=curIds->begin();cellId!=curIds->end();cellId++) + for(int j=0;jgetIJ(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 MEDCouplingFieldDiscretizationGauss::splitIntoSingleGaussDicrPerCellType(std::vector& 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 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 getGaussLocalizationIdsOfOneType(INTERP_KERNEL::NormalizedCellType type) const throw(INTERP_KERNEL::Exception); void getCellIdsHavingGaussLocalization(int locId, std::vector& cellIds) const throw(INTERP_KERNEL::Exception); const MEDCouplingGaussLocalization& getGaussLocalization(int locId) const throw(INTERP_KERNEL::Exception); - std::vector 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(disc->clonePart(partBg,partEnd))); ret->setMesh(m); std::vector 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 ret1; + std::vector 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