* Authors: CDMATH
*/
+#include <iostream>
+#include <stdio.h>
+#include <stdlib.h>
+#include <iterator>
+#include <algorithm>
+#include <cassert>
+
#include "Mesh.hxx"
#include "Node.hxx"
#include "Cell.hxx"
#include "Face.hxx"
-#include "MEDFileMesh.hxx"
-#include "MEDLoader.hxx"
-#include "MEDCouplingUMesh.hxx"
#include "MEDCouplingIMesh.hxx"
+#include "MEDCouplingUMesh.hxx"
#include "MEDCouplingFieldDouble.hxx"
-#include "CdmathException.hxx"
+#include "MEDFileMesh.hxx"
+#include "MEDLoader.hxx"
-#include <iostream>
-#include <stdio.h>
-#include <stdlib.h>
-#include <iterator>
-#include <algorithm>
-#include <string.h>
+#include "CdmathException.hxx"
using namespace MEDCoupling;
using namespace std;
//----------------------------------------------------------------------
{
_mesh=NULL;
+ _meshNotDeleted=false;
_cells=NULL;
_nodes=NULL;
_faces=NULL;
_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);
_faceGroupNames.resize(0);
_faceGroups.resize(0);
+ _faceGroupsIds.resize(0);
_nodeGroupNames.resize(0);
_nodeGroups.resize(0);
+ _nodeGroupsIds.resize(0);
_indexFacePeriodicSet=false;
_name="";
+ _epsilon=1e-6;
}
//----------------------------------------------------------------------
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
return _name;
}
-Mesh::Mesh( const MEDCoupling::MEDCouplingIMesh* mesh )
+Mesh::Mesh( MEDCoupling::MCAuto<const MEDCoupling::MEDCouplingMesh> 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;
+ _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
- double *originPtr = new double[_spaceDim];
- double *dxyzPtr = new double[_spaceDim];
- int *nodeStrctPtr = new int[_spaceDim];
+ MEDCoupling::MEDCouplingStructuredMesh* structuredMesh = dynamic_cast<MEDCoupling::MEDCouplingStructuredMesh*> (_mesh.retn());
+ if(structuredMesh)
+ {
+ _isStructured=true;
+ _nxyz=structuredMesh->getCellGridStructure();
+ }
+ else
+ _isStructured=false;
- 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);
- delete [] originPtr;
- delete [] dxyzPtr;
- delete [] nodeStrctPtr;
- delete [] Box0 ;
setMesh();
}
-Mesh::Mesh( const MEDCoupling::MEDCouplingUMesh* mesh )
+//----------------------------------------------------------------------
+Mesh::Mesh( const std::string filename, const std::string & meshName, int meshLevel)
+//----------------------------------------------------------------------
{
- _spaceDim=mesh->getSpaceDimension();
- _meshDim=mesh->getMeshDimension();
- _isStructured=false;
- double* Box0=new double[2*_spaceDim];
- mesh->getBoundingBox(Box0);
- _name=mesh->getName();
- _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];
- }
-
- _mesh=mesh->deepCopy();
- delete [] Box0 ;
- setMesh();
+ readMeshMed(filename, meshName, meshLevel);
}
//----------------------------------------------------------------------
_spaceDim = mesh.getSpaceDimension() ;
_meshDim = mesh.getMeshDimension() ;
_name=mesh.getName();
- _xMax=mesh.getXMax();
- _yMin=mesh.getYMin();
- _yMax=mesh.getYMax();
- _zMin=mesh.getZMin();
- _zMax=mesh.getZMax();
+ _epsilon=mesh.getComparisonEpsilon();
_isStructured=mesh.isStructured();
if(_isStructured)
- {
_nxyz = mesh.getCellGridStructure() ;
- _dxyz=mesh.getDXYZ();
- }
_numberOfNodes = mesh.getNumberOfNodes();
_numberOfFaces = mesh.getNumberOfFaces();
_numberOfCells = mesh.getNumberOfCells();
_numberOfEdges = mesh.getNumberOfEdges();
_faceGroupNames = mesh.getNameOfFaceGroups() ;
- _faceGroups = mesh.getFaceGroups() ;
+ _faceGroups = mesh.getMEDCouplingFaceGroups() ;
+ _faceGroupsIds = mesh.getFaceGroups() ;
_nodeGroupNames = mesh.getNameOfNodeGroups() ;
- _nodeGroups = mesh.getNodeGroups() ;
+ _nodeGroups = mesh.getMEDCouplingNodeGroups() ;
+ _nodeGroupsIds = 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();
_eltsTypes=mesh.getElementTypes();
_eltsTypesNames=mesh.getElementTypesNames();
- MCAuto<MEDCouplingMesh> m1=mesh.getMEDCouplingMesh()->deepCopy();
- _mesh=m1;
-}
+ MCAuto<MEDCouplingMesh> 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::Mesh( const std::string filename, int meshLevel )
-//----------------------------------------------------------------------
-{
- readMeshMed(filename, meshLevel);
+ _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());
_meshDim=_mesh->getMeshDimension();
_spaceDim=_mesh->getSpaceDimension();
_name=_mesh->getName();
+ _epsilon=1e-6;
_indexFacePeriodicSet=false;
- MEDCoupling::MEDCouplingIMesh* structuredMesh = dynamic_cast<MEDCoupling::MEDCouplingIMesh*> (_mesh.retn());
+ _meshNotDeleted=true;
+
+ MEDCoupling::MEDCouplingStructuredMesh* structuredMesh = dynamic_cast<MEDCoupling::MEDCouplingStructuredMesh*> (_mesh.retn());
if(structuredMesh)
{
_isStructured=true;
- _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;
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<<endl<< "Loaded file "<< filename<<endl;
- cout<<"Mesh name= "<<m->getName()<<", mesh dim="<< _meshDim<< ", space dim="<< _spaceDim<< ", nb cells= "<<getNumberOfCells()<< ", nb nodes= "<<getNumberOfNodes()<<endl;
+ cout<<"Mesh name = "<<m->getName()<<", mesh dim = "<< _meshDim<< ", space dim = "<< _spaceDim<< ", nb cells= "<<getNumberOfCells()<< ", nb nodes= "<<getNumberOfNodes()<<endl;
m->decrRef();
- mu->decrRef();
}
-void
-Mesh::setGroupAtFaceByCoords(double x, double y, double z, double eps, std::string groupName)
+//----------------------------------------------------------------------
+Mesh::Mesh( std::vector<double> points, std::string meshName )
+//----------------------------------------------------------------------
{
- int nbFace=getNumberOfFaces();
- bool flag=false;
- for (int iface=0;iface<nbFace;iface++)
- {
- double FX=_faces[iface].x();
- double FY=_faces[iface].y();
- double FZ=_faces[iface].z();
- if (abs(FX-x)<eps && abs(FY-y)<eps && abs(FZ-z)<eps)
- {
- _faces[iface].setGroupName(groupName);
- std::vector< int > nodesID= _faces[iface].getNodesId();
- int nbNodes = _faces[iface].getNumberOfNodes();
- for(int inode=0 ; inode<nbNodes ; inode++)
- _nodes[nodesID[inode]].setGroupName(groupName);
-
- flag=true;
- }
- }
- if (flag)
+ int nx=points.size();
+
+ if(nx<2)
+ throw CdmathException("Mesh::Mesh( vector<double> points, string meshName) : nx < 2, vector should contain at least two values");
+ int i=0;
+ while( i<nx-1 && points[i+1]>points[i] )
+ i++;
+ if( i!=nx-1 )
{
- _faceGroupNames.insert(_faceGroupNames.begin(),groupName);
- _nodeGroupNames.insert(_nodeGroupNames.begin(),groupName);
- //To do : update _faceGroups and _nodeGroups
+ //cout<< points << endl;
+ throw CdmathException("Mesh::Mesh( vector<double> points, string meshName) : vector values should be sorted");
}
-}
-
-void
-Mesh::setGroupAtPlan(double value, int direction, double eps, std::string groupName)
-{
- int nbFace=getNumberOfFaces();
- bool flag=false;
- for (int iface=0;iface<nbFace;iface++)
- {
- double cord=_faces[iface].getBarryCenter()[direction];
- if (abs(cord-value)<eps)
- {
- _faces[iface].setGroupName(groupName);
- std::vector< int > nodesID= _faces[iface].getNodesId();
- int nbNodes = _faces[iface].getNumberOfNodes();
- for(int inode=0 ; inode<nbNodes ; inode++)
- {
- _nodes[nodesID[inode]].setGroupName(groupName);
- }
-
- flag=true;
- }
- }
- if (flag)
+
+ _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; i<nx- 1 ; i++)
{
- _faceGroupNames.insert(_faceGroupNames.begin(),groupName);
- _nodeGroupNames.insert(_nodeGroupNames.begin(),groupName);
- //To do : update _faceGroups, _nodeGroups
+ nodal_con[0]=i;
+ nodal_con[1]=i+1;
+ mesh1d->insertNextCell(INTERP_KERNEL::NORM_SEG2, 2, nodal_con);
+ coords[i+1]=points[i + 1];
}
-}
+ mesh1d->finishInsertingCells();
-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<nbNodes ; inode++)
- {
- std::vector<int>::const_iterator it = std::find(_boundaryNodeIds.begin(),_boundaryNodeIds.end(),nodesID[inode]);
- if( it != _boundaryNodeIds.end() )
- _boundaryNodeIds.push_back(nodesID[inode]);
- }
- }
-}
+ DataArrayDouble * coords_arr = DataArrayDouble::New();
+ coords_arr->alloc(nx,1);
+ std::copy(coords,coords+nx,coords_arr->getPointer());
+ mesh1d->setCoords(coords_arr);
-std::map<int,int>
-Mesh::getIndexFacePeriodic( void ) const
-{
- return _indexFacePeriodicMap;
-}
+ delete [] coords;
+ coords_arr->decrRef();
-void
-Mesh::setPeriodicFaces(bool check_groups, bool use_central_inversion)
-{
- if(_indexFacePeriodicSet)
- return;
-
- double eps=1.E-10;
+ _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;
- 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)
- {
- 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)<eps || abs(x-xi)<eps )// 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)<eps ) // Case of a central inversion
- && fabs(my_face.getMeasure() - face_i.getMeasure())<eps
- && fabs(my_face.getXN() + face_i.getXN())<eps
- && fabs(my_face.getYN() + face_i.getYN())<eps
- && fabs(my_face.getZN() + face_i.getZN())<eps )
- {
- 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[_boundaryFaceIds[iface]];
- double xi=face_i.x();
- double yi=face_i.y();
- double zi=face_i.z();
- if ( ((abs(y-yi)<eps && abs(x-xi)<eps) || (abs(x-xi)<eps && abs(z-zi)<eps) || (abs(y-yi)<eps && abs(z-zi)<eps))// 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)<eps )// Case of a central inversion
- && fabs(my_face.getMeasure() - face_i.getMeasure())<eps
- && fabs(my_face.getXN() + face_i.getXN())<eps
- && fabs(my_face.getYN() + face_i.getYN())<eps
- && fabs(my_face.getZN() + face_i.getZN())<eps )
- {
- 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;
- }
- _indexFacePeriodicSet=true;
+ setMesh();
}
-int
-Mesh::getIndexFacePeriodic(int indexFace, bool check_groups, bool use_central_inversion)
+//----------------------------------------------------------------------
+Mesh::Mesh( double xmin, double xmax, int nx, std::string meshName )
+//----------------------------------------------------------------------
{
- if (!_faces[indexFace].isBorder())
- {
- cout<<"Pb with indexFace= "<<indexFace<<endl;
- throw CdmathException("Mesh::getIndexFacePeriodic: not a border face" );
- }
-
- if(!_indexFacePeriodicSet)
- setPeriodicFaces(check_groups, use_central_inversion);
+ 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");
- std::map<int,int>::const_iterator it = _indexFacePeriodicMap.find(indexFace);
- if( it != _indexFacePeriodicMap.end() )
- return it->second;
- else
- {
- cout<<"Pb with indexFace= "<<indexFace<<endl;
- throw CdmathException("Mesh::getIndexFacePeriodic: not a periodic face" );
- }
-}
+ double dx = (xmax - xmin)/nx ;
-bool
-Mesh::isBorderNode(int nodeid) const
-{
- return _nodes[nodeid].isBorder();
-}
+ _spaceDim = 1 ;
+ _meshDim = 1 ;
+ _name=meshName;
+ _epsilon=1e-6;
+ _indexFacePeriodicSet=false;
-bool
-Mesh::isBorderFace(int faceid) const
-{
- return _faces[faceid].isBorder();
-}
+ _nxyz.resize(_spaceDim);
+ _nxyz[0]=nx;
-std::vector< int >
-Mesh::getBoundaryFaceIds() const
-{
- return _boundaryFaceIds;
-}
+ double originPtr[_spaceDim];
+ double dxyzPtr[_spaceDim];
+ mcIdType nodeStrctPtr[_spaceDim];
-std::vector< int >
-Mesh::getBoundaryNodeIds() const
-{
- return _boundaryNodeIds;
-}
+ originPtr[0]=xmin;
+ nodeStrctPtr[0]=nx+1;
+ dxyzPtr[0]=dx;
-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;
+ _isStructured = true;
-std::vector< INTERP_KERNEL::NormalizedCellType >
-Mesh::getElementTypes() const
-{
- return _eltsTypes;
+ setMesh();
}
-std::vector< string >
-Mesh::getElementTypesNames() const
+//----------------------------------------------------------------------
+Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, int split_to_triangles_policy, std::string meshName)
+//----------------------------------------------------------------------
{
- 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 " <<endl;
- cout<< _eltsTypes[i]<<endl;
- throw CdmathException("Mesh::getElementTypes : recognised cell med types are NORM_POINT1, NORM_SEG2, NORM_TRI3, NORM_QUAD4, NORM_TETRA4, NORM_PYRA5, NORM_PENTA6, NORM_HEXA8, NORM_POLYGON, NORM_POLYHED");
- }
- }
- return result;
-}
+ 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");
-void
-Mesh::setGroups( const MEDFileUMesh* medmesh, MEDCouplingUMesh* mu)
-{
- //Searching for face groups
- vector<string> faceGroups=medmesh->getGroupsNames() ;
+ double dx = (xmax - xmin)/nx ;
+ double dy = (ymax - ymin)/ny ;
- for (unsigned int i=0;i<faceGroups.size();i++ )
- {
- string groupName=faceGroups[i];
- vector<int> 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<int>::iterator it = find(nonEmptyGrp.begin(), nonEmptyGrp.end(), -1);
- if (it != nonEmptyGrp.end())
- {
- cout<<"Boundary face group named "<< groupName << " found"<<endl;
- MEDCouplingUMesh *m=medmesh->getGroup(-1,groupName.c_str());
- mu->unPolyze();
- _faceGroups.insert(_faceGroups.begin(),m);
- _faceGroupNames.insert(_faceGroupNames.begin(),groupName);
- DataArrayDouble *baryCell = m->computeIsoBarycenterOfNodesPerCell();//computeCellCenterOfMass() ;
- const double *coorBary=baryCell->getConstPointer();
+ _spaceDim = 2 ;
+ _meshDim = 2 ;
+ _name=meshName;
+ _epsilon=1e-6;
+ _indexFacePeriodicSet=false;
+ _nxyz.resize(_spaceDim);
+ _nxyz[0]=nx;
+ _nxyz[1]=ny;
- int nbCellsSubMesh=m->getNumberOfCells();
- for (int ic(0), k(0); ic<nbCellsSubMesh; ic++, k+=_spaceDim)
- {
- vector<double> coorBaryXyz(3,0);
- for (int d=0; d<_spaceDim; d++)
- coorBaryXyz[d] = coorBary[k+d];
- Point p1(coorBaryXyz[0],coorBaryXyz[1],coorBaryXyz[2]) ;
+ double originPtr[_spaceDim];
+ double dxyzPtr[_spaceDim];
+ mcIdType nodeStrctPtr[_spaceDim];
- int flag=0;
- /* double min=INFINITY;
- Point p3;
- int closestFaceNb;*/
- for (int iface=0;iface<_numberOfFaces;iface++ )
- {
- Point p2=_faces[iface].getBarryCenter();
- /*if(groupName=="Wall" and p1.distance(p2)<min)
- {
- min=p1.distance(p2);
- p3=p2;
- closestFaceNb=iface;
- }*/
- if(p1.distance(p2)<1.E-10)
- {
- _faces[iface].setGroupName(groupName);
- flag=1;
- break;
- }
- }
- /* if(groupName=="Wall" and min>1.E-10)
- {
- cout.precision(15);
- cout<<"Subcell number "<< ic <<" Min distance to Wall = "<<min <<" p= "<<p1[0]<<" , "<<p1[1]<<" , "<<p1[2]<<endl;
- cout<<" attained at face "<< closestFaceNb << " p_min= "<<p3[0]<<" , "<<p3[1]<<" , "<<p3[2]<<endl;
- }*/
- if (flag==0)
- throw CdmathException("No face belonging to group " + groupName + " found");
- }
- baryCell->decrRef();
- //m->decrRef();
- }
- }
+ originPtr[0]=xmin;
+ originPtr[1]=ymin;
+ nodeStrctPtr[0]=nx+1;
+ nodeStrctPtr[1]=ny+1;
+ dxyzPtr[0]=dx;
+ dxyzPtr[1]=dy;
- //Searching for node groups
- vector<string> nodeGroups=medmesh->getGroupsOnSpecifiedLev(1) ;
+ _mesh=MEDCouplingIMesh::New(meshName,
+ _spaceDim,
+ nodeStrctPtr,
+ nodeStrctPtr+_spaceDim,
+ originPtr,
+ originPtr+_spaceDim,
+ dxyzPtr,
+ dxyzPtr+_spaceDim);
+ _meshNotDeleted=true;
+ _isStructured = true;
- for (unsigned int i=0;i<nodeGroups.size();i++ )
- {
- string groupName=nodeGroups[i];
- DataArrayInt * nodeGroup=medmesh->getNodeGroupArr( groupName );
- const int *nodeids=nodeGroup->getConstPointer();
+ if(split_to_triangles_policy==0 || split_to_triangles_policy==1)
+ {
+ _mesh=_mesh->buildUnstructured();//simplexize is not available for structured meshes
+ _mesh->simplexize(split_to_triangles_policy);
+ _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");
+ }
+
+ setMesh();
+}
- if(nodeids!=NULL)
- {
- cout<<"Boundary node group named "<< groupName << " found"<<endl;
+//----------------------------------------------------------------------
+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(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");
+ 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");
+ 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");
- _nodeGroups.insert(_nodeGroups.begin(),nodeGroup);
- _nodeGroupNames.insert(_nodeGroupNames.begin(),groupName);
+ _spaceDim = 3;
+ _meshDim = 3;
+ _name=meshName;
+ _epsilon=1e-6;
- int nbNodesSubMesh=nodeGroup->getNumberOfTuples();//nodeGroup->getNbOfElems();
+ double dx = (xmax - xmin)/nx ;
+ double dy = (ymax - ymin)/ny ;
+ double dz = (zmax - zmin)/nz ;
- DataArrayDouble *coo = mu->getCoords() ;
- const double *cood=coo->getConstPointer();
+ _nxyz.resize(_spaceDim);
+ _nxyz[0]=nx;
+ _nxyz[1]=ny;
+ _nxyz[2]=nz;
- for (int ic(0); ic<nbNodesSubMesh; ic++)
- {
- vector<double> 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]) ;
+ double originPtr[_spaceDim];
+ double dxyzPtr[_spaceDim];
+ mcIdType nodeStrctPtr[_spaceDim];
- int flag=0;
- for (int inode=0;inode<_numberOfNodes;inode++ )
- {
- Point p2=_nodes[inode].getPoint();
- if(p1.distance(p2)<1.E-10)
- {
- _nodes[inode].setGroupName(groupName);
- flag=1;
- break;
- }
- }
- if (flag==0)
- throw CdmathException("No node belonging to group " + groupName + " found");
- }
- }
- }
+ 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;
+
+ _mesh=MEDCouplingIMesh::New(meshName,
+ _spaceDim,
+ nodeStrctPtr,
+ nodeStrctPtr+_spaceDim,
+ originPtr,
+ originPtr+_spaceDim,
+ dxyzPtr,
+ dxyzPtr+_spaceDim);
+ _meshNotDeleted=true;
+ _isStructured = true;
+
+ if( split_to_tetrahedra_policy == 0 )
+ {
+ _mesh=_mesh->buildUnstructured();//simplexize is not available for structured meshes
+ _mesh->simplexize(INTERP_KERNEL::PLANAR_FACE_5);
+ _isStructured = false;
+ }
+ 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();
}
//----------------------------------------------------------------------
Mesh::setMesh( void )
//----------------------------------------------------------------------
{
- DataArrayInt *desc = DataArrayInt::New();
- DataArrayInt *descI = DataArrayInt::New();
- DataArrayInt *revDesc = DataArrayInt::New();
- DataArrayInt *revDescI = DataArrayInt::New();
+ /* 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 int *tmp = desc->getConstPointer();
- const int *tmpI=descI->getConstPointer();
-
- const int *tmpA =revDesc->getConstPointer();
- const int *tmpAI=revDescI->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 int *work=tmp+tmpI[id];//corresponds to buildDescendingConnectivity
+ const mcIdType *tmpA =revDesc->getConstPointer();//Lists the cells surrounding each face
+ const mcIdType *tmpAI=revDescI->getConstPointer();
//Test du type d'éléments contenus dans le maillage afin d'éviter les éléments contenant des points de gauss
_eltsTypes=mu->getAllGeoTypesSorted();
}
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
- DataArrayInt *revNode =DataArrayInt::New();
- DataArrayInt *revNodeI=DataArrayInt::New();
+ DataArrayIdType *revNode =DataArrayIdType::New();
+ DataArrayIdType *revNodeI=DataArrayIdType::New();
mu->getReverseNodalConnectivity(revNode,revNodeI) ;
- const int *tmpN =revNode->getConstPointer();
- const int *tmpNI=revNodeI->getConstPointer();
+ const mcIdType *tmpN =revNode->getConstPointer();//Used to know which cells surround a given node
+ const mcIdType *tmpNI=revNodeI->getConstPointer();
- DataArrayInt *revCell =DataArrayInt::New();
- DataArrayInt *revCellI=DataArrayInt::New();
- mu2->getReverseNodalConnectivity(revCell,revCellI) ;
- const int *tmpC =revCell->getConstPointer();
- const int *tmpCI=revCellI->getConstPointer();
+ DataArrayIdType *revCell =DataArrayIdType::New();
+ DataArrayIdType *revCellI=DataArrayIdType::New();
+ mu2->getReverseNodalConnectivity(revCell,revCellI);
+ const mcIdType *tmpC =revCell->getConstPointer();//Used to know which faces surround a given node
+ const mcIdType *tmpCI=revCellI->getConstPointer();
- const DataArrayInt *nodal = mu2->getNodalConnectivity() ;
- const DataArrayInt *nodalI = mu2->getNodalConnectivityIndex() ;
- const int *tmpNE =nodal->getConstPointer();
- const int *tmpNEI=nodalI->getConstPointer();
+ const DataArrayIdType *nodal = mu2->getNodalConnectivity() ;
+ const DataArrayIdType *nodalI = mu2->getNodalConnectivityIndex() ;
+ 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] ;
+ _cells = std::shared_ptr<Cell>(new Cell[_numberOfCells], std::default_delete<Cell[]>()) ;
- _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 = std::shared_ptr<Node>(new Node[_numberOfNodes], std::default_delete<Node[]>()) ;//This array may be resized if isolated nodes are found
_numberOfFaces = mu2->getNumberOfCells();
- _faces = new Face[_numberOfFaces] ;
+ _faces = std::shared_ptr<Face>(new Face[_numberOfFaces], std::default_delete<Face[]>()) ;
_indexFacePeriodicSet=false;
//Definition used if _meshDim =3 to determine the edges
- DataArrayInt *desc2 =DataArrayInt::New();
- DataArrayInt *descI2=DataArrayInt::New();
- DataArrayInt *revDesc2 =DataArrayInt::New();
- DataArrayInt *revDescI2=DataArrayInt::New();
- DataArrayInt *revNode2 =DataArrayInt::New();
- DataArrayInt *revNodeI2=DataArrayInt::New();
- const int *tmpN2 ;
- const int *tmpNI2;
+ DataArrayIdType *desc2 =DataArrayIdType::New();
+ DataArrayIdType *descI2=DataArrayIdType::New();
+ DataArrayIdType *revDesc2 =DataArrayIdType::New();
+ DataArrayIdType *revDescI2=DataArrayIdType::New();
+ DataArrayIdType *revNode2 =DataArrayIdType::New();
+ DataArrayIdType *revNodeI2=DataArrayIdType::New();
+ const mcIdType *tmpN2 ;
+ const mcIdType *tmpNI2;
MEDCouplingUMesh* mu3;
if (_meshDim == 1)
}
// _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 int *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) ) ;
-
- std::vector<int> nodeIdsOfCell ;
+ 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<mcIdType> nodeIdsOfCell ;
mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
- for( int el=0;el<nbVertices;el++ )
- {
- ci.addNodeId(el,nodeIdsOfCell[el]) ;
- ci.addFaceId(el,nodeIdsOfCell[el]) ;
- }
-
- //This assumes that in a 1D mesh the cell numbers form a consecutive sequence
- //going either from left to right (xn=-1) or right to left (xn=1)
- double xn = (cood[nodeIdsOfCell[nbVertices-1]] - cood[nodeIdsOfCell[0]] > 0.0) ? -1.0 : 1.0;
-
- for( int el=0;el<nbFaces;el++ )
+
+ mcIdType nbFaces=tmpI[id+1]-tmpI[id];//should be equal to 2
+ assert( nbFaces==2);
+ const mcIdType *work=tmp+tmpI[id];
+
+ /* compute the normal to the face */
+ xn = cood[nodeIdsOfCell[0]*_spaceDim ] - cood[nodeIdsOfCell[nbFaces-1]*_spaceDim ];
+ if(_spaceDim>1)
+ 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;
}
- _cells[id] = ci ;
+
+ Cell ci( nbVertices, nbFaces, surf[id], p ) ;//nbCells=nbFaces=2
+ for( int el=0;el<nbFaces;el++ )
+ {
+ ci.addNodeId(el,nodeIdsOfCell[el]) ;//global node number
+ ci.addNormalVector(el,xn,yn,zn) ;
+ ci.addFaceId(el,work[el]) ;
+ xn = - xn; yn=-yn; zn=-zn;
+ }
+ _cells.get()[id] = ci ;
}
-
- for( int id(0), k(0); id<_numberOfNodes; id++, k+=_spaceDim)
+
+ for( int id(0); id<_numberOfFaces; id++ )
{
- Point p(cood[k], 0.0, 0.0) ;
- const int *workc=tmpN+tmpNI[id];
- int nbCells=tmpNI[id+1]-tmpNI[id];
-
- const int *workf=tmpC+tmpCI[id];
- int nbFaces=tmpCI[id+1]-tmpCI[id];
- const int *workn=tmpA+tmpAI[id];
- int nbNeighbourNodes=tmpAI[id+1]-tmpAI[id];
- Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
-
- for( int el=0;el<nbCells;el++ )
- vi.addCellId(el,workc[el]) ;
- for( int el=0;el<nbFaces;el++ )
- vi.addFaceId(el,workf[el]) ;
- for( int el=0;el<nbNeighbourNodes;el++ )
- vi.addNeighbourNodeId(el,workn[el]) ;
- _nodes[id] = vi ;
+ const mcIdType *workv=tmpNE+tmpNEI[id]+1;
+ mcIdType nbNodes= tmpNEI[id+1]-tmpNEI[id]-1;//Normally equal to 1.
+ assert( nbNodes==1);
+
+ std::vector<double> 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_id<nbNodes;node_id++)//This loop could be deleted since nbNodes=1. Trying to merge with setMesh
+ fi.addNodeId(node_id,workv[node_id]) ;//global node number
+
+ fi.addCellId(0,workc[0]) ;
+ for(int cell_id=1; cell_id<nbCells;cell_id++)
+ {
+ int cell_idx=0;
+ if (workc[cell_id]!=workc[cell_id-1])//For some meshes (bad ones) the same cell can appear several times
+ {
+ fi.addCellId(cell_idx+1,workc[cell_id]) ;
+ cell_idx++;
+ }
+ }
+ if(nbCells==1)
+ _boundaryFaceIds.push_back(id);
+ _faces.get()[id] = fi ;
}
- _boundaryNodeIds.push_back(0);
- _boundaryNodeIds.push_back(_numberOfNodes-1);
-
- for(int id(0), k(0); id<_numberOfFaces; id++, k+=_spaceDim)
+
+ int correctNbNodes=0;
+ for( int id=0;id<_numberOfNodes;id++ )
{
- Point p(cood[k], 0.0, 0.0) ;
- const int *workc=tmpA+tmpAI[id];
- int nbCells=tmpAI[id+1]-tmpAI[id];
-
- const int *workv=tmpNE+tmpNEI[id]+1;
- Face fi( 1, nbCells, 1.0, p, 1.0, 0.0, 0.0) ;
- fi.addNodeId(0,workv[0]) ;
-
- for(int idCell=0; idCell<nbCells; idCell++)
- fi.addCellId(idCell,workc[idCell]) ;
-
- _faces[id] = fi ;
+ const mcIdType *workc=tmpN+tmpNI[id];
+ mcIdType nbCells=tmpNI[id+1]-tmpNI[id];
+
+ if( nbCells>0)//To make sure this is not an isolated node
+ {
+ correctNbNodes++;
+ std::vector<double> 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;el<nbCells;el++ )
+ vi.addCellId(el,workc[el]) ;
+ for( int el=0;el<nbNeighbourNodes;el++ )
+ vi.addNeighbourNodeId(el,workn[el]) ;//global node number
+ for( int el=0;el<nbFaces;el++ )
+ vi.addFaceId(el,workf[el],_faces.get()[workf[el]].isBorder()) ;
+ if(vi.isBorder())
+ _boundaryNodeIds.push_back(id);
+ _nodes.get()[id] = vi ;
+ }
+ }
+ if( _numberOfNodes!=correctNbNodes)//To do : reduce the size of pointer _nodes
+ {
+ cout<<"Found isolated nodes : correctNbNodes= "<<correctNbNodes<<", _numberOfNodes= "<<_numberOfNodes<<endl;
+ _numberOfNodes = correctNbNodes;
+ //memcpy(_nodes,mesh.getNodes(),correctNbNodes*sizeof(*mesh.getNodes())) ;
}
- _boundaryFaceIds.push_back(0);
- _boundaryFaceIds.push_back(_numberOfFaces-1);
}
- else if(_spaceDim==2 || _spaceDim==3)
+ else if(_meshDim==2 || _meshDim==3)
{
- DataArrayDouble *barySeg = mu2->computeIsoBarycenterOfNodesPerCell();//computeCellCenterOfMass() ;
+ DataArrayDouble *barySeg = mu2->computeIsoBarycenterOfNodesPerCell();//computeCellCenterOfMass() ;//Used as face center
const double *coorBarySeg=barySeg->getConstPointer();
- MEDCouplingFieldDouble* fieldn;
+ MEDCouplingFieldDouble* fieldl=mu2->getMeasureField(true);
+ DataArrayDouble *longueur = fieldl->getArray();
+ const double *lon=longueur->getConstPointer();//The lenght/area of each face
+
+ MEDCouplingFieldDouble* fieldn;//The normal to each face
DataArrayDouble *normal;
const double *tmpNormal;
/*Building mesh cells */
for(int id(0), k(0); id<_numberOfCells; id++, k+=_spaceDim)
{
- const int *work=tmp+tmpI[id];
- int nbFaces=tmpI[id+1]-tmpI[id];
+ const mcIdType *work=tmp+tmpI[id];
+ mcIdType nbFaces=tmpI[id+1]-tmpI[id];
- int nbVertices=mu->getNumberOfNodesInCell(id) ;
+ mcIdType nbVertices=mu->getNumberOfNodesInCell(id) ;
vector<double> coorBaryXyz(3,0);
for (int d=0; d<_spaceDim; d++)
Cell ci( nbVertices, nbFaces, surf[id], p ) ;
/* Filling cell nodes */
- std::vector<int> nodeIdsOfCell ;
+ std::vector<mcIdType> nodeIdsOfCell ;
mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
for( int el=0;el<nbVertices;el++ )
ci.addNodeId(el,nodeIdsOfCell[el]) ;
if(_spaceDim==_meshDim)//use the normal field generated by buildOrthogonalField()
for( int el=0;el<nbFaces;el++ )
{
- int faceIndex=(abs(work[el])-1);//=work[el] since Fortran type numbering was used, and negative sign means anticlockwise numbering
+ mcIdType faceIndex=(abs(work[el])-1);//=work[el] since Fortran type numbering was used, and negative sign means anticlockwise numbering
vector<double> xyzn(3,0);//Outer normal to the cell
if (work[el]<0)
for (int d=0; d<_spaceDim; d++)
ci.addFaceId(el,faceIndex) ;
}
else//build normals associated to the couple (cell id, face el)
- {
+ {//Case _meshDim=1 should be moved upper since we are in the 2D/3D branch
if(_meshDim==1)//we know in this case there are only two faces around the cell id, each face is composed of a single node
{//work[0]= first face global number, work[1]= second face global number
- int indexFace0=abs(work[0])-1;//=work[0] since Fortran type numbering was used, and negative sign means anticlockwise numbering
- int indexFace1=abs(work[1])-1;//=work[1] since Fortran type numbering was used, and negative sign means anticlockwise numbering
- int idNodeA=(tmpNE+tmpNEI[indexFace0]+1)[0];//global number of the first face node work[0]=(abs(work[0])-1)
- int idNodeB=(tmpNE+tmpNEI[indexFace1]+1)[0];//global number of the second face node work[1]=(abs(work[1])-1)
+ mcIdType indexFace0=abs(work[0])-1;//=work[0] since Fortran type numbering was used, and negative sign means anticlockwise numbering
+ mcIdType indexFace1=abs(work[1])-1;//=work[1] since Fortran type numbering was used, and negative sign means anticlockwise numbering
+ mcIdType idNodeA=(tmpNE+tmpNEI[indexFace0]+1)[0];//global number of the first face node work[0]=(abs(work[0])-1)
+ mcIdType idNodeB=(tmpNE+tmpNEI[indexFace1]+1)[0];//global number of the second face node work[1]=(abs(work[1])-1)
Vector vecAB(3);
for(int i=0;i<_spaceDim;i++)
vecAB[i]=coo->getIJ(idNodeB,i) - coo->getIJ(idNodeA,i);
for( int el=0;el<nbFaces;el++ )
{
int faceIndex=(abs(work[el])-1);//=work[el] since Fortran type numbering was used, and negative sign means anticlockwise numbering
- const int *workv=tmpNE+tmpNEI[faceIndex]+1;
- int nbNodes= tmpNEI[faceIndex+1]-tmpNEI[faceIndex]-1;
+ const mcIdType *workv=tmpNE+tmpNEI[faceIndex]+1;
+ mcIdType nbNodes= tmpNEI[faceIndex+1]-tmpNEI[faceIndex]-1;
if(nbNodes!=2)//We want to compute the normal to a straight line, not a curved interface composed of more thant 2 points
{
cout<<"Mesh name "<< mu->getName()<< " space dim= "<< _spaceDim <<" mesh dim= "<< _meshDim <<endl;
throw CdmathException("Mesh::setMesh number of nodes around a face should be 2");
}
- int idNodeA=workv[0];
- int idNodeB=workv[1];
+ mcIdType idNodeA=workv[0];
+ mcIdType idNodeB=workv[1];
vector<double> nodeA(_spaceDim), nodeB(_spaceDim), nodeP(_spaceDim);
for(int i=0;i<_spaceDim;i++)
{
}
}
}
- _cells[id] = ci ;
+ _cells.get()[id] = ci ;
}
- MEDCouplingFieldDouble* fieldl=mu2->getMeasureField(true);
- DataArrayDouble *longueur = fieldl->getArray();
- const double *lon=longueur->getConstPointer();
-
if(_spaceDim!=_meshDim)
{
/* Since spaceDim!=meshDim, don't build normal to faces */
/*Building mesh faces */
for(int id(0), k(0); id<_numberOfFaces; id++, k+=_spaceDim)
{
- vector<double> coorBarySegXyz(3,0);
+ vector<double> coorBarySegXyz(3);
for (int d=0; d<_spaceDim; d++)
coorBarySegXyz[d] = coorBarySeg[k+d];
Point p(coorBarySegXyz[0],coorBarySegXyz[1],coorBarySegXyz[2]) ;
- const int *workc=tmpA+tmpAI[id];
- int nbCells=tmpAI[id+1]-tmpAI[id];
+ const mcIdType *workc=tmpA+tmpAI[id];
+ mcIdType nbCells=tmpAI[id+1]-tmpAI[id];
if (nbCells>2 && _spaceDim==_meshDim)
{
if (nbCells==1)
_boundaryFaceIds.push_back(id);
- const int *workv=tmpNE+tmpNEI[id]+1;
- int nbNodes= tmpNEI[id+1]-tmpNEI[id]-1;
+ const mcIdType *workv=tmpNE+tmpNEI[id]+1;
+ mcIdType nbNodes= tmpNEI[id+1]-tmpNEI[id]-1;
Face fi;
if(_spaceDim==_meshDim)//Euclidean flat mesh geometry
}
}
- _faces[id] = fi ;
+ _faces.get()[id] = fi ;
}
/*Building mesh nodes, should be done after building mesh faces in order to detect boundary nodes*/
+ int correctNbNodes=0;
for(int id(0), k(0); id<_numberOfNodes; id++, k+=_spaceDim)
{
- vector<double> coorP(3,0);
- for (int d=0; d<_spaceDim; d++)
- coorP[d] = cood[k+d];
- Point p(coorP[0],coorP[1],coorP[2]) ;
-
- const int *workc=tmpN+tmpNI[id];
- int nbCells=tmpNI[id+1]-tmpNI[id];
- const int *workf=tmpC+tmpCI[id];
- int nbFaces=tmpCI[id+1]-tmpCI[id];
- const int *workn;
- int nbNeighbourNodes;
- if (_meshDim == 1)
- {
- workn=tmpA+tmpAI[id];
- nbNeighbourNodes=tmpAI[id+1]-tmpAI[id];
- }
- else if (_meshDim == 2)
- {
- workn=tmpC+tmpCI[id];
- nbNeighbourNodes=tmpCI[id+1]-tmpCI[id];
- }
- else//_meshDim == 3
- {
- workn=tmpN2+tmpNI2[id];
- nbNeighbourNodes=tmpNI2[id+1]-tmpNI2[id];
- }
- Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
-
- for( int el=0;el<nbCells;el++ )
- vi.addCellId(el,workc[el]) ;
- for( int el=0;el<nbNeighbourNodes;el++ )
- vi.addNeighbourNodeId(el,workn[el]) ;
- //Detection of border nodes
- bool isBorder=false;
- for( int el=0;el<nbFaces;el++ )
- {
- vi.addFaceId(el,workf[el],_faces[workf[el]].isBorder()) ;
- isBorder= isBorder || _faces[workf[el]].isBorder();
- }
- if(isBorder)
- _boundaryNodeIds.push_back(id);
- _nodes[id] = vi ;
- }
+ const mcIdType *workc=tmpN+tmpNI[id];
+ mcIdType nbCells=tmpNI[id+1]-tmpNI[id];
+ if( nbCells>0)//To make sure this is not an isolated node
+ {
+ correctNbNodes++;
+ vector<double> coorP(3);
+ for (int d=0; d<_spaceDim; d++)
+ coorP[d] = cood[k+d];
+ Point p(coorP[0],coorP[1],coorP[2]) ;
+
+ const mcIdType *workf=tmpC+tmpCI[id];
+ mcIdType nbFaces=tmpCI[id+1]-tmpCI[id];
+ const mcIdType *workn;
+ mcIdType nbNeighbourNodes;
+ if (_meshDim == 1)
+ {
+ workn=tmpA+tmpAI[id];
+ nbNeighbourNodes=tmpAI[id+1]-tmpAI[id];
+ }
+ else if (_meshDim == 2)
+ {
+ workn=tmpC+tmpCI[id];
+ nbNeighbourNodes=tmpCI[id+1]-tmpCI[id];
+ }
+ else//_meshDim == 3
+ {
+ workn=tmpN2+tmpNI2[id];
+ nbNeighbourNodes=tmpNI2[id+1]-tmpNI2[id];
+ }
+ Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
+
+ for( int el=0;el<nbCells;el++ )
+ vi.addCellId(el,workc[el]) ;
+ for( int el=0;el<nbNeighbourNodes;el++ )
+ vi.addNeighbourNodeId(el,workn[el]) ;
+ //Detection of border nodes
+ for( int el=0;el<nbFaces;el++ )
+ vi.addFaceId(el,workf[el],_faces.get()[workf[el]].isBorder()) ;
+ if(vi.isBorder())
+ _boundaryNodeIds.push_back(id);
+ _nodes.get()[id] = vi ;
+ }
+ }
+ if( _numberOfNodes!=correctNbNodes)//To do : reduce the size of pointer _nodes
+ {
+ cout<<"Found isolated nodes : correctNbNodes= "<<correctNbNodes<<", _numberOfNodes= "<<_numberOfNodes<<endl;
+ _numberOfNodes = correctNbNodes;
+ }
+
if(_spaceDim==_meshDim)
fieldn->decrRef();
fieldl->decrRef();
else
throw CdmathException("Mesh::setMesh space dimension should be 1, 2 or 3");
- //definition of the bounding box for unstructured meshes
- if(!_isStructured)//Structured case is treated in function readMeshMed
- {
- double Box0[2*_spaceDim];
- coo->getMinMaxPerComponent(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];
- }
- }
-
//Set boundary groups
- _faceGroupNames.push_back("Boundary");
+ _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();
revDescI2->decrRef();
mu3->decrRef();
}
-
+
return mu;
}
-//----------------------------------------------------------------------
-Mesh::Mesh( double xmin, double xmax, int nx, std::string meshName )
-//----------------------------------------------------------------------
+void
+Mesh::setGroupAtFaceByCoords(double x, double y, double z, double eps, std::string groupName, bool isBoundaryGroup)
{
- 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];
- int *nodeStrctPtr = new int[_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;
-
- DataArrayInt *desc=DataArrayInt::New();
- DataArrayInt *descI=DataArrayInt::New();
- DataArrayInt *revDesc=DataArrayInt::New();
- DataArrayInt *revDescI=DataArrayInt::New();
- MEDCouplingUMesh* mu=_mesh->buildUnstructured();
- MEDCouplingUMesh *mu2=mu->buildDescendingConnectivity(desc,descI,revDesc,revDescI);
-
- const int *tmp=desc->getConstPointer();
- const int *tmpI=descI->getConstPointer();
-
- const int *tmpA =revDesc->getConstPointer();
- const int *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<int> nodeIdsOfCell ;
- mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
- for( int el=0;el<nbVertices;el++ )
+ 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)<eps && abs(FY-y)<eps && abs(FZ-z)<eps)
+ {
+ faceIds.insert(faceIds.end(),_boundaryFaceIds[i]);
+ _faces.get()[_boundaryFaceIds[i]].setGroupName(groupName);
+ }
+ }
+ else
+ for (int iface=0;iface<_numberOfFaces;iface++)
{
- ci.addNodeId(el,nodeIdsOfCell[el]) ;
- ci.addFaceId(el,nodeIdsOfCell[el]) ;
+ Fi=_faces.get()[iface];
+ FX=Fi.x();
+ FY=Fi.y();
+ FZ=Fi.z();
+ if (abs(FX-x)<eps && abs(FY-y)<eps && abs(FZ-z)<eps)
+ {
+ faceIds.insert(faceIds.end(),iface);
+ _faces.get()[iface].setGroupName(groupName);
+ }
}
- double xn = (cood[nodeIdsOfCell[nbVertices-1]] - cood[nodeIdsOfCell[0]] > 0.0) ? -1.0 : 1.0;
-
- int nbFaces=tmpI[id+1]-tmpI[id];
- const int *work=tmp+tmpI[id];
-
- for( int el=0;el<nbFaces;el++ )
+ if (faceIds.size()>0)
+ {
+ std::vector< std::string >::iterator it = std::find(_faceGroupNames.begin(), _faceGroupNames.end(), groupName);
+ if(it == _faceGroupNames.end())//No group named groupName
{
- ci.addNormalVector(el,xn,0.0,0.0) ;
- ci.addFaceId(el,work[el]) ;
- xn = - xn;
+ _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;
}
-
- _cells[id] = ci ;
-
- comp=comp+2;
- }
-
- //Suppress the following since tmpN=tmpA
- DataArrayInt *revNode=DataArrayInt::New();
- DataArrayInt *revNodeI=DataArrayInt::New();
- mu->getReverseNodalConnectivity(revNode,revNodeI) ;
- const int *tmpN=revNode->getConstPointer();
- const int *tmpNI=revNodeI->getConstPointer();
-
- for( int id=0;id<_numberOfNodes;id++ )
- {
- std::vector<double> coo ;
- mu->getCoordinatesOfNode(id,coo);
- Point p(coo[0],0.0,0.0) ;
- const int *workc=tmpN+tmpNI[id];
- int nbCells=tmpNI[id+1]-tmpNI[id];
- int nbFaces=1;
- const int *workn=tmpA+tmpAI[id];
- int nbNeighbourNodes=tmpAI[id+1]-tmpAI[id];
-
- Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
- for( int el=0;el<nbCells;el++ )
- vi.addCellId(el,workc[el]) ;
- for( int el=0;el<nbFaces;el++ )
- vi.addFaceId(el,id) ;
- for( int el=0;el<nbNeighbourNodes;el++ )
- vi.addNeighbourNodeId(el,workn[el]) ;
- _nodes[id] = vi ;
-
-
- int nbVertices=1;
- Face fi( nbVertices, nbCells, 1.0, p, 1., 0., 0. ) ;
- for( int el=0;el<nbVertices;el++ )
- fi.addNodeId(el,id) ;
-
- for( int el=0;el<nbCells;el++ )
- fi.addCellId(el,workc[el]) ;
- _faces[id] = fi ;
}
-
- //Set boundary groups
- _faceGroupNames.push_back("Boundary");
- _nodeGroupNames.push_back("Boundary");
- _boundaryFaceIds.push_back(0);
- _boundaryFaceIds.push_back(_numberOfFaces-1);
- _boundaryNodeIds.push_back(0);
- _boundaryNodeIds.push_back(_numberOfNodes-1);
-
- fieldl->decrRef();
- baryCell->decrRef();
- desc->decrRef();
- descI->decrRef();
- revDesc->decrRef();
- revDescI->decrRef();
- revNode->decrRef();
- revNodeI->decrRef();
- mu2->decrRef();
- mu->decrRef();
}
-//----------------------------------------------------------------------
-Mesh::Mesh( std::vector<double> points, std::string meshName )
-//----------------------------------------------------------------------
+void
+Mesh::setGroupAtNodeByCoords(double x, double y, double z, double eps, std::string groupName, bool isBoundaryGroup)
{
- int nx=points.size();
-
- if(nx<2)
- throw CdmathException("Mesh::Mesh( vector<double> points, string meshName) : nx < 2, vector should contain at least two values");
- int i=0;
- while( i<nx-1 && points[i+1]>points[i] )
- i++;
- if( i!=nx-1 )
- {
- //cout<< points << endl;
- throw CdmathException("Mesh::Mesh( vector<double> points, string meshName) : vector values should be sorted");
- }
-
- _spaceDim = 1 ;
- _meshDim = 1 ;
- _name=meshName;
- _xMin=points[0];
- _xMax=points[nx-1];
- _yMin=0.;
- _yMax=0.;
- _zMin=0.;
- _zMax=0.;
+ 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)<eps && abs(NY-y)<eps && abs(NZ-z)<eps)
+ {
+ nodeIds.insert(nodeIds.end(),_boundaryNodeIds[i]);
+ _nodes.get()[_boundaryNodeIds[i]].setGroupName(groupName);
+ }
+ }
+ else
+ for (int inode=0;inode<_numberOfNodes;inode++)
+ {
+ NX=_nodes.get()[inode].x();
+ NY=_nodes.get()[inode].y();
+ NZ=_nodes.get()[inode].z();
+ if (abs(NX-x)<eps && abs(NY-y)<eps && abs(NZ-z)<eps)
+ {
+ nodeIds.insert(nodeIds.end(),inode);
+ _nodes.get()[inode].setGroupName(groupName);
+ }
+ }
- _isStructured = false;
- _indexFacePeriodicSet=false;
-
- MEDCouplingUMesh * mesh1d = MEDCouplingUMesh::New(meshName, 1);
- mesh1d->allocateCells(nx - 1);
- double * coords = new double[nx];
- int * nodal_con = new int[2];
- coords[0]=points[0];
- for(int i=0; i<nx- 1 ; i++)
+ if (nodeIds.size()>0)
{
- nodal_con[0]=i;
- nodal_con[1]=i+1;
- mesh1d->insertNextCell(INTERP_KERNEL::NORM_SEG2, 2, nodal_con);
- coords[i+1]=points[i + 1];
+ 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;
+ }
}
- 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, nodal_con;
- coords_arr->decrRef();
-
- _mesh=mesh1d;
-
- DataArrayInt *desc=DataArrayInt::New();
- DataArrayInt *descI=DataArrayInt::New();
- DataArrayInt *revDesc=DataArrayInt::New();
- DataArrayInt *revDescI=DataArrayInt::New();
- MEDCouplingUMesh* mu=_mesh->buildUnstructured();
- MEDCouplingUMesh *mu2=mu->buildDescendingConnectivity(desc,descI,revDesc,revDescI);
-
- DataArrayDouble *baryCell = mu->computeCellCenterOfMass() ;
- const double *coorBary=baryCell->getConstPointer();
-
- _numberOfCells = _mesh->getNumberOfCells() ;
- _cells = new Cell[_numberOfCells] ;
-
- _numberOfEdges = _numberOfCells;
-
- _eltsTypes=mu->getAllGeoTypesSorted();
+}
- MEDCouplingFieldDouble* fieldl=mu->getMeasureField(true);
- DataArrayDouble *longueur = fieldl->getArray();
- const double *lon=longueur->getConstPointer();
+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)<eps)
+ {
+ faceIds.insert(faceIds.end(),_boundaryFaceIds[i]);
+ _faces.get()[_boundaryFaceIds[i]].setGroupName(groupName);
+ }
+ }
+ else
+ for (int iface=0;iface<_numberOfFaces;iface++)
+ {
+ cord=_faces.get()[iface].getBarryCenter()[direction];
+ if (abs(cord-value)<eps)
+ {
+ faceIds.insert(faceIds.end(),iface);
+ _faces.get()[iface].setGroupName(groupName);
+ }
+ }
- 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 ) ;
+ /* Construction of the node group */
+ if(isBoundaryGroup)
+ for(int i=0; i<_boundaryNodeIds.size(); i++)
+ {
+ cord=_nodes.get()[_boundaryNodeIds[i]].getPoint()[direction];
+ if (abs(cord-value)<eps)
+ {
+ nodeIds.insert(nodeIds.end(),_boundaryNodeIds[i]);
+ _nodes.get()[_boundaryNodeIds[i]].setGroupName(groupName);
+ }
+ }
+ else
+ for (int inode=0;inode<_numberOfNodes;inode++)
+ {
+ cord=_nodes.get()[inode].getPoint()[direction];
+ if (abs(cord-value)<eps)
+ {
+ nodeIds.insert(nodeIds.end(),inode);
+ _nodes.get()[inode].setGroupName(groupName);
+ }
+ }
- std::vector<int> nodeIdsOfCell ;
- mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
- for( int el=0;el<nbVertices;el++ )
+ if (faceIds.size()>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
{
- ci.addNodeId(el,nodeIdsOfCell[el]) ;
- ci.addFaceId(el,nodeIdsOfCell[el]) ;
+ 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;
}
- _cells[id] = ci ;
- comp=comp+2;
}
+ 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<nbNodes ; inode++)
+ {
+ std::vector<int>::const_iterator it = std::find(_boundaryNodeIds.begin(),_boundaryNodeIds.end(),nodesID[inode]);
+ if( it != _boundaryNodeIds.end() )
+ _boundaryNodeIds.push_back(nodesID[inode]);
+ }
+ }
+}
- //Suppress the following since tmpN=tmpA
- DataArrayInt *revNode=DataArrayInt::New();
- DataArrayInt *revNodeI=DataArrayInt::New();
- mu->getReverseNodalConnectivity(revNode,revNodeI) ;
- const int *tmpN=revNode->getConstPointer();
- const int *tmpNI=revNodeI->getConstPointer();
-
- _numberOfNodes = mu->getNumberOfNodes() ;
- _nodes = new Node[_numberOfNodes] ;
- _numberOfFaces = _numberOfNodes;
- _faces = new Face[_numberOfFaces] ;
-
- for( int id=0;id<_numberOfNodes;id++ )
- {
- std::vector<double> coo ;
- mu->getCoordinatesOfNode(id,coo);
- Point p(coo[0],0.0,0.0) ;
- const int *workc=tmpN+tmpNI[id];
- int nbCells=tmpNI[id+1]-tmpNI[id];
- int nbFaces=1;
- const int *workn=tmpN+tmpNI[id];
- int 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;el<nbCells;el++ )
- vi.addCellId(el,workc[el]) ;
- for( int el=0;el<nbFaces;el++ )
- vi.addFaceId(el,id) ;
- for( int el=0;el<nbNeighbourNodes;el++ )
- vi.addNeighbourNodeId(el,workn[el]) ;
- _nodes[id] = vi ;
-
- for( int el=0;el<nbVertices;el++ )
- fi.addNodeId(el,id) ;
-
- for( int el=0;el<nbCells;el++ )
- fi.addCellId(el,workc[el]) ;
- _faces[id] = fi ;
- }
- //Set boundary groups
- _faceGroupNames.push_back("Boundary");
- _nodeGroupNames.push_back("Boundary");
- _boundaryFaceIds.push_back(0);
- _boundaryFaceIds.push_back(_numberOfFaces-1);
- _boundaryNodeIds.push_back(0);
- _boundaryNodeIds.push_back(_numberOfNodes-1);
-
- fieldl->decrRef();
- baryCell->decrRef();
- desc->decrRef();
- descI->decrRef();
- revDesc->decrRef();
- revDescI->decrRef();
- revNode->decrRef();
- revNodeI->decrRef();
- mu2->decrRef();
- mu->decrRef();
+std::map<int,int>
+Mesh::getIndexFacePeriodic( void ) const
+{
+ return _indexFacePeriodicMap;
}
-//----------------------------------------------------------------------
-Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, int split_to_triangles_policy, std::string meshName)
-//----------------------------------------------------------------------
+
+void
+Mesh::setPeriodicFaces(bool check_groups, bool use_central_inversion)
{
- 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");
+ if(_indexFacePeriodicSet)
+ return;
+
+ for (int indexFace=0;indexFace<_boundaryFaceIds.size() ; indexFace++)
+ {
+ Face my_face=_faces.get()[_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)
+ {
+ 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;
+ }
+ _indexFacePeriodicSet=true;
+}
- _xMin=xmin;
- _xMax=xmax;
- _yMin=ymin;
- _yMax=ymax;
- _zMin=0.;
- _zMax=0.;
+int
+Mesh::getIndexFacePeriodic(int indexFace, bool check_groups, bool use_central_inversion)
+{
+ if (!_faces.get()[indexFace].isBorder())
+ {
+ cout<<"Pb with indexFace= "<<indexFace<<endl;
+ throw CdmathException("Mesh::getIndexFacePeriodic: not a border face" );
+ }
+
+ if(!_indexFacePeriodicSet)
+ setPeriodicFaces(check_groups, use_central_inversion);
+ std::map<int,int>::const_iterator it = _indexFacePeriodicMap.find(indexFace);
+ if( it != _indexFacePeriodicMap.end() )
+ return it->second;
+ else
+ {
+ cout<<"Pb with indexFace= "<<indexFace<<endl;
+ throw CdmathException("Mesh::getIndexFacePeriodic: not a periodic face" );
+ }
+}
- double dx = (xmax - xmin)/nx ;
- double dy = (ymax - ymin)/ny ;
+bool
+Mesh::isBorderNode(int nodeid) const
+{
+ return _nodes.get()[nodeid].isBorder();
+}
- _spaceDim = 2 ;
- _meshDim = 2 ;
- _name=meshName;
- _indexFacePeriodicSet=false;
- _isStructured = true;
- _nxyz.resize(_spaceDim);
- _nxyz[0]=nx;
- _nxyz[1]=ny;
+bool
+Mesh::isBorderFace(int faceid) const
+{
+ 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];
- int *nodeStrctPtr = new int[_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);
+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);
- }
- 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 " <<endl;
+ cout<< _eltsTypes[i]<<endl;
+ throw CdmathException("Mesh::getElementTypes : recognised cell med types are NORM_POINT1, NORM_SEG2, NORM_TRI3, NORM_QUAD4, NORM_TETRA4, NORM_PYRA5, NORM_PENTA6, NORM_HEXA8, NORM_POLYGON, NORM_POLYHED");
}
-
- delete [] originPtr;
- delete [] dxyzPtr;
- delete [] nodeStrctPtr;
-
- setMesh();
+ }
+ return result;
}
-//----------------------------------------------------------------------
-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)
-//----------------------------------------------------------------------
+void
+Mesh::setFaceGroups( const MEDFileUMesh* medmesh, MEDCouplingUMesh* mu)
{
- 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");
- 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");
- 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;
- _indexFacePeriodicSet=false;
- _xMin=xmin;
- _xMax=xmax;
- _yMin=ymin;
- _yMax=ymax;
- _zMin=zmin;
- _zMax=zmax;
+ int nbCellsSubMesh;
+ bool foundFace;
+
+ /* Searching for face groups */
+ vector<string> faceGroups=medmesh->getGroupsNames() ;
- double dx = (xmax - xmin)/nx ;
- double dy = (ymax - ymin)/ny ;
- double dz = (zmax - zmin)/nz ;
+ for (unsigned int i=0;i<faceGroups.size();i++ )
+ {
+ string groupName=faceGroups[i];
+ vector<mcIdType> 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<mcIdType>::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<int>(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<nbCellsSubMesh; ic++, k+=_spaceDim)
+ {
+ vector<double> coorBaryXyz(3,0);
+ for (int d=0; d<_spaceDim; d++)
+ coorBaryXyz[d] = coorBary[k+d];
+ Point p1(coorBaryXyz[0],coorBaryXyz[1],coorBaryXyz[2]) ;
- _isStructured = true;
- _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<string> nodeGroups=medmesh->getGroupsOnSpecifiedLev(1) ;
- _nxyz.resize(_spaceDim);
- _nxyz[0]=nx;
- _nxyz[1]=ny;
- _nxyz[2]=nz;
+ for (unsigned int i=0;i<nodeGroups.size();i++ )
+ {
+ string groupName=nodeGroups[i];
+ DataArrayIdType * nodeGroup=medmesh->getNodeGroupArr( groupName );
+ const mcIdType *nodeids=nodeGroup->getConstPointer();
- double *originPtr = new double[_spaceDim];
- double *dxyzPtr = new double[_spaceDim];
- int *nodeStrctPtr = new int[_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);
+ 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);
- }
- else if( split_to_tetrahedra_policy == 1 )
- {
- _mesh=_mesh->buildUnstructured();
- _mesh->simplexize(INTERP_KERNEL::PLANAR_FACE_6);
- }
- 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<int>(nbNodesSubMesh));//Vector of boundary faces
+ /* Node identification */
+ for (int ic(0); ic<nbNodesSubMesh; ic++)
+ {
+ vector<double> 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 );
+ }
+ }
+ }
}
//----------------------------------------------------------------------
return _meshDim ;
}
-std::vector<double>
-Mesh::getDXYZ() const
-{
- if(!_isStructured)
- throw CdmathException("std::vector<double> Mesh::getDXYZ() : dx,dy and dz are defined only for structured meshes !");
-
- return _dxyz;
-}
-
-std::vector<int>
+std::vector<mcIdType>
Mesh::getCellGridStructure() const
{
if(!_isStructured)
{
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];
}
{
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];
}
double
Mesh::getXMin( void ) const
//----------------------------------------------------------------------
-{
- return _xMin ;
+{
+ double Box0[2*_meshDim];
+ _mesh->getBoundingBox(Box0);
+
+ return Box0[0] ;
}
//----------------------------------------------------------------------
Mesh::getXMax( void ) const
//----------------------------------------------------------------------
{
- return _xMax ;
+ double Box0[2*_meshDim];
+ _mesh->getBoundingBox(Box0);
+
+ return Box0[1] ;
}
//----------------------------------------------------------------------
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] ;
}
//----------------------------------------------------------------------
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] ;
}
//----------------------------------------------------------------------
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] ;
}
//----------------------------------------------------------------------
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] ;
}
//----------------------------------------------------------------------
}
//----------------------------------------------------------------------
-Face*
+std::shared_ptr<Face>
Mesh::getFaces ( void ) const
//----------------------------------------------------------------------
{
}
//----------------------------------------------------------------------
-Cell*
+std::shared_ptr<Cell>
Mesh::getCells ( void ) const
//----------------------------------------------------------------------
{
Mesh::getCell ( int i ) const
//----------------------------------------------------------------------
{
- return _cells[i] ;
+ return _cells.get()[i] ;
}
//----------------------------------------------------------------------
Mesh::getFace ( int i ) const
//----------------------------------------------------------------------
{
- return _faces[i] ;
+ return _faces.get()[i] ;
}
//----------------------------------------------------------------------
Mesh::getNode ( int i ) const
//----------------------------------------------------------------------
{
- return _nodes[i] ;
+ return _nodes.get()[i] ;
}
//----------------------------------------------------------------------
-Node*
+std::shared_ptr<Node>
Mesh::getNodes ( void ) const
//----------------------------------------------------------------------
{
return _faceGroupNames;
}
-vector<MEDCoupling::MEDCouplingUMesh *>
+vector< std::vector<int> >
Mesh::getFaceGroups( void ) const
+{
+ return _faceGroupsIds;
+}
+
+vector<MEDCoupling::MEDCouplingUMesh *>
+Mesh::getMEDCouplingFaceGroups( void ) const
{
return _faceGroups;
}
return _nodeGroupNames;
}
-vector<MEDCoupling::DataArrayInt *>
+vector< std::vector<int> >
Mesh::getNodeGroups( void ) const
+{
+ return _nodeGroupsIds;
+}
+
+
+vector<MEDCoupling::DataArrayIdType *>
+Mesh::getMEDCouplingNodeGroups( void ) const
{
return _nodeGroups;
}
_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() ;
- _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() ;
+ _faceGroupsIds = mesh.getFaceGroups() ;
+ _faceGroups = mesh.getMEDCouplingFaceGroups() ;
_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);
+ _nodeGroupsIds = mesh.getNodeGroups() ;
+ _nodeGroups = mesh.getMEDCouplingNodeGroups() ;
- 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)
_eltsTypes=mesh.getElementTypes();
_eltsTypesNames=mesh.getElementTypesNames();
- MCAuto<MEDCouplingMesh> 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;
}
}
//----------------------------------------------------------------------
double
-Mesh::minRatioVolSurf()
+Mesh::minRatioVolSurf() const
{
double dx_min = 1e30;
for(int i=0; i<_numberOfCells; i++)
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"<<endl;
+ cout<<"Known face group names are " ;
+ for(int i=0; i<_faceGroupNames.size(); i++)
+ cout<< _faceGroupNames[i]<<" ";
+ cout<< endl;
+ throw CdmathException("Required face group does not exist");
+ }
+ else if ( count_occurences <= nth_occurence)//Too many groups found
+ {
+ cout<<"Mesh::getBoundaryGroupMesh Error : "<<count_occurences<<" groups have name " << groupName<<", but you asked fo occurencer number "<<nth_occurence<<"which is too large"<<endl;
+ cout<<"Call function getBoundaryGroupMesh ( string groupName, int nth_group_match ) with nth_group_match between 0 and "<<count_occurences-1<<" to discriminate them "<<endl ;
+ throw CdmathException("Several face groups have the same name but you asked for an occurence that does not exsit");
+ }
+ else if( count_occurences >1 )//Wrning several groups found
+ cout<<"Warning : "<<count_occurences<<" groups have name " << groupName<<". Searching occurence number "<<nth_occurence<<endl;
+
+ //search occurence of group name in known group name list
+ std::vector<std::string>::const_iterator it = _faceGroupNames.begin();
+ for (int i = 0; i<nth_occurence+1; i++)
+ it = std::find(it,_faceGroupNames.end(),groupName);
+
+ return Mesh(_faceGroups[it - _faceGroupNames.begin()]);
+}
+
int
Mesh::getMaxNbNeighbours(EntityType type) const
{
- double result=0;
+ int result=0;
if (type==CELLS)
{
+ int nbNeib;//local number of neighbours
for(int i=0; i<_numberOfCells; i++)
- if(result < _cells[i].getNumberOfFaces())
- result=_cells[i].getNumberOfFaces();
+ {
+ Cell Ci = _cells.get()[i];
+ //Careful with mesh with junctions : do not just take Ci.getNumberOfFaces()
+ nbNeib=0;
+ for(int j=0; j<Ci.getNumberOfFaces(); j++)
+ nbNeib+=_faces.get()[Ci.getFacesId()[j]].getNumberOfCells()-1;//Without junction this would be +=1
+
+ if(result < nbNeib)
+ result=nbNeib;
+ }
}
else if(type==NODES)
{
for(int i=0; i<_numberOfNodes; i++)
- if(result < _nodes[i].getNumberOfEdges())
- result=_nodes[i].getNumberOfEdges();
+ if(result < _nodes.get()[i].getNumberOfEdges())
+ result=_nodes.get()[i].getNumberOfEdges();
}
else
throw CdmathException("Mesh::getMaxNbNeighbours : entity type is not accepted. Should be CELLS or NODES");
Mesh::writeVTK ( const std::string fileName ) const
//----------------------------------------------------------------------
{
+ if( !_meshNotDeleted )
+ throw CdmathException("Mesh::writeVTK : Cannot save mesh : no MEDCouplingMesh loaded (may be deleted)");
+
string fname=fileName+".vtu";
_mesh->writeVTK(fname.c_str()) ;
}
//----------------------------------------------------------------------
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 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;
//meshMEDFile.setMeshAtLevel(0,mu);
//for(int i=0; i< _faceGroups.size(); i++)
//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<std::string>::const_iterator it = std::find(_faceGroupNames.begin(),_faceGroupNames.end(),groupName);
+ if( it == _faceGroupNames.end() )
{
cout<<"Mesh::getGroupFaceIds Error : face group " << groupName << " does not exist"<<endl;
throw CdmathException("Required face group does not exist");
}
else
- {
- std::vector< int > 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<std::string>::const_iterator it = std::find(_nodeGroupNames.begin(),_nodeGroupNames.end(),groupName);
+ if( it == _nodeGroupNames.end() )
{
cout<<"Mesh::getGroupNodeIds Error : node group " << groupName << " does not exist"<<endl;
throw CdmathException("Required node group does not exist");
}
else
- {
- std::vector< int > 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
{
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
getNode(nodeIds[i]).setGroupName(groupName);
}
+void Mesh::deleteMEDCouplingMesh()
+{
+ if(_meshNotDeleted)
+ {
+ (_mesh.retn())->decrRef();
+ _meshNotDeleted=false;
+ }
+ else
+ throw CdmathException("Mesh::deleteMEDCouplingMesh() : mesh is not loaded");
+};
