From ac61d12c1547ae4de8d651218ecf63defafd3b86 Mon Sep 17 00:00:00 2001 From: michael Date: Fri, 10 Dec 2021 14:20:33 +0100 Subject: [PATCH] Start using profile to deal with 1) boundary conditions and 2) mesh deletion in memory --- CDMATH/mesh/inc/Field.hxx | 15 +++- CDMATH/mesh/src/Field.cxx | 163 +++++++++++++++++++++++++------------- 2 files changed, 122 insertions(+), 56 deletions(-) diff --git a/CDMATH/mesh/inc/Field.hxx b/CDMATH/mesh/inc/Field.hxx index 542a688..489588c 100755 --- a/CDMATH/mesh/inc/Field.hxx +++ b/CDMATH/mesh/inc/Field.hxx @@ -12,6 +12,7 @@ namespace MEDCoupling { class MEDCouplingFieldDouble; class DataArrayDouble; + class MEDFileField1TS; } #include "DoubleTab.hxx" @@ -72,8 +73,7 @@ class Field */ Field( const std::string filename, EntityType fieldType, const std::string & fieldName = "", - int iteration = -1, int order = -1, int meshLevel=0, - int numberOfComponents=1, double time=0.0); + int iteration = -1, int order = -1, int meshLevel=0); /** * constructor with data @@ -251,6 +251,16 @@ class Field void setFieldByDataArrayDouble ( const MEDCoupling::DataArrayDouble* array ); + /** + * \brief Delete the medcoupling mesh to save memory space + */ + void deleteMEDCouplingUMesh(); + + /** + * \brief Returns true iff an unstructured mesh has been loaded + */ + bool meshNotDeleted() const {return _mesh.meshNotDeleted();} + Vector getNormEuclidean( void ) const ; double max( int component=0 ) const ; @@ -357,6 +367,7 @@ class Field int _numberOfComponents; double _time; std::string _fieldName; + MEDCoupling::MEDFileField1TS *_ff ;//To save the field when the mesh has been deleted private: diff --git a/CDMATH/mesh/src/Field.cxx b/CDMATH/mesh/src/Field.cxx index 268af97..3252115 100755 --- a/CDMATH/mesh/src/Field.cxx +++ b/CDMATH/mesh/src/Field.cxx @@ -10,6 +10,7 @@ #include "Face.hxx" #include "Field.hxx" #include "MEDFileMesh.hxx" +#include "MEDFileField1TS.hxx" #include "MEDLoader.hxx" #include "MEDCouplingUMesh.hxx" #include "MEDCouplingFieldDouble.hxx" @@ -29,6 +30,7 @@ Field::Field( EntityType typeField ) //---------------------------------------------------------------------- { _field=NULL; + _ff=NULL; _typeField=typeField; _numberOfComponents=0; } @@ -45,6 +47,7 @@ Field::~Field( void ) Field::Field(const std::string fieldName, EntityType type, const Mesh& mesh, int numberOfComponents, double time) { _field = NULL; + _ff=NULL; _mesh=Mesh(mesh); _typeField=type; _numberOfComponents=numberOfComponents; @@ -89,19 +92,16 @@ void Field::buildFieldMemoryStructure() _field->setArray(array); _field->setTime(_time,0,0); array->decrRef(); - mu->decrRef(); } Field::Field( const std::string filename, EntityType type, const std::string & fieldName, - int iteration, int order, int meshLevel, - int numberOfComponents, double time) + int iteration, int order, int meshLevel) { _field = NULL; + _ff=NULL; _mesh=Mesh(filename + ".med", meshLevel); _typeField=type; - _numberOfComponents=numberOfComponents; - _time=time; _fieldName=fieldName; readFieldMed(filename, type, fieldName, iteration, order); @@ -111,6 +111,7 @@ Field::Field(const std::string meshFileName, EntityType type, const std::vector< const std::string & fieldName, int meshLevel, double time ) { _field = NULL; + _ff=NULL; _mesh=Mesh(meshFileName + ".med", meshLevel); _typeField=type; _numberOfComponents=Vconstant.size(); @@ -124,11 +125,12 @@ Field::Field(const std::string meshFileName, EntityType type, const std::vector< for (int ielem=0 ; ielemgetArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()]=Vconstant[jcomp]; + _field->getArray()->getPointer()[jcomp+ielem*nbcomp]=Vconstant[jcomp]; } Field::Field(const Mesh& M, EntityType type, const Vector Vconstant, const std::string & fieldName, double time) { _field = NULL; + _ff=NULL; _mesh=Mesh(M); _typeField=type; _numberOfComponents=Vconstant.size(); @@ -142,11 +144,12 @@ Field::Field(const Mesh& M, EntityType type, const Vector Vconstant, const std:: for (int ielem=0 ; ielemgetArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()]=Vconstant[jcomp]; + _field->getArray()->getPointer()[jcomp+ielem*nbcomp]=Vconstant[jcomp]; } Field::Field(const Mesh& M, EntityType type, const vector Vconstant, const std::string & fieldName, double time) { _field = NULL; + _ff=NULL; _mesh=Mesh(M); _typeField=type; _numberOfComponents=Vconstant.size(); @@ -160,7 +163,7 @@ Field::Field(const Mesh& M, EntityType type, const vector Vconstant, con for (int ielem=0 ; ielemgetArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()]=Vconstant[jcomp]; + _field->getArray()->getPointer()[jcomp+ielem*nbcomp]=Vconstant[jcomp]; } Field::Field( int nDim, const vector Vconstant, EntityType type, double xmin, double xmax, int nx, string leftSide, string rightSide, @@ -169,6 +172,7 @@ Field::Field( int nDim, const vector Vconstant, EntityType type, const std::string & fieldName, double time,double epsilon) { _field = NULL; + _ff=NULL; _typeField=type; _numberOfComponents=Vconstant.size(); _time=time; @@ -206,7 +210,7 @@ Field::Field( int nDim, const vector Vconstant, EntityType type, for (int ielem=0 ; ielemgetArray()->getPointer()[jcomp+ielem*_field->getNumberOfComponents()]=Vconstant[jcomp]; + _field->getArray()->getPointer()[jcomp+ielem*nbcomp]=Vconstant[jcomp]; } Field::Field(const Mesh M, const Vector VV_Left, const Vector VV_Right, double disc_pos, EntityType type, int direction, const std::string & fieldName, double time) @@ -215,6 +219,7 @@ Field::Field(const Mesh M, const Vector VV_Left, const Vector VV_Right, double d throw CdmathException( "Field::Field: Vectors VV_Left and VV_Right have different sizes"); _field = NULL; + _ff=NULL; _mesh=Mesh(M); _typeField=type; _numberOfComponents=VV_Left.getNumberOfRows(); @@ -240,9 +245,9 @@ Field::Field(const Mesh M, const Vector VV_Left, const Vector VV_Right, double d for (int i=0; i< nbcomp; i++) if (component_value< disc_pos ) - _field->getArray()->getPointer()[j+i*_field->getNumberOfComponents()] = VV_Left[i]; + _field->getArray()->getPointer()[j+i*nbcomp] = VV_Left[i]; else - _field->getArray()->getPointer()[j+i*_field->getNumberOfComponents()] = VV_Right[i]; + _field->getArray()->getPointer()[j+i*nbcomp] = VV_Right[i]; } } Field::Field( int nDim, const vector VV_Left, vector VV_Right, @@ -292,6 +297,7 @@ Field::Field(const Mesh M, const Vector Vin, const Vector Vout, double radius, } _field = NULL; + _ff=NULL; _mesh=Mesh(M); _typeField=type; _numberOfComponents=Vout.size(); @@ -318,10 +324,10 @@ Field::Field(const Mesh M, const Vector Vin, const Vector Vout, double radius, } if((currentPoint-Center).norm()getArray()->getPointer()[j+i*_field->getNumberOfComponents()]=Vin[j]; + _field->getArray()->getPointer()[j+i*nbcomp]=Vin[j]; else for(int j=0;jgetArray()->getPointer()[j+i*_field->getNumberOfComponents()]=Vout[j]; + _field->getArray()->getPointer()[j+i*nbcomp]=Vout[j]; } } @@ -344,23 +350,33 @@ Field::readFieldMed( const std::string & fileNameRadical, size_t iField = 0; std::string attributedFieldName; _field = NULL; + _ff=NULL; // Get the name of the right field that we will attribute to the Field. if (fieldName == "") { if (fieldNames.size() > 0) + { + cout<<"Warning : No field name imposed, taking the first field name found : "<< fieldNames[0]<<" in file "<0 ) + { attributedFieldName = fieldName; + + if( mycount> 1 ) + cout<<"Warning : " << mycount << " fields are associated to the name " << fieldName <<" in file "< meshNames - = MEDCoupling::GetMeshNamesOnField(completeFileName, attributedFieldName); + std::vector meshNames = MEDCoupling::GetMeshNamesOnField(completeFileName, attributedFieldName); if (meshNames.size() == 0) { std::ostringstream message; - message << "No mesh associated to " << fieldName - << " in file " << completeFileName; + message << "No mesh associated to " << fieldName<< " in file " << completeFileName; throw CdmathException(message.str().c_str()); } + else + if( meshNames.size() > 1 ) + { + cout<<"Warning : " << meshNames.size() << " meshes are associated to field named " << fieldName <<" in file "< ( - MEDCoupling::ReadFieldCell( completeFileName, - attributedMeshName, 0, - attributedFieldName, iteration, order) ); - break; - case NODES: - _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > ( - MEDCoupling::ReadFieldNode( completeFileName, - attributedMeshName, 0, - attributedFieldName, iteration, order) ); - break; - case FACES: - _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > ( - MEDCoupling::ReadFieldCell( completeFileName, - attributedMeshName, -1, - attributedFieldName, iteration, order) ); - break; + switch (type) + { + case CELLS: + _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > ( + MEDCoupling::ReadFieldCell( completeFileName, + attributedMeshName, 0, + attributedFieldName, iteration, order) ); + break; + case NODES: + _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > ( + MEDCoupling::ReadFieldNode( completeFileName, + attributedMeshName, 0, + attributedFieldName, iteration, order) ); + break; + case FACES: + _field = dynamic_cast< MEDCoupling::MEDCouplingFieldDouble * > ( + MEDCoupling::ReadFieldCell( completeFileName, + attributedMeshName, -1, + attributedFieldName, iteration, order) ); + break; } + + //Read and store the number of components + _numberOfComponents = _field->getNumberOfComponents() ; + _time = _field->getTime(iteration, order); + + cout<<"Found field " << fieldName << " in file " << completeFileName <getArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()]+f(ielem, jcomp); + fres(ielem, jcomp)=_field->getArray()->getConstPointer()[jcomp+ielem*nbComp]+f(ielem, jcomp); return fres; } @@ -840,7 +870,7 @@ Field::operator- ( const Field& f ) const int nbElem=f.getNumberOfElements(); for (int ielem=0 ; ielemgetArray()->getConstPointer()[jcomp+ielem*_field->getNumberOfComponents()]-f(ielem, jcomp); + fres(ielem, jcomp)=_field->getArray()->getConstPointer()[jcomp+ielem*nbComp]-f(ielem, jcomp); return fres; } @@ -914,7 +944,7 @@ Field::operator*= ( double s ) int nbElem=getNumberOfElements(); for (int i=0 ; igetArray()->getPointer()[i+j*_field->getNumberOfComponents()]*=s; + _field->getArray()->getPointer()[i+j*nbComp]*=s; return *this; } @@ -927,7 +957,7 @@ Field::operator/= ( double s ) int nbElem=getNumberOfElements(); for (int i=0 ; igetArray()->getPointer()[i+j*_field->getNumberOfComponents()]/=s; + _field->getArray()->getPointer()[i+j*nbComp]/=s; return *this; } @@ -940,7 +970,7 @@ Field::operator-= ( double s ) int nbElem=getNumberOfElements(); for (int i=0 ; igetArray()->getPointer()[i+j*_field->getNumberOfComponents()]-=s; + _field->getArray()->getPointer()[i+j*nbComp]-=s; return *this; } @@ -953,7 +983,7 @@ Field::operator+= ( double s ) int nbElem=getNumberOfElements(); for (int i=0 ; igetArray()->getPointer()[i+j*_field->getNumberOfComponents()]+=s; + _field->getArray()->getPointer()[i+j*nbComp]+=s; return *this; } @@ -962,6 +992,9 @@ void Field::writeVTK (std::string fileName, bool fromScratch) const //---------------------------------------------------------------------- { + if( !_mesh.isStructured() && !_mesh.meshNotDeleted() ) + throw CdmathException("Field::writeVTK : Cannot save field in VTK format : unstructured mesh with no MEDCouplingUMesh loaded. Use med format."); + string fname=fileName+".pvd"; int iter,order; double time=_field->getTime(iter,order); @@ -1099,10 +1132,27 @@ Field::writeMED ( const std::string fileName, bool fromScratch) const //---------------------------------------------------------------------- { string fname=fileName+".med"; - if (fromScratch) - MEDCoupling::WriteField(fname.c_str(),_field,fromScratch); - else - MEDCoupling::WriteFieldUsingAlreadyWrittenMesh(fname.c_str(),_field); + + if(_mesh.isStructured() || _mesh.meshNotDeleted()) + if (fromScratch) + MEDCoupling::WriteField(fname.c_str(),_field,fromScratch); + else + MEDCoupling::WriteFieldUsingAlreadyWrittenMesh(fname.c_str(),_field); + else//The mesh has ben deleted, use _ff instead of _field to save the values + { + //MEDFileUMesh * meshMEDFile = MEDFileUMesh::New(); + //meshMEDFile->setMeshAtLevel(0,_field->getMesh()->buildUnstructured()); + //meshMEDFile->write(fname.c_str(), fromScratch); + //MEDCoupling::WriteUMesh(fname.c_str(),_field->getMesh()->buildUnstructured(),fromScratch); + //MEDCoupling::WriteMesh(fname.c_str(),_field->getMesh(),fromScratch); + //MEDCoupling::MEDCouplingUMesh* fmesh = dynamic_cast (_field->getMesh()->deepCopy()); + //cout<<" checkConsecutiveCellTypes : "<< fmesh->checkConsecutiveCellTypes() <advancedRepr() <getMesh()->buildUnstructured()->checkConsecutiveCellTypes()<setFieldNoProfileSBT( _field ); + _ff->write(fname.c_str(), fromScratch); + } } Field @@ -1173,3 +1223,8 @@ std::ostream& operator<<(std::ostream& out, const Field& field ) out<< field.getField().retn()->getArray()->repr(); return out; } + +void Field::deleteMEDCouplingUMesh() +{ + return _mesh.deleteMEDCouplingUMesh(); +} -- 2.39.2