From fda0c19aa95c91b3f56ce4576675bae62a604bf0 Mon Sep 17 00:00:00 2001 From: michael Date: Fri, 10 Dec 2021 14:18:22 +0100 Subject: [PATCH] Improved Mesh class, especially when meshDim=1 (graphs) --- CDMATH/mesh/inc/Mesh.hxx | 61 +- CDMATH/mesh/src/Mesh.cxx | 1129 ++++++++++++++++++++------------------ 2 files changed, 644 insertions(+), 546 deletions(-) diff --git a/CDMATH/mesh/inc/Mesh.hxx b/CDMATH/mesh/inc/Mesh.hxx index a7ea9c9..99ce471 100644 --- a/CDMATH/mesh/inc/Mesh.hxx +++ b/CDMATH/mesh/inc/Mesh.hxx @@ -106,7 +106,7 @@ public: //---------------------------------------------------------------- * @param nx : Number of cells in x direction * @param ny : Number of cells in y direction * @param nz : Number of cells in z direction - * @param split_to_tetrahedra_policy : each cuboid will be split into 5 tetrahedra if value is INTERP_KERNEL::PLANAR_FACE_5 or 6 tetrahedra if the value is INTERP_KERNEL::PLANAR_FACE_6 + * @param split_to_tetrahedra_policy : each cuboid will be split into 5 tetrahedra if value is 0 or 6 tetrahedra if the value is 1 */ 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") ; @@ -256,10 +256,15 @@ public: //---------------------------------------------------------------- MEDCoupling::MCAuto getMEDCouplingMesh ( void ) const ; /** - * \brief computes the skin surrounding the mesh + * \brief return the skin surrounding the mesh */ Mesh getBoundaryMesh ( void ) const ; + /** + * \brief return a group surrounding the mesh + */ + Mesh getBoundaryGroupMesh ( std::string groupName, int nth_group_match = 0 ) const ; + /** * \brief return the list of face group names * return _faceGroupNames @@ -303,11 +308,11 @@ public: //---------------------------------------------------------------- /** * @return list of face group Ids */ - std::vector< int > getGroupFaceIds(std::string groupName) const; + std::vector< int > getFaceGroupIds(std::string groupName, bool isBoundaryGroup=true) const; /** * @return list of node group Ids * */ - std::vector< int > getGroupNodeIds(std::string groupName) const; + std::vector< int > getNodeGroupIds(std::string groupName, bool isBoundaryGroup=true) const; /** * \brief write mesh in the VTK format @@ -319,9 +324,11 @@ public: //---------------------------------------------------------------- */ void writeMED ( const std::string fileName ) const ; - void setGroupAtPlan(double value, int direction, double eps, std::string groupName) ; + void setGroupAtPlan(double value, int direction, double eps, std::string groupName, bool isBoundaryGroup=true) ; - void setGroupAtFaceByCoords(double x, double y, double z, double eps, std::string groupName) ; + void setGroupAtFaceByCoords(double x, double y, double z, double eps, std::string groupName, bool isBoundaryGroup=true) ; + + void setGroupAtNodeByCoords(double x, double y, double z, double eps, std::string groupName, bool isBoundaryGroup=true) ; void setFaceGroupByIds(std::vector< int > faceIds, std::string groupName) ; @@ -349,9 +356,9 @@ public: //---------------------------------------------------------------- double getComparisonEpsilon() const {return _epsilon;}; void setComparisonEpsilon(double epsilon){ _epsilon=epsilon;}; // Quick comparison of two meshes to see if they are identical with high probability (three cells are compared) - void checkFastEquivalWith( Mesh m) const { return getMEDCouplingMesh()->checkFastEquivalWith(m.getMEDCouplingMesh(),1e-6);}; + void checkFastEquivalWith( Mesh m) const { return getMEDCouplingMesh()->checkFastEquivalWith(m.getMEDCouplingMesh(),_epsilon);}; // Deep comparison of two meshes to see if they are identical Except for their names and units - bool isEqualWithoutConsideringStr( Mesh m) const { return getMEDCouplingMesh()->isEqualWithoutConsideringStr(m.getMEDCouplingMesh(),1e-6);}; + bool isEqualWithoutConsideringStr( Mesh m) const { return getMEDCouplingMesh()->isEqualWithoutConsideringStr(m.getMEDCouplingMesh(),_epsilon);}; std::vector< std::string > getElementTypesNames() const ; /** @@ -364,12 +371,29 @@ public: //---------------------------------------------------------------- */ int getMaxNbNeighbours(EntityType type) const; + /** + * \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 _meshNotDeleted;} + private: //---------------------------------------------------------------- MEDCoupling::MEDCouplingUMesh* setMesh( void ) ; + void setGroups( const MEDCoupling::MEDFileUMesh* medmesh, MEDCoupling::MEDCouplingUMesh* mu) ;//Read all face and node group + void addNewFaceGroup( const MEDCoupling::MEDCouplingUMesh *m);//adds one face group in the vectors _faceGroups, _faceGroupNames and _faceGroupIds + + /* + * The MEDCoupling mesh + */ + MEDCoupling::MCAuto _mesh; - void setGroups( const MEDCoupling::MEDFileUMesh* medmesh, MEDCoupling::MEDCouplingUMesh* mu) ; - + bool _meshNotDeleted; + std::string _name; /** @@ -456,11 +480,20 @@ private: //---------------------------------------------------------------- * The list of node groups. */ std::vector _nodeGroups; + /* - * The mesh MEDCoupling + * The list of face id in each face groups. + */ + std::vector< std::vector > _faceGroupsIds; + + /* + * The list of node id in each node groups. + */ + std::vector< std::vector > _nodeGroupsIds; + + /* + * Elements types (SEG2, TRI3, QUAD4, HEXA6 ...) */ - MEDCoupling::MCAuto _mesh; - 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; @@ -476,7 +509,7 @@ private: //---------------------------------------------------------------- /* List of boundary nodes*/ std::vector< int > _boundaryNodeIds; /* Boundary mesh */ - const MEDCoupling::MEDCouplingUMesh * _boundaryMesh; + MEDCoupling::MEDCouplingUMesh * _boundaryMesh; double _epsilon; }; diff --git a/CDMATH/mesh/src/Mesh.cxx b/CDMATH/mesh/src/Mesh.cxx index 1b1d4fe..44ca040 100644 --- a/CDMATH/mesh/src/Mesh.cxx +++ b/CDMATH/mesh/src/Mesh.cxx @@ -20,6 +20,7 @@ #include #include #include +#include using namespace MEDCoupling; using namespace std; @@ -29,6 +30,7 @@ Mesh::Mesh( void ) //---------------------------------------------------------------------- { _mesh=NULL; + _meshNotDeleted=false; _cells=NULL; _nodes=NULL; _faces=NULL; @@ -47,12 +49,18 @@ Mesh::Mesh( void ) _zMax=0.; _nxyz.resize(0); _dxyz.resize(0.); + _boundaryMesh=NULL; + _boundaryFaceIds.resize(0); + _boundaryNodeIds.resize(0); _faceGroupNames.resize(0); _faceGroups.resize(0); + _faceGroupsIds.resize(0); _nodeGroupNames.resize(0); _nodeGroups.resize(0); + _nodeGroupsIds.resize(0); _indexFacePeriodicSet=false; _name=""; + _epsilon=1e-6; } //---------------------------------------------------------------------- @@ -62,9 +70,12 @@ Mesh::~Mesh( void ) delete [] _cells; delete [] _nodes; delete [] _faces; + //for(int i=0; i< _faceGroups.size(); i++) // _faceGroups[i]->decrRef(); // _nodeGroups[i]->decrRef(); + if( _meshNotDeleted) + (_mesh.retn())->decrRef(); } std::string @@ -77,12 +88,12 @@ Mesh::Mesh( const MEDCoupling::MEDCouplingIMesh* mesh ) { _spaceDim=mesh->getSpaceDimension(); _meshDim=mesh->getMeshDimension(); - _isStructured=true; _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]; @@ -116,10 +127,14 @@ Mesh::Mesh( const MEDCoupling::MEDCouplingIMesh* mesh ) originPtr+_spaceDim, dxyzPtr, dxyzPtr+_spaceDim); + _meshNotDeleted=true; + _isStructured=true; + delete [] originPtr; delete [] dxyzPtr; delete [] nodeStrctPtr; delete [] Box0 ; + setMesh(); } @@ -127,10 +142,10 @@ Mesh::Mesh( const MEDCoupling::MEDCouplingUMesh* mesh ) { _spaceDim=mesh->getSpaceDimension(); _meshDim=mesh->getMeshDimension(); - _isStructured=false; double* Box0=new double[2*_spaceDim]; mesh->getBoundingBox(Box0); _name=mesh->getName(); + _epsilon=1e-6; _indexFacePeriodicSet=false; _xMin=Box0[0]; @@ -147,10 +162,21 @@ Mesh::Mesh( const MEDCoupling::MEDCouplingUMesh* mesh ) } _mesh=mesh->deepCopy(); + _mesh=mesh->buildUnstructured(); + _meshNotDeleted=true; + _isStructured=false; delete [] Box0 ; + setMesh(); } +//---------------------------------------------------------------------- +Mesh::Mesh( const std::string filename, int meshLevel ) +//---------------------------------------------------------------------- +{ + readMeshMed(filename, meshLevel); +} + //---------------------------------------------------------------------- Mesh::Mesh( const Mesh& mesh ) //---------------------------------------------------------------------- @@ -158,6 +184,7 @@ Mesh::Mesh( const Mesh& mesh ) _spaceDim = mesh.getSpaceDimension() ; _meshDim = mesh.getMeshDimension() ; _name=mesh.getName(); + _epsilon=mesh.getComparisonEpsilon(); _xMax=mesh.getXMax(); _yMin=mesh.getYMin(); _yMax=mesh.getYMax(); @@ -208,13 +235,7 @@ Mesh::Mesh( const Mesh& mesh ) MCAuto m1=mesh.getMEDCouplingMesh()->deepCopy(); _mesh=m1; -} - -//---------------------------------------------------------------------- -Mesh::Mesh( const std::string filename, int meshLevel ) -//---------------------------------------------------------------------- -{ - readMeshMed(filename, meshLevel); + _meshNotDeleted=mesh.meshNotDeleted(); } //---------------------------------------------------------------------- @@ -229,11 +250,13 @@ Mesh::readMeshMed( const std::string filename, const int meshLevel) _meshDim=_mesh->getMeshDimension(); _spaceDim=_mesh->getSpaceDimension(); _name=_mesh->getName(); + _epsilon=1e-6; _indexFacePeriodicSet=false; MEDCoupling::MEDCouplingIMesh* structuredMesh = dynamic_cast (_mesh.retn()); if(structuredMesh) { _isStructured=true; + _meshNotDeleted=false; _dxyz=structuredMesh->getDXYZ(); _nxyz=structuredMesh->getCellGridStructure(); double* Box0=new double[2*_spaceDim]; @@ -253,73 +276,218 @@ Mesh::readMeshMed( const std::string filename, const int meshLevel) } } else + { _isStructured=false; + _meshNotDeleted=true; + } MEDCouplingUMesh* mu = setMesh(); setGroups(m, mu); cout<getGroup(-1,groupName.c_str()); - mu->unPolyze(); - _faceGroups.insert(_faceGroups.begin(),m); - _faceGroupNames.insert(_faceGroupNames.begin(),groupName); + 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(); - - int nbCellsSubMesh=m->getNumberOfCells(); + /* Face identification */ for (int ic(0), k(0); ic coorBaryXyz(3,0); @@ -570,41 +742,27 @@ Mesh::setGroups( const MEDFileUMesh* medmesh, MEDCouplingUMesh* mu) coorBaryXyz[d] = coorBary[k+d]; Point p1(coorBaryXyz[0],coorBaryXyz[1],coorBaryXyz[2]) ; - int flag=0; - /* double min=INFINITY; - Point p3; - int closestFaceNb;*/ + foundFace=false; for (int iface=0;iface<_numberOfFaces;iface++ ) { Point p2=_faces[iface].getBarryCenter(); - /*if(groupName=="Wall" and p1.distance(p2)1.E-10) - { - cout.precision(15); - cout<<"Subcell number "<< ic <<" Min distance to Wall = "<decrRef(); //m->decrRef(); } } - //Searching for node groups + /* Searching for node groups */ vector nodeGroups=medmesh->getGroupsOnSpecifiedLev(1) ; for (unsigned int i=0;igetNumberOfTuples();//nodeGroup->getNbOfElems(); + 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); @@ -632,19 +790,20 @@ Mesh::setGroups( const MEDFileUMesh* medmesh, MEDCouplingUMesh* mu) coorP[d] = cood[nodeids[ic]*_spaceDim+d]; Point p1(coorP[0],coorP[1],coorP[2]) ; - int flag=0; + foundNode=false; for (int inode=0;inode<_numberOfNodes;inode++ ) { Point p2=_nodes[inode].getPoint(); - if(p1.distance(p2)<1.E-10) + if(p1.distance(p2)<_epsilon) { _nodes[inode].setGroupName(groupName); - flag=1; + _nodeGroupsIds[_nodeGroupsIds.size()-1][ic]=inode; + foundNode=true; break; } } - if (flag==0) - throw CdmathException("No node belonging to group " + groupName + " found"); + if (not foundNode) + throw CdmathException("No node found for group " + groupName ); } } } @@ -655,26 +814,34 @@ MEDCouplingUMesh* Mesh::setMesh( void ) //---------------------------------------------------------------------- { + /* This is the main function translating medcouplingumesh info into Mesh class to be used when designing numerical methods + * We need the level 0 mesh to extract the cell-node connectvity + * We need the level -1 mesh to extract the cell-face and face-node connectivities (use o build descending connectivity) + * Be careful : the nodes in the medcoupling mesh are not necessarily all conected to a cell/face. + * Mesh class discard isolated nodes, hence the number of nodes in Mesh class can be lower than the number of nodes in medcouplingumesh. + */ + DataArrayIdType *desc = DataArrayIdType::New(); DataArrayIdType *descI = DataArrayIdType::New(); DataArrayIdType *revDesc = DataArrayIdType::New(); DataArrayIdType *revDescI = DataArrayIdType::New(); MEDCouplingUMesh* mu = _mesh->buildUnstructured(); - + MEDCouplingUMesh* mu2;//mesh of dimension N-1 containing the cell interfaces->cell/face connectivity + mu->unPolyze(); - /* Save the boundary mesh for later use*/ - _boundaryMesh = mu->computeSkin(); + mu->sortCellsInMEDFileFrmt( ); - MEDCouplingUMesh* mu2=mu->buildDescendingConnectivity2(desc,descI,revDesc,revDescI);//mesh of dimension N-1 containing the cell interfaces - - const mcIdType *tmp = desc->getConstPointer(); + if(_meshDim<2) + mu2=mu->buildDescendingConnectivity(desc,descI,revDesc,revDescI); + else + mu2=mu->buildDescendingConnectivity2(desc,descI,revDesc,revDescI); + + const mcIdType *tmp = desc->getConstPointer();//Lists the faces surrounding each cell const mcIdType *tmpI=descI->getConstPointer(); - const mcIdType *tmpA =revDesc->getConstPointer(); + const mcIdType *tmpA =revDesc->getConstPointer();//Lists the cells surrounding each face const mcIdType *tmpAI=revDescI->getConstPointer(); - //const int *work=tmp+tmpI[id];//corresponds to buildDescendingConnectivity - //Test du type d'éléments contenus dans le maillage afin d'éviter les éléments contenant des points de gauss _eltsTypes=mu->getAllGeoTypesSorted(); for(int i=0; i<_eltsTypes.size();i++) @@ -694,37 +861,37 @@ Mesh::setMesh( void ) } DataArrayDouble *baryCell = mu->computeCellCenterOfMass() ; - const double *coorBary=baryCell->getConstPointer(); + const double *coorBary=baryCell->getConstPointer();//Used for cell center coordinates MEDCouplingFieldDouble* fields=mu->getMeasureField(true); DataArrayDouble *surface = fields->getArray(); - const double *surf=surface->getConstPointer(); + const double *surf=surface->getConstPointer();//Used for cell lenght/surface/volume DataArrayDouble *coo = mu->getCoords() ; - const double *cood=coo->getConstPointer(); + const double *cood=coo->getConstPointer();//Used for nodes coordinates DataArrayIdType *revNode =DataArrayIdType::New(); DataArrayIdType *revNodeI=DataArrayIdType::New(); mu->getReverseNodalConnectivity(revNode,revNodeI) ; - const mcIdType *tmpN =revNode->getConstPointer(); + const mcIdType *tmpN =revNode->getConstPointer();//Used to know which cells surround a given node const mcIdType *tmpNI=revNodeI->getConstPointer(); DataArrayIdType *revCell =DataArrayIdType::New(); DataArrayIdType *revCellI=DataArrayIdType::New(); - mu2->getReverseNodalConnectivity(revCell,revCellI) ; - const mcIdType *tmpC =revCell->getConstPointer(); + mu2->getReverseNodalConnectivity(revCell,revCellI); + const mcIdType *tmpC =revCell->getConstPointer();//Used to know which faces surround a given node const mcIdType *tmpCI=revCellI->getConstPointer(); const DataArrayIdType *nodal = mu2->getNodalConnectivity() ; const DataArrayIdType *nodalI = mu2->getNodalConnectivityIndex() ; - const mcIdType *tmpNE =nodal->getConstPointer(); + const mcIdType *tmpNE =nodal->getConstPointer();//Used to know which nodes surround a given face const mcIdType *tmpNEI=nodalI->getConstPointer(); _numberOfCells = mu->getNumberOfCells() ; _cells = new Cell[_numberOfCells] ; - _numberOfNodes = mu->getNumberOfNodes() ; - _nodes = new Node[_numberOfNodes] ; + _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 _numberOfFaces = mu2->getNumberOfCells(); _faces = new Face[_numberOfFaces] ; @@ -756,86 +923,134 @@ Mesh::setMesh( void ) } // _cells, _nodes and _faces initialization: - if (_spaceDim == 1) + if (_meshDim == 1) { + double xn, yn=0., zn=0.;//Components of the normal vector at a cell interface + double norm; for( int id=0;id<_numberOfCells;id++ ) { - const mcIdType *work=tmp+tmpI[id]; - int nbFaces=tmpI[id+1]-tmpI[id]; - - int nbVertices=mu->getNumberOfNodesInCell(id) ; - - Cell ci( nbVertices, nbFaces, surf[id], Point(coorBary[id], 0.0, 0.0) ) ; - + Point p(0.0,0.0,0.0) ; + for(int idim=0; idim<_spaceDim; idim++) + p[idim]=coorBary[id*_spaceDim+idim]; + + mcIdType nbVertices=mu->getNumberOfNodesInCell(id) ;//should be equal to 2 + assert( nbVertices==2); std::vector nodeIdsOfCell ; mu->getNodeIdsOfCell(id,nodeIdsOfCell) ; - for( int el=0;el 0.0) ? -1.0 : 1.0; - - for( int el=0;el1) + yn = cood[nodeIdsOfCell[0]*_spaceDim+1] - cood[nodeIdsOfCell[nbFaces-1]*_spaceDim+1]; + if(_spaceDim>2) + zn = cood[nodeIdsOfCell[0]*_spaceDim+2] - cood[nodeIdsOfCell[nbFaces-1]*_spaceDim+2]; + norm = sqrt(xn*xn+yn*yn+zn*zn); + if(norm<_epsilon) + throw CdmathException("!!! Mesh::setMesh Normal vector has norm 0 !!!"); + else { - ci.addNormalVector(el,xn,0.0,0.0) ; - int indexFace=abs(work[el])-1; - ci.addFaceId(el,indexFace) ; - xn = - xn; + xn /= norm; + yn /= norm; + zn /= norm; + } + + Cell ci( nbVertices, nbFaces, surf[id], p ) ;//nbCells=nbFaces=2 + for( int el=0;el coo(0) ; + mu2->getCoordinatesOfNode(workv[0],coo); + Point p(0,0.0,0.0) ; + for(int idim=0; idim<_spaceDim; idim++) + p[idim]=coo[idim]; + + const mcIdType *workc=tmpA+tmpAI[id]; + mcIdType nbCells=tmpAI[id+1]-tmpAI[id]; + 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_id0)//To make sure this is not an isolated node + { + correctNbNodes++; + std::vector coo(0) ; + mu->getCoordinatesOfNode(id,coo); + Point p(0,0.0,0.0) ; + for(int idim=0; idim<_spaceDim; idim++) + p[idim]=coo[idim]; + + const mcIdType *workf=tmpC+tmpCI[id]; + mcIdType nbFaces=tmpCI[id+1]-tmpCI[id]; + assert( nbFaces==1); + + const mcIdType *workn=tmpN+tmpNI[id]; + mcIdType nbNeighbourNodes=tmpNI[id+1]-tmpNI[id]; + + Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ; + for( int el=0;el1 ) + { //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(); @@ -1112,185 +1337,10 @@ Mesh::setMesh( void ) revDescI2->decrRef(); mu3->decrRef(); } - + return mu; } -//---------------------------------------------------------------------- -Mesh::Mesh( double xmin, double xmax, int nx, std::string meshName ) -//---------------------------------------------------------------------- -{ - 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; - _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); - delete [] originPtr; - delete [] dxyzPtr; - delete [] nodeStrctPtr; - - DataArrayIdType *desc=DataArrayIdType::New(); - DataArrayIdType *descI=DataArrayIdType::New(); - DataArrayIdType *revDesc=DataArrayIdType::New(); - DataArrayIdType *revDescI=DataArrayIdType::New(); - MEDCouplingUMesh* mu=_mesh->buildUnstructured(); - MEDCouplingUMesh *mu2=mu->buildDescendingConnectivity(desc,descI,revDesc,revDescI); - - const mcIdType *tmp=desc->getConstPointer(); - const mcIdType *tmpI=descI->getConstPointer(); - - const mcIdType *tmpA =revDesc->getConstPointer(); - const mcIdType *tmpAI=revDescI->getConstPointer(); - - _eltsTypes=mu->getAllGeoTypesSorted(); - - DataArrayDouble *baryCell = mu->computeCellCenterOfMass() ; - const double *coorBary=baryCell->getConstPointer(); - - _numberOfCells = _mesh->getNumberOfCells() ; - _cells = new Cell[_numberOfCells] ; - - _numberOfNodes = mu->getNumberOfNodes() ; - _nodes = new Node[_numberOfNodes] ; - - _numberOfFaces = _numberOfNodes; - _faces = new Face[_numberOfFaces] ; - - _numberOfEdges = _numberOfCells; - - MEDCouplingFieldDouble* fieldl=mu->getMeasureField(true); - DataArrayDouble *longueur = fieldl->getArray(); - const double *lon=longueur->getConstPointer(); - - DataArrayDouble *coo = mu->getCoords() ; - const double *cood=coo->getConstPointer(); - - int comp=0; - for( int id=0;id<_numberOfCells;id++ ) - { - int nbVertices=mu->getNumberOfNodesInCell(id) ; - Point p(coorBary[id],0.0,0.0) ; - Cell ci( nbVertices, nbVertices, lon[id], p ) ; - - std::vector nodeIdsOfCell ; - mu->getNodeIdsOfCell(id,nodeIdsOfCell) ; - for( int el=0;el 0.0) ? -1.0 : 1.0; - - mcIdType nbFaces=tmpI[id+1]-tmpI[id]; - const mcIdType *work=tmp+tmpI[id]; - - for( int el=0;elgetReverseNodalConnectivity(revNode,revNodeI) ; - const mcIdType *tmpN=revNode->getConstPointer(); - const mcIdType *tmpNI=revNodeI->getConstPointer(); - - for( int id=0;id<_numberOfNodes;id++ ) - { - std::vector coo ; - mu->getCoordinatesOfNode(id,coo); - Point p(coo[0],0.0,0.0) ; - const mcIdType *workc=tmpN+tmpNI[id]; - mcIdType nbCells=tmpNI[id+1]-tmpNI[id]; - mcIdType nbFaces=1; - const mcIdType *workn=tmpA+tmpAI[id]; - mcIdType nbNeighbourNodes=tmpAI[id+1]-tmpAI[id]; - - Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ; - for( int el=0;eldecrRef(); - baryCell->decrRef(); - desc->decrRef(); - descI->decrRef(); - revDesc->decrRef(); - revDescI->decrRef(); - revNode->decrRef(); - revNodeI->decrRef(); - mu2->decrRef(); - mu->decrRef(); -} - //---------------------------------------------------------------------- Mesh::Mesh( std::vector points, std::string meshName ) //---------------------------------------------------------------------- @@ -1311,6 +1361,7 @@ Mesh::Mesh( std::vector points, std::string meshName ) _spaceDim = 1 ; _meshDim = 1 ; _name=meshName; + _epsilon=1e-6; _xMin=points[0]; _xMax=points[nx-1]; _yMin=0.; @@ -1343,110 +1394,67 @@ Mesh::Mesh( std::vector points, std::string meshName ) delete [] coords, nodal_con; coords_arr->decrRef(); - _mesh=mesh1d; - - DataArrayIdType *desc=DataArrayIdType::New(); - DataArrayIdType *descI=DataArrayIdType::New(); - DataArrayIdType *revDesc=DataArrayIdType::New(); - DataArrayIdType *revDescI=DataArrayIdType::New(); - MEDCouplingUMesh* mu=_mesh->buildUnstructured(); - MEDCouplingUMesh *mu2=mu->buildDescendingConnectivity(desc,descI,revDesc,revDescI); - - DataArrayDouble *baryCell = mu->computeCellCenterOfMass() ; - const double *coorBary=baryCell->getConstPointer(); + _mesh=mesh1d->buildUnstructured();//To enable writeMED. Because we declared the mesh as unstructured, we decide to build the unstructured data (not mandatory) + _meshNotDeleted=true; - _numberOfCells = _mesh->getNumberOfCells() ; - _cells = new Cell[_numberOfCells] ; + setMesh(); +} - _numberOfEdges = _numberOfCells; +//---------------------------------------------------------------------- +Mesh::Mesh( double xmin, double xmax, int nx, std::string meshName ) +//---------------------------------------------------------------------- +{ + 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"); - _eltsTypes=mu->getAllGeoTypesSorted(); + double dx = (xmax - xmin)/nx ; - MEDCouplingFieldDouble* fieldl=mu->getMeasureField(true); - DataArrayDouble *longueur = fieldl->getArray(); - const double *lon=longueur->getConstPointer(); + _spaceDim = 1 ; + _meshDim = 1 ; + _name=meshName; + _epsilon=1e-6; + _isStructured = true; + _indexFacePeriodicSet=false; - int comp=0; - for( int id=0;id<_numberOfCells;id++ ) - { - int nbVertices=mu->getNumberOfNodesInCell(id) ; - Point p(coorBary[id],0.0,0.0) ; - Cell ci( nbVertices, nbVertices, lon[id], p ) ; + _xMin=xmin; + _xMax=xmax; + _yMin=0.; + _yMax=0.; + _zMin=0.; + _zMax=0.; - std::vector nodeIdsOfCell ; - mu->getNodeIdsOfCell(id,nodeIdsOfCell) ; - for( int el=0;elgetReverseNodalConnectivity(revNode,revNodeI) ; - const mcIdType *tmpN=revNode->getConstPointer(); - const mcIdType *tmpNI=revNodeI->getConstPointer(); + originPtr[0]=xmin; + nodeStrctPtr[0]=nx+1; + dxyzPtr[0]=dx; - _numberOfNodes = mu->getNumberOfNodes() ; - _nodes = new Node[_numberOfNodes] ; - _numberOfFaces = _numberOfNodes; - _faces = new Face[_numberOfFaces] ; + _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 - for( int id=0;id<_numberOfNodes;id++ ) - { - std::vector coo ; - mu->getCoordinatesOfNode(id,coo); - Point p(coo[0],0.0,0.0) ; - const mcIdType *workc=tmpN+tmpNI[id]; - mcIdType nbCells=tmpNI[id+1]-tmpNI[id]; - mcIdType nbFaces=1; - const mcIdType *workn=tmpN+tmpNI[id]; - mcIdType nbNeighbourNodes=tmpNI[id+1]-tmpNI[id]; - Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ; - int nbVertices=1; - /* provisoire !!!!!!!!!!!!*/ - // Point pf(0.0,0.0,0.0) ; - Face fi( nbVertices, nbCells, 0.0, p, 0., 0., 0. ) ; - - for( int el=0;eldecrRef(); - baryCell->decrRef(); - desc->decrRef(); - descI->decrRef(); - revDesc->decrRef(); - revDescI->decrRef(); - revNode->decrRef(); - revNodeI->decrRef(); - mu2->decrRef(); - mu->decrRef(); + setMesh(); } + //---------------------------------------------------------------------- Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, int split_to_triangles_policy, std::string meshName) //---------------------------------------------------------------------- @@ -1472,8 +1480,8 @@ Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, _spaceDim = 2 ; _meshDim = 2 ; _name=meshName; + _epsilon=1e-6; _indexFacePeriodicSet=false; - _isStructured = true; _nxyz.resize(_spaceDim); _nxyz[0]=nx; _nxyz[1]=ny; @@ -1501,11 +1509,15 @@ Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, 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; 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) { @@ -1536,7 +1548,7 @@ Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, _spaceDim = 3; _meshDim = 3; _name=meshName; - _indexFacePeriodicSet=false; + _epsilon=1e-6; _xMin=xmin; _xMax=xmax; _yMin=ymin; @@ -1548,7 +1560,6 @@ Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double dy = (ymax - ymin)/ny ; double dz = (zmax - zmin)/nz ; - _isStructured = true; _dxyz.resize(_spaceDim); _dxyz[0]=dx; _dxyz[1]=dy; @@ -1581,16 +1592,22 @@ Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, 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; 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 ) { @@ -1844,12 +1861,15 @@ Mesh::operator= ( const Mesh& mesh ) _spaceDim = mesh.getSpaceDimension() ; _meshDim = mesh.getMeshDimension() ; _name = mesh.getName(); + _epsilon=mesh.getComparisonEpsilon(); _numberOfNodes = mesh.getNumberOfNodes(); _numberOfFaces = mesh.getNumberOfFaces(); _numberOfCells = mesh.getNumberOfCells(); _numberOfEdges = mesh.getNumberOfEdges(); _isStructured = mesh.isStructured(); + _meshNotDeleted = mesh.meshNotDeleted(); + if(_isStructured) { _nxyz = mesh.getCellGridStructure() ; @@ -1946,16 +1966,57 @@ Mesh::getBoundaryMesh ( void ) const return Mesh(_boundaryMesh); } +Mesh +Mesh::getBoundaryGroupMesh ( std::string groupName, int nth_occurence ) const +{ + //count occurences of groupName in known group name list + int count_occurences = std::count (_faceGroupNames.begin(),_faceGroupNames.end(),groupName); + + if( count_occurences ==0 )//No group found + { + cout<<"Mesh::getBoundaryGroupMesh Error : face group " << groupName << " does not exist"<1 )//Wrning several groups found + cout<<"Warning : "<::const_iterator it = _faceGroupNames.begin(); + for (int i = 0; iwriteVTK(fname.c_str()) ; } @@ -1982,10 +2046,14 @@ void Mesh::writeMED ( const std::string fileName ) const //---------------------------------------------------------------------- { + if( !_meshNotDeleted ) + throw CdmathException("Mesh::writeMED : Cannot save mesh : no MEDCouplingUMesh loaded"); + string fname=fileName+".med"; MEDCouplingUMesh* mu=_mesh->buildUnstructured(); MEDCoupling::WriteUMesh(fname.c_str(),mu,true); + /* Try to save mesh groups */ //MEDFileUMesh meshMEDFile; //meshMEDFile.setMeshAtLevel(0,mu); //for(int i=0; i< _faceGroups.size(); i++) @@ -1994,51 +2062,34 @@ Mesh::writeMED ( const std::string fileName ) const //MEDCoupling::meshMEDFile.write(fname.c_str(),2) ; //else //MEDCoupling::meshMEDFile.write(fname.c_str(),1) ; - - - mu->decrRef(); } std::vector< int > -Mesh::getGroupFaceIds(std::string groupName) const +Mesh::getFaceGroupIds(std::string groupName, bool isBoundaryGroup) const { - if( std::find(_faceGroupNames.begin(),_faceGroupNames.end(),groupName) == _faceGroupNames.end() ) + std::vector::const_iterator it = std::find(_faceGroupNames.begin(),_faceGroupNames.end(),groupName); + if( it == _faceGroupNames.end() ) { cout<<"Mesh::getGroupFaceIds Error : face group " << groupName << " does not exist"< result(0); - for(int i=0; i<_boundaryFaceIds.size(); i++) - { - vector< std::string > groupNames = getFace(_boundaryFaceIds[i]).getGroupNames(); - if( std::find(groupNames.begin(), groupNames.end(),groupName) != groupNames.end() ) - result.push_back(_boundaryFaceIds[i]); - } - return result; - } + { + return _faceGroupsIds[it-_faceGroupNames.begin()]; + } } std::vector< int > -Mesh::getGroupNodeIds(std::string groupName) const +Mesh::getNodeGroupIds(std::string groupName, bool isBoundaryGroup) const { - if( std::find(_nodeGroupNames.begin(),_nodeGroupNames.end(),groupName) == _nodeGroupNames.end() ) + std::vector::const_iterator it = std::find(_nodeGroupNames.begin(),_nodeGroupNames.end(),groupName); + if( it == _nodeGroupNames.end() ) { cout<<"Mesh::getGroupNodeIds Error : node group " << groupName << " does not exist"< result(0); - for(int i=0; i<_boundaryFaceIds.size(); i++) - { - vector< std::string > groupNames = getNode(_boundaryFaceIds[i]).getGroupNames(); - if( std::find(groupNames.begin(), groupNames.end(),groupName) != groupNames.end() ) - result.push_back(_boundaryFaceIds[i]); - } - return result; - } + return _nodeGroupsIds[it-_nodeGroupNames.begin()]; } void @@ -2046,6 +2097,10 @@ Mesh::setFaceGroupByIds(std::vector< int > faceIds, std::string groupName) { for(int i=0; i< faceIds.size(); i++) getFace(faceIds[i]).setGroupName(groupName); + + _faceGroupsIds.insert( _faceGroupsIds.end(),faceIds); + _faceGroupNames.insert(_faceGroupNames.end(), groupName); + _faceGroups.insert( _faceGroups.end(), NULL); } void @@ -2055,3 +2110,13 @@ Mesh::setNodeGroupByIds(std::vector< int > nodeIds, std::string groupName) getNode(nodeIds[i]).setGroupName(groupName); } +void Mesh::deleteMEDCouplingUMesh() +{ + if(_meshNotDeleted) + { + (_mesh.retn())->decrRef(); + _meshNotDeleted=false; + } + else + throw CdmathException("Mesh::deleteMEDCouplingMesh() : mesh is not loaded"); +}; -- 2.30.2