4 * Created on: 22 janv. 2012
13 #include "MEDFileMesh.hxx"
14 #include "MEDLoader.hxx"
15 #include "MEDCouplingUMesh.hxx"
16 #include "MEDCouplingIMesh.hxx"
17 #include "MEDCouplingFieldDouble.hxx"
19 #include "CdmathException.hxx"
28 using namespace MEDCoupling;
31 //----------------------------------------------------------------------
33 //----------------------------------------------------------------------
54 _faceGroupNames.resize(0);
55 _faceGroups.resize(0);
56 _nodeGroupNames.resize(0);
57 _nodeGroups.resize(0);
58 _indexFacePeriodicSet=false;
62 //----------------------------------------------------------------------
64 //----------------------------------------------------------------------
69 //for(int i=0; i< _faceGroups.size(); i++)
70 // _faceGroups[i]->decrRef();
71 // _nodeGroups[i]->decrRef();
75 Mesh::getName( void ) const
80 Mesh::Mesh( const MEDCoupling::MEDCouplingIMesh* mesh )
82 _spaceDim=mesh->getSpaceDimension();
83 _meshDim=mesh->getMeshDimension();
85 _dxyz=mesh->getDXYZ();
86 _nxyz=mesh->getCellGridStructure();
87 double* Box0=new double[2*_spaceDim];
88 mesh->getBoundingBox(Box0);
89 _name=mesh->getName();
90 _indexFacePeriodicSet=false;
105 double *originPtr = new double[_spaceDim];
106 double *dxyzPtr = new double[_spaceDim];
107 int *nodeStrctPtr = new int[_spaceDim];
109 for(int i=0;i<_spaceDim;i++)
111 originPtr[i]=Box0[2*i];
112 nodeStrctPtr[i]=_nxyz[i]+1;
115 _mesh=MEDCouplingIMesh::New(_name,
118 nodeStrctPtr+_spaceDim,
125 delete [] nodeStrctPtr;
130 Mesh::Mesh( const MEDCoupling::MEDCouplingUMesh* mesh )
132 _spaceDim=mesh->getSpaceDimension();
133 _meshDim=mesh->getMeshDimension();
135 double* Box0=new double[2*_spaceDim];
136 mesh->getBoundingBox(Box0);
137 _name=mesh->getName();
138 _indexFacePeriodicSet=false;
153 _mesh=mesh->deepCopy();
158 //----------------------------------------------------------------------
159 Mesh::Mesh( const Mesh& mesh )
160 //----------------------------------------------------------------------
162 _spaceDim = mesh.getSpaceDimension() ;
163 _meshDim = mesh.getMeshDimension() ;
164 _name=mesh.getName();
165 _xMax=mesh.getXMax();
166 _yMin=mesh.getYMin();
167 _yMax=mesh.getYMax();
168 _zMin=mesh.getZMin();
169 _zMax=mesh.getZMax();
170 _isStructured=mesh.isStructured();
173 _nxyz = mesh.getCellGridStructure() ;
174 _dxyz=mesh.getDXYZ();
176 _numberOfNodes = mesh.getNumberOfNodes();
177 _numberOfFaces = mesh.getNumberOfFaces();
178 _numberOfCells = mesh.getNumberOfCells();
179 _numberOfEdges = mesh.getNumberOfEdges();
181 _faceGroupNames = mesh.getNameOfFaceGroups() ;
182 _faceGroups = mesh.getFaceGroups() ;
183 _nodeGroupNames = mesh.getNameOfNodeGroups() ;
184 _nodeGroups = mesh.getNodeGroups() ;
186 _nodes = new Node[_numberOfNodes] ;
187 _faces = new Face[_numberOfFaces] ;
188 _cells = new Cell[_numberOfCells] ;
190 //memcpy(_nodes,mesh.getNodes(),sizeof(*mesh.getNodes())) ;
191 //memcpy(_cells,mesh.getCells(),sizeof(*mesh.getCells())) ;
192 //memcpy(_faces,mesh.getFaces(),sizeof(*mesh.getFaces())) ;
194 for (int i=0;i<_numberOfNodes;i++)
195 _nodes[i]=mesh.getNode(i);
197 for (int i=0;i<_numberOfFaces;i++)
198 _faces[i]=mesh.getFace(i);
200 for (int i=0;i<_numberOfCells;i++)
201 _cells[i]=mesh.getCell(i);
203 _indexFacePeriodicSet= mesh.isIndexFacePeriodicSet();
204 if(_indexFacePeriodicSet)
205 _indexFacePeriodicMap=mesh.getIndexFacePeriodic();
207 _boundaryFaceIds=mesh.getBoundaryFaceIds();
208 _boundaryNodeIds=mesh.getBoundaryNodeIds();
210 _eltsTypes=mesh.getElementTypes();
211 _eltsTypesNames=mesh.getElementTypesNames();
213 MCAuto<MEDCouplingMesh> m1=mesh.getMEDCouplingMesh()->deepCopy();
217 //----------------------------------------------------------------------
218 Mesh::Mesh( const std::string filename, int meshLevel )
219 //----------------------------------------------------------------------
221 readMeshMed(filename, meshLevel);
224 //----------------------------------------------------------------------
226 Mesh::readMeshMed( const std::string filename, const int meshLevel)
227 //----------------------------------------------------------------------
229 MEDFileUMesh *m=MEDFileUMesh::New(filename.c_str());//reads the first mesh encountered in the file, otherwise call New (const char *fileName, const char *mName, int dt=-1, int it=-1)
230 _mesh=m->getMeshAtLevel(meshLevel);
231 _mesh->checkConsistencyLight();
232 _mesh->setName(_mesh->getName());
233 _meshDim=_mesh->getMeshDimension();
234 _spaceDim=_mesh->getSpaceDimension();
235 _name=_mesh->getName();
236 _indexFacePeriodicSet=false;
237 MEDCoupling::MEDCouplingIMesh* structuredMesh = dynamic_cast<MEDCoupling::MEDCouplingIMesh*> (_mesh.retn());
241 _dxyz=structuredMesh->getDXYZ();
242 _nxyz=structuredMesh->getCellGridStructure();
243 double* Box0=new double[2*_spaceDim];
244 structuredMesh->getBoundingBox(Box0);
262 MEDCouplingUMesh* mu = setMesh();
265 cout<<endl<< "Loaded file "<< filename<<endl;
266 cout<<"Mesh name= "<<m->getName()<<", mesh dim="<< _meshDim<< ", space dim="<< _spaceDim<< ", nb cells= "<<getNumberOfCells()<< ", nb nodes= "<<getNumberOfNodes()<<endl;
273 Mesh::setGroupAtFaceByCoords(double x, double y, double z, double eps, std::string groupName)
275 int nbFace=getNumberOfFaces();
277 for (int iface=0;iface<nbFace;iface++)
279 double FX=_faces[iface].x();
280 double FY=_faces[iface].y();
281 double FZ=_faces[iface].z();
282 if (abs(FX-x)<eps && abs(FY-y)<eps && abs(FZ-z)<eps)
284 _faces[iface].setGroupName(groupName);
285 std::vector< int > nodesID= _faces[iface].getNodesId();
286 int nbNodes = _faces[iface].getNumberOfNodes();
287 for(int inode=0 ; inode<nbNodes ; inode++)
288 _nodes[nodesID[inode]].setGroupName(groupName);
295 _faceGroupNames.insert(_faceGroupNames.begin(),groupName);
296 _nodeGroupNames.insert(_nodeGroupNames.begin(),groupName);
297 //To do : update _faceGroups and _nodeGroups
302 Mesh::setGroupAtPlan(double value, int direction, double eps, std::string groupName)
304 int nbFace=getNumberOfFaces();
306 for (int iface=0;iface<nbFace;iface++)
308 double cord=_faces[iface].getBarryCenter()[direction];
309 if (abs(cord-value)<eps)
311 _faces[iface].setGroupName(groupName);
312 std::vector< int > nodesID= _faces[iface].getNodesId();
313 int nbNodes = _faces[iface].getNumberOfNodes();
314 for(int inode=0 ; inode<nbNodes ; inode++)
316 _nodes[nodesID[inode]].setGroupName(groupName);
324 _faceGroupNames.insert(_faceGroupNames.begin(),groupName);
325 _nodeGroupNames.insert(_nodeGroupNames.begin(),groupName);
326 //To do : update _faceGroups, _nodeGroups
331 Mesh::setBoundaryNodesFromFaces()
333 for (int iface=0;iface<_boundaryFaceIds.size();iface++)
335 std::vector< int > nodesID= _faces[_boundaryFaceIds[iface]].getNodesId();
336 int nbNodes = _faces[_boundaryFaceIds[iface]].getNumberOfNodes();
337 for(int inode=0 ; inode<nbNodes ; inode++)
339 std::vector<int>::const_iterator it = std::find(_boundaryNodeIds.begin(),_boundaryNodeIds.end(),nodesID[inode]);
340 if( it != _boundaryNodeIds.end() )
341 _boundaryNodeIds.push_back(nodesID[inode]);
347 Mesh::getIndexFacePeriodic( void ) const
349 return _indexFacePeriodicMap;
353 Mesh::setPeriodicFaces(bool check_groups, bool use_central_inversion)
355 if(_indexFacePeriodicSet)
360 for (int indexFace=0;indexFace<_boundaryFaceIds.size() ; indexFace++)
362 Face my_face=_faces[_boundaryFaceIds[indexFace]];
366 for (int iface=0;iface<_boundaryFaceIds.size() ; iface++)
369 iface_perio=_boundaryFaceIds[iface];
375 double x=my_face.x();
376 double y=my_face.y();
378 for (int iface=0;iface<_boundaryFaceIds.size() ; iface++)
380 Face face_i=_faces[_boundaryFaceIds[iface]];
381 double xi=face_i.x();
382 double yi=face_i.y();
383 if ( (abs(y-yi)<eps || abs(x-xi)<eps )// Case of a square geometry
384 && ( !check_groups || my_face.getGroupName()!=face_i.getGroupName()) //In case groups need to be checked
385 && ( !use_central_inversion || abs(y+yi) + abs(x+xi)<eps ) // Case of a central inversion
386 && fabs(my_face.getMeasure() - face_i.getMeasure())<eps
387 && fabs(my_face.getXN() + face_i.getXN())<eps
388 && fabs(my_face.getYN() + face_i.getYN())<eps
389 && fabs(my_face.getZN() + face_i.getZN())<eps )
391 iface_perio=_boundaryFaceIds[iface];
398 double x=my_face.x();
399 double y=my_face.y();
400 double z=my_face.z();
402 for (int iface=0;iface<_boundaryFaceIds.size() ; iface++)
404 Face face_i=_faces[_boundaryFaceIds[iface]];
405 double xi=face_i.x();
406 double yi=face_i.y();
407 double zi=face_i.z();
408 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
409 && ( !check_groups || my_face.getGroupName()!=face_i.getGroupName()) //In case groups need to be checked
410 && ( !use_central_inversion || abs(y+yi) + abs(x+xi) + abs(z+zi)<eps )// Case of a central inversion
411 && fabs(my_face.getMeasure() - face_i.getMeasure())<eps
412 && fabs(my_face.getXN() + face_i.getXN())<eps
413 && fabs(my_face.getYN() + face_i.getYN())<eps
414 && fabs(my_face.getZN() + face_i.getZN())<eps )
416 iface_perio=_boundaryFaceIds[iface];
422 throw CdmathException("Mesh::setPeriodicFaces: Mesh dimension should be 1, 2 or 3");
425 throw CdmathException("Mesh::setPeriodicFaces: periodic face not found, iface_perio==-1 " );
427 _indexFacePeriodicMap[_boundaryFaceIds[indexFace]]=iface_perio;
429 _indexFacePeriodicSet=true;
433 Mesh::getIndexFacePeriodic(int indexFace, bool check_groups, bool use_central_inversion)
435 if (!_faces[indexFace].isBorder())
437 cout<<"Pb with indexFace= "<<indexFace<<endl;
438 throw CdmathException("Mesh::getIndexFacePeriodic: not a border face" );
441 if(!_indexFacePeriodicSet)
442 setPeriodicFaces(check_groups, use_central_inversion);
444 std::map<int,int>::const_iterator it = _indexFacePeriodicMap.find(indexFace);
445 if( it != _indexFacePeriodicMap.end() )
449 cout<<"Pb with indexFace= "<<indexFace<<endl;
450 throw CdmathException("Mesh::getIndexFacePeriodic: not a periodic face" );
455 Mesh::isBorderNode(int nodeid) const
457 return _nodes[nodeid].isBorder();
461 Mesh::isBorderFace(int faceid) const
463 return _faces[faceid].isBorder();
467 Mesh::getBoundaryFaceIds() const
469 return _boundaryFaceIds;
473 Mesh::getBoundaryNodeIds() const
475 return _boundaryNodeIds;
479 Mesh::isTriangular() const
481 return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_TRI3;
484 Mesh::isTetrahedral() const
486 return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_TETRA4;
489 Mesh::isQuadrangular() const
491 return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_QUAD4;
494 Mesh::isHexahedral() const
496 return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_HEXA8;
499 Mesh::isStructured() const
501 return _isStructured;
504 std::vector< INTERP_KERNEL::NormalizedCellType >
505 Mesh::getElementTypes() const
510 std::vector< string >
511 Mesh::getElementTypesNames() const
513 std::vector< string > result(0);
514 for(int i=0; i< _eltsTypes.size(); i++)
516 if( _eltsTypes[i]==INTERP_KERNEL::NORM_POINT1)
517 result.push_back("Points");
518 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_SEG2)
519 result.push_back("Segments");
520 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_TRI3)
521 result.push_back("Triangles");
522 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_QUAD4)
523 result.push_back("Quadrangles");
524 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_POLYGON)
525 result.push_back("Polygons");
526 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_TETRA4)
527 result.push_back("Tetrahedra");
528 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_PYRA5)
529 result.push_back("Pyramids");
530 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_PENTA6)
531 result.push_back("Pentahedra");
532 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_HEXA8)
533 result.push_back("Hexahedra");
534 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_POLYHED)
535 result.push_back("Polyhedrons");
538 cout<< "Mesh " + _name + " contains an element of type " <<endl;
539 cout<< _eltsTypes[i]<<endl;
540 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");
547 Mesh::setGroups( const MEDFileUMesh* medmesh, MEDCouplingUMesh* mu)
549 //Searching for face groups
550 vector<string> faceGroups=medmesh->getGroupsNames() ;
552 for (unsigned int i=0;i<faceGroups.size();i++ )
554 string groupName=faceGroups[i];
555 vector<int> nonEmptyGrp(medmesh->getGrpNonEmptyLevels(groupName));
556 //We check if the group has a relative dimension equal to -1
557 //before call to the function getGroup(-1,groupName.c_str())
558 vector<int>::iterator it = find(nonEmptyGrp.begin(), nonEmptyGrp.end(), -1);
559 if (it != nonEmptyGrp.end())
561 cout<<"Boundary face group named "<< groupName << " found"<<endl;
562 MEDCouplingUMesh *m=medmesh->getGroup(-1,groupName.c_str());
564 _faceGroups.insert(_faceGroups.begin(),m);
565 _faceGroupNames.insert(_faceGroupNames.begin(),groupName);
566 DataArrayDouble *baryCell = m->computeIsoBarycenterOfNodesPerCell();//computeCellCenterOfMass() ;
567 const double *coorBary=baryCell->getConstPointer();
569 int nbCellsSubMesh=m->getNumberOfCells();
570 for (int ic(0), k(0); ic<nbCellsSubMesh; ic++, k+=_spaceDim)
572 vector<double> coorBaryXyz(3,0);
573 for (int d=0; d<_spaceDim; d++)
574 coorBaryXyz[d] = coorBary[k+d];
575 Point p1(coorBaryXyz[0],coorBaryXyz[1],coorBaryXyz[2]) ;
578 /* double min=INFINITY;
581 for (int iface=0;iface<_numberOfFaces;iface++ )
583 Point p2=_faces[iface].getBarryCenter();
584 /*if(groupName=="Wall" and p1.distance(p2)<min)
590 if(p1.distance(p2)<1.E-10)
592 _faces[iface].setGroupName(groupName);
597 /* if(groupName=="Wall" and min>1.E-10)
600 cout<<"Subcell number "<< ic <<" Min distance to Wall = "<<min <<" p= "<<p1[0]<<" , "<<p1[1]<<" , "<<p1[2]<<endl;
601 cout<<" attained at face "<< closestFaceNb << " p_min= "<<p3[0]<<" , "<<p3[1]<<" , "<<p3[2]<<endl;
604 throw CdmathException("No face belonging to group " + groupName + " found");
611 //Searching for node groups
612 vector<string> nodeGroups=medmesh->getGroupsOnSpecifiedLev(1) ;
614 for (unsigned int i=0;i<nodeGroups.size();i++ )
616 string groupName=nodeGroups[i];
617 DataArrayInt * nodeGroup=medmesh->getNodeGroupArr( groupName );
618 const int *nodeids=nodeGroup->getConstPointer();
622 cout<<"Boundary node group named "<< groupName << " found"<<endl;
624 _nodeGroups.insert(_nodeGroups.begin(),nodeGroup);
625 _nodeGroupNames.insert(_nodeGroupNames.begin(),groupName);
627 int nbNodesSubMesh=nodeGroup->getNumberOfTuples();//nodeGroup->getNbOfElems();
629 DataArrayDouble *coo = mu->getCoords() ;
630 const double *cood=coo->getConstPointer();
632 for (int ic(0); ic<nbNodesSubMesh; ic++)
634 vector<double> coorP(3,0);
635 for (int d=0; d<_spaceDim; d++)
636 coorP[d] = cood[nodeids[ic]*_spaceDim+d];
637 Point p1(coorP[0],coorP[1],coorP[2]) ;
640 for (int inode=0;inode<_numberOfNodes;inode++ )
642 Point p2=_nodes[inode].getPoint();
643 if(p1.distance(p2)<1.E-10)
645 _nodes[inode].setGroupName(groupName);
651 throw CdmathException("No node belonging to group " + groupName + " found");
657 //----------------------------------------------------------------------
659 Mesh::setMesh( void )
660 //----------------------------------------------------------------------
662 DataArrayInt *desc = DataArrayInt::New();
663 DataArrayInt *descI = DataArrayInt::New();
664 DataArrayInt *revDesc = DataArrayInt::New();
665 DataArrayInt *revDescI = DataArrayInt::New();
666 MEDCouplingUMesh* mu = _mesh->buildUnstructured();
669 /* Save the boundary mesh for later use*/
670 _boundaryMesh = mu->computeSkin();
672 MEDCouplingUMesh* mu2=mu->buildDescendingConnectivity2(desc,descI,revDesc,revDescI);//mesh of dimension N-1 containing the cell interfaces
674 const int *tmp = desc->getConstPointer();
675 const int *tmpI=descI->getConstPointer();
677 const int *tmpA =revDesc->getConstPointer();
678 const int *tmpAI=revDescI->getConstPointer();
680 //const int *work=tmp+tmpI[id];//corresponds to buildDescendingConnectivity
682 //Test du type d'éléments contenus dans le maillage afin d'éviter les éléments contenant des points de gauss
683 _eltsTypes=mu->getAllGeoTypesSorted();
684 for(int i=0; i<_eltsTypes.size();i++)
687 _eltsTypes[i]!= INTERP_KERNEL::NORM_POINT1 && _eltsTypes[i]!= INTERP_KERNEL::NORM_SEG2
688 && _eltsTypes[i]!= INTERP_KERNEL::NORM_TRI3 && _eltsTypes[i]!= INTERP_KERNEL::NORM_QUAD4
689 && _eltsTypes[i]!= INTERP_KERNEL::NORM_TETRA4 && _eltsTypes[i]!= INTERP_KERNEL::NORM_PYRA5
690 && _eltsTypes[i]!= INTERP_KERNEL::NORM_PENTA6 && _eltsTypes[i]!= INTERP_KERNEL::NORM_HEXA8
691 && _eltsTypes[i]!= INTERP_KERNEL::NORM_POLYGON&& _eltsTypes[i]!= INTERP_KERNEL::NORM_POLYHED
694 cout<< "Mesh " + mu->getName() + " contains an element of type " <<endl;
695 cout<< _eltsTypes[i]<<endl;
696 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");
700 DataArrayDouble *baryCell = mu->computeCellCenterOfMass() ;
701 const double *coorBary=baryCell->getConstPointer();
703 MEDCouplingFieldDouble* fields=mu->getMeasureField(true);
704 DataArrayDouble *surface = fields->getArray();
705 const double *surf=surface->getConstPointer();
707 DataArrayDouble *coo = mu->getCoords() ;
708 const double *cood=coo->getConstPointer();
710 DataArrayInt *revNode =DataArrayInt::New();
711 DataArrayInt *revNodeI=DataArrayInt::New();
712 mu->getReverseNodalConnectivity(revNode,revNodeI) ;
713 const int *tmpN =revNode->getConstPointer();
714 const int *tmpNI=revNodeI->getConstPointer();
716 DataArrayInt *revCell =DataArrayInt::New();
717 DataArrayInt *revCellI=DataArrayInt::New();
718 mu2->getReverseNodalConnectivity(revCell,revCellI) ;
719 const int *tmpC =revCell->getConstPointer();
720 const int *tmpCI=revCellI->getConstPointer();
722 const DataArrayInt *nodal = mu2->getNodalConnectivity() ;
723 const DataArrayInt *nodalI = mu2->getNodalConnectivityIndex() ;
724 const int *tmpNE =nodal->getConstPointer();
725 const int *tmpNEI=nodalI->getConstPointer();
727 _numberOfCells = mu->getNumberOfCells() ;
728 _cells = new Cell[_numberOfCells] ;
730 _numberOfNodes = mu->getNumberOfNodes() ;
731 _nodes = new Node[_numberOfNodes] ;
733 _numberOfFaces = mu2->getNumberOfCells();
734 _faces = new Face[_numberOfFaces] ;
736 _indexFacePeriodicSet=false;
738 //Definition used if _meshDim =3 to determine the edges
739 DataArrayInt *desc2 =DataArrayInt::New();
740 DataArrayInt *descI2=DataArrayInt::New();
741 DataArrayInt *revDesc2 =DataArrayInt::New();
742 DataArrayInt *revDescI2=DataArrayInt::New();
743 DataArrayInt *revNode2 =DataArrayInt::New();
744 DataArrayInt *revNodeI2=DataArrayInt::New();
747 MEDCouplingUMesh* mu3;
750 _numberOfEdges = mu->getNumberOfCells();
751 else if (_meshDim == 2)
752 _numberOfEdges = mu2->getNumberOfCells();
755 mu3=mu2->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2);//1D mesh of segments
756 _numberOfEdges = mu3->getNumberOfCells();
757 mu3->getReverseNodalConnectivity(revNode2,revNodeI2) ;
758 tmpN2 =revNode2->getConstPointer();
759 tmpNI2=revNodeI2->getConstPointer();
762 // _cells, _nodes and _faces initialization:
765 for( int id=0;id<_numberOfCells;id++ )
767 const int *work=tmp+tmpI[id];
768 int nbFaces=tmpI[id+1]-tmpI[id];
770 int nbVertices=mu->getNumberOfNodesInCell(id) ;
772 Cell ci( nbVertices, nbFaces, surf[id], Point(coorBary[id], 0.0, 0.0) ) ;
774 std::vector<int> nodeIdsOfCell ;
775 mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
776 for( int el=0;el<nbVertices;el++ )
778 ci.addNodeId(el,nodeIdsOfCell[el]) ;
779 ci.addFaceId(el,nodeIdsOfCell[el]) ;
782 //This assumes that in a 1D mesh the cell numbers form a consecutive sequence
783 //going either from left to right (xn=-1) or right to left (xn=1)
784 double xn = (cood[nodeIdsOfCell[nbVertices-1]] - cood[nodeIdsOfCell[0]] > 0.0) ? -1.0 : 1.0;
786 for( int el=0;el<nbFaces;el++ )
788 ci.addNormalVector(el,xn,0.0,0.0) ;
789 int indexFace=abs(work[el])-1;
790 ci.addFaceId(el,indexFace) ;
796 for( int id(0), k(0); id<_numberOfNodes; id++, k+=_spaceDim)
798 Point p(cood[k], 0.0, 0.0) ;
799 const int *workc=tmpN+tmpNI[id];
800 int nbCells=tmpNI[id+1]-tmpNI[id];
802 const int *workf=tmpC+tmpCI[id];
803 int nbFaces=tmpCI[id+1]-tmpCI[id];
804 const int *workn=tmpA+tmpAI[id];
805 int nbNeighbourNodes=tmpAI[id+1]-tmpAI[id];
806 Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
808 for( int el=0;el<nbCells;el++ )
809 vi.addCellId(el,workc[el]) ;
810 for( int el=0;el<nbFaces;el++ )
811 vi.addFaceId(el,workf[el]) ;
812 for( int el=0;el<nbNeighbourNodes;el++ )
813 vi.addNeighbourNodeId(el,workn[el]) ;
816 _boundaryNodeIds.push_back(0);
817 _boundaryNodeIds.push_back(_numberOfNodes-1);
819 for(int id(0), k(0); id<_numberOfFaces; id++, k+=_spaceDim)
821 Point p(cood[k], 0.0, 0.0) ;
822 const int *workc=tmpA+tmpAI[id];
823 int nbCells=tmpAI[id+1]-tmpAI[id];
825 const int *workv=tmpNE+tmpNEI[id]+1;
826 Face fi( 1, nbCells, 1.0, p, 1.0, 0.0, 0.0) ;
827 fi.addNodeId(0,workv[0]) ;
829 for(int idCell=0; idCell<nbCells; idCell++)
830 fi.addCellId(idCell,workc[idCell]) ;
834 _boundaryFaceIds.push_back(0);
835 _boundaryFaceIds.push_back(_numberOfFaces-1);
837 else if(_spaceDim==2 || _spaceDim==3)
839 DataArrayDouble *barySeg = mu2->computeIsoBarycenterOfNodesPerCell();//computeCellCenterOfMass() ;
840 const double *coorBarySeg=barySeg->getConstPointer();
842 MEDCouplingFieldDouble* fieldn;
843 DataArrayDouble *normal;
844 const double *tmpNormal;
846 if(_spaceDim==_meshDim)
847 fieldn = mu2->buildOrthogonalField();//Compute the normal to each cell interface
849 fieldn = mu->buildOrthogonalField();//compute the 3D normal vector to the 2D cell
851 normal = fieldn->getArray();
852 tmpNormal = normal->getConstPointer();
854 /*Building mesh cells */
855 for(int id(0), k(0); id<_numberOfCells; id++, k+=_spaceDim)
857 const int *work=tmp+tmpI[id];
858 int nbFaces=tmpI[id+1]-tmpI[id];
860 int nbVertices=mu->getNumberOfNodesInCell(id) ;
862 vector<double> coorBaryXyz(3,0);
863 for (int d=0; d<_spaceDim; d++)
864 coorBaryXyz[d] = coorBary[k+d];
866 Point p(coorBaryXyz[0],coorBaryXyz[1],coorBaryXyz[2]) ;
867 Cell ci( nbVertices, nbFaces, surf[id], p ) ;
869 /* Filling cell nodes */
870 std::vector<int> nodeIdsOfCell ;
871 mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
872 for( int el=0;el<nbVertices;el++ )
873 ci.addNodeId(el,nodeIdsOfCell[el]) ;
875 /* Filling cell faces */
876 if(_spaceDim==_meshDim)//use the normal field generated by buildOrthogonalField()
877 for( int el=0;el<nbFaces;el++ )
879 int faceIndex=(abs(work[el])-1);//=work[el] since Fortran type numbering was used, and negative sign means anticlockwise numbering
880 vector<double> xyzn(3,0);//Outer normal to the cell
882 for (int d=0; d<_spaceDim; d++)
883 xyzn[d] = -tmpNormal[_spaceDim*faceIndex+d];
885 for (int d=0; d<_spaceDim; d++)
886 xyzn[d] = +tmpNormal[_spaceDim*faceIndex+d];
887 ci.addNormalVector(el,xyzn[0],xyzn[1],xyzn[2]) ;
888 ci.addFaceId(el,faceIndex) ;
890 else//build normals associated to the couple (cell id, face el)
892 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
893 {//work[0]= first face global number, work[1]= second face global number
894 int indexFace0=abs(work[0])-1;//=work[0] since Fortran type numbering was used, and negative sign means anticlockwise numbering
895 int indexFace1=abs(work[1])-1;//=work[1] since Fortran type numbering was used, and negative sign means anticlockwise numbering
896 int idNodeA=(tmpNE+tmpNEI[indexFace0]+1)[0];//global number of the first face node work[0]=(abs(work[0])-1)
897 int idNodeB=(tmpNE+tmpNEI[indexFace1]+1)[0];//global number of the second face node work[1]=(abs(work[1])-1)
899 for(int i=0;i<_spaceDim;i++)
900 vecAB[i]=coo->getIJ(idNodeB,i) - coo->getIJ(idNodeA,i);
902 ci.addNormalVector(0,-vecAB[0],-vecAB[1],-vecAB[2]) ;
903 ci.addNormalVector(1,vecAB[0],vecAB[1],vecAB[2]) ;
904 ci.addFaceId(0,indexFace0) ;
905 ci.addFaceId(1,indexFace1) ;
907 else//_meshDim==2, number of faces around the cell id is variable, each face is composed of two nodes
910 for (int d=0; d<_spaceDim; d++)
911 xyzn[d] = tmpNormal[_spaceDim*id+d];
912 for( int el=0;el<nbFaces;el++ )
914 int faceIndex=(abs(work[el])-1);//=work[el] since Fortran type numbering was used, and negative sign means anticlockwise numbering
915 const int *workv=tmpNE+tmpNEI[faceIndex]+1;
916 int nbNodes= tmpNEI[faceIndex+1]-tmpNEI[faceIndex]-1;
917 if(nbNodes!=2)//We want to compute the normal to a straight line, not a curved interface composed of more thant 2 points
919 cout<<"Mesh name "<< mu->getName()<< " space dim= "<< _spaceDim <<" mesh dim= "<< _meshDim <<endl;
920 cout<<"For cell id "<<id<<" and local face number "<<el<<", the number of nodes is "<< nbNodes<< ", total number of faces is "<< nbFaces <<endl;
921 throw CdmathException("Mesh::setMesh number of nodes around a face should be 2");
924 int idNodeA=workv[0];
925 int idNodeB=workv[1];
926 vector<double> nodeA(_spaceDim), nodeB(_spaceDim), nodeP(_spaceDim);
927 for(int i=0;i<_spaceDim;i++)
929 nodeA[i]=coo->getIJ(idNodeA,i);
930 nodeB[i]=coo->getIJ(idNodeB,i);
931 nodeP[i]=coorBary[_spaceDim*id+i];
933 //Let P be the barycenter of the cell id
934 Vector vecAB(3), vecPA(3);
935 for(int i=0;i<_spaceDim;i++)
937 vecAB[i]=coo->getIJ(idNodeB,i) - coo->getIJ(idNodeA,i);
938 vecPA[i]=coo->getIJ(idNodeA,i) - coorBary[_spaceDim*id+i];
941 Vector normale = xyzn % vecAB;//Normal to the edge
942 normale/=normale.norm();
945 ci.addNormalVector(el,normale[0],normale[1],normale[2]) ;
947 ci.addNormalVector(el,-normale[0],-normale[1],-normale[2]) ;
948 ci.addFaceId(el,faceIndex) ;
955 MEDCouplingFieldDouble* fieldl=mu2->getMeasureField(true);
956 DataArrayDouble *longueur = fieldl->getArray();
957 const double *lon=longueur->getConstPointer();
959 if(_spaceDim!=_meshDim)
961 /* Since spaceDim!=meshDim, don't build normal to faces */
967 /*Building mesh faces */
968 for(int id(0), k(0); id<_numberOfFaces; id++, k+=_spaceDim)
970 vector<double> coorBarySegXyz(3,0);
971 for (int d=0; d<_spaceDim; d++)
972 coorBarySegXyz[d] = coorBarySeg[k+d];
973 Point p(coorBarySegXyz[0],coorBarySegXyz[1],coorBarySegXyz[2]) ;
974 const int *workc=tmpA+tmpAI[id];
975 int nbCells=tmpAI[id+1]-tmpAI[id];
977 if (nbCells>2 && _spaceDim==_meshDim)
979 cout<<"Warning : nbCells>2, numberOfFaces="<<_numberOfFaces<<endl;
980 cout<<"nbCells= "<<nbCells<<", _spaceDim="<<_spaceDim<<", _meshDim="<<_meshDim<<endl;
981 for(int icell=0; icell<nbCells; icell++)
982 cout<<workc[icell]<<", ";
984 throw CdmathException("Wrong mesh : nbCells>2 and spaceDim==meshDim");
987 _boundaryFaceIds.push_back(id);
989 const int *workv=tmpNE+tmpNEI[id]+1;
990 int nbNodes= tmpNEI[id+1]-tmpNEI[id]-1;
993 if(_spaceDim==_meshDim)//Euclidean flat mesh geometry
995 fi=Face( nbNodes, nbCells, lon[id], p, tmpNormal[k], tmpNormal[k+1], 0.0) ;
997 fi=Face( nbNodes, nbCells, lon[id], p, tmpNormal[k], tmpNormal[k+1], tmpNormal[k+2]) ;
998 else//Curved mesh geometry
999 fi=Face( nbNodes, nbCells, lon[id], p, 0.0, 0.0, 0.0) ;//Since spaceDim!=meshDim, normal to face is not defined
1001 for(int node_id=0; node_id<nbNodes;node_id++)
1002 fi.addNodeId(node_id,workv[node_id]) ;
1004 fi.addCellId(0,workc[0]) ;
1005 for(int cell_id=1; cell_id<nbCells;cell_id++)
1008 if (workc[cell_id]!=workc[cell_id-1])//For some meshes (bad ones) the same cell can appear several times
1010 fi.addCellId(cell_idx+1,workc[cell_id]) ;
1018 /*Building mesh nodes, should be done after building mesh faces in order to detect boundary nodes*/
1019 for(int id(0), k(0); id<_numberOfNodes; id++, k+=_spaceDim)
1021 vector<double> coorP(3,0);
1022 for (int d=0; d<_spaceDim; d++)
1023 coorP[d] = cood[k+d];
1024 Point p(coorP[0],coorP[1],coorP[2]) ;
1026 const int *workc=tmpN+tmpNI[id];
1027 int nbCells=tmpNI[id+1]-tmpNI[id];
1028 const int *workf=tmpC+tmpCI[id];
1029 int nbFaces=tmpCI[id+1]-tmpCI[id];
1031 int nbNeighbourNodes;
1034 workn=tmpA+tmpAI[id];
1035 nbNeighbourNodes=tmpAI[id+1]-tmpAI[id];
1037 else if (_meshDim == 2)
1039 workn=tmpC+tmpCI[id];
1040 nbNeighbourNodes=tmpCI[id+1]-tmpCI[id];
1044 workn=tmpN2+tmpNI2[id];
1045 nbNeighbourNodes=tmpNI2[id+1]-tmpNI2[id];
1047 Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
1049 for( int el=0;el<nbCells;el++ )
1050 vi.addCellId(el,workc[el]) ;
1051 for( int el=0;el<nbNeighbourNodes;el++ )
1052 vi.addNeighbourNodeId(el,workn[el]) ;
1053 //Detection of border nodes
1054 bool isBorder=false;
1055 for( int el=0;el<nbFaces;el++ )
1057 vi.addFaceId(el,workf[el],_faces[workf[el]].isBorder()) ;
1058 isBorder= isBorder || _faces[workf[el]].isBorder();
1061 _boundaryNodeIds.push_back(id);
1065 if(_spaceDim==_meshDim)
1071 throw CdmathException("Mesh::setMesh space dimension should be 1, 2 or 3");
1073 //definition of the bounding box for unstructured meshes
1074 if(!_isStructured)//Structured case is treated in function readMeshMed
1076 double Box0[2*_spaceDim];
1077 coo->getMinMaxPerComponent(Box0);
1093 //Set boundary groups
1094 _faceGroupNames.push_back("Boundary");
1095 _nodeGroupNames.push_back("Boundary");
1100 revDescI->decrRef();
1102 baryCell->decrRef();
1105 revNodeI->decrRef();
1107 revCellI->decrRef();
1111 revNode2->decrRef();
1112 revNodeI2->decrRef();
1115 revDesc2->decrRef();
1116 revDescI2->decrRef();
1123 //----------------------------------------------------------------------
1124 Mesh::Mesh( double xmin, double xmax, int nx, std::string meshName )
1125 //----------------------------------------------------------------------
1128 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx) : nx <= 0");
1130 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx) : xmin >= xmax");
1132 double dx = (xmax - xmin)/nx ;
1137 _isStructured = true;
1138 _indexFacePeriodicSet=false;
1147 _dxyz.resize(_spaceDim);
1149 _nxyz.resize(_spaceDim);
1152 double *originPtr = new double[_spaceDim];
1153 double *dxyzPtr = new double[_spaceDim];
1154 int *nodeStrctPtr = new int[_spaceDim];
1157 nodeStrctPtr[0]=nx+1;
1160 _mesh=MEDCouplingIMesh::New(meshName,
1163 nodeStrctPtr+_spaceDim,
1165 originPtr+_spaceDim,
1168 delete [] originPtr;
1170 delete [] nodeStrctPtr;
1172 DataArrayInt *desc=DataArrayInt::New();
1173 DataArrayInt *descI=DataArrayInt::New();
1174 DataArrayInt *revDesc=DataArrayInt::New();
1175 DataArrayInt *revDescI=DataArrayInt::New();
1176 MEDCouplingUMesh* mu=_mesh->buildUnstructured();
1177 MEDCouplingUMesh *mu2=mu->buildDescendingConnectivity(desc,descI,revDesc,revDescI);
1179 const int *tmp=desc->getConstPointer();
1180 const int *tmpI=descI->getConstPointer();
1182 const int *tmpA =revDesc->getConstPointer();
1183 const int *tmpAI=revDescI->getConstPointer();
1185 _eltsTypes=mu->getAllGeoTypesSorted();
1187 DataArrayDouble *baryCell = mu->computeCellCenterOfMass() ;
1188 const double *coorBary=baryCell->getConstPointer();
1190 _numberOfCells = _mesh->getNumberOfCells() ;
1191 _cells = new Cell[_numberOfCells] ;
1193 _numberOfNodes = mu->getNumberOfNodes() ;
1194 _nodes = new Node[_numberOfNodes] ;
1196 _numberOfFaces = _numberOfNodes;
1197 _faces = new Face[_numberOfFaces] ;
1199 _numberOfEdges = _numberOfCells;
1201 MEDCouplingFieldDouble* fieldl=mu->getMeasureField(true);
1202 DataArrayDouble *longueur = fieldl->getArray();
1203 const double *lon=longueur->getConstPointer();
1205 DataArrayDouble *coo = mu->getCoords() ;
1206 const double *cood=coo->getConstPointer();
1209 for( int id=0;id<_numberOfCells;id++ )
1211 int nbVertices=mu->getNumberOfNodesInCell(id) ;
1212 Point p(coorBary[id],0.0,0.0) ;
1213 Cell ci( nbVertices, nbVertices, lon[id], p ) ;
1215 std::vector<int> nodeIdsOfCell ;
1216 mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
1217 for( int el=0;el<nbVertices;el++ )
1219 ci.addNodeId(el,nodeIdsOfCell[el]) ;
1220 ci.addFaceId(el,nodeIdsOfCell[el]) ;
1223 double xn = (cood[nodeIdsOfCell[nbVertices-1]] - cood[nodeIdsOfCell[0]] > 0.0) ? -1.0 : 1.0;
1225 int nbFaces=tmpI[id+1]-tmpI[id];
1226 const int *work=tmp+tmpI[id];
1228 for( int el=0;el<nbFaces;el++ )
1230 ci.addNormalVector(el,xn,0.0,0.0) ;
1231 ci.addFaceId(el,work[el]) ;
1240 //Suppress the following since tmpN=tmpA
1241 DataArrayInt *revNode=DataArrayInt::New();
1242 DataArrayInt *revNodeI=DataArrayInt::New();
1243 mu->getReverseNodalConnectivity(revNode,revNodeI) ;
1244 const int *tmpN=revNode->getConstPointer();
1245 const int *tmpNI=revNodeI->getConstPointer();
1247 for( int id=0;id<_numberOfNodes;id++ )
1249 std::vector<double> coo ;
1250 mu->getCoordinatesOfNode(id,coo);
1251 Point p(coo[0],0.0,0.0) ;
1252 const int *workc=tmpN+tmpNI[id];
1253 int nbCells=tmpNI[id+1]-tmpNI[id];
1255 const int *workn=tmpA+tmpAI[id];
1256 int nbNeighbourNodes=tmpAI[id+1]-tmpAI[id];
1258 Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
1259 for( int el=0;el<nbCells;el++ )
1260 vi.addCellId(el,workc[el]) ;
1261 for( int el=0;el<nbFaces;el++ )
1262 vi.addFaceId(el,id) ;
1263 for( int el=0;el<nbNeighbourNodes;el++ )
1264 vi.addNeighbourNodeId(el,workn[el]) ;
1269 Face fi( nbVertices, nbCells, 1.0, p, 1., 0., 0. ) ;
1270 for( int el=0;el<nbVertices;el++ )
1271 fi.addNodeId(el,id) ;
1273 for( int el=0;el<nbCells;el++ )
1274 fi.addCellId(el,workc[el]) ;
1278 //Set boundary groups
1279 _faceGroupNames.push_back("Boundary");
1280 _nodeGroupNames.push_back("Boundary");
1281 _boundaryFaceIds.push_back(0);
1282 _boundaryFaceIds.push_back(_numberOfFaces-1);
1283 _boundaryNodeIds.push_back(0);
1284 _boundaryNodeIds.push_back(_numberOfNodes-1);
1287 baryCell->decrRef();
1291 revDescI->decrRef();
1293 revNodeI->decrRef();
1298 //----------------------------------------------------------------------
1299 Mesh::Mesh( std::vector<double> points, std::string meshName )
1300 //----------------------------------------------------------------------
1302 int nx=points.size();
1305 throw CdmathException("Mesh::Mesh( vector<double> points, string meshName) : nx < 2, vector should contain at least two values");
1307 while( i<nx-1 && points[i+1]>points[i] )
1311 //cout<< points << endl;
1312 throw CdmathException("Mesh::Mesh( vector<double> points, string meshName) : vector values should be sorted");
1325 _isStructured = false;
1326 _indexFacePeriodicSet=false;
1328 MEDCouplingUMesh * mesh1d = MEDCouplingUMesh::New(meshName, 1);
1329 mesh1d->allocateCells(nx - 1);
1330 double * coords = new double[nx];
1331 int * nodal_con = new int[2];
1332 coords[0]=points[0];
1333 for(int i=0; i<nx- 1 ; i++)
1337 mesh1d->insertNextCell(INTERP_KERNEL::NORM_SEG2, 2, nodal_con);
1338 coords[i+1]=points[i + 1];
1340 mesh1d->finishInsertingCells();
1342 DataArrayDouble * coords_arr = DataArrayDouble::New();
1343 coords_arr->alloc(nx,1);
1344 std::copy(coords,coords+nx,coords_arr->getPointer());
1345 mesh1d->setCoords(coords_arr);
1347 delete [] coords, nodal_con;
1348 coords_arr->decrRef();
1352 DataArrayInt *desc=DataArrayInt::New();
1353 DataArrayInt *descI=DataArrayInt::New();
1354 DataArrayInt *revDesc=DataArrayInt::New();
1355 DataArrayInt *revDescI=DataArrayInt::New();
1356 MEDCouplingUMesh* mu=_mesh->buildUnstructured();
1357 MEDCouplingUMesh *mu2=mu->buildDescendingConnectivity(desc,descI,revDesc,revDescI);
1359 DataArrayDouble *baryCell = mu->computeCellCenterOfMass() ;
1360 const double *coorBary=baryCell->getConstPointer();
1362 _numberOfCells = _mesh->getNumberOfCells() ;
1363 _cells = new Cell[_numberOfCells] ;
1365 _numberOfEdges = _numberOfCells;
1367 _eltsTypes=mu->getAllGeoTypesSorted();
1369 MEDCouplingFieldDouble* fieldl=mu->getMeasureField(true);
1370 DataArrayDouble *longueur = fieldl->getArray();
1371 const double *lon=longueur->getConstPointer();
1374 for( int id=0;id<_numberOfCells;id++ )
1376 int nbVertices=mu->getNumberOfNodesInCell(id) ;
1377 Point p(coorBary[id],0.0,0.0) ;
1378 Cell ci( nbVertices, nbVertices, lon[id], p ) ;
1380 std::vector<int> nodeIdsOfCell ;
1381 mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
1382 for( int el=0;el<nbVertices;el++ )
1384 ci.addNodeId(el,nodeIdsOfCell[el]) ;
1385 ci.addFaceId(el,nodeIdsOfCell[el]) ;
1392 //Suppress the following since tmpN=tmpA
1393 DataArrayInt *revNode=DataArrayInt::New();
1394 DataArrayInt *revNodeI=DataArrayInt::New();
1395 mu->getReverseNodalConnectivity(revNode,revNodeI) ;
1396 const int *tmpN=revNode->getConstPointer();
1397 const int *tmpNI=revNodeI->getConstPointer();
1399 _numberOfNodes = mu->getNumberOfNodes() ;
1400 _nodes = new Node[_numberOfNodes] ;
1401 _numberOfFaces = _numberOfNodes;
1402 _faces = new Face[_numberOfFaces] ;
1404 for( int id=0;id<_numberOfNodes;id++ )
1406 std::vector<double> coo ;
1407 mu->getCoordinatesOfNode(id,coo);
1408 Point p(coo[0],0.0,0.0) ;
1409 const int *workc=tmpN+tmpNI[id];
1410 int nbCells=tmpNI[id+1]-tmpNI[id];
1412 const int *workn=tmpN+tmpNI[id];
1413 int nbNeighbourNodes=tmpNI[id+1]-tmpNI[id];
1414 Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
1416 /* provisoire !!!!!!!!!!!!*/
1417 // Point pf(0.0,0.0,0.0) ;
1418 Face fi( nbVertices, nbCells, 0.0, p, 0., 0., 0. ) ;
1420 for( int el=0;el<nbCells;el++ )
1421 vi.addCellId(el,workc[el]) ;
1422 for( int el=0;el<nbFaces;el++ )
1423 vi.addFaceId(el,id) ;
1424 for( int el=0;el<nbNeighbourNodes;el++ )
1425 vi.addNeighbourNodeId(el,workn[el]) ;
1428 for( int el=0;el<nbVertices;el++ )
1429 fi.addNodeId(el,id) ;
1431 for( int el=0;el<nbCells;el++ )
1432 fi.addCellId(el,workc[el]) ;
1435 //Set boundary groups
1436 _faceGroupNames.push_back("Boundary");
1437 _nodeGroupNames.push_back("Boundary");
1438 _boundaryFaceIds.push_back(0);
1439 _boundaryFaceIds.push_back(_numberOfFaces-1);
1440 _boundaryNodeIds.push_back(0);
1441 _boundaryNodeIds.push_back(_numberOfNodes-1);
1444 baryCell->decrRef();
1448 revDescI->decrRef();
1450 revNodeI->decrRef();
1454 //----------------------------------------------------------------------
1455 Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, int split_to_triangles_policy, std::string meshName)
1456 //----------------------------------------------------------------------
1459 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : nx <= 0 or ny <= 0");
1461 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : xmin >= xmax");
1463 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : ymin >= ymax");
1473 double dx = (xmax - xmin)/nx ;
1474 double dy = (ymax - ymin)/ny ;
1479 _indexFacePeriodicSet=false;
1480 _isStructured = true;
1481 _nxyz.resize(_spaceDim);
1485 _dxyz.resize(_spaceDim);
1489 double *originPtr = new double[_spaceDim];
1490 double *dxyzPtr = new double[_spaceDim];
1491 int *nodeStrctPtr = new int[_spaceDim];
1495 nodeStrctPtr[0]=nx+1;
1496 nodeStrctPtr[1]=ny+1;
1500 _mesh=MEDCouplingIMesh::New(meshName,
1503 nodeStrctPtr+_spaceDim,
1505 originPtr+_spaceDim,
1509 if(split_to_triangles_policy==0 || split_to_triangles_policy==1)
1511 _mesh=_mesh->buildUnstructured();
1512 _mesh->simplexize(split_to_triangles_policy);
1514 else if (split_to_triangles_policy != -1)
1516 cout<< "split_to_triangles_policy = "<< split_to_triangles_policy << endl;
1517 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : Unknown splitting policy");
1520 delete [] originPtr;
1522 delete [] nodeStrctPtr;
1527 //----------------------------------------------------------------------
1528 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)
1529 //----------------------------------------------------------------------
1531 if(nx<=0 || ny<=0 || nz<=0)
1532 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");
1534 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : xmin >= xmax");
1536 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : ymin >= ymax");
1538 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : zmin >= zmax");
1543 _indexFacePeriodicSet=false;
1551 double dx = (xmax - xmin)/nx ;
1552 double dy = (ymax - ymin)/ny ;
1553 double dz = (zmax - zmin)/nz ;
1555 _isStructured = true;
1556 _dxyz.resize(_spaceDim);
1561 _nxyz.resize(_spaceDim);
1566 double *originPtr = new double[_spaceDim];
1567 double *dxyzPtr = new double[_spaceDim];
1568 int *nodeStrctPtr = new int[_spaceDim];
1573 nodeStrctPtr[0]=nx+1;
1574 nodeStrctPtr[1]=ny+1;
1575 nodeStrctPtr[2]=nz+1;
1580 _mesh=MEDCouplingIMesh::New(meshName,
1583 nodeStrctPtr+_spaceDim,
1585 originPtr+_spaceDim,
1589 if( split_to_tetrahedra_policy == 0 )
1591 _mesh=_mesh->buildUnstructured();
1592 _mesh->simplexize(INTERP_KERNEL::PLANAR_FACE_5);
1594 else if( split_to_tetrahedra_policy == 1 )
1596 _mesh=_mesh->buildUnstructured();
1597 _mesh->simplexize(INTERP_KERNEL::PLANAR_FACE_6);
1599 else if ( split_to_tetrahedra_policy != -1 )
1601 cout<< "split_to_tetrahedra_policy = "<< split_to_tetrahedra_policy << endl;
1602 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");
1605 delete [] originPtr;
1607 delete [] nodeStrctPtr;
1612 //----------------------------------------------------------------------
1614 Mesh::getSpaceDimension( void ) const
1615 //----------------------------------------------------------------------
1620 //----------------------------------------------------------------------
1622 Mesh::getMeshDimension( void ) const
1623 //----------------------------------------------------------------------
1629 Mesh::getDXYZ() const
1632 throw CdmathException("std::vector<double> Mesh::getDXYZ() : dx,dy and dz are defined only for structured meshes !");
1638 Mesh::getCellGridStructure() const
1641 throw CdmathException("std::vector<int> Mesh::getCellGridStructure() : nx, ny and nz are defined only for structured meshes !");
1646 //----------------------------------------------------------------------
1648 Mesh::getNx( void ) const
1649 //----------------------------------------------------------------------
1652 throw CdmathException("int Mesh::getNx( void ) : Nx is defined only for structured meshes !");
1657 //----------------------------------------------------------------------
1659 Mesh::getNy( void ) const
1660 //----------------------------------------------------------------------
1663 throw CdmathException("int Mesh::getNy( void ) : Ny is defined only for structured meshes !");
1665 throw CdmathException("int double& Field::operator[ielem] : Ny is not defined in dimension < 2!");
1670 //----------------------------------------------------------------------
1672 Mesh::getNz( void ) const
1673 //----------------------------------------------------------------------
1676 throw CdmathException("int Mesh::getNz( void ) : Nz is defined only for structured meshes !");
1678 throw CdmathException("int Mesh::getNz( void ) : Nz is not defined in dimension < 3!");
1683 //----------------------------------------------------------------------
1685 Mesh::getXMin( void ) const
1686 //----------------------------------------------------------------------
1691 //----------------------------------------------------------------------
1693 Mesh::getXMax( void ) const
1694 //----------------------------------------------------------------------
1699 //----------------------------------------------------------------------
1701 Mesh::getYMin( void ) const
1702 //----------------------------------------------------------------------
1707 //----------------------------------------------------------------------
1709 Mesh::getYMax( void ) const
1710 //----------------------------------------------------------------------
1715 //----------------------------------------------------------------------
1717 Mesh::getZMin( void ) const
1718 //----------------------------------------------------------------------
1723 //----------------------------------------------------------------------
1725 Mesh::getZMax( void ) const
1726 //----------------------------------------------------------------------
1731 //----------------------------------------------------------------------
1732 MCAuto<MEDCouplingMesh>
1733 Mesh::getMEDCouplingMesh( void ) const
1734 //----------------------------------------------------------------------
1739 //----------------------------------------------------------------------
1741 Mesh::getNumberOfNodes ( void ) const
1742 //----------------------------------------------------------------------
1744 return _numberOfNodes ;
1747 //----------------------------------------------------------------------
1749 Mesh::getNumberOfCells ( void ) const
1750 //----------------------------------------------------------------------
1752 return _numberOfCells ;
1755 //----------------------------------------------------------------------
1757 Mesh::getNumberOfFaces ( void ) const
1758 //----------------------------------------------------------------------
1760 return _numberOfFaces ;
1763 //----------------------------------------------------------------------
1765 Mesh::getNumberOfEdges ( void ) const
1766 //----------------------------------------------------------------------
1768 return _numberOfEdges ;
1771 //----------------------------------------------------------------------
1773 Mesh::getFaces ( void ) const
1774 //----------------------------------------------------------------------
1779 //----------------------------------------------------------------------
1781 Mesh::getCells ( void ) const
1782 //----------------------------------------------------------------------
1787 //----------------------------------------------------------------------
1789 Mesh::getCell ( int i ) const
1790 //----------------------------------------------------------------------
1795 //----------------------------------------------------------------------
1797 Mesh::getFace ( int i ) const
1798 //----------------------------------------------------------------------
1803 //----------------------------------------------------------------------
1805 Mesh::getNode ( int i ) const
1806 //----------------------------------------------------------------------
1811 //----------------------------------------------------------------------
1813 Mesh::getNodes ( void ) const
1814 //----------------------------------------------------------------------
1820 Mesh::getNameOfFaceGroups( void ) const
1822 return _faceGroupNames;
1825 vector<MEDCoupling::MEDCouplingUMesh *>
1826 Mesh::getFaceGroups( void ) const
1832 Mesh::getNameOfNodeGroups( void ) const
1834 return _nodeGroupNames;
1837 vector<MEDCoupling::DataArrayInt *>
1838 Mesh::getNodeGroups( void ) const
1843 //----------------------------------------------------------------------
1845 Mesh::operator= ( const Mesh& mesh )
1846 //----------------------------------------------------------------------
1848 _spaceDim = mesh.getSpaceDimension() ;
1849 _meshDim = mesh.getMeshDimension() ;
1850 _name = mesh.getName();
1851 _numberOfNodes = mesh.getNumberOfNodes();
1852 _numberOfFaces = mesh.getNumberOfFaces();
1853 _numberOfCells = mesh.getNumberOfCells();
1854 _numberOfEdges = mesh.getNumberOfEdges();
1856 _isStructured = mesh.isStructured();
1859 _nxyz = mesh.getCellGridStructure() ;
1860 _dxyz=mesh.getDXYZ();
1861 _xMin=mesh.getXMin();
1862 _xMax=mesh.getXMax();
1863 _yMin=mesh.getYMin();
1864 _yMax=mesh.getYMax();
1865 _zMin=mesh.getZMin();
1866 _zMax=mesh.getZMax();
1868 _faceGroupNames = mesh.getNameOfFaceGroups() ;
1869 _faceGroups = mesh.getFaceGroups() ;
1870 _nodeGroupNames = mesh.getNameOfNodeGroups() ;
1871 _nodeGroups = mesh.getNodeGroups() ;
1889 _nodes = new Node[_numberOfNodes] ;
1890 _faces = new Face[_numberOfFaces] ;
1891 _cells = new Cell[_numberOfCells] ;
1893 for (int i=0;i<_numberOfNodes;i++)
1894 _nodes[i]=mesh.getNode(i);
1896 for (int i=0;i<_numberOfFaces;i++)
1897 _faces[i]=mesh.getFace(i);
1899 for (int i=0;i<_numberOfCells;i++)
1900 _cells[i]=mesh.getCell(i);
1902 _indexFacePeriodicSet= mesh.isIndexFacePeriodicSet();
1903 if(_indexFacePeriodicSet)
1904 _indexFacePeriodicMap=mesh.getIndexFacePeriodic();
1906 _boundaryFaceIds=mesh.getBoundaryFaceIds();
1907 _boundaryNodeIds=mesh.getBoundaryNodeIds();
1909 _eltsTypes=mesh.getElementTypes();
1910 _eltsTypesNames=mesh.getElementTypesNames();
1912 MCAuto<MEDCouplingMesh> m1=mesh.getMEDCouplingMesh()->deepCopy();
1917 bool Mesh::isIndexFacePeriodicSet() const
1919 return _indexFacePeriodicSet;
1921 //----------------------------------------------------------------------
1923 Mesh::minRatioVolSurf()
1925 double dx_min = 1e30;
1926 for(int i=0; i<_numberOfCells; i++)
1928 Cell Ci = getCell(i);
1932 for(int k=0; k< Ci.getNumberOfFaces(); k++)
1934 int indexFace=Ci.getFacesId()[k];
1935 Face Fk = getFace(indexFace);
1936 perimeter+=Fk.getMeasure();
1938 dx_min = min(dx_min,Ci.getMeasure()/perimeter);
1941 dx_min = min(dx_min,Ci.getMeasure());
1948 Mesh::getBoundaryMesh ( void ) const
1950 return Mesh(_boundaryMesh);
1954 Mesh::getMaxNbNeighbours(EntityType type) const
1960 for(int i=0; i<_numberOfCells; i++)
1961 if(result < _cells[i].getNumberOfFaces())
1962 result=_cells[i].getNumberOfFaces();
1964 else if(type==NODES)
1966 for(int i=0; i<_numberOfNodes; i++)
1967 if(result < _nodes[i].getNumberOfEdges())
1968 result=_nodes[i].getNumberOfEdges();
1971 throw CdmathException("Mesh::getMaxNbNeighbours : entity type is not accepted. Should be CELLS or NODES");
1975 //----------------------------------------------------------------------
1977 Mesh::writeVTK ( const std::string fileName ) const
1978 //----------------------------------------------------------------------
1980 string fname=fileName+".vtu";
1981 _mesh->writeVTK(fname.c_str()) ;
1984 //----------------------------------------------------------------------
1986 Mesh::writeMED ( const std::string fileName ) const
1987 //----------------------------------------------------------------------
1989 string fname=fileName+".med";
1990 MEDCouplingUMesh* mu=_mesh->buildUnstructured();
1991 MEDCoupling::WriteUMesh(fname.c_str(),mu,true);
1993 //MEDFileUMesh meshMEDFile;
1994 //meshMEDFile.setMeshAtLevel(0,mu);
1995 //for(int i=0; i< _faceGroups.size(); i++)
1996 //meshMEDFile.setMeshAtLevel(-1,_faceGroups[i]);
1998 //MEDCoupling::meshMEDFile.write(fname.c_str(),2) ;
2000 //MEDCoupling::meshMEDFile.write(fname.c_str(),1) ;
2007 Mesh::getGroupFaceIds(std::string groupName) const
2009 if( std::find(_faceGroupNames.begin(),_faceGroupNames.end(),groupName) == _faceGroupNames.end() )
2011 cout<<"Mesh::getGroupFaceIds Error : face group " << groupName << " does not exist"<<endl;
2012 throw CdmathException("Required face group does not exist");
2016 std::vector< int > result(0);
2017 for(int i=0; i<_boundaryFaceIds.size(); i++)
2019 vector< std::string > groupNames = getFace(_boundaryFaceIds[i]).getGroupNames();
2020 if( std::find(groupNames.begin(), groupNames.end(),groupName) != groupNames.end() )
2021 result.push_back(_boundaryFaceIds[i]);
2028 Mesh::getGroupNodeIds(std::string groupName) const
2030 if( std::find(_nodeGroupNames.begin(),_nodeGroupNames.end(),groupName) == _nodeGroupNames.end() )
2032 cout<<"Mesh::getGroupNodeIds Error : node group " << groupName << " does not exist"<<endl;
2033 throw CdmathException("Required node group does not exist");
2037 std::vector< int > result(0);
2038 for(int i=0; i<_boundaryFaceIds.size(); i++)
2040 vector< std::string > groupNames = getNode(_boundaryFaceIds[i]).getGroupNames();
2041 if( std::find(groupNames.begin(), groupNames.end(),groupName) != groupNames.end() )
2042 result.push_back(_boundaryFaceIds[i]);
2049 Mesh::setFaceGroupByIds(std::vector< int > faceIds, std::string groupName)
2051 for(int i=0; i< faceIds.size(); i++)
2052 getFace(faceIds[i]).setGroupName(groupName);
2056 Mesh::setNodeGroupByIds(std::vector< int > nodeIds, std::string groupName)
2058 for(int i=0; i< nodeIds.size(); i++)
2059 getNode(nodeIds[i]).setGroupName(groupName);