From 258aed605783a5c801e860879cdd4190b407a7d3 Mon Sep 17 00:00:00 2001 From: michael Date: Thu, 13 Jan 2022 17:56:21 +0100 Subject: [PATCH] Improved writing of meshes and fields --- CDMATH/mesh/inc/Field.hxx | 8 +- CDMATH/mesh/inc/Mesh.hxx | 187 ++-- CDMATH/mesh/src/Field.cxx | 51 +- CDMATH/mesh/src/Mesh.cxx | 1890 +++++++++++++++++-------------------- 4 files changed, 970 insertions(+), 1166 deletions(-) diff --git a/CDMATH/mesh/inc/Field.hxx b/CDMATH/mesh/inc/Field.hxx index 489588c..fe79b58 100755 --- a/CDMATH/mesh/inc/Field.hxx +++ b/CDMATH/mesh/inc/Field.hxx @@ -15,7 +15,6 @@ namespace MEDCoupling class MEDFileField1TS; } -#include "DoubleTab.hxx" #include "Vector.hxx" #include "Mesh.hxx" @@ -87,7 +86,7 @@ class Field * */ Field(const std::string meshfileName, EntityType fieldType, const std::vector Vconstant,const std::string & fieldName = "", - int meshLevel=0, double time=0.0); + int meshLevel=0, double time=0.0, std::string meshName=""); /** * constructor with data @@ -254,7 +253,7 @@ class Field /** * \brief Delete the medcoupling mesh to save memory space */ - void deleteMEDCouplingUMesh(); + void deleteMEDCouplingMesh(); /** * \brief Returns true iff an unstructured mesh has been loaded @@ -347,7 +346,7 @@ class Field void writeVTK ( const std::string fileName, bool fromScratch=true ) const ; - void writeMED ( const std::string fileName, bool fromScratch=true ) const ; + void writeMED ( const std::string fileName, bool fromScratch=true ) ; void writeCSV ( const std::string fileName ) const ; @@ -367,7 +366,6 @@ 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/inc/Mesh.hxx b/CDMATH/mesh/inc/Mesh.hxx index 99ce471..23eb46e 100644 --- a/CDMATH/mesh/inc/Mesh.hxx +++ b/CDMATH/mesh/inc/Mesh.hxx @@ -8,9 +8,12 @@ #ifndef MESH_HXX_ #define MESH_HXX_ +#include + +#include +#include "NormalizedGeometricTypes" + #include "MEDCouplingUMesh.hxx" -#include "MEDCouplingIMesh.hxx" -#include "MEDCouplingFieldDouble.hxx" /** * Mesh class is defined by @@ -25,19 +28,15 @@ namespace MEDCoupling { +class MEDFileMesh; class MEDFileUMesh; class MEDCouplingMesh; -class MEDCouplingIMesh; -class MEDCouplingUMesh; class DataArrayIdType; } namespace ParaMEDMEM { class DataArrayIdType; } -#include -#include "NormalizedGeometricTypes" - class Node; class Cell; class Face; @@ -62,12 +61,18 @@ public: //---------------------------------------------------------------- */ Mesh ( void ) ; + /** + * \brief constructor with data from a medcoupling mesh + * @param medcoupling mesh + */ + Mesh( MEDCoupling::MCAuto mesh ) ; + /** * \brief constructor with data to load a general unstructured mesh * @param filename name of mesh file * @param meshLevel : relative mesh dimension : 0->cells, 1->Faces etc */ - Mesh ( const std::string filename, int meshLevel=0 ) ; + Mesh ( const std::string filename, const std::string & meshName="" , int meshLevel=0) ; /** * \brief constructor with data for a regular 1D grid @@ -110,15 +115,12 @@ public: //---------------------------------------------------------------- */ Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz, int split_to_tetrahedra_policy=-1, std::string meshName="MESH3D_Regular_Cuboid_Grid") ; - Mesh( const MEDCoupling::MEDCouplingIMesh* mesh ) ; - Mesh( const MEDCoupling::MEDCouplingUMesh* mesh ) ; - /** * \brief constructor with data * @param filename : file name of mesh med file * @param meshLevel : relative mesh dimension : 0->cells, 1->Faces etc */ - void readMeshMed( const std::string filename, int meshLevel=0 ) ; + void readMeshMed( const std::string filename, const std::string & meshName="" , int meshLevel=0) ; /** * \brief constructor by copy @@ -144,29 +146,11 @@ public: //---------------------------------------------------------------- int getSpaceDimension( void ) const ; /** - * \brief Mesh dimension + * \brief return Mesh dimension * @return _meshDim */ int getMeshDimension( void ) const ; - /** - * \brief return The nodes in this mesh - * @return _nodes - */ - Node* getNodes ( void ) const ; - - /** - * \brief return The cells in this mesh - * @return _vertices - */ - Cell* getCells ( void ) const ; - - /** - * \brief return the faces in this mesh - * @return _vertices - */ - Face* getFaces ( void ) const ; - /** * \brief return the number of nodes in this mesh * @return _numberOfNodes @@ -239,12 +223,10 @@ public: //---------------------------------------------------------------- double getZMax( void ) const ; - std::vector getDXYZ() const ;// for structured meshes - std::vector getCellGridStructure() const;// for structured meshes /** - * \brief surcharge operator = + * \brief overload operator = * @param mesh : The Mesh object to be copied */ const Mesh& operator= ( const Mesh& mesh ) ; @@ -267,13 +249,13 @@ public: //---------------------------------------------------------------- /** * \brief return the list of face group names - * return _faceGroupNames + * @return _faceGroupNames */ std::vector getNameOfFaceGroups( void ) const ; /** * \brief return the list of node group names - * return _nodeGroupNames + * @return _nodeGroupNames */ std::vector getNameOfNodeGroups( void ) const ; @@ -289,8 +271,8 @@ public: //---------------------------------------------------------------- */ std::vector getNodeGroups( void ) const ; - /* - * Functions to extract boundary nodes and faces Ids + /** + * \brief Functions to extract boundary nodes and faces Ids */ /** * \brief return the list of boundary faces Ids @@ -302,8 +284,8 @@ public: //---------------------------------------------------------------- * @return _boundaryNodeIds */ std::vector< int > getBoundaryNodeIds() const; - /* - * Functions to extract group nodes and faces ids + /** + * \brief Functions to extract group nodes and faces ids */ /** * @return list of face group Ids @@ -322,7 +304,7 @@ public: //---------------------------------------------------------------- /** * \brief write mesh in the MED format */ - void writeMED ( const std::string fileName ) const ; + void writeMED ( const std::string fileName, bool fromScratch = true ) const; void setGroupAtPlan(double value, int direction, double eps, std::string groupName, bool isBoundaryGroup=true) ; @@ -362,7 +344,7 @@ public: //---------------------------------------------------------------- std::vector< std::string > getElementTypesNames() const ; /** - * \brief Compute the minimum value over all cells of the ratio cell perimeter/cell vaolume + * \brief Compute the minimum value over all cells of the ratio cell perimeter/cell volume */ double minRatioVolSurf() const; @@ -374,7 +356,7 @@ public: //---------------------------------------------------------------- /** * \brief Delete the medcoupling mesh to save memory space */ - void deleteMEDCouplingUMesh(); + void deleteMEDCouplingMesh(); /** * \brief Returns true iff an unstructured mesh has been loaded @@ -384,13 +366,14 @@ public: //---------------------------------------------------------------- private: //---------------------------------------------------------------- MEDCoupling::MEDCouplingUMesh* setMesh( void ) ; - void setGroups( const MEDCoupling::MEDFileUMesh* medmesh, MEDCoupling::MEDCouplingUMesh* mu) ;//Read all face and node group + void setFaceGroups( const MEDCoupling::MEDFileUMesh* medmesh, MEDCoupling::MEDCouplingUMesh* mu) ;//Read all face groups + void setNodeGroups( const MEDCoupling::MEDFileMesh* medmesh, MEDCoupling::MEDCouplingUMesh* mu) ;//Read all node groups void addNewFaceGroup( const MEDCoupling::MEDCouplingUMesh *m);//adds one face group in the vectors _faceGroups, _faceGroupNames and _faceGroupIds - /* - * The MEDCoupling mesh + /** + * \brief The MEDCoupling mesh */ - MEDCoupling::MCAuto _mesh; + MEDCoupling::MCAuto _mesh;// This is either a MEDCouplingUMesh or a MEDCouplingStructuredMesh bool _meshNotDeleted; @@ -406,100 +389,106 @@ private: //---------------------------------------------------------------- */ int _meshDim ; - /* - * Structured mesh parameters + /** + * \brief Signal a structured mesh + */ + bool _isStructured; + /** + * \brief Number of cells in each direction (Structured meshes) */ - - bool _isStructured; - - double _xMin; - - double _xMax; - - double _yMin; - - double _yMax; - - double _zMin; - - double _zMax; - std::vector _nxyz; - std::vector _dxyz; - /* - * The nodes in this mesh. + /** + * \brief The nodes in this mesh. */ - Node *_nodes; + std::shared_ptr _nodes; - /* - * The number of nodes in this mesh. + /** + * \brief The number of nodes in this mesh. */ int _numberOfNodes; - /* - * The faces in this mesh. + /** + * \brief The faces in this mesh. */ - Face *_faces; + std::shared_ptr _faces; - /* - * The numbers of faces in this mesh. + /** + * \brief The numbers of faces in this mesh. */ int _numberOfFaces; - /* - * The cells in this mesh. + /** + * \brief The cells in this mesh. */ - Cell *_cells; + std::shared_ptr _cells; - /* - * The number of cells in this mesh. + /** + * \brief The number of cells in this mesh. */ int _numberOfCells; - /* - * The number of edges in this mesh. + /** + * \brief return The nodes in this mesh + * @return _nodes */ - int _numberOfEdges;//Useful to deduce the number of non zero coefficients in the finite element matrix + std::shared_ptr getNodes ( void ) const ; - /* - * The names of face groups. + /** + * \brief return The cells in this mesh + * @return _vertices + */ + std::shared_ptr getCells ( void ) const ; + + /** + * \brief return the faces in this mesh + * @return _vertices + */ + std::shared_ptr getFaces ( void ) const ; + + /** + * \brief The number of edges in this mesh. + */ + int _numberOfEdges;//Useful to deduce the number of non zero coefficients in a finite element matrix + + /** + * \brief The names of face groups. */ std::vector _faceGroupNames; - /* - * The names of node groups. + /** + * \brief The names of node groups. */ std::vector _nodeGroupNames; - /* - * The list of face groups. + /** + * \brief The list of face groups. */ std::vector _faceGroups; - /* - * The list of node groups. + /** + * \brief The list of node groups. */ std::vector _nodeGroups; - /* - * The list of face id in each face groups. + /** + * \brief The list of face id in each face groups. */ std::vector< std::vector > _faceGroupsIds; - /* - * The list of node id in each node groups. + /** + * \brief The list of node id in each node groups. */ std::vector< std::vector > _nodeGroupsIds; - /* - * Elements types (SEG2, TRI3, QUAD4, HEXA6 ...) + /** + * \brief Elements types (SEG2, TRI3, QUAD4, HEXA6 ...) */ std::vector< INTERP_KERNEL::NormalizedCellType > _eltsTypes;//List of cell types contained in the mesh std::vector< std::string > _eltsTypesNames;//List of cell types contained in the mesh std::vector< INTERP_KERNEL::NormalizedCellType > getElementTypes() const; - /* - * Tools to manage periodic boundary conditions in square/cube geometries + /** + * \brief Tools to manage periodic boundary conditions in square/cube geometries */ bool _indexFacePeriodicSet; std::map _indexFacePeriodicMap; diff --git a/CDMATH/mesh/src/Field.cxx b/CDMATH/mesh/src/Field.cxx index 3252115..2997451 100755 --- a/CDMATH/mesh/src/Field.cxx +++ b/CDMATH/mesh/src/Field.cxx @@ -30,7 +30,6 @@ Field::Field( EntityType typeField ) //---------------------------------------------------------------------- { _field=NULL; - _ff=NULL; _typeField=typeField; _numberOfComponents=0; } @@ -47,7 +46,6 @@ 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; @@ -68,7 +66,7 @@ void Field::buildFieldMemoryStructure() }else if(_typeField==NODES) { _field=MEDCouplingFieldDouble::New(ON_NODES); - array->alloc(_mesh.getNumberOfNodes(),_numberOfComponents); + array->alloc(mu->getNumberOfNodes(),_numberOfComponents); _field->setMesh(mu); }else if(_typeField==FACES) { @@ -99,8 +97,6 @@ Field::Field( const std::string filename, EntityType type, int iteration, int order, int meshLevel) { _field = NULL; - _ff=NULL; - _mesh=Mesh(filename + ".med", meshLevel); _typeField=type; _fieldName=fieldName; @@ -108,11 +104,10 @@ Field::Field( const std::string filename, EntityType type, } Field::Field(const std::string meshFileName, EntityType type, const std::vector Vconstant, - const std::string & fieldName, int meshLevel, double time ) + const std::string & fieldName, int meshLevel, double time, std::string meshName ) { _field = NULL; - _ff=NULL; - _mesh=Mesh(meshFileName + ".med", meshLevel); + _mesh=Mesh(meshFileName + ".med", meshName, meshLevel); _typeField=type; _numberOfComponents=Vconstant.size(); _time=time; @@ -130,7 +125,6 @@ Field::Field(const std::string meshFileName, EntityType type, const std::vector< 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(); @@ -149,7 +143,6 @@ Field::Field(const Mesh& M, EntityType type, const Vector Vconstant, const std:: 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(); @@ -172,7 +165,6 @@ 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; @@ -219,7 +211,6 @@ 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(); @@ -297,7 +288,6 @@ 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(); @@ -350,7 +340,6 @@ 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 == "") { @@ -430,7 +419,8 @@ Field::readFieldMed( const std::string & fileNameRadical, _numberOfComponents = _field->getNumberOfComponents() ; _time = _field->getTime(iteration, order); - cout<<"Found field " << fieldName << " in file " << completeFileName <getMesh()->getName()<< endl; + _mesh=Mesh( completeFileName, _field->getMesh()->getName()); } @@ -992,9 +982,6 @@ 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); @@ -1128,7 +1115,7 @@ Field::writeCSV ( const std::string fileName ) const //---------------------------------------------------------------------- void -Field::writeMED ( const std::string fileName, bool fromScratch) const +Field::writeMED ( const std::string fileName, bool fromScratch) //---------------------------------------------------------------------- { string fname=fileName+".med"; @@ -1138,20 +1125,16 @@ Field::writeMED ( const std::string fileName, bool fromScratch) const 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 + else//The mesh has ben deleted, use a MEDFileField1TS 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); + if ( not fromScratch) + { + MEDCoupling::MCAuto ff=MEDFileField1TS::New();//To save the field when the mesh has been deleted + ff->setFieldNoProfileSBT( _field ); + ff->write(fname.c_str(), fromScratch); + } + else + throw CdmathException("Field::writeMED Error !!! The mesh has been deleted, cannot write field from scratch"); } } @@ -1224,7 +1207,7 @@ std::ostream& operator<<(std::ostream& out, const Field& field ) return out; } -void Field::deleteMEDCouplingUMesh() +void Field::deleteMEDCouplingMesh() { - return _mesh.deleteMEDCouplingUMesh(); + return _mesh.deleteMEDCouplingMesh(); } diff --git a/CDMATH/mesh/src/Mesh.cxx b/CDMATH/mesh/src/Mesh.cxx index 44ca040..8976721 100644 --- a/CDMATH/mesh/src/Mesh.cxx +++ b/CDMATH/mesh/src/Mesh.cxx @@ -5,22 +5,26 @@ * Authors: CDMATH */ +#include +#include +#include +#include +#include +#include + #include "Mesh.hxx" #include "Node.hxx" #include "Cell.hxx" #include "Face.hxx" -#include "CdmathException.hxx" +#include "MEDCouplingIMesh.hxx" +#include "MEDCouplingUMesh.hxx" +#include "MEDCouplingFieldDouble.hxx" + #include "MEDFileMesh.hxx" #include "MEDLoader.hxx" -#include -#include -#include -#include -#include -#include -#include +#include "CdmathException.hxx" using namespace MEDCoupling; using namespace std; @@ -41,14 +45,7 @@ Mesh::Mesh( void ) _numberOfCells = 0; _numberOfEdges = 0; _isStructured=false; - _xMin=0.; - _xMax=0.; - _yMin=0.; - _yMax=0.; - _zMin=0.; - _zMax=0.; _nxyz.resize(0); - _dxyz.resize(0.); _boundaryMesh=NULL; _boundaryFaceIds.resize(0); _boundaryNodeIds.resize(0); @@ -67,10 +64,6 @@ Mesh::Mesh( void ) Mesh::~Mesh( void ) //---------------------------------------------------------------------- { - delete [] _cells; - delete [] _nodes; - delete [] _faces; - //for(int i=0; i< _faceGroups.size(); i++) // _faceGroups[i]->decrRef(); // _nodeGroups[i]->decrRef(); @@ -84,97 +77,34 @@ Mesh::getName( void ) const return _name; } -Mesh::Mesh( const MEDCoupling::MEDCouplingIMesh* mesh ) -{ - _spaceDim=mesh->getSpaceDimension(); - _meshDim=mesh->getMeshDimension(); - _dxyz=mesh->getDXYZ(); - _nxyz=mesh->getCellGridStructure(); - double* Box0=new double[2*_spaceDim]; - mesh->getBoundingBox(Box0); - _name=mesh->getName(); - _epsilon=1e-6; - _indexFacePeriodicSet=false; - - _xMin=Box0[0]; - _xMax=Box0[1]; - if (_spaceDim>=2) - { - _yMin=Box0[2]; - _yMax=Box0[3]; - } - if (_spaceDim>=3) - { - _zMin=Box0[4]; - _zMax=Box0[5]; - } - - double *originPtr = new double[_spaceDim]; - double *dxyzPtr = new double[_spaceDim]; - mcIdType *nodeStrctPtr = new mcIdType[_spaceDim]; - - for(int i=0;i<_spaceDim;i++) - { - originPtr[i]=Box0[2*i]; - nodeStrctPtr[i]=_nxyz[i]+1; - dxyzPtr[i]=_dxyz[i]; - } - _mesh=MEDCouplingIMesh::New(_name, - _spaceDim, - nodeStrctPtr, - nodeStrctPtr+_spaceDim, - originPtr, - originPtr+_spaceDim, - dxyzPtr, - dxyzPtr+_spaceDim); - _meshNotDeleted=true; - _isStructured=true; - - delete [] originPtr; - delete [] dxyzPtr; - delete [] nodeStrctPtr; - delete [] Box0 ; - - setMesh(); -} - -Mesh::Mesh( const MEDCoupling::MEDCouplingUMesh* mesh ) +Mesh::Mesh( MEDCoupling::MCAuto mesh ) { _spaceDim=mesh->getSpaceDimension(); _meshDim=mesh->getMeshDimension(); - double* Box0=new double[2*_spaceDim]; - mesh->getBoundingBox(Box0); _name=mesh->getName(); _epsilon=1e-6; _indexFacePeriodicSet=false; + _meshNotDeleted=true; - _xMin=Box0[0]; - _xMax=Box0[1]; - if (_spaceDim>=2) - { - _yMin=Box0[2]; - _yMax=Box0[3]; - } - if (_spaceDim>=3) - { - _zMin=Box0[4]; - _zMax=Box0[5]; - } + _mesh= mesh->clone(false);//Clone because you will need to buildUnstructured. No deep copy : it is assumed node coordinates and cell connectivity will not change - _mesh=mesh->deepCopy(); - _mesh=mesh->buildUnstructured(); - _meshNotDeleted=true; - _isStructured=false; - delete [] Box0 ; + MEDCoupling::MEDCouplingStructuredMesh* structuredMesh = dynamic_cast (_mesh.retn()); + if(structuredMesh) + { + _isStructured=true; + _nxyz=structuredMesh->getCellGridStructure(); + } + else + _isStructured=false; setMesh(); } //---------------------------------------------------------------------- -Mesh::Mesh( const std::string filename, int meshLevel ) +Mesh::Mesh( const std::string filename, const std::string & meshName, int meshLevel) //---------------------------------------------------------------------- { - readMeshMed(filename, meshLevel); + readMeshMed(filename, meshName, meshLevel); } //---------------------------------------------------------------------- @@ -185,17 +115,9 @@ Mesh::Mesh( const Mesh& mesh ) _meshDim = mesh.getMeshDimension() ; _name=mesh.getName(); _epsilon=mesh.getComparisonEpsilon(); - _xMax=mesh.getXMax(); - _yMin=mesh.getYMin(); - _yMax=mesh.getYMax(); - _zMin=mesh.getZMin(); - _zMax=mesh.getZMax(); _isStructured=mesh.isStructured(); if(_isStructured) - { _nxyz = mesh.getCellGridStructure() ; - _dxyz=mesh.getDXYZ(); - } _numberOfNodes = mesh.getNumberOfNodes(); _numberOfFaces = mesh.getNumberOfFaces(); _numberOfCells = mesh.getNumberOfCells(); @@ -206,23 +128,10 @@ Mesh::Mesh( const Mesh& mesh ) _nodeGroupNames = mesh.getNameOfNodeGroups() ; _nodeGroups = mesh.getNodeGroups() ; - _nodes = new Node[_numberOfNodes] ; - _faces = new Face[_numberOfFaces] ; - _cells = new Cell[_numberOfCells] ; + _nodes = mesh.getNodes() ; + _faces = mesh.getFaces() ; + _cells = mesh.getCells() ; - //memcpy(_nodes,mesh.getNodes(),sizeof(*mesh.getNodes())) ; - //memcpy(_cells,mesh.getCells(),sizeof(*mesh.getCells())) ; - //memcpy(_faces,mesh.getFaces(),sizeof(*mesh.getFaces())) ; - - for (int i=0;i<_numberOfNodes;i++) - _nodes[i]=mesh.getNode(i); - - for (int i=0;i<_numberOfFaces;i++) - _faces[i]=mesh.getFace(i); - - for (int i=0;i<_numberOfCells;i++) - _cells[i]=mesh.getCell(i); - _indexFacePeriodicSet= mesh.isIndexFacePeriodicSet(); if(_indexFacePeriodicSet) _indexFacePeriodicMap=mesh.getIndexFacePeriodic(); @@ -233,17 +142,23 @@ Mesh::Mesh( const Mesh& mesh ) _eltsTypes=mesh.getElementTypes(); _eltsTypesNames=mesh.getElementTypesNames(); - MCAuto m1=mesh.getMEDCouplingMesh()->deepCopy(); + MCAuto m1=mesh.getMEDCouplingMesh()->clone(false);//Clone because you will need to buildUnstructured. No deep copy : it is assumed node coordinates and cell connectivity will not change + _mesh=m1; _meshNotDeleted=mesh.meshNotDeleted(); } //---------------------------------------------------------------------- void -Mesh::readMeshMed( const std::string filename, const int meshLevel) +Mesh::readMeshMed( const std::string filename, const std::string & meshName, int meshLevel) //---------------------------------------------------------------------- { - MEDFileUMesh *m=MEDFileUMesh::New(filename.c_str());//reads the first mesh encountered in the file, otherwise call New (const char *fileName, const char *mName, int dt=-1, int it=-1) + MEDFileMesh *m; + if( meshName == "" ) + m=MEDFileMesh::New(filename.c_str());//reads the first mesh encountered in the file, otherwise call New (const char *fileName, const char *mName, int dt=-1, int it=-1) + else + m=MEDFileMesh::New(filename.c_str(), meshName.c_str());//seeks the mesh named meshName in the file + _mesh=m->getMeshAtLevel(meshLevel); _mesh->checkConsistencyLight(); _mesh->setName(_mesh->getName()); @@ -252,37 +167,22 @@ Mesh::readMeshMed( const std::string filename, const int meshLevel) _name=_mesh->getName(); _epsilon=1e-6; _indexFacePeriodicSet=false; - MEDCoupling::MEDCouplingIMesh* structuredMesh = dynamic_cast (_mesh.retn()); + _meshNotDeleted=true; + + MEDCoupling::MEDCouplingStructuredMesh* structuredMesh = dynamic_cast (_mesh.retn()); if(structuredMesh) { _isStructured=true; - _meshNotDeleted=false; - _dxyz=structuredMesh->getDXYZ(); _nxyz=structuredMesh->getCellGridStructure(); - double* Box0=new double[2*_spaceDim]; - structuredMesh->getBoundingBox(Box0); - - _xMin=Box0[0]; - _xMax=Box0[1]; - if (_spaceDim>=2) - { - _yMin=Box0[2]; - _yMax=Box0[3]; - } - if (_spaceDim>=3) - { - _zMin=Box0[4]; - _zMax=Box0[5]; - } } else - { _isStructured=false; - _meshNotDeleted=true; - } MEDCouplingUMesh* mu = setMesh(); - setGroups(m, mu); + setNodeGroups(m, mu);//Works for both cartesan and unstructured meshes + MEDFileUMesh *umedfile=dynamic_cast< MEDFileUMesh * > (m); + if(umedfile) + setFaceGroups(umedfile, mu);//Works only for unstructured meshes cout< points, string meshName) : nx < 2, vector should contain at least two values"); + int i=0; + while( ipoints[i] ) + i++; + if( i!=nx-1 ) { - std::vector< std::string >::iterator it = std::find(_faceGroupNames.begin(), _faceGroupNames.end(), groupName); - if(it == _faceGroupNames.end())//No group named groupName - { - _faceGroupNames.insert(_faceGroupNames.end(),groupName); - _faceGroupsIds.insert( _faceGroupsIds.end(),faceIds); - _faceGroups.insert( _faceGroups.end(), NULL);//No mesh created. Create one ? - } - else - { - std::vector< int > faceGroupIds = _faceGroupsIds[it-_faceGroupNames.begin()]; - faceGroupIds.insert( faceGroupIds.end(), faceIds.begin(), faceIds.end()); - /* Detect and erase duplicates face ids */ - sort( faceGroupIds.begin(), faceGroupIds.end() ); - faceGroupIds.erase( unique( faceGroupIds.begin(), faceGroupIds.end() ), faceGroupIds.end() ); - _faceGroupsIds[it-_faceGroupNames.begin()] = faceGroupIds; - } - } + //cout<< points << endl; + throw CdmathException("Mesh::Mesh( vector points, string meshName) : vector values should be sorted"); + } + + _spaceDim = 1 ; + _meshDim = 1 ; + _name=meshName; + _epsilon=1e-6; + _indexFacePeriodicSet=false; + + MEDCouplingUMesh * mesh1d = MEDCouplingUMesh::New(meshName, 1); + mesh1d->allocateCells(nx - 1); + double * coords = new double[nx]; + mcIdType nodal_con[2]; + coords[0]=points[0]; + for(int i=0; iinsertNextCell(INTERP_KERNEL::NORM_SEG2, 2, nodal_con); + coords[i+1]=points[i + 1]; + } + mesh1d->finishInsertingCells(); + + DataArrayDouble * coords_arr = DataArrayDouble::New(); + coords_arr->alloc(nx,1); + std::copy(coords,coords+nx,coords_arr->getPointer()); + mesh1d->setCoords(coords_arr); + + delete [] coords; + coords_arr->decrRef(); + + _mesh=mesh1d;//To enable writeMED. Because we declared the mesh as unstructured, we decide to build the unstructured data (not mandatory) + _meshNotDeleted=true; + _isStructured = false; + + setMesh(); } -void -Mesh::setGroupAtNodeByCoords(double x, double y, double z, double eps, std::string groupName, bool isBoundaryGroup) +//---------------------------------------------------------------------- +Mesh::Mesh( double xmin, double xmax, int nx, std::string meshName ) +//---------------------------------------------------------------------- { - std::vector< int > nodeIds(0); - double NX, NY, NZ; - Node Ni; - - /* Construction of the node group */ - if(isBoundaryGroup) - for(int i=0; i<_boundaryNodeIds.size(); i++) - { - Ni=_nodes[_boundaryNodeIds[i]]; - NX=Ni.x(); - NY=Ni.y(); - NZ=Ni.z(); - if (abs(NX-x)=xmax) + throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx) : xmin >= xmax"); - if (nodeIds.size()>0) - { - std::vector< std::string >::iterator it = std::find(_nodeGroupNames.begin(), _nodeGroupNames.end(), groupName); - if(it == _nodeGroupNames.end())//No group named groupName - { - _nodeGroupNames.insert(_nodeGroupNames.end(),groupName); - _nodeGroupsIds.insert( _nodeGroupsIds.end(),nodeIds); - _nodeGroups.insert( _nodeGroups.end(), NULL);//No mesh created. Create one ? - } - else - { - std::vector< int > nodeGroupIds = _nodeGroupsIds[it-_nodeGroupNames.begin()]; - nodeGroupIds.insert( nodeGroupIds.end(), nodeIds.begin(), nodeIds.end()); - /* Detect and erase duplicates node ids */ - sort( nodeGroupIds.begin(), nodeGroupIds.end() ); - nodeGroupIds.erase( unique( nodeGroupIds.begin(), nodeGroupIds.end() ), nodeGroupIds.end() ); - _nodeGroupsIds[it-_nodeGroupNames.begin()] = nodeGroupIds; - } - } + double dx = (xmax - xmin)/nx ; + + _spaceDim = 1 ; + _meshDim = 1 ; + _name=meshName; + _epsilon=1e-6; + _indexFacePeriodicSet=false; + + _nxyz.resize(_spaceDim); + _nxyz[0]=nx; + + double originPtr[_spaceDim]; + double dxyzPtr[_spaceDim]; + mcIdType nodeStrctPtr[_spaceDim]; + + originPtr[0]=xmin; + nodeStrctPtr[0]=nx+1; + dxyzPtr[0]=dx; + + _mesh=MEDCouplingIMesh::New(meshName, + _spaceDim, + nodeStrctPtr, + nodeStrctPtr+_spaceDim, + originPtr, + originPtr+_spaceDim, + dxyzPtr, + dxyzPtr+_spaceDim); + _meshNotDeleted=true; + _isStructured = true; + + setMesh(); } -void -Mesh::setGroupAtPlan(double value, int direction, double eps, std::string groupName, bool isBoundaryGroup) +//---------------------------------------------------------------------- +Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, int split_to_triangles_policy, std::string meshName) +//---------------------------------------------------------------------- { - std::vector< int > faceIds(0), nodeIds(0); - double cord; - - /* Construction of the face group */ - if(isBoundaryGroup) - for(int i=0; i<_boundaryFaceIds.size(); i++) - { - cord=_faces[_boundaryFaceIds[i]].getBarryCenter()[direction]; - if (abs(cord-value)=xmax) + throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : xmin >= xmax"); + if(ymin>=ymax) + throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : ymin >= ymax"); - /* Construction of the node group */ - if(isBoundaryGroup) - for(int i=0; i<_boundaryNodeIds.size(); i++) - { - cord=_nodes[_boundaryNodeIds[i]].getPoint()[direction]; - if (abs(cord-value)0) - { - std::vector< std::string >::iterator it = std::find(_faceGroupNames.begin(), _faceGroupNames.end(), groupName); - if(it == _faceGroupNames.end())//No group named groupName - { - _faceGroupNames.insert(_faceGroupNames.end(),groupName); - _faceGroupsIds.insert( _faceGroupsIds.end(),faceIds); - _faceGroups.insert( _faceGroups.end(), NULL);//No mesh created. Create one ? - } - else - { - std::vector< int > faceGroupIds = _faceGroupsIds[it-_faceGroupNames.begin()]; - faceGroupIds.insert( faceGroupIds.end(), faceIds.begin(), faceIds.end()); - /* Detect and erase duplicates face ids */ - sort( faceGroupIds.begin(), faceGroupIds.end() ); - faceGroupIds.erase( unique( faceGroupIds.begin(), faceGroupIds.end() ), faceGroupIds.end() ); - _faceGroupsIds[it-_faceGroupNames.begin()] = faceGroupIds; - } - } - if (nodeIds.size()>0) - { - std::vector< std::string >::iterator it = std::find(_nodeGroupNames.begin(), _nodeGroupNames.end(), groupName); - if(it == _nodeGroupNames.end())//No group named groupName - { - _nodeGroupNames.insert(_nodeGroupNames.end(),groupName); - _nodeGroupsIds.insert( _nodeGroupsIds.end(),nodeIds); - _nodeGroups.insert( _nodeGroups.end(), NULL);//No mesh created. Create one ? - } - else - { - std::vector< int > nodeGroupIds = _nodeGroupsIds[it-_nodeGroupNames.begin()]; - nodeGroupIds.insert( nodeGroupIds.end(), nodeIds.begin(), nodeIds.end()); - /* Detect and erase duplicates node ids */ - sort( nodeGroupIds.begin(), nodeGroupIds.end() ); - nodeGroupIds.erase( unique( nodeGroupIds.begin(), nodeGroupIds.end() ), nodeGroupIds.end() ); - _nodeGroupsIds[it-_nodeGroupNames.begin()] = nodeGroupIds; - } - } -} + _spaceDim = 2 ; + _meshDim = 2 ; + _name=meshName; + _epsilon=1e-6; + _indexFacePeriodicSet=false; + _nxyz.resize(_spaceDim); + _nxyz[0]=nx; + _nxyz[1]=ny; -void -Mesh::setBoundaryNodesFromFaces() -{ - for (int iface=0;iface<_boundaryFaceIds.size();iface++) - { - std::vector< int > nodesID= _faces[_boundaryFaceIds[iface]].getNodesId(); - int nbNodes = _faces[_boundaryFaceIds[iface]].getNumberOfNodes(); - for(int inode=0 ; inode::const_iterator it = std::find(_boundaryNodeIds.begin(),_boundaryNodeIds.end(),nodesID[inode]); - if( it != _boundaryNodeIds.end() ) - _boundaryNodeIds.push_back(nodesID[inode]); - } - } -} + double originPtr[_spaceDim]; + double dxyzPtr[_spaceDim]; + mcIdType nodeStrctPtr[_spaceDim]; -std::map -Mesh::getIndexFacePeriodic( void ) const -{ - return _indexFacePeriodicMap; -} + originPtr[0]=xmin; + originPtr[1]=ymin; + nodeStrctPtr[0]=nx+1; + nodeStrctPtr[1]=ny+1; + dxyzPtr[0]=dx; + dxyzPtr[1]=dy; -void -Mesh::setPeriodicFaces(bool check_groups, bool use_central_inversion) -{ - if(_indexFacePeriodicSet) - return; - - for (int indexFace=0;indexFace<_boundaryFaceIds.size() ; indexFace++) - { - Face my_face=_faces[_boundaryFaceIds[indexFace]]; - int iface_perio=-1; - if(_meshDim==1) - { - for (int iface=0;iface<_boundaryFaceIds.size() ; iface++) - if(iface!=indexFace) - { - iface_perio=_boundaryFaceIds[iface]; - break; - } - } - else if(_meshDim==2) + _mesh=MEDCouplingIMesh::New(meshName, + _spaceDim, + nodeStrctPtr, + nodeStrctPtr+_spaceDim, + originPtr, + originPtr+_spaceDim, + dxyzPtr, + dxyzPtr+_spaceDim); + _meshNotDeleted=true; + _isStructured = true; + + if(split_to_triangles_policy==0 || split_to_triangles_policy==1) { - double x=my_face.x(); - double y=my_face.y(); - - for (int iface=0;iface<_boundaryFaceIds.size() ; iface++) - { - Face face_i=_faces[_boundaryFaceIds[iface]]; - double xi=face_i.x(); - double yi=face_i.y(); - if ( (abs(y-yi)<_epsilon || abs(x-xi)<_epsilon )// Case of a square geometry - && ( !check_groups || my_face.getGroupName()!=face_i.getGroupName()) //In case groups need to be checked - && ( !use_central_inversion || abs(y+yi) + abs(x+xi)<_epsilon ) // Case of a central inversion - && fabs(my_face.getMeasure() - face_i.getMeasure())<_epsilon - && fabs(my_face.getXN() + face_i.getXN())<_epsilon - && fabs(my_face.getYN() + face_i.getYN())<_epsilon - && fabs(my_face.getZN() + face_i.getZN())<_epsilon ) - { - iface_perio=_boundaryFaceIds[iface]; - break; - } - } + _mesh=_mesh->buildUnstructured();//simplexize is not available for structured meshes + _mesh->simplexize(split_to_triangles_policy); + _isStructured = false; } - else if(_meshDim==3) + else if (split_to_triangles_policy != -1) { - double x=my_face.x(); - double y=my_face.y(); - double z=my_face.z(); - - for (int iface=0;iface<_boundaryFaceIds.size() ; iface++) - { - Face face_i=_faces[_boundaryFaceIds[iface]]; - double xi=face_i.x(); - double yi=face_i.y(); - double zi=face_i.z(); - if ( ((abs(y-yi)<_epsilon && abs(x-xi)<_epsilon) || (abs(x-xi)<_epsilon && abs(z-zi)<_epsilon) || (abs(y-yi)<_epsilon && abs(z-zi)<_epsilon))// Case of a cube geometry - && ( !check_groups || my_face.getGroupName()!=face_i.getGroupName()) //In case groups need to be checked - && ( !use_central_inversion || abs(y+yi) + abs(x+xi) + abs(z+zi)<_epsilon )// Case of a central inversion - && fabs(my_face.getMeasure() - face_i.getMeasure())<_epsilon - && fabs(my_face.getXN() + face_i.getXN())<_epsilon - && fabs(my_face.getYN() + face_i.getYN())<_epsilon - && fabs(my_face.getZN() + face_i.getZN())<_epsilon ) - { - iface_perio=_boundaryFaceIds[iface]; - break; - } - } + cout<< "split_to_triangles_policy = "<< split_to_triangles_policy << endl; + throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : Unknown splitting policy"); } - else - throw CdmathException("Mesh::setPeriodicFaces: Mesh dimension should be 1, 2 or 3"); - - if (iface_perio==-1) - throw CdmathException("Mesh::setPeriodicFaces: periodic face not found, iface_perio==-1 " ); - else - _indexFacePeriodicMap[_boundaryFaceIds[indexFace]]=iface_perio; - } - _indexFacePeriodicSet=true; + + setMesh(); } -int -Mesh::getIndexFacePeriodic(int indexFace, bool check_groups, bool use_central_inversion) +//---------------------------------------------------------------------- +Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz, int split_to_tetrahedra_policy, std::string meshName) +//---------------------------------------------------------------------- { - if (!_faces[indexFace].isBorder()) - { - cout<<"Pb with indexFace= "<::const_iterator it = _indexFacePeriodicMap.find(indexFace); - if( it != _indexFacePeriodicMap.end() ) - return it->second; - else - { - cout<<"Pb with indexFace= "<=xmax) + throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : xmin >= xmax"); + if(ymin>=ymax) + throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : ymin >= ymax"); + if(zmin>=zmax) + throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : zmin >= zmax"); -bool -Mesh::isBorderNode(int nodeid) const -{ - return _nodes[nodeid].isBorder(); -} + _spaceDim = 3; + _meshDim = 3; + _name=meshName; + _epsilon=1e-6; -bool -Mesh::isBorderFace(int faceid) const -{ - return _faces[faceid].isBorder(); -} + double dx = (xmax - xmin)/nx ; + double dy = (ymax - ymin)/ny ; + double dz = (zmax - zmin)/nz ; -std::vector< int > -Mesh::getBoundaryFaceIds() const -{ - return _boundaryFaceIds; -} + _nxyz.resize(_spaceDim); + _nxyz[0]=nx; + _nxyz[1]=ny; + _nxyz[2]=nz; -std::vector< int > -Mesh::getBoundaryNodeIds() const -{ - return _boundaryNodeIds; -} + double originPtr[_spaceDim]; + double dxyzPtr[_spaceDim]; + mcIdType nodeStrctPtr[_spaceDim]; -bool -Mesh::isTriangular() const -{ - return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_TRI3; -} -bool -Mesh::isTetrahedral() const -{ - return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_TETRA4; -} -bool -Mesh::isQuadrangular() const -{ - return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_QUAD4; -} -bool -Mesh::isHexahedral() const -{ - return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_HEXA8; -} -bool -Mesh::isStructured() const -{ - return _isStructured; -} + originPtr[0]=xmin; + originPtr[1]=ymin; + originPtr[2]=zmin; + nodeStrctPtr[0]=nx+1; + nodeStrctPtr[1]=ny+1; + nodeStrctPtr[2]=nz+1; + dxyzPtr[0]=dx; + dxyzPtr[1]=dy; + dxyzPtr[2]=dz; -std::vector< INTERP_KERNEL::NormalizedCellType > -Mesh::getElementTypes() const -{ - return _eltsTypes; -} + _mesh=MEDCouplingIMesh::New(meshName, + _spaceDim, + nodeStrctPtr, + nodeStrctPtr+_spaceDim, + originPtr, + originPtr+_spaceDim, + dxyzPtr, + dxyzPtr+_spaceDim); + _meshNotDeleted=true; + _isStructured = true; -std::vector< string > -Mesh::getElementTypesNames() const -{ - std::vector< string > result(0); - for(int i=0; i< _eltsTypes.size(); i++) - { - if( _eltsTypes[i]==INTERP_KERNEL::NORM_POINT1) - result.push_back("Points"); - else if( _eltsTypes[i]==INTERP_KERNEL::NORM_SEG2) - result.push_back("Segments"); - else if( _eltsTypes[i]==INTERP_KERNEL::NORM_TRI3) - result.push_back("Triangles"); - else if( _eltsTypes[i]==INTERP_KERNEL::NORM_QUAD4) - result.push_back("Quadrangles"); - else if( _eltsTypes[i]==INTERP_KERNEL::NORM_POLYGON) - result.push_back("Polygons"); - else if( _eltsTypes[i]==INTERP_KERNEL::NORM_TETRA4) - result.push_back("Tetrahedra"); - else if( _eltsTypes[i]==INTERP_KERNEL::NORM_PYRA5) - result.push_back("Pyramids"); - else if( _eltsTypes[i]==INTERP_KERNEL::NORM_PENTA6) - result.push_back("Pentahedra"); - else if( _eltsTypes[i]==INTERP_KERNEL::NORM_HEXA8) - result.push_back("Hexahedra"); - else if( _eltsTypes[i]==INTERP_KERNEL::NORM_POLYHED) - result.push_back("Polyhedrons"); - else - { - cout<< "Mesh " + _name + " contains an element of type " <buildUnstructured();//simplexize is not available for structured meshes + _mesh->simplexize(INTERP_KERNEL::PLANAR_FACE_5); + _isStructured = false; } - } - return result; -} - -void -Mesh::setGroups( const MEDFileUMesh* medmesh, MEDCouplingUMesh* mu) -{ - int nbCellsSubMesh, nbNodesSubMesh; - bool foundFace, foundNode; - - /* Searching for face groups */ - vector faceGroups=medmesh->getGroupsNames() ; - - for (unsigned int i=0;i nonEmptyGrp(medmesh->getGrpNonEmptyLevels(groupName)); - //We check if the group has a relative dimension equal to -1 - //before call to the function getGroup(-1,groupName.c_str()) - vector::iterator it = find(nonEmptyGrp.begin(), nonEmptyGrp.end(), -1); - if (it != nonEmptyGrp.end()) - { - MEDCouplingUMesh *m=medmesh->getGroup(-1,groupName.c_str()); - m->unPolyze(); - m->sortCellsInMEDFileFrmt( ); - nbCellsSubMesh=m->getNumberOfCells(); - - _faceGroups.insert(_faceGroups.end(),m);//Vector of group meshes - _faceGroupNames.insert(_faceGroupNames.end(),groupName);//Vector of group names - _faceGroupsIds.insert(_faceGroupsIds.end(),std::vector(nbCellsSubMesh));//Vector of group face Ids. The filling of the face ids takes place below. - - DataArrayDouble *baryCell = m->computeIsoBarycenterOfNodesPerCell();//computeCellCenterOfMass() ; - const double *coorBary=baryCell->getConstPointer(); - /* Face identification */ - for (int ic(0), k(0); ic coorBaryXyz(3,0); - for (int d=0; d<_spaceDim; d++) - coorBaryXyz[d] = coorBary[k+d]; - Point p1(coorBaryXyz[0],coorBaryXyz[1],coorBaryXyz[2]) ; - - foundFace=false; - for (int iface=0;iface<_numberOfFaces;iface++ ) - { - Point p2=_faces[iface].getBarryCenter(); - if(p1.distance(p2)<_epsilon) - { - _faces[iface].setGroupName(groupName); - _faceGroupsIds[_faceGroupsIds.size()-1][ic]=iface; - foundFace=true; - break; - } - } - if (not foundFace) - throw CdmathException("No face found for group " + groupName ); - } - baryCell->decrRef(); - //m->decrRef(); - } - } - - /* Searching for node groups */ - vector nodeGroups=medmesh->getGroupsOnSpecifiedLev(1) ; - - for (unsigned int i=0;igetNodeGroupArr( groupName ); - const mcIdType *nodeids=nodeGroup->getConstPointer(); - - if(nodeids!=NULL) - { - _nodeGroups.insert(_nodeGroups.end(),nodeGroup); - _nodeGroupNames.insert(_nodeGroupNames.end(),groupName); - - nbNodesSubMesh=nodeGroup->getNumberOfTuples();//nodeGroup->getNbOfElems(); - - DataArrayDouble *coo = mu->getCoords() ; - const double *cood=coo->getConstPointer(); - - _nodeGroupsIds.insert(_nodeGroupsIds.end(),std::vector(nbNodesSubMesh));//Vector of boundary faces - /* Node identification */ - for (int ic(0); ic coorP(3,0); - for (int d=0; d<_spaceDim; d++) - coorP[d] = cood[nodeids[ic]*_spaceDim+d]; - Point p1(coorP[0],coorP[1],coorP[2]) ; - - foundNode=false; - for (int inode=0;inode<_numberOfNodes;inode++ ) - { - Point p2=_nodes[inode].getPoint(); - if(p1.distance(p2)<_epsilon) - { - _nodes[inode].setGroupName(groupName); - _nodeGroupsIds[_nodeGroupsIds.size()-1][ic]=inode; - foundNode=true; - break; - } - } - if (not foundNode) - throw CdmathException("No node found for group " + groupName ); - } - } - } + else if( split_to_tetrahedra_policy == 1 ) + { + _mesh=_mesh->buildUnstructured();//simplexize is not available for structured meshes + _mesh->simplexize(INTERP_KERNEL::PLANAR_FACE_6); + _isStructured = false; + } + else if ( split_to_tetrahedra_policy != -1 ) + { + cout<< "split_to_tetrahedra_policy = "<< split_to_tetrahedra_policy << endl; + throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : splitting policy value should be 0 or 1"); + } + + setMesh(); } //---------------------------------------------------------------------- @@ -888,13 +496,13 @@ Mesh::setMesh( void ) const mcIdType *tmpNEI=nodalI->getConstPointer(); _numberOfCells = mu->getNumberOfCells() ; - _cells = new Cell[_numberOfCells] ; + _cells = std::shared_ptr(new Cell[_numberOfCells], std::default_delete()) ; _numberOfNodes = mu->getNumberOfNodes() ;//This number may include isolated nodes that will not be loaded. The number will be updated during nodes constructions - _nodes = new Node[_numberOfNodes] ;//This array may be resized if isolated nodes are found + _nodes = std::shared_ptr(new Node[_numberOfNodes], std::default_delete()) ;//This array may be resized if isolated nodes are found _numberOfFaces = mu2->getNumberOfCells(); - _faces = new Face[_numberOfFaces] ; + _faces = std::shared_ptr(new Face[_numberOfFaces], std::default_delete()) ; _indexFacePeriodicSet=false; @@ -966,7 +574,7 @@ Mesh::setMesh( void ) ci.addFaceId(el,work[el]) ; xn = - xn; yn=-yn; zn=-zn; } - _cells[id] = ci ; + _cells.get()[id] = ci ; } for( int id(0); id<_numberOfFaces; id++ ) @@ -986,7 +594,7 @@ Mesh::setMesh( void ) assert( nbCells>0);//To make sure our face is not located on an isolated node Face fi( nbNodes, nbCells, 1.0, p, 1., 0., 0. ) ; - for(int node_id=0; node_id1 ) + { //Set face boundary group + _boundaryMesh = mu->computeSkin(); + _faceGroups.push_back(_boundaryMesh); + } + else + _faceGroups.push_back(NULL); + _nodeGroups.push_back(NULL); + + desc->decrRef(); + descI->decrRef(); + revDesc->decrRef(); + revDescI->decrRef(); + mu2->decrRef(); + baryCell->decrRef(); + fields->decrRef(); + revNode->decrRef(); + revNodeI->decrRef(); + revCell->decrRef(); + revCellI->decrRef(); + + if (_meshDim == 3) + { + revNode2->decrRef(); + revNodeI2->decrRef(); + desc2->decrRef(); + descI2->decrRef(); + revDesc2->decrRef(); + revDescI2->decrRef(); + mu3->decrRef(); + } + + return mu; +} + +void +Mesh::setGroupAtFaceByCoords(double x, double y, double z, double eps, std::string groupName, bool isBoundaryGroup) +{ + std::vector< int > faceIds(0); + double FX, FY, FZ; + Face Fi; + + /* Construction of the face group */ + if(isBoundaryGroup) + for(int i=0; i<_boundaryFaceIds.size(); i++) + { + Fi=_faces.get()[_boundaryFaceIds[i]]; + FX=Fi.x(); + FY=Fi.y(); + FZ=Fi.z(); + if (abs(FX-x)0) + { + std::vector< std::string >::iterator it = std::find(_faceGroupNames.begin(), _faceGroupNames.end(), groupName); + if(it == _faceGroupNames.end())//No group named groupName + { + _faceGroupNames.insert(_faceGroupNames.end(),groupName); + _faceGroupsIds.insert( _faceGroupsIds.end(),faceIds); + _faceGroups.insert( _faceGroups.end(), NULL);//No mesh created. Create one ? + } + else + { + std::vector< int > faceGroupIds = _faceGroupsIds[it-_faceGroupNames.begin()]; + faceGroupIds.insert( faceGroupIds.end(), faceIds.begin(), faceIds.end()); + /* Detect and erase duplicates face ids */ + sort( faceGroupIds.begin(), faceGroupIds.end() ); + faceGroupIds.erase( unique( faceGroupIds.begin(), faceGroupIds.end() ), faceGroupIds.end() ); + _faceGroupsIds[it-_faceGroupNames.begin()] = faceGroupIds; + } + } +} + +void +Mesh::setGroupAtNodeByCoords(double x, double y, double z, double eps, std::string groupName, bool isBoundaryGroup) +{ + std::vector< int > nodeIds(0); + double NX, NY, NZ; + Node Ni; + + /* Construction of the node group */ + if(isBoundaryGroup) + for(int i=0; i<_boundaryNodeIds.size(); i++) + { + Ni=_nodes.get()[_boundaryNodeIds[i]]; + NX=Ni.x(); + NY=Ni.y(); + NZ=Ni.z(); + if (abs(NX-x)0) + { + std::vector< std::string >::iterator it = std::find(_nodeGroupNames.begin(), _nodeGroupNames.end(), groupName); + if(it == _nodeGroupNames.end())//No group named groupName + { + _nodeGroupNames.insert(_nodeGroupNames.end(),groupName); + _nodeGroupsIds.insert( _nodeGroupsIds.end(),nodeIds); + _nodeGroups.insert( _nodeGroups.end(), NULL);//No mesh created. Create one ? + } + else + { + std::vector< int > nodeGroupIds = _nodeGroupsIds[it-_nodeGroupNames.begin()]; + nodeGroupIds.insert( nodeGroupIds.end(), nodeIds.begin(), nodeIds.end()); + /* Detect and erase duplicates node ids */ + sort( nodeGroupIds.begin(), nodeGroupIds.end() ); + nodeGroupIds.erase( unique( nodeGroupIds.begin(), nodeGroupIds.end() ), nodeGroupIds.end() ); + _nodeGroupsIds[it-_nodeGroupNames.begin()] = nodeGroupIds; + } + } +} + +void +Mesh::setGroupAtPlan(double value, int direction, double eps, std::string groupName, bool isBoundaryGroup) +{ + std::vector< int > faceIds(0), nodeIds(0); + double cord; + + /* Construction of the face group */ + if(isBoundaryGroup) + for(int i=0; i<_boundaryFaceIds.size(); i++) + { + cord=_faces.get()[_boundaryFaceIds[i]].getBarryCenter()[direction]; + if (abs(cord-value)0) + { + std::vector< std::string >::iterator it = std::find(_faceGroupNames.begin(), _faceGroupNames.end(), groupName); + if(it == _faceGroupNames.end())//No group named groupName + { + _faceGroupNames.insert(_faceGroupNames.end(),groupName); + _faceGroupsIds.insert( _faceGroupsIds.end(),faceIds); + _faceGroups.insert( _faceGroups.end(), NULL);//No mesh created. Create one ? + } + else + {cout<<"_faceGroupNames.size()="<<_faceGroupNames.size()<<", _faceGroupsIds.size()="<<_faceGroupsIds.size()< faceGroupIds = _faceGroupsIds[it-_faceGroupNames.begin()]; + faceGroupIds.insert( faceGroupIds.end(), faceIds.begin(), faceIds.end()); + /* Detect and erase duplicates face ids */ + sort( faceGroupIds.begin(), faceGroupIds.end() ); + faceGroupIds.erase( unique( faceGroupIds.begin(), faceGroupIds.end() ), faceGroupIds.end() ); + _faceGroupsIds[it-_faceGroupNames.begin()] = faceGroupIds; + } + } + if (nodeIds.size()>0) + { + std::vector< std::string >::iterator it = std::find(_nodeGroupNames.begin(), _nodeGroupNames.end(), groupName); + if(it == _nodeGroupNames.end())//No group named groupName + { + _nodeGroupNames.insert(_nodeGroupNames.end(),groupName); + _nodeGroupsIds.insert( _nodeGroupsIds.end(),nodeIds); + _nodeGroups.insert( _nodeGroups.end(), NULL);//No mesh created. Create one ? + } + else + { + std::vector< int > nodeGroupIds = _nodeGroupsIds[it-_nodeGroupNames.begin()]; + nodeGroupIds.insert( nodeGroupIds.end(), nodeIds.begin(), nodeIds.end()); + /* Detect and erase duplicates node ids */ + sort( nodeGroupIds.begin(), nodeGroupIds.end() ); + nodeGroupIds.erase( unique( nodeGroupIds.begin(), nodeGroupIds.end() ), nodeGroupIds.end() ); + _nodeGroupsIds[it-_nodeGroupNames.begin()] = nodeGroupIds; + } + } +} + +void +Mesh::setBoundaryNodesFromFaces() +{ + for (int iface=0;iface<_boundaryFaceIds.size();iface++) + { + std::vector< int > nodesID= _faces.get()[_boundaryFaceIds[iface]].getNodesId(); + int nbNodes = _faces.get()[_boundaryFaceIds[iface]].getNumberOfNodes(); + for(int inode=0 ; inode::const_iterator it = std::find(_boundaryNodeIds.begin(),_boundaryNodeIds.end(),nodesID[inode]); + if( it != _boundaryNodeIds.end() ) + _boundaryNodeIds.push_back(nodesID[inode]); + } + } +} + +std::map +Mesh::getIndexFacePeriodic( void ) const +{ + return _indexFacePeriodicMap; +} + +void +Mesh::setPeriodicFaces(bool check_groups, bool use_central_inversion) +{ + if(_indexFacePeriodicSet) + return; + + for (int indexFace=0;indexFace<_boundaryFaceIds.size() ; indexFace++) { - double Box0[2*_spaceDim]; - coo->getMinMaxPerComponent(Box0); - - _xMin=Box0[0]; - _xMax=Box0[1]; - if (_spaceDim>=2) + Face my_face=_faces.get()[_boundaryFaceIds[indexFace]]; + int iface_perio=-1; + if(_meshDim==1) { - _yMin=Box0[2]; - _yMax=Box0[3]; + for (int iface=0;iface<_boundaryFaceIds.size() ; iface++) + if(iface!=indexFace) + { + iface_perio=_boundaryFaceIds[iface]; + break; + } } - if (_spaceDim>=3) + else if(_meshDim==2) { - _zMin=Box0[4]; - _zMax=Box0[5]; + double x=my_face.x(); + double y=my_face.y(); + + for (int iface=0;iface<_boundaryFaceIds.size() ; iface++) + { + Face face_i=_faces.get()[_boundaryFaceIds[iface]]; + double xi=face_i.x(); + double yi=face_i.y(); + if ( (abs(y-yi)<_epsilon || abs(x-xi)<_epsilon )// Case of a square geometry + && ( !check_groups || my_face.getGroupName()!=face_i.getGroupName()) //In case groups need to be checked + && ( !use_central_inversion || abs(y+yi) + abs(x+xi)<_epsilon ) // Case of a central inversion + && fabs(my_face.getMeasure() - face_i.getMeasure())<_epsilon + && fabs(my_face.getXN() + face_i.getXN())<_epsilon + && fabs(my_face.getYN() + face_i.getYN())<_epsilon + && fabs(my_face.getZN() + face_i.getZN())<_epsilon ) + { + iface_perio=_boundaryFaceIds[iface]; + break; + } + } } + else if(_meshDim==3) + { + double x=my_face.x(); + double y=my_face.y(); + double z=my_face.z(); + + for (int iface=0;iface<_boundaryFaceIds.size() ; iface++) + { + Face face_i=_faces.get()[_boundaryFaceIds[iface]]; + double xi=face_i.x(); + double yi=face_i.y(); + double zi=face_i.z(); + if ( ((abs(y-yi)<_epsilon && abs(x-xi)<_epsilon) || (abs(x-xi)<_epsilon && abs(z-zi)<_epsilon) || (abs(y-yi)<_epsilon && abs(z-zi)<_epsilon))// Case of a cube geometry + && ( !check_groups || my_face.getGroupName()!=face_i.getGroupName()) //In case groups need to be checked + && ( !use_central_inversion || abs(y+yi) + abs(x+xi) + abs(z+zi)<_epsilon )// Case of a central inversion + && fabs(my_face.getMeasure() - face_i.getMeasure())<_epsilon + && fabs(my_face.getXN() + face_i.getXN())<_epsilon + && fabs(my_face.getYN() + face_i.getYN())<_epsilon + && fabs(my_face.getZN() + face_i.getZN())<_epsilon ) + { + iface_perio=_boundaryFaceIds[iface]; + break; + } + } + } + else + throw CdmathException("Mesh::setPeriodicFaces: Mesh dimension should be 1, 2 or 3"); + + if (iface_perio==-1) + throw CdmathException("Mesh::setPeriodicFaces: periodic face not found, iface_perio==-1 " ); + else + _indexFacePeriodicMap[_boundaryFaceIds[indexFace]]=iface_perio; } - - //Set boundary groups - _faceGroupNames.push_back("Boundary"); - _nodeGroupNames.push_back("Boundary"); - _faceGroupsIds.push_back(_boundaryFaceIds); - _nodeGroupsIds.push_back(_boundaryNodeIds); - if( _meshDim>1 ) - { //Set face boundary group - _boundaryMesh = mu->computeSkin(); - _faceGroups.push_back(_boundaryMesh); - } - else - _faceGroups.push_back(NULL); - _nodeGroups.push_back(NULL); - - desc->decrRef(); - descI->decrRef(); - revDesc->decrRef(); - revDescI->decrRef(); - mu2->decrRef(); - baryCell->decrRef(); - fields->decrRef(); - revNode->decrRef(); - revNodeI->decrRef(); - revCell->decrRef(); - revCellI->decrRef(); - - if (_meshDim == 3) - { - revNode2->decrRef(); - revNodeI2->decrRef(); - desc2->decrRef(); - descI2->decrRef(); - revDesc2->decrRef(); - revDescI2->decrRef(); - mu3->decrRef(); - } - - return mu; + _indexFacePeriodicSet=true; } -//---------------------------------------------------------------------- -Mesh::Mesh( std::vector points, std::string meshName ) -//---------------------------------------------------------------------- +int +Mesh::getIndexFacePeriodic(int indexFace, bool check_groups, bool use_central_inversion) { - int nx=points.size(); - - if(nx<2) - throw CdmathException("Mesh::Mesh( vector points, string meshName) : nx < 2, vector should contain at least two values"); - int i=0; - while( ipoints[i] ) - i++; - if( i!=nx-1 ) - { - //cout<< points << endl; - throw CdmathException("Mesh::Mesh( vector points, string meshName) : vector values should be sorted"); - } - - _spaceDim = 1 ; - _meshDim = 1 ; - _name=meshName; - _epsilon=1e-6; - _xMin=points[0]; - _xMax=points[nx-1]; - _yMin=0.; - _yMax=0.; - _zMin=0.; - _zMax=0.; + if (!_faces.get()[indexFace].isBorder()) + { + cout<<"Pb with indexFace= "<allocateCells(nx - 1); - double * coords = new double[nx]; - mcIdType * nodal_con = new mcIdType[2]; - coords[0]=points[0]; - for(int i=0; i::const_iterator it = _indexFacePeriodicMap.find(indexFace); + if( it != _indexFacePeriodicMap.end() ) + return it->second; + else { - nodal_con[0]=i; - nodal_con[1]=i+1; - mesh1d->insertNextCell(INTERP_KERNEL::NORM_SEG2, 2, nodal_con); - coords[i+1]=points[i + 1]; + cout<<"Pb with indexFace= "<finishInsertingCells(); - - DataArrayDouble * coords_arr = DataArrayDouble::New(); - coords_arr->alloc(nx,1); - std::copy(coords,coords+nx,coords_arr->getPointer()); - mesh1d->setCoords(coords_arr); - - delete [] coords, nodal_con; - coords_arr->decrRef(); - - _mesh=mesh1d->buildUnstructured();//To enable writeMED. Because we declared the mesh as unstructured, we decide to build the unstructured data (not mandatory) - _meshNotDeleted=true; - - setMesh(); } -//---------------------------------------------------------------------- -Mesh::Mesh( double xmin, double xmax, int nx, std::string meshName ) -//---------------------------------------------------------------------- +bool +Mesh::isBorderNode(int nodeid) const { - if(nx<=0) - throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx) : nx <= 0"); - if(xmin>=xmax) - throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx) : xmin >= xmax"); - - double dx = (xmax - xmin)/nx ; - - _spaceDim = 1 ; - _meshDim = 1 ; - _name=meshName; - _epsilon=1e-6; - _isStructured = true; - _indexFacePeriodicSet=false; - - _xMin=xmin; - _xMax=xmax; - _yMin=0.; - _yMax=0.; - _zMin=0.; - _zMax=0.; - - _dxyz.resize(_spaceDim); - _dxyz[0]=dx; - _nxyz.resize(_spaceDim); - _nxyz[0]=nx; - - double *originPtr = new double[_spaceDim]; - double *dxyzPtr = new double[_spaceDim]; - mcIdType *nodeStrctPtr = new mcIdType[_spaceDim]; - - originPtr[0]=xmin; - nodeStrctPtr[0]=nx+1; - dxyzPtr[0]=dx; - - _mesh=MEDCouplingIMesh::New(meshName, - _spaceDim, - nodeStrctPtr, - nodeStrctPtr+_spaceDim, - originPtr, - originPtr+_spaceDim, - dxyzPtr, - dxyzPtr+_spaceDim); - _meshNotDeleted=true;//Because the mesh is structured cartesian : no data in memory. No nodes and cell coordinates stored - - delete [] originPtr; - delete [] dxyzPtr; - delete [] nodeStrctPtr; - - setMesh(); + return _nodes.get()[nodeid].isBorder(); } -//---------------------------------------------------------------------- -Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, int split_to_triangles_policy, std::string meshName) -//---------------------------------------------------------------------- +bool +Mesh::isBorderFace(int faceid) const { - if(nx<=0 || ny<=0) - throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : nx <= 0 or ny <= 0"); - if(xmin>=xmax) - throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : xmin >= xmax"); - if(ymin>=ymax) - throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : ymin >= ymax"); - - _xMin=xmin; - _xMax=xmax; - _yMin=ymin; - _yMax=ymax; - _zMin=0.; - _zMax=0.; - - - double dx = (xmax - xmin)/nx ; - double dy = (ymax - ymin)/ny ; - - _spaceDim = 2 ; - _meshDim = 2 ; - _name=meshName; - _epsilon=1e-6; - _indexFacePeriodicSet=false; - _nxyz.resize(_spaceDim); - _nxyz[0]=nx; - _nxyz[1]=ny; + return _faces.get()[faceid].isBorder(); +} - _dxyz.resize(_spaceDim); - _dxyz[0]=dx; - _dxyz[1]=dy; +std::vector< int > +Mesh::getBoundaryFaceIds() const +{ + return _boundaryFaceIds; +} - double *originPtr = new double[_spaceDim]; - double *dxyzPtr = new double[_spaceDim]; - mcIdType *nodeStrctPtr = new mcIdType[_spaceDim]; +std::vector< int > +Mesh::getBoundaryNodeIds() const +{ + return _boundaryNodeIds; +} - originPtr[0]=xmin; - originPtr[1]=ymin; - nodeStrctPtr[0]=nx+1; - nodeStrctPtr[1]=ny+1; - dxyzPtr[0]=dx; - dxyzPtr[1]=dy; +bool +Mesh::isTriangular() const +{ + return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_TRI3; +} +bool +Mesh::isTetrahedral() const +{ + return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_TETRA4; +} +bool +Mesh::isQuadrangular() const +{ + return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_QUAD4; +} +bool +Mesh::isHexahedral() const +{ + return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_HEXA8; +} +bool +Mesh::isStructured() const +{ + return _isStructured; +} - _mesh=MEDCouplingIMesh::New(meshName, - _spaceDim, - nodeStrctPtr, - nodeStrctPtr+_spaceDim, - originPtr, - originPtr+_spaceDim, - dxyzPtr, - dxyzPtr+_spaceDim); - _meshNotDeleted=true;//Because the mesh is structured cartesian : no data in memory. No nodes and cell coordinates stored - _isStructured = true; +std::vector< INTERP_KERNEL::NormalizedCellType > +Mesh::getElementTypes() const +{ + return _eltsTypes; +} - if(split_to_triangles_policy==0 || split_to_triangles_policy==1) - { - _mesh=_mesh->buildUnstructured(); - _mesh->simplexize(split_to_triangles_policy); - _meshNotDeleted=true;//Now the mesh is unstructured and stored with nodes and cell coordinates - _isStructured = false; - } - else if (split_to_triangles_policy != -1) - { - cout<< "split_to_triangles_policy = "<< split_to_triangles_policy << endl; - throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : Unknown splitting policy"); +std::vector< string > +Mesh::getElementTypesNames() const +{ + std::vector< string > result(0); + for(int i=0; i< _eltsTypes.size(); i++) + { + if( _eltsTypes[i]==INTERP_KERNEL::NORM_POINT1) + result.push_back("Points"); + else if( _eltsTypes[i]==INTERP_KERNEL::NORM_SEG2) + result.push_back("Segments"); + else if( _eltsTypes[i]==INTERP_KERNEL::NORM_TRI3) + result.push_back("Triangles"); + else if( _eltsTypes[i]==INTERP_KERNEL::NORM_QUAD4) + result.push_back("Quadrangles"); + else if( _eltsTypes[i]==INTERP_KERNEL::NORM_POLYGON) + result.push_back("Polygons"); + else if( _eltsTypes[i]==INTERP_KERNEL::NORM_TETRA4) + result.push_back("Tetrahedra"); + else if( _eltsTypes[i]==INTERP_KERNEL::NORM_PYRA5) + result.push_back("Pyramids"); + else if( _eltsTypes[i]==INTERP_KERNEL::NORM_PENTA6) + result.push_back("Pentahedra"); + else if( _eltsTypes[i]==INTERP_KERNEL::NORM_HEXA8) + result.push_back("Hexahedra"); + else if( _eltsTypes[i]==INTERP_KERNEL::NORM_POLYHED) + result.push_back("Polyhedrons"); + else + { + cout<< "Mesh " + _name + " contains an element of type " <=xmax) - throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : xmin >= xmax"); - if(ymin>=ymax) - throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : ymin >= ymax"); - if(zmin>=zmax) - throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : zmin >= zmax"); - - _spaceDim = 3; - _meshDim = 3; - _name=meshName; - _epsilon=1e-6; - _xMin=xmin; - _xMax=xmax; - _yMin=ymin; - _yMax=ymax; - _zMin=zmin; - _zMax=zmax; + int nbCellsSubMesh; + bool foundFace; + + /* Searching for face groups */ + vector faceGroups=medmesh->getGroupsNames() ; - double dx = (xmax - xmin)/nx ; - double dy = (ymax - ymin)/ny ; - double dz = (zmax - zmin)/nz ; + for (unsigned int i=0;i nonEmptyGrp(medmesh->getGrpNonEmptyLevels(groupName)); + //We check if the group has a relative dimension equal to -1 + //before call to the function getGroup(-1,groupName.c_str()) + vector::iterator it = find(nonEmptyGrp.begin(), nonEmptyGrp.end(), -1); + if (it != nonEmptyGrp.end()) + { + MEDCouplingUMesh *m=medmesh->getGroup(-1,groupName.c_str()); + m->unPolyze(); + nbCellsSubMesh=m->getNumberOfCells(); + + _faceGroups.insert(_faceGroups.end(),m);//Vector of group meshes + _faceGroupNames.insert(_faceGroupNames.end(),groupName);//Vector of group names + _faceGroupsIds.insert(_faceGroupsIds.end(),std::vector(nbCellsSubMesh));//Vector of group face Ids. The filling of the face ids takes place below. + + DataArrayDouble *baryCell = m->computeIsoBarycenterOfNodesPerCell();//computeCellCenterOfMass() ; + const double *coorBary=baryCell->getConstPointer(); + // Face identification + for (int ic(0), k(0); ic coorBaryXyz(3,0); + for (int d=0; d<_spaceDim; d++) + coorBaryXyz[d] = coorBary[k+d]; + Point p1(coorBaryXyz[0],coorBaryXyz[1],coorBaryXyz[2]) ; - _dxyz.resize(_spaceDim); - _dxyz[0]=dx; - _dxyz[1]=dy; - _dxyz[2]=dz; + foundFace=false; + for (int iface=0;iface<_numberOfFaces;iface++ ) + { + Point p2=_faces.get()[iface].getBarryCenter(); + if(p1.distance(p2)<_epsilon) + { + _faces.get()[iface].setGroupName(groupName); + _faceGroupsIds[_faceGroupsIds.size()-1][ic]=iface; + foundFace=true; + break; + } + } + if (not foundFace) + throw CdmathException("No face found for group " + groupName ); + } + baryCell->decrRef(); + //m->decrRef(); + } + } +} +void +Mesh::setNodeGroups( const MEDFileMesh* medmesh, MEDCouplingUMesh* mu) +{ + int nbNodesSubMesh; + bool foundNode; + + /* Searching for node groups */ + vector nodeGroups=medmesh->getGroupsOnSpecifiedLev(1) ; - _nxyz.resize(_spaceDim); - _nxyz[0]=nx; - _nxyz[1]=ny; - _nxyz[2]=nz; + for (unsigned int i=0;igetNodeGroupArr( groupName ); + const mcIdType *nodeids=nodeGroup->getConstPointer(); - double *originPtr = new double[_spaceDim]; - double *dxyzPtr = new double[_spaceDim]; - mcIdType *nodeStrctPtr = new mcIdType[_spaceDim]; + if(nodeids!=NULL) + { + _nodeGroups.insert(_nodeGroups.end(),nodeGroup); + _nodeGroupNames.insert(_nodeGroupNames.end(),groupName); - originPtr[0]=xmin; - originPtr[1]=ymin; - originPtr[2]=zmin; - nodeStrctPtr[0]=nx+1; - nodeStrctPtr[1]=ny+1; - nodeStrctPtr[2]=nz+1; - dxyzPtr[0]=dx; - dxyzPtr[1]=dy; - dxyzPtr[2]=dz; + nbNodesSubMesh=nodeGroup->getNumberOfTuples();//nodeGroup->getNbOfElems(); - _mesh=MEDCouplingIMesh::New(meshName, - _spaceDim, - nodeStrctPtr, - nodeStrctPtr+_spaceDim, - originPtr, - originPtr+_spaceDim, - dxyzPtr, - dxyzPtr+_spaceDim); - _meshNotDeleted=true;//Because the mesh is structured cartesian : no data in memory. Nno nodes and cell coordinates stored - _isStructured = true; + DataArrayDouble *coo = mu->getCoords() ; + const double *cood=coo->getConstPointer(); - if( split_to_tetrahedra_policy == 0 ) - { - _mesh=_mesh->buildUnstructured(); - _mesh->simplexize(INTERP_KERNEL::PLANAR_FACE_5); - _meshNotDeleted=true;//Now the mesh is unstructured and stored with nodes and cell coordinates - _isStructured = false; - } - else if( split_to_tetrahedra_policy == 1 ) - { - _mesh=_mesh->buildUnstructured(); - _mesh->simplexize(INTERP_KERNEL::PLANAR_FACE_6); - _meshNotDeleted=true;//Now the mesh is unstructured and stored with nodes and cell coordinates - _isStructured = false; - } - else if ( split_to_tetrahedra_policy != -1 ) - { - cout<< "split_to_tetrahedra_policy = "<< split_to_tetrahedra_policy << endl; - throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : splitting policy value should be 0 or 1"); - } + _nodeGroupsIds.insert(_nodeGroupsIds.end(),std::vector(nbNodesSubMesh));//Vector of boundary faces + /* Node identification */ + for (int ic(0); ic coorP(3,0); + for (int d=0; d<_spaceDim; d++) + coorP[d] = cood[nodeids[ic]*_spaceDim+d]; + Point p1(coorP[0],coorP[1],coorP[2]) ; - delete [] originPtr; - delete [] dxyzPtr; - delete [] nodeStrctPtr; - - setMesh(); + foundNode=false; + for (int inode=0;inode<_numberOfNodes;inode++ ) + { + Point p2=_nodes.get()[inode].getPoint(); + if(p1.distance(p2)<_epsilon) + { + _nodes.get()[inode].setGroupName(groupName); + _nodeGroupsIds[_nodeGroupsIds.size()-1][ic]=inode; + foundNode=true; + break; + } + } + if (not foundNode) + throw CdmathException("No node found for group " + groupName ); + } + } + } } //---------------------------------------------------------------------- @@ -1638,15 +1476,6 @@ Mesh::getMeshDimension( void ) const return _meshDim ; } -std::vector -Mesh::getDXYZ() const -{ - if(!_isStructured) - throw CdmathException("std::vector Mesh::getDXYZ() : dx,dy and dz are defined only for structured meshes !"); - - return _dxyz; -} - std::vector Mesh::getCellGridStructure() const { @@ -1674,8 +1503,8 @@ Mesh::getNy( void ) const { if(!_isStructured) throw CdmathException("int Mesh::getNy( void ) : Ny is defined only for structured meshes !"); - if(_spaceDim < 2) - throw CdmathException("int double& Field::operator[ielem] : Ny is not defined in dimension < 2!"); + if(_meshDim < 2) + throw CdmathException("int Mesh::getNy( void ) : Ny is not defined for mesh dimension < 2!"); return _nxyz[1]; } @@ -1687,8 +1516,8 @@ Mesh::getNz( void ) const { if(!_isStructured) throw CdmathException("int Mesh::getNz( void ) : Nz is defined only for structured meshes !"); - if(_spaceDim < 3) - throw CdmathException("int Mesh::getNz( void ) : Nz is not defined in dimension < 3!"); + if(_meshDim < 3) + throw CdmathException("int Mesh::getNz( void ) : Nz is not defined for mesh dimension < 3!"); return _nxyz[2]; } @@ -1697,8 +1526,11 @@ Mesh::getNz( void ) const double Mesh::getXMin( void ) const //---------------------------------------------------------------------- -{ - return _xMin ; +{ + double Box0[2*_meshDim]; + _mesh->getBoundingBox(Box0); + + return Box0[0] ; } //---------------------------------------------------------------------- @@ -1706,7 +1538,10 @@ double Mesh::getXMax( void ) const //---------------------------------------------------------------------- { - return _xMax ; + double Box0[2*_meshDim]; + _mesh->getBoundingBox(Box0); + + return Box0[1] ; } //---------------------------------------------------------------------- @@ -1714,7 +1549,13 @@ double Mesh::getYMin( void ) const //---------------------------------------------------------------------- { - return _yMin ; + if(_meshDim<2) + throw CdmathException("Mesh::getYMin : dimension should be >=2"); + + double Box0[2*_meshDim]; + _mesh->getBoundingBox(Box0); + + return Box0[2] ; } //---------------------------------------------------------------------- @@ -1722,7 +1563,13 @@ double Mesh::getYMax( void ) const //---------------------------------------------------------------------- { - return _yMax ; + if(_meshDim<2) + throw CdmathException("Mesh::getYMax : dimension should be >=2"); + + double Box0[2*_meshDim]; + _mesh->getBoundingBox(Box0); + + return Box0[3] ; } //---------------------------------------------------------------------- @@ -1730,7 +1577,13 @@ double Mesh::getZMin( void ) const //---------------------------------------------------------------------- { - return _zMin ; + if(_meshDim<3) + throw CdmathException("Mesh::getZMin : dimension should be 3"); + + double Box0[2*_meshDim]; + _mesh->getBoundingBox(Box0); + + return Box0[4] ; } //---------------------------------------------------------------------- @@ -1738,7 +1591,13 @@ double Mesh::getZMax( void ) const //---------------------------------------------------------------------- { - return _zMax ; + if(_meshDim<3) + throw CdmathException("Mesh::getZMax : dimension should be 3"); + + double Box0[2*_meshDim]; + _mesh->getBoundingBox(Box0); + + return Box0[5] ; } //---------------------------------------------------------------------- @@ -1782,7 +1641,7 @@ Mesh::getNumberOfEdges ( void ) const } //---------------------------------------------------------------------- -Face* +std::shared_ptr Mesh::getFaces ( void ) const //---------------------------------------------------------------------- { @@ -1790,7 +1649,7 @@ Mesh::getFaces ( void ) const } //---------------------------------------------------------------------- -Cell* +std::shared_ptr Mesh::getCells ( void ) const //---------------------------------------------------------------------- { @@ -1802,7 +1661,7 @@ Cell& Mesh::getCell ( int i ) const //---------------------------------------------------------------------- { - return _cells[i] ; + return _cells.get()[i] ; } //---------------------------------------------------------------------- @@ -1810,7 +1669,7 @@ Face& Mesh::getFace ( int i ) const //---------------------------------------------------------------------- { - return _faces[i] ; + return _faces.get()[i] ; } //---------------------------------------------------------------------- @@ -1818,11 +1677,11 @@ Node& Mesh::getNode ( int i ) const //---------------------------------------------------------------------- { - return _nodes[i] ; + return _nodes.get()[i] ; } //---------------------------------------------------------------------- -Node* +std::shared_ptr Mesh::getNodes ( void ) const //---------------------------------------------------------------------- { @@ -1871,49 +1730,16 @@ Mesh::operator= ( const Mesh& mesh ) _meshNotDeleted = mesh.meshNotDeleted(); if(_isStructured) - { _nxyz = mesh.getCellGridStructure() ; - _dxyz=mesh.getDXYZ(); - _xMin=mesh.getXMin(); - _xMax=mesh.getXMax(); - _yMin=mesh.getYMin(); - _yMax=mesh.getYMax(); - _zMin=mesh.getZMin(); - _zMax=mesh.getZMax(); - } + _faceGroupNames = mesh.getNameOfFaceGroups() ; _faceGroups = mesh.getFaceGroups() ; _nodeGroupNames = mesh.getNameOfNodeGroups() ; _nodeGroups = mesh.getNodeGroups() ; - if (_nodes) - { - delete [] _nodes ; - _nodes=NULL; - } - if (_faces) - { - delete [] _faces ; - _faces=NULL; - } - if (_cells) - { - delete [] _cells ; - _cells=NULL; - } - - _nodes = new Node[_numberOfNodes] ; - _faces = new Face[_numberOfFaces] ; - _cells = new Cell[_numberOfCells] ; - - for (int i=0;i<_numberOfNodes;i++) - _nodes[i]=mesh.getNode(i); - - for (int i=0;i<_numberOfFaces;i++) - _faces[i]=mesh.getFace(i); - - for (int i=0;i<_numberOfCells;i++) - _cells[i]=mesh.getCell(i); + _nodes = mesh.getNodes() ; + _faces = mesh.getFaces() ; + _cells = mesh.getCells() ; _indexFacePeriodicSet= mesh.isIndexFacePeriodicSet(); if(_indexFacePeriodicSet) @@ -1925,8 +1751,8 @@ Mesh::operator= ( const Mesh& mesh ) _eltsTypes=mesh.getElementTypes(); _eltsTypesNames=mesh.getElementTypesNames(); - MCAuto m1=mesh.getMEDCouplingMesh()->deepCopy(); - _mesh=m1; + _mesh=mesh.getMEDCouplingMesh()->clone(false);//Clone because you will need to buildUnstructured. No deep copy : it is assumed node coordinates and cell connectivity will not change + return *this; } @@ -2008,11 +1834,11 @@ Mesh::getMaxNbNeighbours(EntityType type) const int nbNeib;//local number of neighbours for(int i=0; i<_numberOfCells; i++) { - Cell Ci = _cells[i]; + Cell Ci = _cells.get()[i]; //Careful with mesh with junctions : do not just take Ci.getNumberOfFaces() nbNeib=0; for(int j=0; jwriteVTK(fname.c_str()) ; @@ -2043,15 +1869,23 @@ Mesh::writeVTK ( const std::string fileName ) const //---------------------------------------------------------------------- void -Mesh::writeMED ( const std::string fileName ) const +Mesh::writeMED ( const std::string fileName, bool fromScratch ) const //---------------------------------------------------------------------- { if( !_meshNotDeleted ) - throw CdmathException("Mesh::writeMED : Cannot save mesh : no MEDCouplingUMesh loaded"); - + throw CdmathException("Mesh::writeMED : Cannot save mesh : no MEDCouplingMesh loaded (may be deleted)"); + string fname=fileName+".med"; - MEDCouplingUMesh* mu=_mesh->buildUnstructured(); - MEDCoupling::WriteUMesh(fname.c_str(),mu,true); + if(_isStructured)//Check if we have a medcouplingimesh that can't be written directly + { + const MEDCoupling::MEDCouplingIMesh* iMesh = dynamic_cast< const MEDCoupling::MEDCouplingIMesh* > ((const MEDCoupling::MEDCouplingMesh*) _mesh); + if(iMesh)//medcouplingimesh : Use convertToCartesian in order to write mesh + MEDCoupling::WriteMesh(fname.c_str(),iMesh->convertToCartesian(), fromScratch); + else//medcouplingcmesh : save directly + MEDCoupling::WriteMesh(fname.c_str(),_mesh, fromScratch); + } + else + MEDCoupling::WriteMesh(fname.c_str(),_mesh, fromScratch); /* Try to save mesh groups */ //MEDFileUMesh meshMEDFile; @@ -2110,7 +1944,7 @@ Mesh::setNodeGroupByIds(std::vector< int > nodeIds, std::string groupName) getNode(nodeIds[i]).setGroupName(groupName); } -void Mesh::deleteMEDCouplingUMesh() +void Mesh::deleteMEDCouplingMesh() { if(_meshNotDeleted) { -- 2.39.2