4 * Created on: 22 janv. 2012
20 #include "MEDCouplingIMesh.hxx"
21 #include "MEDCouplingUMesh.hxx"
22 #include "MEDCouplingFieldDouble.hxx"
24 #include "MEDFileMesh.hxx"
25 #include "MEDLoader.hxx"
27 #include "CdmathException.hxx"
29 using namespace MEDCoupling;
32 //----------------------------------------------------------------------
34 //----------------------------------------------------------------------
37 _meshNotDeleted=false;
50 _boundaryFaceIds.resize(0);
51 _boundaryNodeIds.resize(0);
52 _faceGroupNames.resize(0);
53 _faceGroups.resize(0);
54 _faceGroupsIds.resize(0);
55 _nodeGroupNames.resize(0);
56 _nodeGroups.resize(0);
57 _nodeGroupsIds.resize(0);
58 _indexFacePeriodicSet=false;
63 //----------------------------------------------------------------------
65 //----------------------------------------------------------------------
67 //for(int i=0; i< _faceGroups.size(); i++)
68 // _faceGroups[i]->decrRef();
69 // _nodeGroups[i]->decrRef();
71 (_mesh.retn())->decrRef();
75 Mesh::getName( void ) const
80 Mesh::Mesh( MEDCoupling::MCAuto<const MEDCoupling::MEDCouplingMesh> mesh )
82 _spaceDim=mesh->getSpaceDimension();
83 _meshDim=mesh->getMeshDimension();
84 _name=mesh->getName();
86 _indexFacePeriodicSet=false;
89 _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
91 MEDCoupling::MEDCouplingStructuredMesh* structuredMesh = dynamic_cast<MEDCoupling::MEDCouplingStructuredMesh*> (_mesh.retn());
95 _nxyz=structuredMesh->getCellGridStructure();
103 //----------------------------------------------------------------------
104 Mesh::Mesh( const std::string filename, const std::string & meshName, int meshLevel)
105 //----------------------------------------------------------------------
107 readMeshMed(filename, meshName, meshLevel);
110 //----------------------------------------------------------------------
111 Mesh::Mesh( const Mesh& mesh )
112 //----------------------------------------------------------------------
114 _spaceDim = mesh.getSpaceDimension() ;
115 _meshDim = mesh.getMeshDimension() ;
116 _name=mesh.getName();
117 _epsilon=mesh.getComparisonEpsilon();
118 _isStructured=mesh.isStructured();
120 _nxyz = mesh.getCellGridStructure() ;
121 _numberOfNodes = mesh.getNumberOfNodes();
122 _numberOfFaces = mesh.getNumberOfFaces();
123 _numberOfCells = mesh.getNumberOfCells();
124 _numberOfEdges = mesh.getNumberOfEdges();
126 _faceGroupNames = mesh.getNameOfFaceGroups() ;
127 _faceGroups = mesh.getFaceGroups() ;
128 _nodeGroupNames = mesh.getNameOfNodeGroups() ;
129 _nodeGroups = mesh.getNodeGroups() ;
131 _nodes = mesh.getNodes() ;
132 _faces = mesh.getFaces() ;
133 _cells = mesh.getCells() ;
135 _indexFacePeriodicSet= mesh.isIndexFacePeriodicSet();
136 if(_indexFacePeriodicSet)
137 _indexFacePeriodicMap=mesh.getIndexFacePeriodic();
139 _boundaryFaceIds=mesh.getBoundaryFaceIds();
140 _boundaryNodeIds=mesh.getBoundaryNodeIds();
142 _eltsTypes=mesh.getElementTypes();
143 _eltsTypesNames=mesh.getElementTypesNames();
145 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
148 _meshNotDeleted=mesh.meshNotDeleted();
151 //----------------------------------------------------------------------
153 Mesh::readMeshMed( const std::string filename, const std::string & meshName, int meshLevel)
154 //----------------------------------------------------------------------
158 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)
160 m=MEDFileMesh::New(filename.c_str(), meshName.c_str());//seeks the mesh named meshName in the file
162 _mesh=m->getMeshAtLevel(meshLevel);
163 _mesh->checkConsistencyLight();
164 _mesh->setName(_mesh->getName());
165 _meshDim=_mesh->getMeshDimension();
166 _spaceDim=_mesh->getSpaceDimension();
167 _name=_mesh->getName();
169 _indexFacePeriodicSet=false;
170 _meshNotDeleted=true;
172 MEDCoupling::MEDCouplingStructuredMesh* structuredMesh = dynamic_cast<MEDCoupling::MEDCouplingStructuredMesh*> (_mesh.retn());
176 _nxyz=structuredMesh->getCellGridStructure();
181 MEDCouplingUMesh* mu = setMesh();
182 setNodeGroups(m, mu);//Works for both cartesan and unstructured meshes
183 MEDFileUMesh *umedfile=dynamic_cast< MEDFileUMesh * > (m);
185 setFaceGroups(umedfile, mu);//Works only for unstructured meshes
187 cout<<endl<< "Loaded file "<< filename<<endl;
188 cout<<"Mesh name = "<<m->getName()<<", mesh dim = "<< _meshDim<< ", space dim = "<< _spaceDim<< ", nb cells= "<<getNumberOfCells()<< ", nb nodes= "<<getNumberOfNodes()<<endl;
193 //----------------------------------------------------------------------
194 Mesh::Mesh( std::vector<double> points, std::string meshName )
195 //----------------------------------------------------------------------
197 int nx=points.size();
200 throw CdmathException("Mesh::Mesh( vector<double> points, string meshName) : nx < 2, vector should contain at least two values");
202 while( i<nx-1 && points[i+1]>points[i] )
206 //cout<< points << endl;
207 throw CdmathException("Mesh::Mesh( vector<double> points, string meshName) : vector values should be sorted");
214 _indexFacePeriodicSet=false;
216 MEDCouplingUMesh * mesh1d = MEDCouplingUMesh::New(meshName, 1);
217 mesh1d->allocateCells(nx - 1);
218 double * coords = new double[nx];
219 mcIdType nodal_con[2];
221 for(int i=0; i<nx- 1 ; i++)
225 mesh1d->insertNextCell(INTERP_KERNEL::NORM_SEG2, 2, nodal_con);
226 coords[i+1]=points[i + 1];
228 mesh1d->finishInsertingCells();
230 DataArrayDouble * coords_arr = DataArrayDouble::New();
231 coords_arr->alloc(nx,1);
232 std::copy(coords,coords+nx,coords_arr->getPointer());
233 mesh1d->setCoords(coords_arr);
236 coords_arr->decrRef();
238 _mesh=mesh1d;//To enable writeMED. Because we declared the mesh as unstructured, we decide to build the unstructured data (not mandatory)
239 _meshNotDeleted=true;
240 _isStructured = false;
245 //----------------------------------------------------------------------
246 Mesh::Mesh( double xmin, double xmax, int nx, std::string meshName )
247 //----------------------------------------------------------------------
250 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx) : nx <= 0");
252 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx) : xmin >= xmax");
254 double dx = (xmax - xmin)/nx ;
260 _indexFacePeriodicSet=false;
262 _nxyz.resize(_spaceDim);
265 double originPtr[_spaceDim];
266 double dxyzPtr[_spaceDim];
267 mcIdType nodeStrctPtr[_spaceDim];
270 nodeStrctPtr[0]=nx+1;
273 _mesh=MEDCouplingIMesh::New(meshName,
276 nodeStrctPtr+_spaceDim,
281 _meshNotDeleted=true;
282 _isStructured = true;
287 //----------------------------------------------------------------------
288 Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, int split_to_triangles_policy, std::string meshName)
289 //----------------------------------------------------------------------
292 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : nx <= 0 or ny <= 0");
294 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : xmin >= xmax");
296 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : ymin >= ymax");
298 double dx = (xmax - xmin)/nx ;
299 double dy = (ymax - ymin)/ny ;
305 _indexFacePeriodicSet=false;
306 _nxyz.resize(_spaceDim);
310 double originPtr[_spaceDim];
311 double dxyzPtr[_spaceDim];
312 mcIdType nodeStrctPtr[_spaceDim];
316 nodeStrctPtr[0]=nx+1;
317 nodeStrctPtr[1]=ny+1;
321 _mesh=MEDCouplingIMesh::New(meshName,
324 nodeStrctPtr+_spaceDim,
329 _meshNotDeleted=true;
330 _isStructured = true;
332 if(split_to_triangles_policy==0 || split_to_triangles_policy==1)
334 _mesh=_mesh->buildUnstructured();//simplexize is not available for structured meshes
335 _mesh->simplexize(split_to_triangles_policy);
336 _isStructured = false;
338 else if (split_to_triangles_policy != -1)
340 cout<< "split_to_triangles_policy = "<< split_to_triangles_policy << endl;
341 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : Unknown splitting policy");
347 //----------------------------------------------------------------------
348 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)
349 //----------------------------------------------------------------------
351 if(nx<=0 || ny<=0 || nz<=0)
352 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");
354 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : xmin >= xmax");
356 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : ymin >= ymax");
358 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : zmin >= zmax");
365 double dx = (xmax - xmin)/nx ;
366 double dy = (ymax - ymin)/ny ;
367 double dz = (zmax - zmin)/nz ;
369 _nxyz.resize(_spaceDim);
374 double originPtr[_spaceDim];
375 double dxyzPtr[_spaceDim];
376 mcIdType nodeStrctPtr[_spaceDim];
381 nodeStrctPtr[0]=nx+1;
382 nodeStrctPtr[1]=ny+1;
383 nodeStrctPtr[2]=nz+1;
388 _mesh=MEDCouplingIMesh::New(meshName,
391 nodeStrctPtr+_spaceDim,
396 _meshNotDeleted=true;
397 _isStructured = true;
399 if( split_to_tetrahedra_policy == 0 )
401 _mesh=_mesh->buildUnstructured();//simplexize is not available for structured meshes
402 _mesh->simplexize(INTERP_KERNEL::PLANAR_FACE_5);
403 _isStructured = false;
405 else if( split_to_tetrahedra_policy == 1 )
407 _mesh=_mesh->buildUnstructured();//simplexize is not available for structured meshes
408 _mesh->simplexize(INTERP_KERNEL::PLANAR_FACE_6);
409 _isStructured = false;
411 else if ( split_to_tetrahedra_policy != -1 )
413 cout<< "split_to_tetrahedra_policy = "<< split_to_tetrahedra_policy << endl;
414 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");
420 //----------------------------------------------------------------------
422 Mesh::setMesh( void )
423 //----------------------------------------------------------------------
425 /* This is the main function translating medcouplingumesh info into Mesh class to be used when designing numerical methods
426 * We need the level 0 mesh to extract the cell-node connectvity
427 * We need the level -1 mesh to extract the cell-face and face-node connectivities (use o build descending connectivity)
428 * Be careful : the nodes in the medcoupling mesh are not necessarily all conected to a cell/face.
429 * Mesh class discard isolated nodes, hence the number of nodes in Mesh class can be lower than the number of nodes in medcouplingumesh.
432 DataArrayIdType *desc = DataArrayIdType::New();
433 DataArrayIdType *descI = DataArrayIdType::New();
434 DataArrayIdType *revDesc = DataArrayIdType::New();
435 DataArrayIdType *revDescI = DataArrayIdType::New();
436 MEDCouplingUMesh* mu = _mesh->buildUnstructured();
437 MEDCouplingUMesh* mu2;//mesh of dimension N-1 containing the cell interfaces->cell/face connectivity
440 mu->sortCellsInMEDFileFrmt( );
443 mu2=mu->buildDescendingConnectivity(desc,descI,revDesc,revDescI);
445 mu2=mu->buildDescendingConnectivity2(desc,descI,revDesc,revDescI);
447 const mcIdType *tmp = desc->getConstPointer();//Lists the faces surrounding each cell
448 const mcIdType *tmpI=descI->getConstPointer();
450 const mcIdType *tmpA =revDesc->getConstPointer();//Lists the cells surrounding each face
451 const mcIdType *tmpAI=revDescI->getConstPointer();
453 //Test du type d'éléments contenus dans le maillage afin d'éviter les éléments contenant des points de gauss
454 _eltsTypes=mu->getAllGeoTypesSorted();
455 for(int i=0; i<_eltsTypes.size();i++)
458 _eltsTypes[i]!= INTERP_KERNEL::NORM_POINT1 && _eltsTypes[i]!= INTERP_KERNEL::NORM_SEG2
459 && _eltsTypes[i]!= INTERP_KERNEL::NORM_TRI3 && _eltsTypes[i]!= INTERP_KERNEL::NORM_QUAD4
460 && _eltsTypes[i]!= INTERP_KERNEL::NORM_TETRA4 && _eltsTypes[i]!= INTERP_KERNEL::NORM_PYRA5
461 && _eltsTypes[i]!= INTERP_KERNEL::NORM_PENTA6 && _eltsTypes[i]!= INTERP_KERNEL::NORM_HEXA8
462 && _eltsTypes[i]!= INTERP_KERNEL::NORM_POLYGON&& _eltsTypes[i]!= INTERP_KERNEL::NORM_POLYHED
465 cout<< "Mesh " + mu->getName() + " contains an element of type " <<endl;
466 cout<< _eltsTypes[i]<<endl;
467 throw CdmathException("Mesh::setMesh : in order to avoid gauss points, mesh should contain elements of type NORM_POINT1, NORM_SEG2, NORM_TRI3, NORM_QUAD4, NORM_TETRA4, NORM_PYRA5, NORM_PENTA6, NORM_HEXA8, NORM_POLYGON, NORM_POLYHED");
471 DataArrayDouble *baryCell = mu->computeCellCenterOfMass() ;
472 const double *coorBary=baryCell->getConstPointer();//Used for cell center coordinates
474 MEDCouplingFieldDouble* fields=mu->getMeasureField(true);
475 DataArrayDouble *surface = fields->getArray();
476 const double *surf=surface->getConstPointer();//Used for cell lenght/surface/volume
478 DataArrayDouble *coo = mu->getCoords() ;
479 const double *cood=coo->getConstPointer();//Used for nodes coordinates
481 DataArrayIdType *revNode =DataArrayIdType::New();
482 DataArrayIdType *revNodeI=DataArrayIdType::New();
483 mu->getReverseNodalConnectivity(revNode,revNodeI) ;
484 const mcIdType *tmpN =revNode->getConstPointer();//Used to know which cells surround a given node
485 const mcIdType *tmpNI=revNodeI->getConstPointer();
487 DataArrayIdType *revCell =DataArrayIdType::New();
488 DataArrayIdType *revCellI=DataArrayIdType::New();
489 mu2->getReverseNodalConnectivity(revCell,revCellI);
490 const mcIdType *tmpC =revCell->getConstPointer();//Used to know which faces surround a given node
491 const mcIdType *tmpCI=revCellI->getConstPointer();
493 const DataArrayIdType *nodal = mu2->getNodalConnectivity() ;
494 const DataArrayIdType *nodalI = mu2->getNodalConnectivityIndex() ;
495 const mcIdType *tmpNE =nodal->getConstPointer();//Used to know which nodes surround a given face
496 const mcIdType *tmpNEI=nodalI->getConstPointer();
498 _numberOfCells = mu->getNumberOfCells() ;
499 _cells = std::shared_ptr<Cell>(new Cell[_numberOfCells], std::default_delete<Cell[]>()) ;
501 _numberOfNodes = mu->getNumberOfNodes() ;//This number may include isolated nodes that will not be loaded. The number will be updated during nodes constructions
502 _nodes = std::shared_ptr<Node>(new Node[_numberOfNodes], std::default_delete<Node[]>()) ;//This array may be resized if isolated nodes are found
504 _numberOfFaces = mu2->getNumberOfCells();
505 _faces = std::shared_ptr<Face>(new Face[_numberOfFaces], std::default_delete<Face[]>()) ;
507 _indexFacePeriodicSet=false;
509 //Definition used if _meshDim =3 to determine the edges
510 DataArrayIdType *desc2 =DataArrayIdType::New();
511 DataArrayIdType *descI2=DataArrayIdType::New();
512 DataArrayIdType *revDesc2 =DataArrayIdType::New();
513 DataArrayIdType *revDescI2=DataArrayIdType::New();
514 DataArrayIdType *revNode2 =DataArrayIdType::New();
515 DataArrayIdType *revNodeI2=DataArrayIdType::New();
516 const mcIdType *tmpN2 ;
517 const mcIdType *tmpNI2;
518 MEDCouplingUMesh* mu3;
521 _numberOfEdges = mu->getNumberOfCells();
522 else if (_meshDim == 2)
523 _numberOfEdges = mu2->getNumberOfCells();
526 mu3=mu2->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2);//1D mesh of segments
527 _numberOfEdges = mu3->getNumberOfCells();
528 mu3->getReverseNodalConnectivity(revNode2,revNodeI2) ;
529 tmpN2 =revNode2->getConstPointer();
530 tmpNI2=revNodeI2->getConstPointer();
533 // _cells, _nodes and _faces initialization:
536 double xn, yn=0., zn=0.;//Components of the normal vector at a cell interface
538 for( int id=0;id<_numberOfCells;id++ )
540 Point p(0.0,0.0,0.0) ;
541 for(int idim=0; idim<_spaceDim; idim++)
542 p[idim]=coorBary[id*_spaceDim+idim];
544 mcIdType nbVertices=mu->getNumberOfNodesInCell(id) ;//should be equal to 2
545 assert( nbVertices==2);
546 std::vector<mcIdType> nodeIdsOfCell ;
547 mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
549 mcIdType nbFaces=tmpI[id+1]-tmpI[id];//should be equal to 2
551 const mcIdType *work=tmp+tmpI[id];
553 /* compute the normal to the face */
554 xn = cood[nodeIdsOfCell[0]*_spaceDim ] - cood[nodeIdsOfCell[nbFaces-1]*_spaceDim ];
556 yn = cood[nodeIdsOfCell[0]*_spaceDim+1] - cood[nodeIdsOfCell[nbFaces-1]*_spaceDim+1];
558 zn = cood[nodeIdsOfCell[0]*_spaceDim+2] - cood[nodeIdsOfCell[nbFaces-1]*_spaceDim+2];
559 norm = sqrt(xn*xn+yn*yn+zn*zn);
561 throw CdmathException("!!! Mesh::setMesh Normal vector has norm 0 !!!");
569 Cell ci( nbVertices, nbFaces, surf[id], p ) ;//nbCells=nbFaces=2
570 for( int el=0;el<nbFaces;el++ )
572 ci.addNodeId(el,nodeIdsOfCell[el]) ;//global node number
573 ci.addNormalVector(el,xn,yn,zn) ;
574 ci.addFaceId(el,work[el]) ;
575 xn = - xn; yn=-yn; zn=-zn;
577 _cells.get()[id] = ci ;
580 for( int id(0); id<_numberOfFaces; id++ )
582 const mcIdType *workv=tmpNE+tmpNEI[id]+1;
583 mcIdType nbNodes= tmpNEI[id+1]-tmpNEI[id]-1;//Normally equal to 1.
586 std::vector<double> coo(0) ;
587 mu2->getCoordinatesOfNode(workv[0],coo);
589 for(int idim=0; idim<_spaceDim; idim++)
592 const mcIdType *workc=tmpA+tmpAI[id];
593 mcIdType nbCells=tmpAI[id+1]-tmpAI[id];
594 assert( nbCells>0);//To make sure our face is not located on an isolated node
596 Face fi( nbNodes, nbCells, 1.0, p, 1., 0., 0. ) ;
597 for(int node_id=0; node_id<nbNodes;node_id++)//This loop could be deleted since nbNodes=1. Trying to merge with setMesh
598 fi.addNodeId(node_id,workv[node_id]) ;//global node number
600 fi.addCellId(0,workc[0]) ;
601 for(int cell_id=1; cell_id<nbCells;cell_id++)
604 if (workc[cell_id]!=workc[cell_id-1])//For some meshes (bad ones) the same cell can appear several times
606 fi.addCellId(cell_idx+1,workc[cell_id]) ;
611 _boundaryFaceIds.push_back(id);
612 _faces.get()[id] = fi ;
615 int correctNbNodes=0;
616 for( int id=0;id<_numberOfNodes;id++ )
618 const mcIdType *workc=tmpN+tmpNI[id];
619 mcIdType nbCells=tmpNI[id+1]-tmpNI[id];
621 if( nbCells>0)//To make sure this is not an isolated node
624 std::vector<double> coo(0) ;
625 mu->getCoordinatesOfNode(id,coo);
627 for(int idim=0; idim<_spaceDim; idim++)
630 const mcIdType *workf=tmpC+tmpCI[id];
631 mcIdType nbFaces=tmpCI[id+1]-tmpCI[id];
634 const mcIdType *workn=tmpN+tmpNI[id];
635 mcIdType nbNeighbourNodes=tmpNI[id+1]-tmpNI[id];
637 Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
638 for( int el=0;el<nbCells;el++ )
639 vi.addCellId(el,workc[el]) ;
640 for( int el=0;el<nbNeighbourNodes;el++ )
641 vi.addNeighbourNodeId(el,workn[el]) ;//global node number
642 for( int el=0;el<nbFaces;el++ )
643 vi.addFaceId(el,workf[el],_faces.get()[workf[el]].isBorder()) ;
645 _boundaryNodeIds.push_back(id);
646 _nodes.get()[id] = vi ;
649 if( _numberOfNodes!=correctNbNodes)//To do : reduce the size of pointer _nodes
651 cout<<"Found isolated nodes : correctNbNodes= "<<correctNbNodes<<", _numberOfNodes= "<<_numberOfNodes<<endl;
652 _numberOfNodes = correctNbNodes;
653 //memcpy(_nodes,mesh.getNodes(),correctNbNodes*sizeof(*mesh.getNodes())) ;
656 else if(_meshDim==2 || _meshDim==3)
658 DataArrayDouble *barySeg = mu2->computeIsoBarycenterOfNodesPerCell();//computeCellCenterOfMass() ;//Used as face center
659 const double *coorBarySeg=barySeg->getConstPointer();
661 MEDCouplingFieldDouble* fieldl=mu2->getMeasureField(true);
662 DataArrayDouble *longueur = fieldl->getArray();
663 const double *lon=longueur->getConstPointer();//The lenght/area of each face
665 MEDCouplingFieldDouble* fieldn;//The normal to each face
666 DataArrayDouble *normal;
667 const double *tmpNormal;
669 if(_spaceDim==_meshDim)
670 fieldn = mu2->buildOrthogonalField();//Compute the normal to each cell interface
672 fieldn = mu->buildOrthogonalField();//compute the 3D normal vector to the 2D cell
674 normal = fieldn->getArray();
675 tmpNormal = normal->getConstPointer();
677 /*Building mesh cells */
678 for(int id(0), k(0); id<_numberOfCells; id++, k+=_spaceDim)
680 const mcIdType *work=tmp+tmpI[id];
681 mcIdType nbFaces=tmpI[id+1]-tmpI[id];
683 mcIdType nbVertices=mu->getNumberOfNodesInCell(id) ;
685 vector<double> coorBaryXyz(3,0);
686 for (int d=0; d<_spaceDim; d++)
687 coorBaryXyz[d] = coorBary[k+d];
689 Point p(coorBaryXyz[0],coorBaryXyz[1],coorBaryXyz[2]) ;
690 Cell ci( nbVertices, nbFaces, surf[id], p ) ;
692 /* Filling cell nodes */
693 std::vector<mcIdType> nodeIdsOfCell ;
694 mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
695 for( int el=0;el<nbVertices;el++ )
696 ci.addNodeId(el,nodeIdsOfCell[el]) ;
698 /* Filling cell faces */
699 if(_spaceDim==_meshDim)//use the normal field generated by buildOrthogonalField()
700 for( int el=0;el<nbFaces;el++ )
702 mcIdType faceIndex=(abs(work[el])-1);//=work[el] since Fortran type numbering was used, and negative sign means anticlockwise numbering
703 vector<double> xyzn(3,0);//Outer normal to the cell
705 for (int d=0; d<_spaceDim; d++)
706 xyzn[d] = -tmpNormal[_spaceDim*faceIndex+d];
708 for (int d=0; d<_spaceDim; d++)
709 xyzn[d] = +tmpNormal[_spaceDim*faceIndex+d];
710 ci.addNormalVector(el,xyzn[0],xyzn[1],xyzn[2]) ;
711 ci.addFaceId(el,faceIndex) ;
713 else//build normals associated to the couple (cell id, face el)
714 {//Case _meshDim=1 should be moved upper since we are in the 2D/3D branch
715 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
716 {//work[0]= first face global number, work[1]= second face global number
717 mcIdType indexFace0=abs(work[0])-1;//=work[0] since Fortran type numbering was used, and negative sign means anticlockwise numbering
718 mcIdType indexFace1=abs(work[1])-1;//=work[1] since Fortran type numbering was used, and negative sign means anticlockwise numbering
719 mcIdType idNodeA=(tmpNE+tmpNEI[indexFace0]+1)[0];//global number of the first face node work[0]=(abs(work[0])-1)
720 mcIdType idNodeB=(tmpNE+tmpNEI[indexFace1]+1)[0];//global number of the second face node work[1]=(abs(work[1])-1)
722 for(int i=0;i<_spaceDim;i++)
723 vecAB[i]=coo->getIJ(idNodeB,i) - coo->getIJ(idNodeA,i);
725 ci.addNormalVector(0,-vecAB[0],-vecAB[1],-vecAB[2]) ;
726 ci.addNormalVector(1,vecAB[0],vecAB[1],vecAB[2]) ;
727 ci.addFaceId(0,indexFace0) ;
728 ci.addFaceId(1,indexFace1) ;
730 else//_meshDim==2, number of faces around the cell id is variable, each face is composed of two nodes
733 for (int d=0; d<_spaceDim; d++)
734 xyzn[d] = tmpNormal[_spaceDim*id+d];
735 for( int el=0;el<nbFaces;el++ )
737 int faceIndex=(abs(work[el])-1);//=work[el] since Fortran type numbering was used, and negative sign means anticlockwise numbering
738 const mcIdType *workv=tmpNE+tmpNEI[faceIndex]+1;
739 mcIdType nbNodes= tmpNEI[faceIndex+1]-tmpNEI[faceIndex]-1;
740 if(nbNodes!=2)//We want to compute the normal to a straight line, not a curved interface composed of more thant 2 points
742 cout<<"Mesh name "<< mu->getName()<< " space dim= "<< _spaceDim <<" mesh dim= "<< _meshDim <<endl;
743 cout<<"For cell id "<<id<<" and local face number "<<el<<", the number of nodes is "<< nbNodes<< ", total number of faces is "<< nbFaces <<endl;
744 throw CdmathException("Mesh::setMesh number of nodes around a face should be 2");
747 mcIdType idNodeA=workv[0];
748 mcIdType idNodeB=workv[1];
749 vector<double> nodeA(_spaceDim), nodeB(_spaceDim), nodeP(_spaceDim);
750 for(int i=0;i<_spaceDim;i++)
752 nodeA[i]=coo->getIJ(idNodeA,i);
753 nodeB[i]=coo->getIJ(idNodeB,i);
754 nodeP[i]=coorBary[_spaceDim*id+i];
756 //Let P be the barycenter of the cell id
757 Vector vecAB(3), vecPA(3);
758 for(int i=0;i<_spaceDim;i++)
760 vecAB[i]=coo->getIJ(idNodeB,i) - coo->getIJ(idNodeA,i);
761 vecPA[i]=coo->getIJ(idNodeA,i) - coorBary[_spaceDim*id+i];
764 Vector normale = xyzn % vecAB;//Normal to the edge
765 normale/=normale.norm();
768 ci.addNormalVector(el,normale[0],normale[1],normale[2]) ;
770 ci.addNormalVector(el,-normale[0],-normale[1],-normale[2]) ;
771 ci.addFaceId(el,faceIndex) ;
775 _cells.get()[id] = ci ;
778 if(_spaceDim!=_meshDim)
780 /* Since spaceDim!=meshDim, don't build normal to faces */
786 /*Building mesh faces */
787 for(int id(0), k(0); id<_numberOfFaces; id++, k+=_spaceDim)
789 vector<double> coorBarySegXyz(3);
790 for (int d=0; d<_spaceDim; d++)
791 coorBarySegXyz[d] = coorBarySeg[k+d];
792 Point p(coorBarySegXyz[0],coorBarySegXyz[1],coorBarySegXyz[2]) ;
793 const mcIdType *workc=tmpA+tmpAI[id];
794 mcIdType nbCells=tmpAI[id+1]-tmpAI[id];
796 if (nbCells>2 && _spaceDim==_meshDim)
798 cout<<"Warning : nbCells>2, numberOfFaces="<<_numberOfFaces<<endl;
799 cout<<"nbCells= "<<nbCells<<", _spaceDim="<<_spaceDim<<", _meshDim="<<_meshDim<<endl;
800 for(int icell=0; icell<nbCells; icell++)
801 cout<<workc[icell]<<", ";
803 throw CdmathException("Wrong mesh : nbCells>2 and spaceDim==meshDim");
806 _boundaryFaceIds.push_back(id);
808 const mcIdType *workv=tmpNE+tmpNEI[id]+1;
809 mcIdType nbNodes= tmpNEI[id+1]-tmpNEI[id]-1;
812 if(_spaceDim==_meshDim)//Euclidean flat mesh geometry
814 fi=Face( nbNodes, nbCells, lon[id], p, tmpNormal[k], tmpNormal[k+1], 0.0) ;
816 fi=Face( nbNodes, nbCells, lon[id], p, tmpNormal[k], tmpNormal[k+1], tmpNormal[k+2]) ;
817 else//Curved mesh geometry
818 fi=Face( nbNodes, nbCells, lon[id], p, 0.0, 0.0, 0.0) ;//Since spaceDim!=meshDim, normal to face is not defined
820 for(int node_id=0; node_id<nbNodes;node_id++)
821 fi.addNodeId(node_id,workv[node_id]) ;
823 fi.addCellId(0,workc[0]) ;
824 for(int cell_id=1; cell_id<nbCells;cell_id++)
827 if (workc[cell_id]!=workc[cell_id-1])//For some meshes (bad ones) the same cell can appear several times
829 fi.addCellId(cell_idx+1,workc[cell_id]) ;
834 _faces.get()[id] = fi ;
837 /*Building mesh nodes, should be done after building mesh faces in order to detect boundary nodes*/
838 int correctNbNodes=0;
839 for(int id(0), k(0); id<_numberOfNodes; id++, k+=_spaceDim)
841 const mcIdType *workc=tmpN+tmpNI[id];
842 mcIdType nbCells=tmpNI[id+1]-tmpNI[id];
844 if( nbCells>0)//To make sure this is not an isolated node
847 vector<double> coorP(3);
848 for (int d=0; d<_spaceDim; d++)
849 coorP[d] = cood[k+d];
850 Point p(coorP[0],coorP[1],coorP[2]) ;
852 const mcIdType *workf=tmpC+tmpCI[id];
853 mcIdType nbFaces=tmpCI[id+1]-tmpCI[id];
854 const mcIdType *workn;
855 mcIdType nbNeighbourNodes;
858 workn=tmpA+tmpAI[id];
859 nbNeighbourNodes=tmpAI[id+1]-tmpAI[id];
861 else if (_meshDim == 2)
863 workn=tmpC+tmpCI[id];
864 nbNeighbourNodes=tmpCI[id+1]-tmpCI[id];
868 workn=tmpN2+tmpNI2[id];
869 nbNeighbourNodes=tmpNI2[id+1]-tmpNI2[id];
871 Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
873 for( int el=0;el<nbCells;el++ )
874 vi.addCellId(el,workc[el]) ;
875 for( int el=0;el<nbNeighbourNodes;el++ )
876 vi.addNeighbourNodeId(el,workn[el]) ;
877 //Detection of border nodes
878 for( int el=0;el<nbFaces;el++ )
879 vi.addFaceId(el,workf[el],_faces.get()[workf[el]].isBorder()) ;
881 _boundaryNodeIds.push_back(id);
882 _nodes.get()[id] = vi ;
885 if( _numberOfNodes!=correctNbNodes)//To do : reduce the size of pointer _nodes
887 cout<<"Found isolated nodes : correctNbNodes= "<<correctNbNodes<<", _numberOfNodes= "<<_numberOfNodes<<endl;
888 _numberOfNodes = correctNbNodes;
891 if(_spaceDim==_meshDim)
897 throw CdmathException("Mesh::setMesh space dimension should be 1, 2 or 3");
899 //Set boundary groups
900 _faceGroupNames.push_back("Boundary");
901 _nodeGroupNames.push_back("Boundary");
902 _faceGroupsIds.push_back(_boundaryFaceIds);
903 _nodeGroupsIds.push_back(_boundaryNodeIds);
905 { //Set face boundary group
906 _boundaryMesh = mu->computeSkin();
907 _faceGroups.push_back(_boundaryMesh);
910 _faceGroups.push_back(NULL);
911 _nodeGroups.push_back(NULL);
928 revNodeI2->decrRef();
932 revDescI2->decrRef();
940 Mesh::setGroupAtFaceByCoords(double x, double y, double z, double eps, std::string groupName, bool isBoundaryGroup)
942 std::vector< int > faceIds(0);
946 /* Construction of the face group */
948 for(int i=0; i<_boundaryFaceIds.size(); i++)
950 Fi=_faces.get()[_boundaryFaceIds[i]];
954 if (abs(FX-x)<eps && abs(FY-y)<eps && abs(FZ-z)<eps)
956 faceIds.insert(faceIds.end(),_boundaryFaceIds[i]);
957 _faces.get()[_boundaryFaceIds[i]].setGroupName(groupName);
961 for (int iface=0;iface<_numberOfFaces;iface++)
963 Fi=_faces.get()[iface];
967 if (abs(FX-x)<eps && abs(FY-y)<eps && abs(FZ-z)<eps)
969 faceIds.insert(faceIds.end(),iface);
970 _faces.get()[iface].setGroupName(groupName);
974 if (faceIds.size()>0)
976 std::vector< std::string >::iterator it = std::find(_faceGroupNames.begin(), _faceGroupNames.end(), groupName);
977 if(it == _faceGroupNames.end())//No group named groupName
979 _faceGroupNames.insert(_faceGroupNames.end(),groupName);
980 _faceGroupsIds.insert( _faceGroupsIds.end(),faceIds);
981 _faceGroups.insert( _faceGroups.end(), NULL);//No mesh created. Create one ?
985 std::vector< int > faceGroupIds = _faceGroupsIds[it-_faceGroupNames.begin()];
986 faceGroupIds.insert( faceGroupIds.end(), faceIds.begin(), faceIds.end());
987 /* Detect and erase duplicates face ids */
988 sort( faceGroupIds.begin(), faceGroupIds.end() );
989 faceGroupIds.erase( unique( faceGroupIds.begin(), faceGroupIds.end() ), faceGroupIds.end() );
990 _faceGroupsIds[it-_faceGroupNames.begin()] = faceGroupIds;
996 Mesh::setGroupAtNodeByCoords(double x, double y, double z, double eps, std::string groupName, bool isBoundaryGroup)
998 std::vector< int > nodeIds(0);
1002 /* Construction of the node group */
1004 for(int i=0; i<_boundaryNodeIds.size(); i++)
1006 Ni=_nodes.get()[_boundaryNodeIds[i]];
1010 if (abs(NX-x)<eps && abs(NY-y)<eps && abs(NZ-z)<eps)
1012 nodeIds.insert(nodeIds.end(),_boundaryNodeIds[i]);
1013 _nodes.get()[_boundaryNodeIds[i]].setGroupName(groupName);
1017 for (int inode=0;inode<_numberOfNodes;inode++)
1019 NX=_nodes.get()[inode].x();
1020 NY=_nodes.get()[inode].y();
1021 NZ=_nodes.get()[inode].z();
1022 if (abs(NX-x)<eps && abs(NY-y)<eps && abs(NZ-z)<eps)
1024 nodeIds.insert(nodeIds.end(),inode);
1025 _nodes.get()[inode].setGroupName(groupName);
1029 if (nodeIds.size()>0)
1031 std::vector< std::string >::iterator it = std::find(_nodeGroupNames.begin(), _nodeGroupNames.end(), groupName);
1032 if(it == _nodeGroupNames.end())//No group named groupName
1034 _nodeGroupNames.insert(_nodeGroupNames.end(),groupName);
1035 _nodeGroupsIds.insert( _nodeGroupsIds.end(),nodeIds);
1036 _nodeGroups.insert( _nodeGroups.end(), NULL);//No mesh created. Create one ?
1040 std::vector< int > nodeGroupIds = _nodeGroupsIds[it-_nodeGroupNames.begin()];
1041 nodeGroupIds.insert( nodeGroupIds.end(), nodeIds.begin(), nodeIds.end());
1042 /* Detect and erase duplicates node ids */
1043 sort( nodeGroupIds.begin(), nodeGroupIds.end() );
1044 nodeGroupIds.erase( unique( nodeGroupIds.begin(), nodeGroupIds.end() ), nodeGroupIds.end() );
1045 _nodeGroupsIds[it-_nodeGroupNames.begin()] = nodeGroupIds;
1051 Mesh::setGroupAtPlan(double value, int direction, double eps, std::string groupName, bool isBoundaryGroup)
1053 std::vector< int > faceIds(0), nodeIds(0);
1056 /* Construction of the face group */
1058 for(int i=0; i<_boundaryFaceIds.size(); i++)
1060 cord=_faces.get()[_boundaryFaceIds[i]].getBarryCenter()[direction];
1061 if (abs(cord-value)<eps)
1063 faceIds.insert(faceIds.end(),_boundaryFaceIds[i]);
1064 _faces.get()[_boundaryFaceIds[i]].setGroupName(groupName);
1068 for (int iface=0;iface<_numberOfFaces;iface++)
1070 cord=_faces.get()[iface].getBarryCenter()[direction];
1071 if (abs(cord-value)<eps)
1073 faceIds.insert(faceIds.end(),iface);
1074 _faces.get()[iface].setGroupName(groupName);
1078 /* Construction of the node group */
1080 for(int i=0; i<_boundaryNodeIds.size(); i++)
1082 cord=_nodes.get()[_boundaryNodeIds[i]].getPoint()[direction];
1083 if (abs(cord-value)<eps)
1085 nodeIds.insert(nodeIds.end(),_boundaryNodeIds[i]);
1086 _nodes.get()[_boundaryNodeIds[i]].setGroupName(groupName);
1090 for (int inode=0;inode<_numberOfNodes;inode++)
1092 cord=_nodes.get()[inode].getPoint()[direction];
1093 if (abs(cord-value)<eps)
1095 nodeIds.insert(nodeIds.end(),inode);
1096 _nodes.get()[inode].setGroupName(groupName);
1100 if (faceIds.size()>0)
1102 std::vector< std::string >::iterator it = std::find(_faceGroupNames.begin(), _faceGroupNames.end(), groupName);
1103 if(it == _faceGroupNames.end())//No group named groupName
1105 _faceGroupNames.insert(_faceGroupNames.end(),groupName);
1106 _faceGroupsIds.insert( _faceGroupsIds.end(),faceIds);
1107 _faceGroups.insert( _faceGroups.end(), NULL);//No mesh created. Create one ?
1110 {cout<<"_faceGroupNames.size()="<<_faceGroupNames.size()<<", _faceGroupsIds.size()="<<_faceGroupsIds.size()<<endl;
1111 std::vector< int > faceGroupIds = _faceGroupsIds[it-_faceGroupNames.begin()];
1112 faceGroupIds.insert( faceGroupIds.end(), faceIds.begin(), faceIds.end());
1113 /* Detect and erase duplicates face ids */
1114 sort( faceGroupIds.begin(), faceGroupIds.end() );
1115 faceGroupIds.erase( unique( faceGroupIds.begin(), faceGroupIds.end() ), faceGroupIds.end() );
1116 _faceGroupsIds[it-_faceGroupNames.begin()] = faceGroupIds;
1119 if (nodeIds.size()>0)
1121 std::vector< std::string >::iterator it = std::find(_nodeGroupNames.begin(), _nodeGroupNames.end(), groupName);
1122 if(it == _nodeGroupNames.end())//No group named groupName
1124 _nodeGroupNames.insert(_nodeGroupNames.end(),groupName);
1125 _nodeGroupsIds.insert( _nodeGroupsIds.end(),nodeIds);
1126 _nodeGroups.insert( _nodeGroups.end(), NULL);//No mesh created. Create one ?
1130 std::vector< int > nodeGroupIds = _nodeGroupsIds[it-_nodeGroupNames.begin()];
1131 nodeGroupIds.insert( nodeGroupIds.end(), nodeIds.begin(), nodeIds.end());
1132 /* Detect and erase duplicates node ids */
1133 sort( nodeGroupIds.begin(), nodeGroupIds.end() );
1134 nodeGroupIds.erase( unique( nodeGroupIds.begin(), nodeGroupIds.end() ), nodeGroupIds.end() );
1135 _nodeGroupsIds[it-_nodeGroupNames.begin()] = nodeGroupIds;
1141 Mesh::setBoundaryNodesFromFaces()
1143 for (int iface=0;iface<_boundaryFaceIds.size();iface++)
1145 std::vector< int > nodesID= _faces.get()[_boundaryFaceIds[iface]].getNodesId();
1146 int nbNodes = _faces.get()[_boundaryFaceIds[iface]].getNumberOfNodes();
1147 for(int inode=0 ; inode<nbNodes ; inode++)
1149 std::vector<int>::const_iterator it = std::find(_boundaryNodeIds.begin(),_boundaryNodeIds.end(),nodesID[inode]);
1150 if( it != _boundaryNodeIds.end() )
1151 _boundaryNodeIds.push_back(nodesID[inode]);
1157 Mesh::getIndexFacePeriodic( void ) const
1159 return _indexFacePeriodicMap;
1163 Mesh::setPeriodicFaces(bool check_groups, bool use_central_inversion)
1165 if(_indexFacePeriodicSet)
1168 for (int indexFace=0;indexFace<_boundaryFaceIds.size() ; indexFace++)
1170 Face my_face=_faces.get()[_boundaryFaceIds[indexFace]];
1174 for (int iface=0;iface<_boundaryFaceIds.size() ; iface++)
1175 if(iface!=indexFace)
1177 iface_perio=_boundaryFaceIds[iface];
1181 else if(_meshDim==2)
1183 double x=my_face.x();
1184 double y=my_face.y();
1186 for (int iface=0;iface<_boundaryFaceIds.size() ; iface++)
1188 Face face_i=_faces.get()[_boundaryFaceIds[iface]];
1189 double xi=face_i.x();
1190 double yi=face_i.y();
1191 if ( (abs(y-yi)<_epsilon || abs(x-xi)<_epsilon )// Case of a square geometry
1192 && ( !check_groups || my_face.getGroupName()!=face_i.getGroupName()) //In case groups need to be checked
1193 && ( !use_central_inversion || abs(y+yi) + abs(x+xi)<_epsilon ) // Case of a central inversion
1194 && fabs(my_face.getMeasure() - face_i.getMeasure())<_epsilon
1195 && fabs(my_face.getXN() + face_i.getXN())<_epsilon
1196 && fabs(my_face.getYN() + face_i.getYN())<_epsilon
1197 && fabs(my_face.getZN() + face_i.getZN())<_epsilon )
1199 iface_perio=_boundaryFaceIds[iface];
1204 else if(_meshDim==3)
1206 double x=my_face.x();
1207 double y=my_face.y();
1208 double z=my_face.z();
1210 for (int iface=0;iface<_boundaryFaceIds.size() ; iface++)
1212 Face face_i=_faces.get()[_boundaryFaceIds[iface]];
1213 double xi=face_i.x();
1214 double yi=face_i.y();
1215 double zi=face_i.z();
1216 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
1217 && ( !check_groups || my_face.getGroupName()!=face_i.getGroupName()) //In case groups need to be checked
1218 && ( !use_central_inversion || abs(y+yi) + abs(x+xi) + abs(z+zi)<_epsilon )// Case of a central inversion
1219 && fabs(my_face.getMeasure() - face_i.getMeasure())<_epsilon
1220 && fabs(my_face.getXN() + face_i.getXN())<_epsilon
1221 && fabs(my_face.getYN() + face_i.getYN())<_epsilon
1222 && fabs(my_face.getZN() + face_i.getZN())<_epsilon )
1224 iface_perio=_boundaryFaceIds[iface];
1230 throw CdmathException("Mesh::setPeriodicFaces: Mesh dimension should be 1, 2 or 3");
1232 if (iface_perio==-1)
1233 throw CdmathException("Mesh::setPeriodicFaces: periodic face not found, iface_perio==-1 " );
1235 _indexFacePeriodicMap[_boundaryFaceIds[indexFace]]=iface_perio;
1237 _indexFacePeriodicSet=true;
1241 Mesh::getIndexFacePeriodic(int indexFace, bool check_groups, bool use_central_inversion)
1243 if (!_faces.get()[indexFace].isBorder())
1245 cout<<"Pb with indexFace= "<<indexFace<<endl;
1246 throw CdmathException("Mesh::getIndexFacePeriodic: not a border face" );
1249 if(!_indexFacePeriodicSet)
1250 setPeriodicFaces(check_groups, use_central_inversion);
1252 std::map<int,int>::const_iterator it = _indexFacePeriodicMap.find(indexFace);
1253 if( it != _indexFacePeriodicMap.end() )
1257 cout<<"Pb with indexFace= "<<indexFace<<endl;
1258 throw CdmathException("Mesh::getIndexFacePeriodic: not a periodic face" );
1263 Mesh::isBorderNode(int nodeid) const
1265 return _nodes.get()[nodeid].isBorder();
1269 Mesh::isBorderFace(int faceid) const
1271 return _faces.get()[faceid].isBorder();
1275 Mesh::getBoundaryFaceIds() const
1277 return _boundaryFaceIds;
1281 Mesh::getBoundaryNodeIds() const
1283 return _boundaryNodeIds;
1287 Mesh::isTriangular() const
1289 return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_TRI3;
1292 Mesh::isTetrahedral() const
1294 return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_TETRA4;
1297 Mesh::isQuadrangular() const
1299 return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_QUAD4;
1302 Mesh::isHexahedral() const
1304 return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_HEXA8;
1307 Mesh::isStructured() const
1309 return _isStructured;
1312 std::vector< INTERP_KERNEL::NormalizedCellType >
1313 Mesh::getElementTypes() const
1318 std::vector< string >
1319 Mesh::getElementTypesNames() const
1321 std::vector< string > result(0);
1322 for(int i=0; i< _eltsTypes.size(); i++)
1324 if( _eltsTypes[i]==INTERP_KERNEL::NORM_POINT1)
1325 result.push_back("Points");
1326 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_SEG2)
1327 result.push_back("Segments");
1328 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_TRI3)
1329 result.push_back("Triangles");
1330 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_QUAD4)
1331 result.push_back("Quadrangles");
1332 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_POLYGON)
1333 result.push_back("Polygons");
1334 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_TETRA4)
1335 result.push_back("Tetrahedra");
1336 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_PYRA5)
1337 result.push_back("Pyramids");
1338 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_PENTA6)
1339 result.push_back("Pentahedra");
1340 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_HEXA8)
1341 result.push_back("Hexahedra");
1342 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_POLYHED)
1343 result.push_back("Polyhedrons");
1346 cout<< "Mesh " + _name + " contains an element of type " <<endl;
1347 cout<< _eltsTypes[i]<<endl;
1348 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");
1355 Mesh::setFaceGroups( const MEDFileUMesh* medmesh, MEDCouplingUMesh* mu)
1360 /* Searching for face groups */
1361 vector<string> faceGroups=medmesh->getGroupsNames() ;
1363 for (unsigned int i=0;i<faceGroups.size();i++ )
1365 string groupName=faceGroups[i];
1366 vector<mcIdType> nonEmptyGrp(medmesh->getGrpNonEmptyLevels(groupName));
1367 //We check if the group has a relative dimension equal to -1
1368 //before call to the function getGroup(-1,groupName.c_str())
1369 vector<mcIdType>::iterator it = find(nonEmptyGrp.begin(), nonEmptyGrp.end(), -1);
1370 if (it != nonEmptyGrp.end())
1372 MEDCouplingUMesh *m=medmesh->getGroup(-1,groupName.c_str());
1374 nbCellsSubMesh=m->getNumberOfCells();
1376 _faceGroups.insert(_faceGroups.end(),m);//Vector of group meshes
1377 _faceGroupNames.insert(_faceGroupNames.end(),groupName);//Vector of group names
1378 _faceGroupsIds.insert(_faceGroupsIds.end(),std::vector<int>(nbCellsSubMesh));//Vector of group face Ids. The filling of the face ids takes place below.
1380 DataArrayDouble *baryCell = m->computeIsoBarycenterOfNodesPerCell();//computeCellCenterOfMass() ;
1381 const double *coorBary=baryCell->getConstPointer();
1382 // Face identification
1383 for (int ic(0), k(0); ic<nbCellsSubMesh; ic++, k+=_spaceDim)
1385 vector<double> coorBaryXyz(3,0);
1386 for (int d=0; d<_spaceDim; d++)
1387 coorBaryXyz[d] = coorBary[k+d];
1388 Point p1(coorBaryXyz[0],coorBaryXyz[1],coorBaryXyz[2]) ;
1391 for (int iface=0;iface<_numberOfFaces;iface++ )
1393 Point p2=_faces.get()[iface].getBarryCenter();
1394 if(p1.distance(p2)<_epsilon)
1396 _faces.get()[iface].setGroupName(groupName);
1397 _faceGroupsIds[_faceGroupsIds.size()-1][ic]=iface;
1403 throw CdmathException("No face found for group " + groupName );
1405 baryCell->decrRef();
1411 Mesh::setNodeGroups( const MEDFileMesh* medmesh, MEDCouplingUMesh* mu)
1416 /* Searching for node groups */
1417 vector<string> nodeGroups=medmesh->getGroupsOnSpecifiedLev(1) ;
1419 for (unsigned int i=0;i<nodeGroups.size();i++ )
1421 string groupName=nodeGroups[i];
1422 DataArrayIdType * nodeGroup=medmesh->getNodeGroupArr( groupName );
1423 const mcIdType *nodeids=nodeGroup->getConstPointer();
1427 _nodeGroups.insert(_nodeGroups.end(),nodeGroup);
1428 _nodeGroupNames.insert(_nodeGroupNames.end(),groupName);
1430 nbNodesSubMesh=nodeGroup->getNumberOfTuples();//nodeGroup->getNbOfElems();
1432 DataArrayDouble *coo = mu->getCoords() ;
1433 const double *cood=coo->getConstPointer();
1435 _nodeGroupsIds.insert(_nodeGroupsIds.end(),std::vector<int>(nbNodesSubMesh));//Vector of boundary faces
1436 /* Node identification */
1437 for (int ic(0); ic<nbNodesSubMesh; ic++)
1439 vector<double> coorP(3,0);
1440 for (int d=0; d<_spaceDim; d++)
1441 coorP[d] = cood[nodeids[ic]*_spaceDim+d];
1442 Point p1(coorP[0],coorP[1],coorP[2]) ;
1445 for (int inode=0;inode<_numberOfNodes;inode++ )
1447 Point p2=_nodes.get()[inode].getPoint();
1448 if(p1.distance(p2)<_epsilon)
1450 _nodes.get()[inode].setGroupName(groupName);
1451 _nodeGroupsIds[_nodeGroupsIds.size()-1][ic]=inode;
1457 throw CdmathException("No node found for group " + groupName );
1463 //----------------------------------------------------------------------
1465 Mesh::getSpaceDimension( void ) const
1466 //----------------------------------------------------------------------
1471 //----------------------------------------------------------------------
1473 Mesh::getMeshDimension( void ) const
1474 //----------------------------------------------------------------------
1479 std::vector<mcIdType>
1480 Mesh::getCellGridStructure() const
1483 throw CdmathException("std::vector<int> Mesh::getCellGridStructure() : nx, ny and nz are defined only for structured meshes !");
1488 //----------------------------------------------------------------------
1490 Mesh::getNx( void ) const
1491 //----------------------------------------------------------------------
1494 throw CdmathException("int Mesh::getNx( void ) : Nx is defined only for structured meshes !");
1499 //----------------------------------------------------------------------
1501 Mesh::getNy( void ) const
1502 //----------------------------------------------------------------------
1505 throw CdmathException("int Mesh::getNy( void ) : Ny is defined only for structured meshes !");
1507 throw CdmathException("int Mesh::getNy( void ) : Ny is not defined for mesh dimension < 2!");
1512 //----------------------------------------------------------------------
1514 Mesh::getNz( void ) const
1515 //----------------------------------------------------------------------
1518 throw CdmathException("int Mesh::getNz( void ) : Nz is defined only for structured meshes !");
1520 throw CdmathException("int Mesh::getNz( void ) : Nz is not defined for mesh dimension < 3!");
1525 //----------------------------------------------------------------------
1527 Mesh::getXMin( void ) const
1528 //----------------------------------------------------------------------
1530 double Box0[2*_meshDim];
1531 _mesh->getBoundingBox(Box0);
1536 //----------------------------------------------------------------------
1538 Mesh::getXMax( void ) const
1539 //----------------------------------------------------------------------
1541 double Box0[2*_meshDim];
1542 _mesh->getBoundingBox(Box0);
1547 //----------------------------------------------------------------------
1549 Mesh::getYMin( void ) const
1550 //----------------------------------------------------------------------
1553 throw CdmathException("Mesh::getYMin : dimension should be >=2");
1555 double Box0[2*_meshDim];
1556 _mesh->getBoundingBox(Box0);
1561 //----------------------------------------------------------------------
1563 Mesh::getYMax( void ) const
1564 //----------------------------------------------------------------------
1567 throw CdmathException("Mesh::getYMax : dimension should be >=2");
1569 double Box0[2*_meshDim];
1570 _mesh->getBoundingBox(Box0);
1575 //----------------------------------------------------------------------
1577 Mesh::getZMin( void ) const
1578 //----------------------------------------------------------------------
1581 throw CdmathException("Mesh::getZMin : dimension should be 3");
1583 double Box0[2*_meshDim];
1584 _mesh->getBoundingBox(Box0);
1589 //----------------------------------------------------------------------
1591 Mesh::getZMax( void ) const
1592 //----------------------------------------------------------------------
1595 throw CdmathException("Mesh::getZMax : dimension should be 3");
1597 double Box0[2*_meshDim];
1598 _mesh->getBoundingBox(Box0);
1603 //----------------------------------------------------------------------
1604 MCAuto<MEDCouplingMesh>
1605 Mesh::getMEDCouplingMesh( void ) const
1606 //----------------------------------------------------------------------
1611 //----------------------------------------------------------------------
1613 Mesh::getNumberOfNodes ( void ) const
1614 //----------------------------------------------------------------------
1616 return _numberOfNodes ;
1619 //----------------------------------------------------------------------
1621 Mesh::getNumberOfCells ( void ) const
1622 //----------------------------------------------------------------------
1624 return _numberOfCells ;
1627 //----------------------------------------------------------------------
1629 Mesh::getNumberOfFaces ( void ) const
1630 //----------------------------------------------------------------------
1632 return _numberOfFaces ;
1635 //----------------------------------------------------------------------
1637 Mesh::getNumberOfEdges ( void ) const
1638 //----------------------------------------------------------------------
1640 return _numberOfEdges ;
1643 //----------------------------------------------------------------------
1644 std::shared_ptr<Face>
1645 Mesh::getFaces ( void ) const
1646 //----------------------------------------------------------------------
1651 //----------------------------------------------------------------------
1652 std::shared_ptr<Cell>
1653 Mesh::getCells ( void ) const
1654 //----------------------------------------------------------------------
1659 //----------------------------------------------------------------------
1661 Mesh::getCell ( int i ) const
1662 //----------------------------------------------------------------------
1664 return _cells.get()[i] ;
1667 //----------------------------------------------------------------------
1669 Mesh::getFace ( int i ) const
1670 //----------------------------------------------------------------------
1672 return _faces.get()[i] ;
1675 //----------------------------------------------------------------------
1677 Mesh::getNode ( int i ) const
1678 //----------------------------------------------------------------------
1680 return _nodes.get()[i] ;
1683 //----------------------------------------------------------------------
1684 std::shared_ptr<Node>
1685 Mesh::getNodes ( void ) const
1686 //----------------------------------------------------------------------
1692 Mesh::getNameOfFaceGroups( void ) const
1694 return _faceGroupNames;
1697 vector<MEDCoupling::MEDCouplingUMesh *>
1698 Mesh::getFaceGroups( void ) const
1704 Mesh::getNameOfNodeGroups( void ) const
1706 return _nodeGroupNames;
1709 vector<MEDCoupling::DataArrayIdType *>
1710 Mesh::getNodeGroups( void ) const
1715 //----------------------------------------------------------------------
1717 Mesh::operator= ( const Mesh& mesh )
1718 //----------------------------------------------------------------------
1720 _spaceDim = mesh.getSpaceDimension() ;
1721 _meshDim = mesh.getMeshDimension() ;
1722 _name = mesh.getName();
1723 _epsilon=mesh.getComparisonEpsilon();
1724 _numberOfNodes = mesh.getNumberOfNodes();
1725 _numberOfFaces = mesh.getNumberOfFaces();
1726 _numberOfCells = mesh.getNumberOfCells();
1727 _numberOfEdges = mesh.getNumberOfEdges();
1729 _isStructured = mesh.isStructured();
1730 _meshNotDeleted = mesh.meshNotDeleted();
1733 _nxyz = mesh.getCellGridStructure() ;
1735 _faceGroupNames = mesh.getNameOfFaceGroups() ;
1736 _faceGroups = mesh.getFaceGroups() ;
1737 _nodeGroupNames = mesh.getNameOfNodeGroups() ;
1738 _nodeGroups = mesh.getNodeGroups() ;
1740 _nodes = mesh.getNodes() ;
1741 _faces = mesh.getFaces() ;
1742 _cells = mesh.getCells() ;
1744 _indexFacePeriodicSet= mesh.isIndexFacePeriodicSet();
1745 if(_indexFacePeriodicSet)
1746 _indexFacePeriodicMap=mesh.getIndexFacePeriodic();
1748 _boundaryFaceIds=mesh.getBoundaryFaceIds();
1749 _boundaryNodeIds=mesh.getBoundaryNodeIds();
1751 _eltsTypes=mesh.getElementTypes();
1752 _eltsTypesNames=mesh.getElementTypesNames();
1754 _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
1759 bool Mesh::isIndexFacePeriodicSet() const
1761 return _indexFacePeriodicSet;
1763 //----------------------------------------------------------------------
1765 Mesh::minRatioVolSurf() const
1767 double dx_min = 1e30;
1768 for(int i=0; i<_numberOfCells; i++)
1770 Cell Ci = getCell(i);
1774 for(int k=0; k< Ci.getNumberOfFaces(); k++)
1776 int indexFace=Ci.getFacesId()[k];
1777 Face Fk = getFace(indexFace);
1778 perimeter+=Fk.getMeasure();
1780 dx_min = min(dx_min,Ci.getMeasure()/perimeter);
1783 dx_min = min(dx_min,Ci.getMeasure());
1790 Mesh::getBoundaryMesh ( void ) const
1792 return Mesh(_boundaryMesh);
1796 Mesh::getBoundaryGroupMesh ( std::string groupName, int nth_occurence ) const
1798 //count occurences of groupName in known group name list
1799 int count_occurences = std::count (_faceGroupNames.begin(),_faceGroupNames.end(),groupName);
1801 if( count_occurences ==0 )//No group found
1803 cout<<"Mesh::getBoundaryGroupMesh Error : face group " << groupName << " does not exist"<<endl;
1804 cout<<"Known face group names are " ;
1805 for(int i=0; i<_faceGroupNames.size(); i++)
1806 cout<< _faceGroupNames[i]<<" ";
1808 throw CdmathException("Required face group does not exist");
1810 else if ( count_occurences <= nth_occurence)//Too many groups found
1812 cout<<"Mesh::getBoundaryGroupMesh Error : "<<count_occurences<<" groups have name " << groupName<<", but you asked fo occurencer number "<<nth_occurence<<"which is too large"<<endl;
1813 cout<<"Call function getBoundaryGroupMesh ( string groupName, int nth_group_match ) with nth_group_match between 0 and "<<count_occurences-1<<" to discriminate them "<<endl ;
1814 throw CdmathException("Several face groups have the same name but you asked for an occurence that does not exsit");
1816 else if( count_occurences >1 )//Wrning several groups found
1817 cout<<"Warning : "<<count_occurences<<" groups have name " << groupName<<". Searching occurence number "<<nth_occurence<<endl;
1819 //search occurence of group name in known group name list
1820 std::vector<std::string>::const_iterator it = _faceGroupNames.begin();
1821 for (int i = 0; i<nth_occurence+1; i++)
1822 it = std::find(it,_faceGroupNames.end(),groupName);
1824 return Mesh(_faceGroups[it - _faceGroupNames.begin()]);
1828 Mesh::getMaxNbNeighbours(EntityType type) const
1834 int nbNeib;//local number of neighbours
1835 for(int i=0; i<_numberOfCells; i++)
1837 Cell Ci = _cells.get()[i];
1838 //Careful with mesh with junctions : do not just take Ci.getNumberOfFaces()
1840 for(int j=0; j<Ci.getNumberOfFaces(); j++)
1841 nbNeib+=_faces.get()[Ci.getFacesId()[j]].getNumberOfCells()-1;//Without junction this would be +=1
1847 else if(type==NODES)
1849 for(int i=0; i<_numberOfNodes; i++)
1850 if(result < _nodes.get()[i].getNumberOfEdges())
1851 result=_nodes.get()[i].getNumberOfEdges();
1854 throw CdmathException("Mesh::getMaxNbNeighbours : entity type is not accepted. Should be CELLS or NODES");
1858 //----------------------------------------------------------------------
1860 Mesh::writeVTK ( const std::string fileName ) const
1861 //----------------------------------------------------------------------
1863 if( !_meshNotDeleted )
1864 throw CdmathException("Mesh::writeVTK : Cannot save mesh : no MEDCouplingMesh loaded (may be deleted)");
1866 string fname=fileName+".vtu";
1867 _mesh->writeVTK(fname.c_str()) ;
1870 //----------------------------------------------------------------------
1872 Mesh::writeMED ( const std::string fileName, bool fromScratch ) const
1873 //----------------------------------------------------------------------
1875 if( !_meshNotDeleted )
1876 throw CdmathException("Mesh::writeMED : Cannot save mesh : no MEDCouplingMesh loaded (may be deleted)");
1878 string fname=fileName+".med";
1879 if(_isStructured)//Check if we have a medcouplingimesh that can't be written directly
1881 const MEDCoupling::MEDCouplingIMesh* iMesh = dynamic_cast< const MEDCoupling::MEDCouplingIMesh* > ((const MEDCoupling::MEDCouplingMesh*) _mesh);
1882 if(iMesh)//medcouplingimesh : Use convertToCartesian in order to write mesh
1883 MEDCoupling::WriteMesh(fname.c_str(),iMesh->convertToCartesian(), fromScratch);
1884 else//medcouplingcmesh : save directly
1885 MEDCoupling::WriteMesh(fname.c_str(),_mesh, fromScratch);
1888 MEDCoupling::WriteMesh(fname.c_str(),_mesh, fromScratch);
1890 /* Try to save mesh groups */
1891 //MEDFileUMesh meshMEDFile;
1892 //meshMEDFile.setMeshAtLevel(0,mu);
1893 //for(int i=0; i< _faceGroups.size(); i++)
1894 //meshMEDFile.setMeshAtLevel(-1,_faceGroups[i]);
1896 //MEDCoupling::meshMEDFile.write(fname.c_str(),2) ;
1898 //MEDCoupling::meshMEDFile.write(fname.c_str(),1) ;
1902 Mesh::getFaceGroupIds(std::string groupName, bool isBoundaryGroup) const
1904 std::vector<std::string>::const_iterator it = std::find(_faceGroupNames.begin(),_faceGroupNames.end(),groupName);
1905 if( it == _faceGroupNames.end() )
1907 cout<<"Mesh::getGroupFaceIds Error : face group " << groupName << " does not exist"<<endl;
1908 throw CdmathException("Required face group does not exist");
1912 return _faceGroupsIds[it-_faceGroupNames.begin()];
1917 Mesh::getNodeGroupIds(std::string groupName, bool isBoundaryGroup) const
1919 std::vector<std::string>::const_iterator it = std::find(_nodeGroupNames.begin(),_nodeGroupNames.end(),groupName);
1920 if( it == _nodeGroupNames.end() )
1922 cout<<"Mesh::getGroupNodeIds Error : node group " << groupName << " does not exist"<<endl;
1923 throw CdmathException("Required node group does not exist");
1926 return _nodeGroupsIds[it-_nodeGroupNames.begin()];
1930 Mesh::setFaceGroupByIds(std::vector< int > faceIds, std::string groupName)
1932 for(int i=0; i< faceIds.size(); i++)
1933 getFace(faceIds[i]).setGroupName(groupName);
1935 _faceGroupsIds.insert( _faceGroupsIds.end(),faceIds);
1936 _faceGroupNames.insert(_faceGroupNames.end(), groupName);
1937 _faceGroups.insert( _faceGroups.end(), NULL);
1941 Mesh::setNodeGroupByIds(std::vector< int > nodeIds, std::string groupName)
1943 for(int i=0; i< nodeIds.size(); i++)
1944 getNode(nodeIds[i]).setGroupName(groupName);
1947 void Mesh::deleteMEDCouplingMesh()
1951 (_mesh.retn())->decrRef();
1952 _meshNotDeleted=false;
1955 throw CdmathException("Mesh::deleteMEDCouplingMesh() : mesh is not loaded");