From 50a33d75cc8045222df1febc22379554a1ad6828 Mon Sep 17 00:00:00 2001 From: ageay Date: Thu, 24 Feb 2011 09:00:11 +0000 Subject: [PATCH] Gauss Points localization computation in MEDCoupling. --- src/MEDCoupling/MEDCouplingField.cxx | 13 ++ src/MEDCoupling/MEDCouplingField.hxx | 2 + .../MEDCouplingFieldDiscretization.cxx | 151 +++++++++++++++++- .../MEDCouplingFieldDiscretization.hxx | 4 + src/MEDCoupling/MEDCouplingMemArray.cxx | 58 +++++++ src/MEDCoupling/MEDCouplingMemArray.hxx | 3 + src/MEDCoupling/Makefile.am | 2 +- 7 files changed, 229 insertions(+), 4 deletions(-) diff --git a/src/MEDCoupling/MEDCouplingField.cxx b/src/MEDCoupling/MEDCouplingField.cxx index c1879a752..82d13c906 100644 --- a/src/MEDCoupling/MEDCouplingField.cxx +++ b/src/MEDCoupling/MEDCouplingField.cxx @@ -104,6 +104,19 @@ void MEDCouplingField::setNature(NatureOfField nat) throw(INTERP_KERNEL::Excepti _nature=nat; } +/*! + * This method returns is case of success an instance of DataArrayDouble the user is in reponsability to deal with. + * If 'this->_mesh' is not set an exception will be thrown. + * For a field on node the array of coords will be returned. For a field on cell a ParaMEDMEM::DataArrayDouble instance + * containing the barycenter of cells will be returned. And for a field on gauss point the explicit position of gauss points. + */ +DataArrayDouble *MEDCouplingField::getLocalizationOfDiscr() const throw(INTERP_KERNEL::Exception) +{ + if(!_mesh) + throw INTERP_KERNEL::Exception("MEDCouplingField::getLocalizationOfDiscr : No mesh set !"); + return _type->getLocalizationOfDiscValues(_mesh); +} + /*! * This method retrieves the measure field of 'this'. If no '_mesh' is defined an exception will be thrown. * Warning the retrieved field life cycle is the responsability of caller. diff --git a/src/MEDCoupling/MEDCouplingField.hxx b/src/MEDCoupling/MEDCouplingField.hxx index 09244c9cb..840c8a6e7 100644 --- a/src/MEDCoupling/MEDCouplingField.hxx +++ b/src/MEDCoupling/MEDCouplingField.hxx @@ -33,6 +33,7 @@ namespace ParaMEDMEM { class DataArrayInt; + class DataArrayDouble; class MEDCouplingMesh; class MEDCouplingFieldDouble; class MEDCouplingFieldDiscretization; @@ -55,6 +56,7 @@ namespace ParaMEDMEM TypeOfField getTypeOfField() const; NatureOfField getNature() const { return _nature; } virtual void setNature(NatureOfField nat) throw(INTERP_KERNEL::Exception); + DataArrayDouble *getLocalizationOfDiscr() const throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *buildMeasureField(bool isAbs) const throw(INTERP_KERNEL::Exception); MEDCouplingMesh *buildSubMeshData(const int *start, const int *end, DataArrayInt *&di) const; MEDCouplingFieldDiscretization *getDiscretization() const { return _type; } diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx index 92ba48099..898b4d2b3 100644 --- a/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx +++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.cxx @@ -19,13 +19,17 @@ #include "MEDCouplingFieldDiscretization.hxx" #include "MEDCouplingCMesh.hxx" +#include "MEDCouplingUMesh.hxx" #include "MEDCouplingFieldDouble.hxx" -#include "CellModel.hxx" +#include "MEDCouplingAutoRefCountObjectPtr.hxx" +#include "CellModel.hxx" #include "InterpolationUtils.hxx" #include "InterpKernelAutoPtr.hxx" +#include "InterpKernelGaussCoords.hxx" #include +#include #include #include #include @@ -42,6 +46,8 @@ const char MEDCouplingFieldDiscretizationP1::REPR[]="P1"; const TypeOfField MEDCouplingFieldDiscretizationP1::TYPE=ON_NODES; +const int MEDCouplingFieldDiscretizationPerCell::DFT_INVALID_LOCID_VALUE=-1; + const char MEDCouplingFieldDiscretizationGauss::REPR[]="GAUSS"; const TypeOfField MEDCouplingFieldDiscretizationGauss::TYPE=ON_GAUSS_PT; @@ -660,10 +666,19 @@ void MEDCouplingFieldDiscretizationPerCell::buildDiscrPerCellIfNecessary(const M int nbTuples=m->getNumberOfCells(); _discr_per_cell->alloc(nbTuples,1); int *ptr=_discr_per_cell->getPointer(); - std::fill(ptr,ptr+nbTuples,-1); + std::fill(ptr,ptr+nbTuples,DFT_INVALID_LOCID_VALUE); } } +void MEDCouplingFieldDiscretizationPerCell::checkNoOrphanCells() const throw(INTERP_KERNEL::Exception) +{ + if(!_discr_per_cell) + throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell::checkNoOrphanCells : no discretization defined !"); + MEDCouplingAutoRefCountObjectPtr test=_discr_per_cell->getIdsEqual(DFT_INVALID_LOCID_VALUE); + if(test->getNumberOfTuples()!=0) + throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationPerCell::checkNoOrphanCells : presence of orphan cells !"); +} + MEDCouplingFieldDiscretizationGauss::MEDCouplingFieldDiscretizationGauss() { } @@ -761,7 +776,44 @@ void MEDCouplingFieldDiscretizationGauss::renumberArraysForCell(const MEDCouplin DataArrayDouble *MEDCouplingFieldDiscretizationGauss::getLocalizationOfDiscValues(const MEDCouplingMesh *mesh) const { - throw INTERP_KERNEL::Exception("Not implemented yet !"); + checkNoOrphanCells(); + MEDCouplingAutoRefCountObjectPtr umesh=mesh->buildUnstructured();//in general do nothing + int nbOfTuples=getNumberOfTuples(mesh); + DataArrayDouble *ret=DataArrayDouble::New(); + int spaceDim=mesh->getSpaceDimension(); + ret->alloc(nbOfTuples,spaceDim); + std::vector< std::vector > locIds; + std::vector parts=splitIntoSingleGaussDicrPerCellType(locIds); + std::vector< MEDCouplingAutoRefCountObjectPtr > parts2(parts.size()); + std::copy(parts.begin(),parts.end(),parts2.begin()); + MEDCouplingAutoRefCountObjectPtr offsets=buildNbOfGaussPointPerCellField(); + offsets->computeOffsets(); + const int *ptrOffsets=offsets->getConstPointer(); + const double *coords=umesh->getCoords()->getConstPointer(); + const int *connI=umesh->getNodalConnectivityIndex()->getConstPointer(); + const int *conn=umesh->getNodalConnectivity()->getConstPointer(); + double *valsToFill=ret->getPointer(); + for(std::size_t i=0;i::const_iterator it=locIds[i].begin();it!=locIds[i].end();it++) + { + const MEDCouplingGaussLocalization& cli=_loc[*it];//curLocInfo + INTERP_KERNEL::NormalizedCellType typ=cli.getType(); + const std::vector& wg=cli.getWeights(); + calculator.addGaussInfo(typ,INTERP_KERNEL::CellModel::getCellModel(typ).getDimension(), + &cli.getGaussCoords()[0],wg.size(),&cli.getRefCoords()[0], + INTERP_KERNEL::CellModel::getCellModel(typ).getNumberOfNodes()); + } + int nbt=parts2[i]->getNumberOfTuples(); + for(const int *w=parts2[i]->getConstPointer();w!=parts2[i]->getConstPointer()+nbt;w++) + { + const MEDCouplingGaussLocalization& cli=_loc[*w]; + calculator.calculateCoords(cli.getType(),coords,spaceDim,conn+connI[*w]+1,valsToFill+spaceDim*(ptrOffsets[*w])); + } + } + ret->copyStringInfoFrom(*umesh->getCoords()); + return ret; } void MEDCouplingFieldDiscretizationGauss::computeMeshRestrictionFromTupleIds(const MEDCouplingMesh *mesh, const int *partBg, const int *partEnd, @@ -921,6 +973,13 @@ void MEDCouplingFieldDiscretizationGauss::renumberValuesOnCellsR(const MEDCoupli void MEDCouplingFieldDiscretizationGauss::setGaussLocalizationOnType(const MEDCouplingMesh *m, INTERP_KERNEL::NormalizedCellType type, const std::vector& refCoo, const std::vector& gsCoo, const std::vector& wg) throw(INTERP_KERNEL::Exception) { + const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::getCellModel(type); + if((int)cm.getDimension()!=m->getMeshDimension()) + { + std::ostringstream oss; oss << "MEDCouplingFieldDiscretizationGauss::setGaussLocalizationOnType : mismatch of dimensions ! MeshDim==" << m->getMeshDimension(); + oss << " whereas Type '" << cm.getRepr() << "' has dimension " << cm.getDimension() << " !"; + throw INTERP_KERNEL::Exception(oss.str().c_str()); + } buildDiscrPerCellIfNecessary(m); int id=_loc.size(); MEDCouplingGaussLocalization elt(type,refCoo,gsCoo,wg); @@ -1038,6 +1097,30 @@ int MEDCouplingFieldDiscretizationGauss::getOffsetOfCell(int cellId) const throw return ret; } +/*! + * This method do the assumption that there is no orphan cell. If there is an exception is thrown. + * This method makes the assumption too that '_discr_per_cell' is defined. If not an exception is thrown. + * This method returns a newly created array with number of tuples equals to '_discr_per_cell->getNumberOfTuples' and number of components equal to 1. + * The i_th tuple in returned array is the number of gauss point if the corresponding cell. + */ +DataArrayInt *MEDCouplingFieldDiscretizationGauss::buildNbOfGaussPointPerCellField() const throw(INTERP_KERNEL::Exception) +{ + if(!_discr_per_cell) + throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::buildNbOfGaussPointPerCellField : no discretization array set !"); + int nbOfTuples=_discr_per_cell->getNumberOfTuples(); + MEDCouplingAutoRefCountObjectPtr ret=DataArrayInt::New(); + const int *w=_discr_per_cell->getConstPointer(); + ret->alloc(nbOfTuples,1); + int *valsToFill=ret->getPointer(); + for(int i=0;iincrRef(); + return ret; +} + /*! * This method makes the assumption that _discr_per_cell is set. * This method reduces as much as possible number size of _loc. @@ -1073,6 +1156,68 @@ void MEDCouplingFieldDiscretizationGauss::zipGaussLocalizations() _loc=tmpLoc; } +/*! + * This method is usefull 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< std::vector >& locIds) const throw(INTERP_KERNEL::Exception) +{ + if(!_discr_per_cell) + throw INTERP_KERNEL::Exception("MEDCouplingFieldDiscretizationGauss::splitIntoSingleGaussDicrPerCellType : no descretization set !"); + locIds.clear(); + std::vector ret; + const int *discrPerCell=_discr_per_cell->getConstPointer(); + MEDCouplingAutoRefCountObjectPtr ret2=_discr_per_cell->getIdsNotEqual(-1); + int nbOfTuplesSet=ret2->getNumberOfTuples(); + std::list idsRemaining(ret2->getConstPointer(),ret2->getConstPointer()+nbOfTuplesSet); + std::list::iterator it=idsRemaining.begin(); + while(it!=idsRemaining.end()) + { + std::vector ids; + std::set curLocIds; + std::set curCellTypes; + while(it!=idsRemaining.end()) + { + int curDiscrId=discrPerCell[*it]; + INTERP_KERNEL::NormalizedCellType typ=_loc[curDiscrId].getType(); + if(curCellTypes.find(typ)!=curCellTypes.end()) + { + if(curLocIds.find(curDiscrId)!=curLocIds.end()) + { + curLocIds.insert(curDiscrId); + curCellTypes.insert(typ); + ids.push_back(*it); + it=idsRemaining.erase(it); + } + else + it++; + } + else + { + curLocIds.insert(curDiscrId); + curCellTypes.insert(typ); + ids.push_back(*it); + it=idsRemaining.erase(it); + } + } + it=idsRemaining.begin(); + ret.resize(ret.size()+1); + DataArrayInt *part=DataArrayInt::New(); + part->alloc(ids.size(),1); + std::copy(ids.begin(),ids.end(),part->getPointer()); + ret.back()=part; + locIds.resize(locIds.size()+1); + locIds.back().insert(locIds.back().end(),curLocIds.begin(),curLocIds.end()); + } + return ret; +} + MEDCouplingFieldDiscretizationGaussNE::MEDCouplingFieldDiscretizationGaussNE() { } diff --git a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx index 80a3f49ae..2f543fe96 100644 --- a/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx +++ b/src/MEDCoupling/MEDCouplingFieldDiscretization.hxx @@ -164,10 +164,12 @@ namespace ParaMEDMEM bool isEqual(const MEDCouplingFieldDiscretization *other, double eps) 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: DataArrayInt *_discr_per_cell; + static const int DFT_INVALID_LOCID_VALUE; }; class MEDCOUPLING_EXPORT MEDCouplingFieldDiscretizationGauss : public MEDCouplingFieldDiscretizationPerCell @@ -211,6 +213,8 @@ namespace ParaMEDMEM int getGaussLocalizationIdOfOneType(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< std::vector >& locIds) const throw(INTERP_KERNEL::Exception); + DataArrayInt *buildNbOfGaussPointPerCellField() const throw(INTERP_KERNEL::Exception); protected: MEDCouplingFieldDiscretizationGauss(const MEDCouplingFieldDiscretizationGauss& other); void zipGaussLocalizations(); diff --git a/src/MEDCoupling/MEDCouplingMemArray.cxx b/src/MEDCoupling/MEDCouplingMemArray.cxx index 4e0fd0f4f..bb57d0dea 100644 --- a/src/MEDCoupling/MEDCouplingMemArray.cxx +++ b/src/MEDCoupling/MEDCouplingMemArray.cxx @@ -2344,6 +2344,22 @@ DataArrayInt *DataArrayInt::getIdsEqual(int val) const throw(INTERP_KERNEL::Exce return ret; } +DataArrayInt *DataArrayInt::getIdsNotEqual(int val) const throw(INTERP_KERNEL::Exception) +{ + if(getNumberOfComponents()!=1) + throw INTERP_KERNEL::Exception("DataArrayInt::getIdsNotEqual : the array must have only one component !"); + const int *cptr=getConstPointer(); + std::vector res; + int nbOfTuples=getNumberOfTuples(); + for(int i=0;ialloc(res.size(),1); + std::copy(res.begin(),res.end(),ret->getPointer()); + return ret; +} + DataArrayInt *DataArrayInt::getIdsEqualList(const std::vector& vals) const throw(INTERP_KERNEL::Exception) { if(getNumberOfComponents()!=1) @@ -2361,6 +2377,23 @@ DataArrayInt *DataArrayInt::getIdsEqualList(const std::vector& vals) const return ret; } +DataArrayInt *DataArrayInt::getIdsNotEqualList(const std::vector& vals) const throw(INTERP_KERNEL::Exception) +{ + if(getNumberOfComponents()!=1) + throw INTERP_KERNEL::Exception("DataArrayInt::getIdsNotEqualList : the array must have only one component !"); + std::set vals2(vals.begin(),vals.end()); + const int *cptr=getConstPointer(); + std::vector res; + int nbOfTuples=getNumberOfTuples(); + for(int i=0;ialloc(res.size(),1); + std::copy(res.begin(),res.end(),ret->getPointer()); + return ret; +} + DataArrayInt *DataArrayInt::Aggregate(const DataArrayInt *a1, const DataArrayInt *a2, int offsetA2) { int nbOfComp=a1->getNumberOfComponents(); @@ -2676,6 +2709,31 @@ DataArrayInt *DataArrayInt::deltaShiftIndex() const throw(INTERP_KERNEL::Excepti return ret; } +/*! + * This method performs the work on itself. This method works on array with number of component equal to one and allocated. If not an exception is thrown. + * This method conserves the number of tuples and number of components (1). No reallocation is done. + * For an array [3,5,1,2,0,8] it becomes [0,3,8,9,11,11]. For each i>0 new[i]=new[i-1]+old[i-1] for i==0 new[i]=0. + * This could be usefull for allToAllV in MPI with contiguous policy. + */ +void DataArrayInt::computeOffsets() throw(INTERP_KERNEL::Exception) +{ + checkAllocated(); + if(getNumberOfComponents()!=1) + throw INTERP_KERNEL::Exception("DataArrayInt::computeOffsets : only single component allowed !"); + int nbOfTuples=getNumberOfTuples(); + if(nbOfTuples==0) + return ; + int *work=getPointer(); + int tmp=work[0]; + work[0]=0; + for(int i=1;i& vals) const throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT DataArrayInt *getIdsNotEqualList(const std::vector& vals) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT int getMaxValue(int& tupleId) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT int getMinValue(int& tupleId) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT void applyLin(int a, int b, int compoId) throw(INTERP_KERNEL::Exception); @@ -298,6 +300,7 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT DataArrayInt *buildUnion(const DataArrayInt *other) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT DataArrayInt *buildIntersection(const DataArrayInt *other) const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT DataArrayInt *deltaShiftIndex() const throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT void computeOffsets() throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT std::set getDifferentValues() const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT void useArray(const int *array, bool ownership, DeallocType type, int nbOfTuple, int nbOfCompo); MEDCOUPLING_EXPORT void writeOnPlace(int id, int element0, const int *others, int sizeOfOthers) { _mem.writeOnPlace(id,element0,others,sizeOfOthers); } diff --git a/src/MEDCoupling/Makefile.am b/src/MEDCoupling/Makefile.am index 37058285b..14146ff61 100644 --- a/src/MEDCoupling/Makefile.am +++ b/src/MEDCoupling/Makefile.am @@ -57,7 +57,7 @@ libmedcoupling_la_LDFLAGS= libmedcoupling_la_CPPFLAGS= @CXXTMPDPTHFLAGS@ \ -I$(srcdir) -I$(srcdir)/../INTERP_KERNEL/Bases \ -I$(srcdir)/../INTERP_KERNEL -I$(srcdir)/../INTERP_KERNEL/Geometric2D \ - -I$(srcdir)/../INTERP_KERNEL/ExprEval + -I$(srcdir)/../INTERP_KERNEL/ExprEval -I$(srcdir)/../INTERP_KERNEL/GaussPoints # the geom2D library is included in the interpkernel one libmedcoupling_la_LIBADD= ../INTERP_KERNEL/libinterpkernel.la -- 2.39.2