From eea86dac41fbf7bae9ee6948e45bc0f6e642377e Mon Sep 17 00:00:00 2001 From: michael Date: Thu, 13 Jan 2022 17:52:33 +0100 Subject: [PATCH] Updated IJK Mesh class --- CDMATH/IJKmesh/inc/IJKCell.hxx | 6 +- CDMATH/IJKmesh/inc/IJKMesh.hxx | 386 +++++++++-------- CDMATH/IJKmesh/src/IJKMesh.cxx | 765 ++++++++++++++++----------------- 3 files changed, 592 insertions(+), 565 deletions(-) mode change 100644 => 100755 CDMATH/IJKmesh/inc/IJKCell.hxx diff --git a/CDMATH/IJKmesh/inc/IJKCell.hxx b/CDMATH/IJKmesh/inc/IJKCell.hxx old mode 100644 new mode 100755 index ec1d7b4..1be92a2 --- a/CDMATH/IJKmesh/inc/IJKCell.hxx +++ b/CDMATH/IJKmesh/inc/IJKCell.hxx @@ -128,7 +128,7 @@ class IJKCell private: //---------------------------------------------------------------- /* - * The (i,j,k) coordinate of this cell. + * The (i,j,k) coordinates of this cell. */ std::vector< int > _IJKCoords; /* @@ -137,14 +137,14 @@ class IJKCell int _cellIndex;// necessaire ??? /* - * The coordinate of barycenter the cell. + * The coordinate of the barycenter the cell. */ Point _point ;// necessaire ??? se déduit des cordonnées IJK /* * The number of cells surrounding this cell. */ - int _numberOfCells ; + int _numberOfCells ;//Dépend si on est une cellule frontière /* * The indices of cells surrounding this Node. diff --git a/CDMATH/IJKmesh/inc/IJKMesh.hxx b/CDMATH/IJKmesh/inc/IJKMesh.hxx index d2182c5..9690f46 100644 --- a/CDMATH/IJKmesh/inc/IJKMesh.hxx +++ b/CDMATH/IJKmesh/inc/IJKMesh.hxx @@ -10,24 +10,35 @@ /** * IJKMesh class is defined by - * - case 1: file name of mesh med file (MEDCouplingIMesh) - * - case 2: 1D cartesian, xmin and xmax and number of cells - * - case 3: 2D cartesian, xmin, xmax, ymin and ymax and numbers of cells in x direction and y direction - * - case 4: 3D cartesian, xmin, xmax, ymin, ymax, zmin and zmax and numbers of cells in x direction, y direction and z direction + * - case 1: file name of mesh med file (MEDCouplingCMesh) + * - case 2: 1D cartesian, xmin and xmax and number of cells (MEDCouplingIMesh) + * - case 3: 2D cartesian, xmin, xmax, ymin and ymax and numbers of cells in x direction and y direction (MEDCouplingIMesh) + * - case 4: 3D cartesian, xmin, xmax, ymin, ymax, zmin and zmax and numbers of cells in x direction, y direction and z direction (MEDCouplingIMesh) * - * Mesh structure - * nx , ny , nz cell structure - * nx+1, ny+1, nz+1 node structure + * Mesh cell and node structures are stored in a single MEDCouplingStructuredMesh _mesh + * nx * ny * nz cell structure + * (nx+1) * (ny+1) * (nz+1) node structure * * The face structure is more tricky because it depends on the dimension - * if dim=1, a single grid : nx+1 - * if dim=2, union of two grids : nx,ny+1 and nx,ny+1 - * if dim=3, union of three grids : nx,ny,nz+1, nx,ny+1,nz and nx+1,ny,nz + * if dim=1, a single grid : + * nx+1 nodes + * if dim=2, union of two grids : + * (nx+1)*ny nodes (faces orthogonal to x-axis), origin (0,dy/2) + * nx*(ny+1) nodes (faces orthogonal to x-axis), origin (dx/2,0) + * if dim=3, union of three grids : + * nx*ny*(nz+1) (faces orthogonal to x-axis), origin (0,dy/2,dz/2) + * nx*(ny+1)*nz (faces orthogonal to y-axis), origin (dx/2,0,dz/2) + * (nx+1)*ny*nz (faces orthogonal to z-axis), origin (dx/2,dy/2,0) + * + * Mesh face structures are stored in a vector of meshes _faceMeshes of size meshDim + * The face centers are located on the nodes of the meshes in _faceMeshes + * The origins of the meshes in _faceMeshes are shifted from the origin of _mesh by dx/2 on x-axis, dy/2 on y-axis and dz/2 on z-axis * * - number of nodes surounding each cell : known constant : 2*_Ndim * - number of faces surounding each cell : known constant : 2 _Ndim * - normal vectors surrounding each cell * - measure of each cell : known constant : dx*dy*dz + * - measures of faces : known constants : dx*dy, dx*dz and dy*dz */ namespace MEDCoupling @@ -48,12 +59,12 @@ class IJKNode; class IJKCell; class IJKFace; -typedef enum +enum EntityType { CELLS = 0, NODES = 1, FACES = 2, - } EntityType; + }; #include #include @@ -64,19 +75,19 @@ class IJKMesh public: //---------------------------------------------------------------- /** - * default constructor + * \brief default constructor */ Mesh ( void ) ; /** - * constructor with data to load a structured MEDCouplingIMesh + * \brief constructor with data to load a structured MEDCouplingIMesh * @param filename name of structured 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 ) ; /** - * constructor with data for a regular 1D grid + * \brief constructor with data for a regular 1D grid * @param xmin : minimum x * @param xmax : maximum x * @param nx : Number of cells in x direction @@ -84,7 +95,7 @@ public: //---------------------------------------------------------------- Mesh( double xmin, double xmax, int nx, std::string meshName="MESH1D_Regular_Grid" ) ; /** - * constructor with data for a regular 2D grid + * \brief constructor with data for a regular 2D grid * @param xmin : minimum x * @param xmax : maximum x * @param ymin : minimum y @@ -95,7 +106,7 @@ public: //---------------------------------------------------------------- Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, std::string meshName="MESH2D_Regular_Rectangle_Grid") ; /** - * constructor with data for a regular 3D grid + * \brief constructor with data for a regular 3D grid * @param xmin : minimum x * @param xmax : maximum x * @param ymin : minimum y @@ -108,156 +119,220 @@ public: //---------------------------------------------------------------- */ Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz, std::string meshName="MESH3D_Regular_Cuboid_Grid") ; - Mesh( const MEDCoupling::MEDCouplingIMesh* mesh ) ; + Mesh( const MEDCoupling::MEDCouplingMesh* mesh ) ; + /** + * \brief Computes the global cell number from its IJK position + * @param int i : cell index along x-axis + * @param int j : cell index along y-axis + * @param int k : cell index along z-axis + * @return global cell number + */ + int getCellNumber(int i, int j, int k) const { return _mesh->getCellIdFromPos ( i, j, k); }; + /** + * \brief Computes the global node number from its IJK position + * @param int i : node index along x-axis + * @param int j : node index along y-axis + * @param int k : node index along z-axis + * @return global node number + */ + int getNodeNumber(int i, int j, int k) const { return _mesh->getNodeIdFromPos ( i, j, k); }; + /** + * \brief Computes the global face number from its IJK position + * @param int face grid number corresponding to the direction of the face normal : 0->x, 1->y, 2->z + * @param int i : node index along x-axis + * @param int j : node index along y-axis + * @param int k : node index along z-axis + * @return global node number + */ + int getFaceNumber(int face_grid_number, int i, int j, int k) const { return _faceMeshes[face_grid_number]->getNodeIdFromPos ( i, j, k); }; + /** + * \brief Computes the IJK position of a cell from its index number + * @param int cellId : global cell index + * @return vector of i,j,k indices of the cell + */ + std::vector< int > getIJKCellCoordinates(int cellId) const { return _mesh->getLocationFromCellId(cellId); }; + /** + * \brief Computes the IJK position of a node from its index number + * @param int nodeId : global node index + * @return vector of i,j,k indices of the node + */ + std::vector< int > getIJKNodeCoordinates(int nodeId) const { return _mesh->getLocationFromNodeId(nodeId); }; + /** + * \brief Computes the IJK position of a face from its index number + * @param int face grid number corresponding to the direction of the face normal : 0->x, 1->y, 2->z + * @param int faceId : global face index + * @return vector of i,j,k indices of the node + */ + std::vector< int > getIJKFaceCoordinates( int face_grid_number, int faceId) const { return _faceMeshes[face_grid_number]->getLocationFromNodeId(faceId); }; + /** + * \brief Computes the indices of nodes surrounding a given cell + * @param int cellId : global cell index + * @return vector of indices of the nodes surrounding the cell + */ + std::vector< mcIdType > getNodeIdsOfCell(mcIdType cellId) const { std::vector< mcIdType > conn; _mesh_>getNodeIdsOfCell(mcIdType cellId, conn) ; return conn; }; + /** + * \brief Computes the coordinates of a node + * @param int nodeId : global node index + * @return vector of components of the node coordinates + */ + std::vector< double > getNodeCoordinates (mcIdType nodeId) const { std::vector< double > coo; _mesh_>getCoordinatesOfNode(mcIdType nodeId, coo) ; return coo; }; + /** + * \brief Computes the coordinates of a face center of mass + * @param int face grid number corresponding to the direction of the face normal : 0->x, 1->y, 2->z + * @param int faceId : global face index + * @return vector of components of the face coordinates + */ + std::vector< double > getFaceCoordinates ( int face_grid_number, mcIdType faceId) const { std::vector< double > coo; _faceMeshes[face_grid_number]>getCoordinatesOfNode(mcIdType faceId, coo) ; return coo; }; + /** + * \brief Computes the isobarycenter of a cell + * @param int cellId : cell number + */ + std::vector< double > getCellCenterCoordinates (mcIdType cellId) const ; + + double getMeasureOfAnyCell () const; + /** - * constructor with data + * \brief constructor with data * @param filename : file name of structured 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 ) ; /** - * constructor by copy + * \brief constructor by copy * @param mesh : The Mesh object to be copied */ Mesh ( const IJKMesh & mesh ) ; /** - * destructor + * \brief destructor */ ~Mesh( void ) ; /** - * return mesh name + * \brief return mesh name * @return _name */ - std::string getName( void ) const ; + std::string getName( void ) const { return _mesh->getName (); }; /** - * return Space dimension - * @return _spaceDim + * \brief return Space dimension + * @return spaceDim */ - int getSpaceDimension( void ) const ; + int getSpaceDimension( void ) const { return _mesh->getSpaceDimension (); }; /** - * return Mesh dimension - * @return _meshDim + * \brief return Mesh dimension + * @return meshDim */ - int getMeshDimension( void ) const ; + int getMeshDimension( void ) const { return _mesh->getMeshDimension (); }; /** - * return the number of nodes in this mesh + * \brief return the number of nodes in this mesh * @return _numberOfNodes */ int getNumberOfNodes ( void ) const ; /** - * return the number of faces in this mesh + * \brief return the number of faces in this mesh * @return _numberOfFaces */ int getNumberOfFaces ( void ) const ; /** - * return the number of cells in this mesh + * \brief return the number of cells in this mesh * @return _numberOfCells */ int getNumberOfCells ( void ) const ; /** - * return The cell i in this mesh - * @return _cells[i] - */ - IJKCell& getCell ( int i ) const ; - - /** - * return The face i in this mesh - * @return _faces[i] - */ - IJKFace& getFace ( int i ) const ; - - /** - * return The node i in this mesh - * @return _nodes[i] + * \brief return the number of edges in this mesh + * @return _numberOfEdges */ - IJKNode& getNode ( int i ) const ; + int getNumberOfEdges ( void ) const ; /** - * return number of cell in x direction (structured mesh) - * return _nX + * \brief return number of cell in x direction (structured mesh) + * @return _nX */ int getNx( void ) const ; /** - * return number of cell in y direction (structured mesh) - * return _nY + * \brief return number of cell in y direction (structured mesh) + * @return _nY */ int getNy( void ) const ; /** - * return number of cell in z direction (structured mesh) - * return _nZ + * \brief return number of cell in z direction (structured mesh) + * @return _nZ */ int getNz( void ) const ; - double getXMin( void ) const ; - - double getXMax( void ) const ; - - double getYMin( void ) const ; - - double getYMax( void ) const ; - - double getZMin( void ) const ; - - double getZMax( void ) const ; - - std::vector getDXYZ() const ; - std::vector getCellGridStructure() const; std::vector getNodeGridStructure() const; + std::vector< vector< int > > getFaceGridStructures() const; + /** - * surcharge operator = + * \brief overload operator = * @param mesh : The Mesh object to be copied */ const Mesh& operator= ( const Mesh& mesh ) ; /** - * return the mesh MEDCoupling - * return _mesh + * \brief return the MEDCouplingStructuredMesh + * @return _mesh */ - MEDCoupling::MCAuto getMEDCouplingIMesh ( void ) const ; + MEDCoupling::MCAuto getMEDCouplingStructuredMesh ( void ) const ; /** - * return the list of face group names - * return _faaceGroupNames + * \brief return the MEDCouplingStructuredMesh + * @return _faceMesh + */ + std::vector< MEDCoupling::MCAuto > getFaceMeshes const; + + /** + * \brief return the list of face group names + * @return _faceGroupNames */ std::vector getNameOfFaceGroups( void ) const ; /** - * return the list of node group names - * return _nodeGroupNames + * \brief return the list of node group names + * @return _nodeGroupNames */ std::vector getNameOfNodeGroups( void ) const ; /** - * write mesh in the VTK format + * \brief write the cell mesh in the VTK format */ void writeVTK ( const std::string fileName ) const ; /** - * write mesh in the MED format + * \brief write the cell and face meshes in the VTK format */ - void writeMED ( const std::string fileName ) const ; + void writeVTKAll meshes ( const std::string fileName ) const ; + + /** + * \brief write the cell mesh in the MED format + */ + void writeMED ( const std::string fileName, bool fromScratch = true ) const ; + + /** + * \brief write the cell and face meshes in the MED format + */ + void writeMEDAllMeshes ( const std::string fileName, bool fromScratch = true ) const ; /* * Functions to manage periodic boundary condition in square/cubic geometries */ - void setPeriodicFaces() ; + void setPeriodicFaces(bool check_groups= false, bool use_central_inversion=false) ; int getIndexFacePeriodic(int indexFace, bool check_groups= false, bool use_central_inversion=false); - void setBoundaryNodes(); + void setBoundaryNodesFromFaces(); + std::map getIndexFacePeriodic( void ) const; bool isIndexFacePeriodicSet() const ; bool isBorderNode(int nodeid) const ; @@ -266,133 +341,100 @@ public: //---------------------------------------------------------------- bool isQuadrangular() const ; bool isHexahedral() const ; bool isStructured() const ; - std::string getElementTypes() const ; - + + // epsilon used in mesh comparisons + 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 _mesh()->checkFastEquivalWith(m.getMEDCouplingStructuredMesh(),_epsilon);}; + // Deep comparison of two meshes to see if they are identical Except for their names and units + bool isEqualWithoutConsideringStr( Mesh m) const { return _mesh->isEqualWithoutConsideringStr(m.getMEDCouplingStructuredMesh(),_epsilon);}; + + std::vector< std::string > getElementTypesNames() const ; /** - * 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(); + double minRatioVolSurf() const{ return _cellMeasure / *max_element(begin(_faceMeasures), end(_faceMeasures));}; /** - * Return the maximum number of neighbours around an element (cells around a cell or nodes around a node) + * \brief Compute the maximum number of neighbours around an element (cells around a cell or nodes around a node) */ int getMaxNbNeighbours(EntityType type) const; /** - * return the measure of all cells (length in 1D, surface in 2D or volume in 3D) - * @return _measureOfCells + * return the measure of a cell (length in 1D, surface in 2D or volume in 3D) + * @return _cellMeasure */ - double getCellsMeasure ( void ) const ; + double getCellMeasure ( ) const { return _cellMeasure;}; + + /** + * return the measure of a cell (length in 1D, surface in 2D or volume in 3D) + * @return _faceMeasures + */ + std::vector< double > getFaceMeasures ( ) const { return _faceMeasures;}; /** * return normal vectors around each cell */ - std::vector< Vector > getNormalVectors (void) const ; + std::vector< Vector > getNormalVectors (void) const { return _faceNormals;}; private: //---------------------------------------------------------------- - std::string _name; - - /** - * Space dimension - */ - int _spaceDim ; - - /** - * Mesh dimension + /** + * \brief The cell mesh + * Question : can _mesh be const since no buildUnstructured is applied? */ - int _meshDim ; - - /* - * Structured mesh parameters - */ - - double _xMin; - - double _xMax; - - double _yMin; - - double _yMax; - - double _zMin; - - double _zMax; - - std::vector _nxyz;//Number of cells in each direction - - std::vector _dxyz;//lenght depth and height of each cell - - /* - * The number of nodes in this mesh. + MEDCoupling::MCAuto _mesh;//This is either a MEDCouplingIMesh (creation from scratch) or a MEDCouplingCMesh (loaded from a file) + /** + * \brief The face meshes */ - int _numberOfNodes; - - /* - * The numbers of faces in this mesh. + std::vector< MEDCoupling::MCAuto > _faceMeshes;//These are MEDCouplingIMesh + /** + * \brief Generate the face meshes from the cell mesh */ - int _numberOfFaces; + void setFaceMeshes(); - /* - * The number of cells in this mesh. - */ - int _numberOfCells; + double _cellMeasure; + std::vector< double > _faceMeasures; + std::vector< std::vector< double > > _faceNormals; - /* - * The names of face groups. + /* Boundary data */ + /** + * \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 mesh MEDCouplingIMesh - */ - MEDCoupling::MCAuto _mesh; - std::vector< INTERP_KERNEL::NormalizedCellType > _eltsTypes;//List of cell types contained in the mesh + + /** + * \brief List of boundary faces + */ + std::vector< int > _boundaryFaceIds; + /** + * \brief List of boundary nodes + */ + std::vector< int > _boundaryNodeIds; - /* - * 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; - /* List of boundary faces*/ - std::vector< int > _boundaryFaceIds; - /* List of boundary nodes*/ - std::vector< int > _boundaryNodeIds; - - /* - * The number of nodes surrounding each cell. - */ - int _numberOfNodesSurroundingCells ; - - /* - * The number of faces surounding each cell. - */ - int _numberOfFacesSurroundingCells ; - - /* - * The length or surface or volume of each cell. - */ - double _measureOfCells ; - - /* - * The normal vectors surrounding each cell. - */ - std::vector< Vector > _normalVectorsAroundCells ; + double _epsilon; }; #endif /* IJKMESH_HXX_ */ diff --git a/CDMATH/IJKmesh/src/IJKMesh.cxx b/CDMATH/IJKmesh/src/IJKMesh.cxx index 5559c2d..1a6bdeb 100644 --- a/CDMATH/IJKmesh/src/IJKMesh.cxx +++ b/CDMATH/IJKmesh/src/IJKMesh.cxx @@ -5,130 +5,127 @@ * Authors: CDMATH */ +#include +#include +#include +#include +#include + #include "IJKMesh.hxx" #include "IJKNode.hxx" #include "IJKCell.hxx" #include "IJKFace.hxx" +#include "MEDCouplingIMesh.hxx" +#include "MEDCouplingUMesh.hxx" #include "MEDFileCMesh.hxx" #include "MEDLoader.hxx" -#include "MEDCouplingIMesh.hxx" #include "CdmathException.hxx" -#include -#include -#include -#include -#include - using namespace MEDCoupling; using namespace std; //---------------------------------------------------------------------- -Mesh::Mesh( void ) +IJKMesh::IJKMesh( void ) //---------------------------------------------------------------------- { _mesh=NULL; - _spaceDim = 0 ; - _meshDim = 0 ; - _numberOfNodes = 0; - _numberOfFaces = 0; - _numberOfCells = 0; - _xMin=0.; - _xMax=0.; - _yMin=0.; - _yMax=0.; - _zMin=0.; - _zMax=0.; - _nxyz.resize(0); - _dxyz.resize(0.); + _faceMeshes=std::vector< MEDCoupling::MCAuto >(0) + _measureField=NULL; _faceGroupNames.resize(0); _faceGroups.resize(0); _nodeGroupNames.resize(0); _nodeGroups.resize(0); _indexFacePeriodicSet=false; _name=""; + _epsilon=1e-6; } //---------------------------------------------------------------------- -Mesh::~Mesh( void ) +IJKMesh::~IJKMesh( void ) //---------------------------------------------------------------------- { + _measureField->decrRef(); } std::string -Mesh::getName( void ) const +IJKMesh::getName( void ) const { return _name; } -Mesh::Mesh( const MEDCoupling::MEDCouplingIMesh* mesh ) +IJKMesh::IJKMesh( const MEDCoupling::MEDCouplingStructuredMesh* 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]; - } + _measureField = mesh->getMeasureField(true); + + //_mesh=mesh if _mesh is declared const + _mesh=mesh->clone(false);//No deep copy : it is assumed node coordinates and cell connectivity will not change + setFaceMeshes(); +} + +void +IJKMesh::setFaceMeshes() +{ + int meshDim = _mesh->getMeshDimension(); + + _faceMeshes=std::vector< MEDCoupling::MCAuto >( meshDim ) + + std:vector< int > nodeStr = _mesh->getNodeGridStructure(); + std:vector< int > dxyz = _mesh->getDXYZ(); + std:vector< double > origin = _mesh->getOrigin(); - double *originPtr = new double[_spaceDim]; - double *dxyzPtr = new double[_spaceDim]; - int *nodeStrctPtr = new int[_spaceDim]; + double originPtr[meshDim]; + double dxyzPtr[meshDim]; + mcIdType nodeStrctPtr[meshDim]; - for(int i=0;i<_spaceDim;i++) + /* Prepare the creation of face meshes, and the filling of face normals, face measures and cell measure */ + for(int i=0; i > (meshDim, std:vector< double >(meshDim,0)); + /* Creation of face meshes, and filling of face normals, face measures and cell measure */ + for(int i=0; igetName()+"_faces_"+std::to_string(i), + _mesh->getMeshDimension(), + nodeStrctPtr, + nodeStrctPtr+meshDim, + originPtr, + originPtr+meshDim, + dxyzPtr, + dxyzPtr+meshDim); + for(int j=0; jgetMeasureField(true); + _faceMeshes=m.getFaceMeshes(); _faceGroupNames = m.getNameOfFaceGroups() ; _faceGroups = m.getFaceGroups() ; _nodeGroupNames = m.getNameOfNodeGroups() ; @@ -137,206 +134,101 @@ Mesh::Mesh( const IJKMesh& m ) if(_indexFacePeriodicSet) _indexFacePeriodicMap=m.getIndexFacePeriodic(); - MCAuto m1=m.getMEDCouplingIMesh()->deepCopy(); - _mesh=m1; + //_mesh=m if _mesh is declared const + _mesh=m.getMEDCouplingIMesh()->clone(false); } //---------------------------------------------------------------------- -Mesh::Mesh( const std::string filename, int meshLevel ) +IJKMesh::IJKMesh( const std::string filename, const std::string & meshName, int meshLevel ) //---------------------------------------------------------------------- { - readMeshMed(filename, meshLevel); + readMeshMed(filename, meshName, meshLevel); } //---------------------------------------------------------------------- void -Mesh::readMeshMed( const std::string filename, const int meshLevel) +IJKMesh::readMeshMed( const std::string filename, const std::string & meshName , int meshLevel ) //---------------------------------------------------------------------- { - MEDFileCMesh *m=MEDFileCMesh::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) + MEDFileCMesh * m;//Here we would like a MEDFileStructuredMesh but that class does not exist + + if( meshName == "" ) + m=MEDFileCMesh::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=MEDFileCMesh::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()); - _meshDim=_mesh->getMeshDimension(); - _spaceDim=_mesh->getSpaceDimension(); - _name=_mesh->getName(); + _epsilon=1e-6; _indexFacePeriodicSet=false; - MEDCoupling::MEDCouplingIMesh* structuredMesh = dynamic_cast (_mesh.retn()); - if(structuredMesh) - { - _dxyz=structuredMesh->getDXYZ(); - _nxyz=structuredMesh->getCellGridStructure(); - double* Box0=new double[2*_spaceDim]; - structuredMesh->getBoundingBox(Box0); - - _xMin=Box0[0]; - _xMax=Box0[1]; - std::cout<<"nx= "<<_nxyz[0]; - if (_spaceDim>=2) - { - _yMin=Box0[2]; - _yMax=Box0[3]; - std::cout<<", "<<"ny= "<<_nxyz[1]; - } - if (_spaceDim>=3) - { - _zMin=Box0[4]; - _zMax=Box0[5]; - std::cout<<", "<<"nz= "<<_nxyz[2]; - } - } - else - throw CdmathException("Mesh::readMeshMed med file does not contain a structured MedcouplingIMesh mesh"); + _measureField = mesh->getMeasureField(true); + + setFaceMeshes(); cout<=xmax) - throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx) : xmin >= xmax"); + throw CdmathException("IJKMesh::IJKMesh( double xmin, double xmax, int nx) : xmin >= xmax"); double dx = (xmax - xmin)/nx ; - _spaceDim = 1 ; - _meshDim = 1 ; - _name=meshName; + double meshDim = 1 ; + _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]; - int *nodeStrctPtr = new int[_spaceDim]; + double originPtr[meshDim]; + double dxyzPtr[meshDim]; + mcIdType nodeStrctPtr[meshDim]; originPtr[0]=xmin; nodeStrctPtr[0]=nx+1; dxyzPtr[0]=dx; - _mesh=MEDCouplingIMesh::New(meshName, - _spaceDim, + _mesh=MEDCouplingIIJKMesh::New(meshName, + meshDim, nodeStrctPtr, - nodeStrctPtr+_spaceDim, + nodeStrctPtr+meshDim, originPtr, - originPtr+_spaceDim, + originPtr+meshDim, dxyzPtr, - dxyzPtr+_spaceDim); - delete [] originPtr; - delete [] dxyzPtr; - delete [] nodeStrctPtr; - - _numberOfCells = _mesh->getNumberOfCells() ; - - _numberOfNodes = _mesh->getNumberOfNodes() ; - - _numberOfFaces = _numberOfNodes; - + dxyzPtr+meshDim); + _measureField = _mesh->getMeasureField(true); + setFaceMeshes(); } //---------------------------------------------------------------------- -Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, std::string meshName) +IJKMesh::IJKMesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, std::string meshName) //---------------------------------------------------------------------- { 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"); + throw CdmathException("IJKMesh::IJKMesh( 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"); + throw CdmathException("IJKMesh::IJKMesh( 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.; - + throw CdmathException("IJKMesh::IJKMesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : ymin >= ymax"); double dx = (xmax - xmin)/nx ; double dy = (ymax - ymin)/ny ; - _spaceDim = 2 ; - _meshDim = 2 ; - _name=meshName; + double meshDim = 2 ; + _epsilon=1e-6; _indexFacePeriodicSet=false; - _nxyz.resize(_spaceDim); - _nxyz[0]=nx; - _nxyz[1]=ny; - - _dxyz.resize(_spaceDim); - _dxyz[0]=dx; - _dxyz[1]=dy; - double *originPtr = new double[_spaceDim]; - double *dxyzPtr = new double[_spaceDim]; - int *nodeStrctPtr = new int[_spaceDim]; + double originPtr[meshDim]; + double dxyzPtr[meshDim]; + mcIdType nodeStrctPtr[meshDim]; originPtr[0]=xmin; originPtr[1]=ymin; @@ -346,67 +238,40 @@ Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, dxyzPtr[1]=dy; _mesh=MEDCouplingIMesh::New(meshName, - _spaceDim, + meshDim, nodeStrctPtr, - nodeStrctPtr+_spaceDim, + nodeStrctPtr+meshDim, originPtr, - originPtr+_spaceDim, + originPtr+meshDim, dxyzPtr, - dxyzPtr+_spaceDim); - - delete [] originPtr; - delete [] dxyzPtr; - delete [] nodeStrctPtr; - - _numberOfCells = _mesh->getNumberOfCells() ; - - _numberOfNodes = _mesh->getNumberOfNodes() ; - - _numberOfFaces = nx*(ny+1)+ny*(nx+1); - + dxyzPtr+meshDim); + _measureField = _mesh->getMeasureField(true); + setFaceMeshes(); } //---------------------------------------------------------------------- -Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz, std::string meshName) +IJKMesh::IJKMesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz, std::string meshName) //---------------------------------------------------------------------- { if(nx<=0 || ny<=0 || nz<=0) - throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : nx <= 0 or ny <= 0 or nz <= 0"); + throw CdmathException("IJKMesh::IJKMesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : nx <= 0 or ny <= 0 or nz <= 0"); if(xmin>=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"); + throw CdmathException("IJKMesh::IJKMesh( 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"); + throw CdmathException("IJKMesh::IJKMesh( 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"); + throw CdmathException("IJKMesh::IJKMesh( 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; - _indexFacePeriodicSet=false; - _xMin=xmin; - _xMax=xmax; - _yMin=ymin; - _yMax=ymax; - _zMin=zmin; - _zMax=zmax; + double meshDim = 3; + _epsilon=1e-6; double dx = (xmax - xmin)/nx ; double dy = (ymax - ymin)/ny ; double dz = (zmax - zmin)/nz ; - _dxyz.resize(_spaceDim); - _dxyz[0]=dx; - _dxyz[1]=dy; - _dxyz[2]=dz; - - _nxyz.resize(_spaceDim); - _nxyz[0]=nx; - _nxyz[1]=ny; - _nxyz[2]=nz; - - double *originPtr = new double[_spaceDim]; - double *dxyzPtr = new double[_spaceDim]; - int *nodeStrctPtr = new int[_spaceDim]; + double originPtr[meshDim]; + double dxyzPtr[meshDim]; + mcIdType nodeStrctPtr[meshDim]; originPtr[0]=xmin; originPtr[1]=ymin; @@ -419,277 +284,397 @@ Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, dxyzPtr[2]=dz; _mesh=MEDCouplingIMesh::New(meshName, - _spaceDim, + meshDim, nodeStrctPtr, - nodeStrctPtr+_spaceDim, + nodeStrctPtr+meshDim, originPtr, - originPtr+_spaceDim, + originPtr+meshDim, dxyzPtr, - dxyzPtr+_spaceDim); - - delete [] originPtr; - delete [] dxyzPtr; - delete [] nodeStrctPtr; - - _numberOfCells = _mesh->getNumberOfCells() ; - - _numberOfNodes = _mesh->getNumberOfNodes() ; - - _numberOfFaces = nx*ny*(nz+1)+nx*nz*(ny+1)+ny*nz*(nx+1); - + dxyzPtr+meshDim); + _measureField = _mesh->getMeasureField(true); + setFaceMeshes(); } -//---------------------------------------------------------------------- -int -Mesh::getSpaceDimension( void ) const -//---------------------------------------------------------------------- +void +IJKMesh::setPeriodicFaces() { - return _spaceDim ; + _indexFacePeriodicSet=true; } -//---------------------------------------------------------------------- -int -Mesh::getMeshDimension( void ) const -//---------------------------------------------------------------------- +bool +IJKMesh::isBorderNode(int nodeid) const { - return _meshDim ; + return getNode(nodeid).isBorder(); } -std::vector -Mesh::getDXYZ() const +bool +IJKMesh::isTriangular() const { - return _dxyz; + return false; } -std::vector -Mesh::getCellGridStructure() const +bool +IJKMesh::isQuadrangular() const { - return _nxyz; + return _mesh->getMeshDimension()==2; } - -//---------------------------------------------------------------------- -int -Mesh::getNx( void ) const -//---------------------------------------------------------------------- +bool +IJKMesh::isHexahedral() const { - return _nxyz[0]; + return _mesh->getMeshDimension()==3; } - -//---------------------------------------------------------------------- -int -Mesh::getNy( void ) const -//---------------------------------------------------------------------- +bool +IJKMesh::isStructured() const { - if(_spaceDim < 2) - throw CdmathException("int double& Field::operator[ielem] : Ny is not defined in dimension < 2!"); - - return _nxyz[1]; + return true; } -//---------------------------------------------------------------------- -int -Mesh::getNz( void ) const -//---------------------------------------------------------------------- +std::string +IJKMesh::getElementTypes() const { - if(_spaceDim < 3) - throw CdmathException("int Mesh::getNz( void ) : Nz is not defined in dimension < 3!"); - - return _nxyz[2]; + if( _mesh->getMeshDimension()==1 ) + return "Segments "; + else if( _mesh->getMeshDimension()==2 ) + return "Quadrangles "; + else if( _mesh->getMeshDimension()==3 ) + return "Hexahedra "; + else + { + cout<< "Mesh " + _name + " does not have acceptable dimension. Mesh dimension is " << meshDim<getMeshDimension()]; + _mesh->getBoundingBox(Box0); + + return Box0[0] ; } //---------------------------------------------------------------------- double -Mesh::getXMax( void ) const +IJKMesh::getXMax( void ) const //---------------------------------------------------------------------- { - return _xMax ; + double Box0[2*_mesh->getMeshDimension()]; + _mesh->getBoundingBox(Box0); + + return Box0[1] ; } //---------------------------------------------------------------------- double -Mesh::getYMin( void ) const +IJKMesh::getYMin( void ) const //---------------------------------------------------------------------- { - return _yMin ; + if(_mesh->getMeshDimension()<2) + throw CdmathException("IJKMesh::getYMin : dimension should be >=2"); + + double Box0[2*_mesh->getMeshDimension()]; + _mesh->getBoundingBox(Box0); + + return Box0[2] ; } //---------------------------------------------------------------------- double -Mesh::getYMax( void ) const +IJKMesh::getYMax( void ) const //---------------------------------------------------------------------- { - return _yMax ; + if(_mesh->getMeshDimension()<2) + throw CdmathException("IJKMesh::getYMax : dimension should be >=2"); + + double Box0[2*_mesh->getMeshDimension()]; + _mesh->getBoundingBox(Box0); + + return Box0[3] ; } //---------------------------------------------------------------------- double -Mesh::getZMin( void ) const +IJKMesh::getZMin( void ) const //---------------------------------------------------------------------- { - return _zMin ; + if(_mesh->getMeshDimension()<3) + throw CdmathException("IJKMesh::getZMin : dimension should be 3"); + + double Box0[2*_mesh->getMeshDimension()]; + _mesh->getBoundingBox(Box0); + + return Box0[4] ; } //---------------------------------------------------------------------- double -Mesh::getZMax( void ) const +IJKMesh::getZMax( void ) const //---------------------------------------------------------------------- { - return _zMax ; + if(_mesh->getMeshDimension()<3) + throw CdmathException("IJKMesh::getZMax : dimension should be 3"); + + double Box0[2*_mesh->getMeshDimension()]; + _mesh->getBoundingBox(Box0); + + return Box0[5] ; } //---------------------------------------------------------------------- -MCAuto -Mesh::getMEDCouplingIMesh( void ) const +MCAuto +IJKMesh::getMEDCouplingStructuredMesh( void ) const //---------------------------------------------------------------------- { return _mesh ; } +std::vector< MEDCoupling::MCAuto > +IJKMesh::getFaceMeshes() const +{ + return _faceMeshes; +} //---------------------------------------------------------------------- int -Mesh::getNumberOfNodes ( void ) const +IJKMesh::getNumberOfNodes ( void ) const //---------------------------------------------------------------------- { - return _numberOfNodes ; + return _mesh->getNumberOfNodes (); } //---------------------------------------------------------------------- int -Mesh::getNumberOfCells ( void ) const +IJKMesh::getNumberOfCells ( void ) const //---------------------------------------------------------------------- { - return _numberOfCells ; + return _mesh->getNumberOfCells (); } -//---------------------------------------------------------------------- -int -Mesh::getNumberOfFaces ( void ) const -//---------------------------------------------------------------------- +std::vector +IJKMesh::getCellGridStructure() const { - return _numberOfFaces ; + return _mesh->getCellGridStructure(); } -//---------------------------------------------------------------------- -Cell& -Mesh::getCell ( int i ) const -//---------------------------------------------------------------------- +std::vector +IJKMesh::getNodeGridStructure() const { - return _cells[i] ; + return _mesh->getNodeGridStructure(); } -//---------------------------------------------------------------------- -Face& -Mesh::getFace ( int i ) const -//---------------------------------------------------------------------- +std::vector< vector > +IJKMesh::getFaceGridStructures() const { - return _faces[i] ; + std::vector< vector > result(_mesh->getmeshDimension()); + for(int i = 0; i<_mesh->getmeshDimension(); i++) + result[i] = _faceMeshes[i]->getNodeGridStructure(); + + return result; } //---------------------------------------------------------------------- -Node& -Mesh::getNode ( int i ) const +int +IJKMesh::getNumberOfFaces ( void ) const //---------------------------------------------------------------------- { - return _nodes[i] ; + std::vector NxNyNz = _mesh->getCellGridStructure(); + + switch( _mesh->getMeshDimension () ) + { + case 1: + return NxNyNz[0]+1; + case 2: + return NxNyNz[0]*(1+NxNyNz[1]) + NxNyNz[1]*(1 + NxNyNz[0]); + case 2: + return NxNyNz[0]*NxNyNz[1]*(1+NxNyNz[2]) + NxNyNz[0]*NxNyNz[2]*(1+NxNyNz[1]) + NxNyNz[1]*NxNyNz[2]*(1+NxNyNz[0]); + default + throw CdmathException("IJKMesh::getNumberOfFaces space dimension must be between 1 and 3"); + } } +int +IJKMesh::getNumberOfEdges ( void ) const +{ + std::vector NxNyNz = _mesh->getCellGridStructure(); + + switch( _mesh->getMeshDimension () ) + { + case 1: + return NxNyNz[0]; + case 2: + return NxNyNz[0]*(1+NxNyNz[1]) + NxNyNz[1]*(1 + NxNyNz[0]); + case 2: + return NxNyNz[0]*(1+NxNyNz[1])*(1+NxNyNz[2]) + NxNyNz[1]*(1 + NxNyNz[0])*(1+NxNyNz[2]) + NxNyNz[2]*(1 + NxNyNz[0])*(1+NxNyNz[1]); + default + throw CdmathException("IJKMesh::getNumberOfEdges space dimension must be between 1 and 3"); + } +}; + vector -Mesh::getNameOfFaceGroups( void ) const +IJKMesh::getNameOfFaceGroups( void ) const { return _faceGroupNames; } vector -Mesh::getFaceGroups( void ) const +IJKMesh::getFaceGroups( void ) const { return _faceGroups; } vector -Mesh::getNameOfNodeGroups( void ) const +IJKMesh::getNameOfNodeGroups( void ) const { return _nodeGroupNames; } vector -Mesh::getNodeGroups( void ) const +IJKMesh::getNodeGroups( void ) const { return _nodeGroups; } //---------------------------------------------------------------------- const IJKMesh& -Mesh::operator= ( const IJKMesh& mesh ) +IJKMesh::operator= ( const IJKMesh& mesh ) //---------------------------------------------------------------------- { - _spaceDim = mesh.getSpaceDimension() ; - _meshDim = mesh.getMeshDimension() ; - _name = mesh.getName(); - _numberOfNodes = mesh.getNumberOfNodes(); - _numberOfFaces = mesh.getNumberOfFaces(); - _numberOfCells = mesh.getNumberOfCells(); + _epsilon=mesh.getComparisonEpsilon(); _indexFacePeriodicSet= mesh.isIndexFacePeriodicSet(); if(_indexFacePeriodicSet) _indexFacePeriodicMap=mesh.getIndexFacePeriodic(); - _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() ; - MCAuto m1=mesh.getMEDCouplingIMesh()->deepCopy(); - _mesh=m1; + //_mesh=mesh if _mesh is declared const + _mesh=mesh.getMEDCouplingStructuredMesh()->clone(false); + + _faceMeshes=mesh.getFaceMeshes()); return *this; } - -bool Mesh::isIndexFacePeriodicSet() const + +bool IJKMesh::isIndexFacePeriodicSet() const { return _indexFacePeriodicSet; } //---------------------------------------------------------------------- double -Mesh::minRatioVolSurf() +IJKMesh::minRatioVolSurf() { return dx_min; } int -Mesh::getMaxNbNeighbours(EntityType type) const +IJKMesh::getMaxNbNeighbours(EntityType type) const { - return 2*_meshDim; + return 2*_mesh->getMeshDimension(); } //---------------------------------------------------------------------- void -Mesh::writeVTK ( const std::string fileName ) const +IJKMesh::writeVTK ( const std::string fileName ) const //---------------------------------------------------------------------- { string fname=fileName+".vtu"; _mesh->writeVTK(fname.c_str()) ; } +//---------------------------------------------------------------------- +void +IJKMesh::writeVTKAllMeshes ( const std::string fileName ) const +//---------------------------------------------------------------------- +{ + string fname=fileName+".vtu"; + _mesh->writeVTK(fname.c_str()) ; + for(int i=0; i< _mesh->getMeshDimension(); i++) + _faceMeshes[i]->writeVTK(fname.c_str()) ; +} + +//---------------------------------------------------------------------- +void +IJKMesh::writeMED ( const std::string fileName, bool fromScratch ) const +//---------------------------------------------------------------------- +{ + string fname=fileName+".med"; + //Save cell mesh after checking mesh is imesh + const MEDCoupling::MEDCouplingIMesh* iMesh = dynamic_cast< const MEDCoupling::MEDCouplingIMesh* > ((const MEDCoupling::MEDCouplingStructuredMesh*) _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); +} //---------------------------------------------------------------------- void -Mesh::writeMED ( const std::string fileName ) const +IJKMesh::writeMEDAllMeshes ( const std::string fileName, bool fromScratch ) const //---------------------------------------------------------------------- { + //Save cell mesh after checking mesh is imesh string fname=fileName+".med"; - MEDCoupling::WriteCMesh(fname.c_str(),_mesh,true); + const MEDCoupling::MEDCouplingIMesh* iMesh = dynamic_cast< const MEDCoupling::MEDCouplingIMesh* > ((const MEDCoupling::MEDCouplingStructuredMesh*) _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); + + //Save face meshes after checking mesh is imesh + std::vector< const MEDCoupling::MEDCouplingIMesh* > iMeshes(_mesh->getMeshDimension()); + for(int i=0; i< _mesh->getMeshDimension(); i++) + { + const MEDCoupling::MEDCouplingIMesh* iMeshes[i] = dynamic_cast< const MEDCoupling::MEDCouplingIMesh* > ((const MEDCoupling::MEDCouplingStructuredMesh*) _faceMeshes[i]); + if(iMesh)//medcouplingimesh : Use convertToCartesian in order to write mesh + MEDCoupling::WriteMesh(fname.c_str(),iMeshes[i]->convertToCartesian(), fromScratch); + else//medcouplingcmesh : save directly + MEDCoupling::WriteMesh(fname.c_str(), _faceMeshes[i], fromScratch); + } +} + +std::vector< double > +IJKMesh::getCellCenterCoordinates (mcIdType cellId) const +{ + std::vector< double > result(_mesh->getMeshDimension(),0), coo_node; + std::vector< mcIdType > conn=getNodeIdsOfCell(mcIdType cellId); + int nbNodes=conn.size(); + + for (int i = 0; i< nbNodes; i++) + { + coo_node = getNodeCoordinates (conn[i]); + for(int j=0; j< _mesh->getMeshDimension(); j++) + result[j]+=coo_node[j]; + } + for(int j=0; j< _mesh->getMeshDimension(); j++) + result[j]/=nbNodes; + return result; +}; + +//---------------------------------------------------------------------- +int +IJKMesh::getNx( void ) const +//---------------------------------------------------------------------- +{ + return _mesh->getCellGridStructure()[0]; +} - mu->decrRef(); +//---------------------------------------------------------------------- +int +IJKMesh::getNy( void ) const +//---------------------------------------------------------------------- +{ + if(_mesh->getMeshDimension () < 2) + throw CdmathException("int IJKMesh::getNy( void ) : Ny is not defined in dimension < 2!"); + else + return _mesh->getCellGridStructure()[1]; +} + +//---------------------------------------------------------------------- +int +IJKMesh::getNz( void ) const +//---------------------------------------------------------------------- +{ + if(_mesh->getMeshDimension () < 3) + throw CdmathException("int IJKMesh::getNz( void ) : Nz is not defined in dimension < 3!"); + else + return _mesh->getCellGridStructure()[2]; } -- 2.39.2