#include <iterator>
#include <algorithm>
#include <string.h>
+#include <assert.h>
using namespace MEDCoupling;
using namespace std;
//----------------------------------------------------------------------
{
_mesh=NULL;
+ _meshNotDeleted=false;
_cells=NULL;
_nodes=NULL;
_faces=NULL;
_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;
}
//----------------------------------------------------------------------
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
{
_spaceDim=mesh->getSpaceDimension();
_meshDim=mesh->getMeshDimension();
- _isStructured=true;
_dxyz=mesh->getDXYZ();
_nxyz=mesh->getCellGridStructure();
double* Box0=new double[2*_spaceDim];
mesh->getBoundingBox(Box0);
_name=mesh->getName();
+ _epsilon=1e-6;
_indexFacePeriodicSet=false;
_xMin=Box0[0];
originPtr+_spaceDim,
dxyzPtr,
dxyzPtr+_spaceDim);
+ _meshNotDeleted=true;
+ _isStructured=true;
+
delete [] originPtr;
delete [] dxyzPtr;
delete [] nodeStrctPtr;
delete [] Box0 ;
+
setMesh();
}
{
_spaceDim=mesh->getSpaceDimension();
_meshDim=mesh->getMeshDimension();
- _isStructured=false;
double* Box0=new double[2*_spaceDim];
mesh->getBoundingBox(Box0);
_name=mesh->getName();
+ _epsilon=1e-6;
_indexFacePeriodicSet=false;
_xMin=Box0[0];
}
_mesh=mesh->deepCopy();
+ _mesh=mesh->buildUnstructured();
+ _meshNotDeleted=true;
+ _isStructured=false;
delete [] Box0 ;
+
setMesh();
}
+//----------------------------------------------------------------------
+Mesh::Mesh( const std::string filename, int meshLevel )
+//----------------------------------------------------------------------
+{
+ readMeshMed(filename, meshLevel);
+}
+
//----------------------------------------------------------------------
Mesh::Mesh( const Mesh& mesh )
//----------------------------------------------------------------------
_spaceDim = mesh.getSpaceDimension() ;
_meshDim = mesh.getMeshDimension() ;
_name=mesh.getName();
+ _epsilon=mesh.getComparisonEpsilon();
_xMax=mesh.getXMax();
_yMin=mesh.getYMin();
_yMax=mesh.getYMax();
MCAuto<MEDCouplingMesh> m1=mesh.getMEDCouplingMesh()->deepCopy();
_mesh=m1;
-}
-
-//----------------------------------------------------------------------
-Mesh::Mesh( const std::string filename, int meshLevel )
-//----------------------------------------------------------------------
-{
- readMeshMed(filename, meshLevel);
+ _meshNotDeleted=mesh.meshNotDeleted();
}
//----------------------------------------------------------------------
_meshDim=_mesh->getMeshDimension();
_spaceDim=_mesh->getSpaceDimension();
_name=_mesh->getName();
+ _epsilon=1e-6;
_indexFacePeriodicSet=false;
MEDCoupling::MEDCouplingIMesh* structuredMesh = dynamic_cast<MEDCoupling::MEDCouplingIMesh*> (_mesh.retn());
if(structuredMesh)
{
_isStructured=true;
+ _meshNotDeleted=false;
_dxyz=structuredMesh->getDXYZ();
_nxyz=structuredMesh->getCellGridStructure();
double* Box0=new double[2*_spaceDim];
}
}
else
+ {
_isStructured=false;
+ _meshNotDeleted=true;
+ }
MEDCouplingUMesh* mu = setMesh();
setGroups(m, mu);
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::setGroupAtFaceByCoords(double x, double y, double z, double eps, std::string groupName, bool isBoundaryGroup)
{
- 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)
+ 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[_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[_boundaryFaceIds[i]].setGroupName(groupName);
+ }
+ }
+ else
+ for (int iface=0;iface<_numberOfFaces;iface++)
{
- _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);
+ Fi=_faces[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[iface].setGroupName(groupName);
+ }
+ }
- flag=true;
+ 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
+ {
+ std::vector< int > faceGroupIds = _faceGroupsIds[it-_faceGroupNames.begin()];
+ faceGroupIds.insert( faceGroupIds.end(), faceIds.begin(), faceIds.end());
+ /* Detect and erase duplicates face ids */
+ sort( faceGroupIds.begin(), faceGroupIds.end() );
+ faceGroupIds.erase( unique( faceGroupIds.begin(), faceGroupIds.end() ), faceGroupIds.end() );
+ _faceGroupsIds[it-_faceGroupNames.begin()] = faceGroupIds;
}
}
- if (flag)
+}
+
+void
+Mesh::setGroupAtNodeByCoords(double x, double y, double z, double eps, std::string groupName, bool isBoundaryGroup)
+{
+ std::vector< int > nodeIds(0);
+ double NX, NY, NZ;
+ Node Ni;
+
+ /* Construction of the node group */
+ if(isBoundaryGroup)
+ for(int i=0; i<_boundaryNodeIds.size(); i++)
+ {
+ Ni=_nodes[_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[_boundaryNodeIds[i]].setGroupName(groupName);
+ }
+ }
+ else
+ for (int inode=0;inode<_numberOfNodes;inode++)
+ {
+ NX=_nodes[inode].x();
+ NY=_nodes[inode].y();
+ NZ=_nodes[inode].z();
+ if (abs(NX-x)<eps && abs(NY-y)<eps && abs(NZ-z)<eps)
+ {
+ nodeIds.insert(nodeIds.end(),inode);
+ _nodes[inode].setGroupName(groupName);
+ }
+ }
+
+ if (nodeIds.size()>0)
{
- _faceGroupNames.insert(_faceGroupNames.begin(),groupName);
- _nodeGroupNames.insert(_nodeGroupNames.begin(),groupName);
- //To do : update _faceGroups and _nodeGroups
+ std::vector< std::string >::iterator it = std::find(_nodeGroupNames.begin(), _nodeGroupNames.end(), groupName);
+ if(it == _nodeGroupNames.end())//No group named groupName
+ {
+ _nodeGroupNames.insert(_nodeGroupNames.end(),groupName);
+ _nodeGroupsIds.insert( _nodeGroupsIds.end(),nodeIds);
+ _nodeGroups.insert( _nodeGroups.end(), NULL);//No mesh created. Create one ?
+ }
+ else
+ {
+ std::vector< int > nodeGroupIds = _nodeGroupsIds[it-_nodeGroupNames.begin()];
+ nodeGroupIds.insert( nodeGroupIds.end(), nodeIds.begin(), nodeIds.end());
+ /* Detect and erase duplicates node ids */
+ sort( nodeGroupIds.begin(), nodeGroupIds.end() );
+ nodeGroupIds.erase( unique( nodeGroupIds.begin(), nodeGroupIds.end() ), nodeGroupIds.end() );
+ _nodeGroupsIds[it-_nodeGroupNames.begin()] = nodeGroupIds;
+ }
}
}
void
-Mesh::setGroupAtPlan(double value, int direction, double eps, std::string groupName)
+Mesh::setGroupAtPlan(double value, int direction, double eps, std::string groupName, bool isBoundaryGroup)
{
- int nbFace=getNumberOfFaces();
- bool flag=false;
- for (int iface=0;iface<nbFace;iface++)
- {
- double cord=_faces[iface].getBarryCenter()[direction];
- if (abs(cord-value)<eps)
+ std::vector< int > faceIds(0), nodeIds(0);
+ double cord;
+
+ /* Construction of the face group */
+ if(isBoundaryGroup)
+ for(int i=0; i<_boundaryFaceIds.size(); i++)
+ {
+ cord=_faces[_boundaryFaceIds[i]].getBarryCenter()[direction];
+ if (abs(cord-value)<eps)
+ {
+ faceIds.insert(faceIds.end(),_boundaryFaceIds[i]);
+ _faces[_boundaryFaceIds[i]].setGroupName(groupName);
+ }
+ }
+ else
+ for (int iface=0;iface<_numberOfFaces;iface++)
{
- _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);
- }
+ cord=_faces[iface].getBarryCenter()[direction];
+ if (abs(cord-value)<eps)
+ {
+ faceIds.insert(faceIds.end(),iface);
+ _faces[iface].setGroupName(groupName);
+ }
+ }
- flag=true;
+ /* Construction of the node group */
+ if(isBoundaryGroup)
+ for(int i=0; i<_boundaryNodeIds.size(); i++)
+ {
+ cord=_nodes[_boundaryNodeIds[i]].getPoint()[direction];
+ if (abs(cord-value)<eps)
+ {
+ nodeIds.insert(nodeIds.end(),_boundaryNodeIds[i]);
+ _nodes[_boundaryNodeIds[i]].setGroupName(groupName);
+ }
+ }
+ else
+ for (int inode=0;inode<_numberOfNodes;inode++)
+ {
+ cord=_nodes[inode].getPoint()[direction];
+ if (abs(cord-value)<eps)
+ {
+ nodeIds.insert(nodeIds.end(),inode);
+ _nodes[inode].setGroupName(groupName);
+ }
+ }
+
+ 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
+ {
+ std::vector< int > faceGroupIds = _faceGroupsIds[it-_faceGroupNames.begin()];
+ faceGroupIds.insert( faceGroupIds.end(), faceIds.begin(), faceIds.end());
+ /* Detect and erase duplicates face ids */
+ sort( faceGroupIds.begin(), faceGroupIds.end() );
+ faceGroupIds.erase( unique( faceGroupIds.begin(), faceGroupIds.end() ), faceGroupIds.end() );
+ _faceGroupsIds[it-_faceGroupNames.begin()] = faceGroupIds;
}
}
- if (flag)
+ if (nodeIds.size()>0)
{
- _faceGroupNames.insert(_faceGroupNames.begin(),groupName);
- _nodeGroupNames.insert(_nodeGroupNames.begin(),groupName);
- //To do : update _faceGroups, _nodeGroups
+ 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;
+ }
}
}
if(_indexFacePeriodicSet)
return;
- double eps=1.E-10;
-
for (int indexFace=0;indexFace<_boundaryFaceIds.size() ; indexFace++)
{
Face my_face=_faces[_boundaryFaceIds[indexFace]];
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
+ 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)<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 )
+ && ( !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;
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
+ 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)<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 )
+ && ( !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;
void
Mesh::setGroups( const MEDFileUMesh* medmesh, MEDCouplingUMesh* mu)
{
- //Searching for face groups
+ int nbCellsSubMesh, nbNodesSubMesh;
+ bool foundFace, foundNode;
+
+ /* Searching for face groups */
vector<string> faceGroups=medmesh->getGroupsNames() ;
for (unsigned int i=0;i<faceGroups.size();i++ )
vector<mcIdType>::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);
+ m->unPolyze();
+ m->sortCellsInMEDFileFrmt( );
+ nbCellsSubMesh=m->getNumberOfCells();
+
+ _faceGroups.insert(_faceGroups.end(),m);//Vector of group meshes
+ _faceGroupNames.insert(_faceGroupNames.end(),groupName);//Vector of group names
+ _faceGroupsIds.insert(_faceGroupsIds.end(),std::vector<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();
-
- int nbCellsSubMesh=m->getNumberOfCells();
+ /* Face identification */
for (int ic(0), k(0); ic<nbCellsSubMesh; ic++, k+=_spaceDim)
{
vector<double> coorBaryXyz(3,0);
coorBaryXyz[d] = coorBary[k+d];
Point p1(coorBaryXyz[0],coorBaryXyz[1],coorBaryXyz[2]) ;
- int flag=0;
- /* double min=INFINITY;
- Point p3;
- int closestFaceNb;*/
+ foundFace=false;
for (int iface=0;iface<_numberOfFaces;iface++ )
{
Point p2=_faces[iface].getBarryCenter();
- /*if(groupName=="Wall" and p1.distance(p2)<min)
- {
- min=p1.distance(p2);
- p3=p2;
- closestFaceNb=iface;
- }*/
- if(p1.distance(p2)<1.E-10)
+ if(p1.distance(p2)<_epsilon)
{
_faces[iface].setGroupName(groupName);
- flag=1;
+ _faceGroupsIds[_faceGroupsIds.size()-1][ic]=iface;
+ foundFace=true;
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");
+ if (not foundFace)
+ throw CdmathException("No face found for group " + groupName );
}
baryCell->decrRef();
//m->decrRef();
}
}
- //Searching for node groups
+ /* Searching for node groups */
vector<string> nodeGroups=medmesh->getGroupsOnSpecifiedLev(1) ;
for (unsigned int i=0;i<nodeGroups.size();i++ )
if(nodeids!=NULL)
{
- cout<<"Boundary node group named "<< groupName << " found"<<endl;
+ _nodeGroups.insert(_nodeGroups.end(),nodeGroup);
+ _nodeGroupNames.insert(_nodeGroupNames.end(),groupName);
- _nodeGroups.insert(_nodeGroups.begin(),nodeGroup);
- _nodeGroupNames.insert(_nodeGroupNames.begin(),groupName);
-
- int nbNodesSubMesh=nodeGroup->getNumberOfTuples();//nodeGroup->getNbOfElems();
+ nbNodesSubMesh=nodeGroup->getNumberOfTuples();//nodeGroup->getNbOfElems();
DataArrayDouble *coo = mu->getCoords() ;
const double *cood=coo->getConstPointer();
+ _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);
coorP[d] = cood[nodeids[ic]*_spaceDim+d];
Point p1(coorP[0],coorP[1],coorP[2]) ;
- int flag=0;
+ foundNode=false;
for (int inode=0;inode<_numberOfNodes;inode++ )
{
Point p2=_nodes[inode].getPoint();
- if(p1.distance(p2)<1.E-10)
+ if(p1.distance(p2)<_epsilon)
{
_nodes[inode].setGroupName(groupName);
- flag=1;
+ _nodeGroupsIds[_nodeGroupsIds.size()-1][ic]=inode;
+ foundNode=true;
break;
}
}
- if (flag==0)
- throw CdmathException("No node belonging to group " + groupName + " found");
+ if (not foundNode)
+ throw CdmathException("No node found for group " + groupName );
}
}
}
Mesh::setMesh( void )
//----------------------------------------------------------------------
{
+ /* This is the main function translating medcouplingumesh info into Mesh class to be used when designing numerical methods
+ * We need the level 0 mesh to extract the cell-node connectvity
+ * We need the level -1 mesh to extract the cell-face and face-node connectivities (use o build descending connectivity)
+ * Be careful : the nodes in the medcoupling mesh are not necessarily all conected to a cell/face.
+ * Mesh class discard isolated nodes, hence the number of nodes in Mesh class can be lower than the number of nodes in medcouplingumesh.
+ */
+
DataArrayIdType *desc = DataArrayIdType::New();
DataArrayIdType *descI = DataArrayIdType::New();
DataArrayIdType *revDesc = DataArrayIdType::New();
DataArrayIdType *revDescI = DataArrayIdType::New();
MEDCouplingUMesh* mu = _mesh->buildUnstructured();
-
+ MEDCouplingUMesh* mu2;//mesh of dimension N-1 containing the cell interfaces->cell/face connectivity
+
mu->unPolyze();
- /* Save the boundary mesh for later use*/
- _boundaryMesh = mu->computeSkin();
+ mu->sortCellsInMEDFileFrmt( );
- MEDCouplingUMesh* mu2=mu->buildDescendingConnectivity2(desc,descI,revDesc,revDescI);//mesh of dimension N-1 containing the cell interfaces
-
- const mcIdType *tmp = desc->getConstPointer();
+ if(_meshDim<2)
+ mu2=mu->buildDescendingConnectivity(desc,descI,revDesc,revDescI);
+ else
+ mu2=mu->buildDescendingConnectivity2(desc,descI,revDesc,revDescI);
+
+ const mcIdType *tmp = desc->getConstPointer();//Lists the faces surrounding each cell
const mcIdType *tmpI=descI->getConstPointer();
- const mcIdType *tmpA =revDesc->getConstPointer();
+ const mcIdType *tmpA =revDesc->getConstPointer();//Lists the cells surrounding each face
const mcIdType *tmpAI=revDescI->getConstPointer();
- //const int *work=tmp+tmpI[id];//corresponds to buildDescendingConnectivity
-
//Test du type d'éléments contenus dans le maillage afin d'éviter les éléments contenant des points de gauss
_eltsTypes=mu->getAllGeoTypesSorted();
for(int i=0; i<_eltsTypes.size();i++)
}
DataArrayDouble *baryCell = mu->computeCellCenterOfMass() ;
- const double *coorBary=baryCell->getConstPointer();
+ const double *coorBary=baryCell->getConstPointer();//Used for cell center coordinates
MEDCouplingFieldDouble* fields=mu->getMeasureField(true);
DataArrayDouble *surface = fields->getArray();
- const double *surf=surface->getConstPointer();
+ const double *surf=surface->getConstPointer();//Used for cell lenght/surface/volume
DataArrayDouble *coo = mu->getCoords() ;
- const double *cood=coo->getConstPointer();
+ const double *cood=coo->getConstPointer();//Used for nodes coordinates
DataArrayIdType *revNode =DataArrayIdType::New();
DataArrayIdType *revNodeI=DataArrayIdType::New();
mu->getReverseNodalConnectivity(revNode,revNodeI) ;
- const mcIdType *tmpN =revNode->getConstPointer();
+ const mcIdType *tmpN =revNode->getConstPointer();//Used to know which cells surround a given node
const mcIdType *tmpNI=revNodeI->getConstPointer();
DataArrayIdType *revCell =DataArrayIdType::New();
DataArrayIdType *revCellI=DataArrayIdType::New();
- mu2->getReverseNodalConnectivity(revCell,revCellI) ;
- const mcIdType *tmpC =revCell->getConstPointer();
+ mu2->getReverseNodalConnectivity(revCell,revCellI);
+ const mcIdType *tmpC =revCell->getConstPointer();//Used to know which faces surround a given node
const mcIdType *tmpCI=revCellI->getConstPointer();
const DataArrayIdType *nodal = mu2->getNodalConnectivity() ;
const DataArrayIdType *nodalI = mu2->getNodalConnectivityIndex() ;
- const mcIdType *tmpNE =nodal->getConstPointer();
+ const mcIdType *tmpNE =nodal->getConstPointer();//Used to know which nodes surround a given face
const mcIdType *tmpNEI=nodalI->getConstPointer();
_numberOfCells = mu->getNumberOfCells() ;
_cells = new Cell[_numberOfCells] ;
- _numberOfNodes = mu->getNumberOfNodes() ;
- _nodes = new Node[_numberOfNodes] ;
+ _numberOfNodes = mu->getNumberOfNodes() ;//This number may include isolated nodes that will not be loaded. The number will be updated during nodes constructions
+ _nodes = new Node[_numberOfNodes] ;//This array may be resized if isolated nodes are found
_numberOfFaces = mu2->getNumberOfCells();
_faces = new Face[_numberOfFaces] ;
}
// _cells, _nodes and _faces initialization:
- if (_spaceDim == 1)
+ if (_meshDim == 1)
{
+ double xn, yn=0., zn=0.;//Components of the normal vector at a cell interface
+ double norm;
for( int id=0;id<_numberOfCells;id++ )
{
- const mcIdType *work=tmp+tmpI[id];
- int nbFaces=tmpI[id+1]-tmpI[id];
-
- int nbVertices=mu->getNumberOfNodesInCell(id) ;
-
- Cell ci( nbVertices, nbFaces, surf[id], Point(coorBary[id], 0.0, 0.0) ) ;
-
+ Point p(0.0,0.0,0.0) ;
+ for(int idim=0; idim<_spaceDim; idim++)
+ p[idim]=coorBary[id*_spaceDim+idim];
+
+ mcIdType nbVertices=mu->getNumberOfNodesInCell(id) ;//should be equal to 2
+ assert( nbVertices==2);
std::vector<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;
+ }
+
+ 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[id] = ci ;
}
-
- for( int id(0), k(0); id<_numberOfNodes; id++, k+=_spaceDim)
- {
- Point p(cood[k], 0.0, 0.0) ;
- const mcIdType *workc=tmpN+tmpNI[id];
- mcIdType nbCells=tmpNI[id+1]-tmpNI[id];
-
- const mcIdType *workf=tmpC+tmpCI[id];
- mcIdType nbFaces=tmpCI[id+1]-tmpCI[id];
- const mcIdType *workn=tmpA+tmpAI[id];
- mcIdType 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 ;
- }
- _boundaryNodeIds.push_back(0);
- _boundaryNodeIds.push_back(_numberOfNodes-1);
-
- for(int id(0), k(0); id<_numberOfFaces; id++, k+=_spaceDim)
+
+ for( int id(0); id<_numberOfFaces; id++ )
{
- Point p(cood[k], 0.0, 0.0) ;
- const mcIdType *workc=tmpA+tmpAI[id];
- mcIdType nbCells=tmpAI[id+1]-tmpAI[id];
-
const mcIdType *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]) ;
-
+ 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 b 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[id] = fi ;
}
- _boundaryFaceIds.push_back(0);
- _boundaryFaceIds.push_back(_numberOfFaces-1);
+
+ int correctNbNodes=0;
+ for( int id=0;id<_numberOfNodes;id++ )
+ {
+ 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[workf[el]].isBorder()) ;
+ if(vi.isBorder())
+ _boundaryNodeIds.push_back(id);
+ _nodes[id] = vi ;
+ }
+ }
+ if( _numberOfNodes!=correctNbNodes)
+ cout<<"Found isolated nodes : correctNbNodes= "<<correctNbNodes<<", _numberOfNodes= "<<_numberOfNodes<<endl;
}
- 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;
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
mcIdType indexFace0=abs(work[0])-1;//=work[0] since Fortran type numbering was used, and negative sign means anticlockwise numbering
_cells[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]) ;
}
/*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 mcIdType *workc=tmpN+tmpNI[id];
mcIdType nbCells=tmpNI[id+1]-tmpNI[id];
- 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
- 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 ;
+
+ 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[workf[el]].isBorder()) ;
+ if(vi.isBorder())
+ _boundaryNodeIds.push_back(id);
+ _nodes[id] = vi ;
+ }
}
+ if( _numberOfNodes!=correctNbNodes)
+ cout<<"Found isolated nodes : correctNbNodes= "<<correctNbNodes<<", _numberOfNodes= "<<_numberOfNodes<<endl;
if(_spaceDim==_meshDim)
fieldn->decrRef();
_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 )
-//----------------------------------------------------------------------
-{
- if(nx<=0)
- throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx) : nx <= 0");
- if(xmin>=xmax)
- throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx) : xmin >= xmax");
-
- double dx = (xmax - xmin)/nx ;
-
- _spaceDim = 1 ;
- _meshDim = 1 ;
- _name=meshName;
- _isStructured = true;
- _indexFacePeriodicSet=false;
-
- _xMin=xmin;
- _xMax=xmax;
- _yMin=0.;
- _yMax=0.;
- _zMin=0.;
- _zMax=0.;
-
- _dxyz.resize(_spaceDim);
- _dxyz[0]=dx;
- _nxyz.resize(_spaceDim);
- _nxyz[0]=nx;
-
- double *originPtr = new double[_spaceDim];
- double *dxyzPtr = new double[_spaceDim];
- mcIdType *nodeStrctPtr = new mcIdType[_spaceDim];
-
- originPtr[0]=xmin;
- nodeStrctPtr[0]=nx+1;
- dxyzPtr[0]=dx;
-
- _mesh=MEDCouplingIMesh::New(meshName,
- _spaceDim,
- nodeStrctPtr,
- nodeStrctPtr+_spaceDim,
- originPtr,
- originPtr+_spaceDim,
- dxyzPtr,
- dxyzPtr+_spaceDim);
- delete [] originPtr;
- delete [] dxyzPtr;
- delete [] nodeStrctPtr;
-
- DataArrayIdType *desc=DataArrayIdType::New();
- DataArrayIdType *descI=DataArrayIdType::New();
- DataArrayIdType *revDesc=DataArrayIdType::New();
- DataArrayIdType *revDescI=DataArrayIdType::New();
- MEDCouplingUMesh* mu=_mesh->buildUnstructured();
- MEDCouplingUMesh *mu2=mu->buildDescendingConnectivity(desc,descI,revDesc,revDescI);
-
- const mcIdType *tmp=desc->getConstPointer();
- const mcIdType *tmpI=descI->getConstPointer();
-
- const mcIdType *tmpA =revDesc->getConstPointer();
- const mcIdType *tmpAI=revDescI->getConstPointer();
-
- _eltsTypes=mu->getAllGeoTypesSorted();
-
- DataArrayDouble *baryCell = mu->computeCellCenterOfMass() ;
- const double *coorBary=baryCell->getConstPointer();
-
- _numberOfCells = _mesh->getNumberOfCells() ;
- _cells = new Cell[_numberOfCells] ;
-
- _numberOfNodes = mu->getNumberOfNodes() ;
- _nodes = new Node[_numberOfNodes] ;
-
- _numberOfFaces = _numberOfNodes;
- _faces = new Face[_numberOfFaces] ;
-
- _numberOfEdges = _numberOfCells;
-
- MEDCouplingFieldDouble* fieldl=mu->getMeasureField(true);
- DataArrayDouble *longueur = fieldl->getArray();
- const double *lon=longueur->getConstPointer();
-
- DataArrayDouble *coo = mu->getCoords() ;
- const double *cood=coo->getConstPointer();
-
- int comp=0;
- for( int id=0;id<_numberOfCells;id++ )
- {
- int nbVertices=mu->getNumberOfNodesInCell(id) ;
- Point p(coorBary[id],0.0,0.0) ;
- Cell ci( nbVertices, nbVertices, lon[id], p ) ;
-
- std::vector<mcIdType> nodeIdsOfCell ;
- mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
- for( int el=0;el<nbVertices;el++ )
- {
- ci.addNodeId(el,nodeIdsOfCell[el]) ;
- ci.addFaceId(el,nodeIdsOfCell[el]) ;
- }
-
- double xn = (cood[nodeIdsOfCell[nbVertices-1]] - cood[nodeIdsOfCell[0]] > 0.0) ? -1.0 : 1.0;
-
- mcIdType nbFaces=tmpI[id+1]-tmpI[id];
- const mcIdType *work=tmp+tmpI[id];
-
- for( int el=0;el<nbFaces;el++ )
- {
- ci.addNormalVector(el,xn,0.0,0.0) ;
- ci.addFaceId(el,work[el]) ;
- xn = - xn;
- }
-
- _cells[id] = ci ;
-
- comp=comp+2;
- }
-
- //Suppress the following since tmpN=tmpA
- DataArrayIdType *revNode=DataArrayIdType::New();
- DataArrayIdType *revNodeI=DataArrayIdType::New();
- mu->getReverseNodalConnectivity(revNode,revNodeI) ;
- const mcIdType *tmpN=revNode->getConstPointer();
- const mcIdType *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 mcIdType *workc=tmpN+tmpNI[id];
- mcIdType nbCells=tmpNI[id+1]-tmpNI[id];
- mcIdType nbFaces=1;
- const mcIdType *workn=tmpA+tmpAI[id];
- mcIdType nbNeighbourNodes=tmpAI[id+1]-tmpAI[id];
-
- Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
- for( int el=0;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 )
//----------------------------------------------------------------------
_spaceDim = 1 ;
_meshDim = 1 ;
_name=meshName;
+ _epsilon=1e-6;
_xMin=points[0];
_xMax=points[nx-1];
_yMin=0.;
delete [] coords, nodal_con;
coords_arr->decrRef();
- _mesh=mesh1d;
-
- DataArrayIdType *desc=DataArrayIdType::New();
- DataArrayIdType *descI=DataArrayIdType::New();
- DataArrayIdType *revDesc=DataArrayIdType::New();
- DataArrayIdType *revDescI=DataArrayIdType::New();
- MEDCouplingUMesh* mu=_mesh->buildUnstructured();
- MEDCouplingUMesh *mu2=mu->buildDescendingConnectivity(desc,descI,revDesc,revDescI);
-
- DataArrayDouble *baryCell = mu->computeCellCenterOfMass() ;
- const double *coorBary=baryCell->getConstPointer();
+ _mesh=mesh1d->buildUnstructured();//To enable writeMED. Because we declared the mesh as unstructured, we decide to build the unstructured data (not mandatory)
+ _meshNotDeleted=true;
- _numberOfCells = _mesh->getNumberOfCells() ;
- _cells = new Cell[_numberOfCells] ;
+ setMesh();
+}
- _numberOfEdges = _numberOfCells;
+//----------------------------------------------------------------------
+Mesh::Mesh( double xmin, double xmax, int nx, std::string meshName )
+//----------------------------------------------------------------------
+{
+ if(nx<=0)
+ throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx) : nx <= 0");
+ if(xmin>=xmax)
+ throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx) : xmin >= xmax");
- _eltsTypes=mu->getAllGeoTypesSorted();
+ double dx = (xmax - xmin)/nx ;
- MEDCouplingFieldDouble* fieldl=mu->getMeasureField(true);
- DataArrayDouble *longueur = fieldl->getArray();
- const double *lon=longueur->getConstPointer();
+ _spaceDim = 1 ;
+ _meshDim = 1 ;
+ _name=meshName;
+ _epsilon=1e-6;
+ _isStructured = true;
+ _indexFacePeriodicSet=false;
- int comp=0;
- for( int id=0;id<_numberOfCells;id++ )
- {
- int nbVertices=mu->getNumberOfNodesInCell(id) ;
- Point p(coorBary[id],0.0,0.0) ;
- Cell ci( nbVertices, nbVertices, lon[id], p ) ;
+ _xMin=xmin;
+ _xMax=xmax;
+ _yMin=0.;
+ _yMax=0.;
+ _zMin=0.;
+ _zMax=0.;
- std::vector<mcIdType> nodeIdsOfCell ;
- mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
- for( int el=0;el<nbVertices;el++ )
- {
- ci.addNodeId(el,nodeIdsOfCell[el]) ;
- ci.addFaceId(el,nodeIdsOfCell[el]) ;
- }
- _cells[id] = ci ;
- comp=comp+2;
- }
+ _dxyz.resize(_spaceDim);
+ _dxyz[0]=dx;
+ _nxyz.resize(_spaceDim);
+ _nxyz[0]=nx;
+ double *originPtr = new double[_spaceDim];
+ double *dxyzPtr = new double[_spaceDim];
+ mcIdType *nodeStrctPtr = new mcIdType[_spaceDim];
- //Suppress the following since tmpN=tmpA
- DataArrayIdType *revNode=DataArrayIdType::New();
- DataArrayIdType *revNodeI=DataArrayIdType::New();
- mu->getReverseNodalConnectivity(revNode,revNodeI) ;
- const mcIdType *tmpN=revNode->getConstPointer();
- const mcIdType *tmpNI=revNodeI->getConstPointer();
+ originPtr[0]=xmin;
+ nodeStrctPtr[0]=nx+1;
+ dxyzPtr[0]=dx;
- _numberOfNodes = mu->getNumberOfNodes() ;
- _nodes = new Node[_numberOfNodes] ;
- _numberOfFaces = _numberOfNodes;
- _faces = new Face[_numberOfFaces] ;
+ _mesh=MEDCouplingIMesh::New(meshName,
+ _spaceDim,
+ nodeStrctPtr,
+ nodeStrctPtr+_spaceDim,
+ originPtr,
+ originPtr+_spaceDim,
+ dxyzPtr,
+ dxyzPtr+_spaceDim);
+ _meshNotDeleted=true;//Because the mesh is structured cartesian : no data in memory. No nodes and cell coordinates stored
- for( int id=0;id<_numberOfNodes;id++ )
- {
- std::vector<double> coo ;
- mu->getCoordinatesOfNode(id,coo);
- Point p(coo[0],0.0,0.0) ;
- const mcIdType *workc=tmpN+tmpNI[id];
- mcIdType nbCells=tmpNI[id+1]-tmpNI[id];
- mcIdType nbFaces=1;
- const mcIdType *workn=tmpN+tmpNI[id];
- mcIdType nbNeighbourNodes=tmpNI[id+1]-tmpNI[id];
- Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
- int nbVertices=1;
- /* provisoire !!!!!!!!!!!!*/
- // Point pf(0.0,0.0,0.0) ;
- Face fi( nbVertices, nbCells, 0.0, p, 0., 0., 0. ) ;
-
- for( int el=0;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);
+ delete [] originPtr;
+ delete [] dxyzPtr;
+ delete [] nodeStrctPtr;
- fieldl->decrRef();
- baryCell->decrRef();
- desc->decrRef();
- descI->decrRef();
- revDesc->decrRef();
- revDescI->decrRef();
- revNode->decrRef();
- revNodeI->decrRef();
- mu2->decrRef();
- mu->decrRef();
+ setMesh();
}
+
//----------------------------------------------------------------------
Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, int split_to_triangles_policy, std::string meshName)
//----------------------------------------------------------------------
_spaceDim = 2 ;
_meshDim = 2 ;
_name=meshName;
+ _epsilon=1e-6;
_indexFacePeriodicSet=false;
- _isStructured = true;
_nxyz.resize(_spaceDim);
_nxyz[0]=nx;
_nxyz[1]=ny;
originPtr+_spaceDim,
dxyzPtr,
dxyzPtr+_spaceDim);
+ _meshNotDeleted=true;//Because the mesh is structured cartesian : no data in memory. No nodes and cell coordinates stored
+ _isStructured = true;
if(split_to_triangles_policy==0 || split_to_triangles_policy==1)
{
_mesh=_mesh->buildUnstructured();
_mesh->simplexize(split_to_triangles_policy);
+ _meshNotDeleted=true;//Now the mesh is unstructured and stored with nodes and cell coordinates
+ _isStructured = false;
}
else if (split_to_triangles_policy != -1)
{
_spaceDim = 3;
_meshDim = 3;
_name=meshName;
- _indexFacePeriodicSet=false;
+ _epsilon=1e-6;
_xMin=xmin;
_xMax=xmax;
_yMin=ymin;
double dy = (ymax - ymin)/ny ;
double dz = (zmax - zmin)/nz ;
- _isStructured = true;
_dxyz.resize(_spaceDim);
_dxyz[0]=dx;
_dxyz[1]=dy;
originPtr+_spaceDim,
dxyzPtr,
dxyzPtr+_spaceDim);
+ _meshNotDeleted=true;//Because the mesh is structured cartesian : no data in memory. Nno nodes and cell coordinates stored
+ _isStructured = true;
if( split_to_tetrahedra_policy == 0 )
{
_mesh=_mesh->buildUnstructured();
_mesh->simplexize(INTERP_KERNEL::PLANAR_FACE_5);
+ _meshNotDeleted=true;//Now the mesh is unstructured and stored with nodes and cell coordinates
+ _isStructured = false;
}
else if( split_to_tetrahedra_policy == 1 )
{
_mesh=_mesh->buildUnstructured();
_mesh->simplexize(INTERP_KERNEL::PLANAR_FACE_6);
+ _meshNotDeleted=true;//Now the mesh is unstructured and stored with nodes and cell coordinates
+ _isStructured = false;
}
else if ( split_to_tetrahedra_policy != -1 )
{
_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() ;
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[i];
+ //Careful with mesh with junctions : do not just take Ci.getNumberOfFaces()
+ nbNeib=0;
+ for(int j=0; j<Ci.getNumberOfFaces(); j++)
+ nbNeib+=_faces[Ci.getFacesId()[j]].getNumberOfCells()-1;//Without junction this would be +=1
+
+ if(result < nbNeib)
+ result=nbNeib;
+ }
}
else if(type==NODES)
{
Mesh::writeVTK ( const std::string fileName ) const
//----------------------------------------------------------------------
{
+ if( !_isStructured && !_meshNotDeleted )
+ throw CdmathException("Mesh::writeVTK : Cannot save mesh : no MEDCouplingUMesh loaded");
+
string fname=fileName+".vtu";
_mesh->writeVTK(fname.c_str()) ;
}
Mesh::writeMED ( const std::string fileName ) const
//----------------------------------------------------------------------
{
+ if( !_meshNotDeleted )
+ throw CdmathException("Mesh::writeMED : Cannot save mesh : no MEDCouplingUMesh loaded");
+
string fname=fileName+".med";
MEDCouplingUMesh* mu=_mesh->buildUnstructured();
MEDCoupling::WriteUMesh(fname.c_str(),mu,true);
+ /* Try to save mesh groups */
//MEDFileUMesh meshMEDFile;
//meshMEDFile.setMeshAtLevel(0,mu);
//for(int i=0; i< _faceGroups.size(); i++)
//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::deleteMEDCouplingUMesh()
+{
+ if(_meshNotDeleted)
+ {
+ (_mesh.retn())->decrRef();
+ _meshNotDeleted=false;
+ }
+ else
+ throw CdmathException("Mesh::deleteMEDCouplingMesh() : mesh is not loaded");
+};
