4 * Created on: 22 janv. 2012
13 #include "CdmathException.hxx"
14 #include "MEDFileMesh.hxx"
15 #include "MEDLoader.hxx"
24 using namespace MEDCoupling;
27 //----------------------------------------------------------------------
29 //----------------------------------------------------------------------
50 _faceGroupNames.resize(0);
51 _faceGroups.resize(0);
52 _nodeGroupNames.resize(0);
53 _nodeGroups.resize(0);
54 _indexFacePeriodicSet=false;
58 //----------------------------------------------------------------------
60 //----------------------------------------------------------------------
65 //for(int i=0; i< _faceGroups.size(); i++)
66 // _faceGroups[i]->decrRef();
67 // _nodeGroups[i]->decrRef();
71 Mesh::getName( void ) const
76 Mesh::Mesh( const MEDCoupling::MEDCouplingIMesh* mesh )
78 _spaceDim=mesh->getSpaceDimension();
79 _meshDim=mesh->getMeshDimension();
81 _dxyz=mesh->getDXYZ();
82 _nxyz=mesh->getCellGridStructure();
83 double* Box0=new double[2*_spaceDim];
84 mesh->getBoundingBox(Box0);
85 _name=mesh->getName();
86 _indexFacePeriodicSet=false;
101 double *originPtr = new double[_spaceDim];
102 double *dxyzPtr = new double[_spaceDim];
103 mcIdType *nodeStrctPtr = new mcIdType[_spaceDim];
105 for(int i=0;i<_spaceDim;i++)
107 originPtr[i]=Box0[2*i];
108 nodeStrctPtr[i]=_nxyz[i]+1;
111 _mesh=MEDCouplingIMesh::New(_name,
114 nodeStrctPtr+_spaceDim,
121 delete [] nodeStrctPtr;
126 Mesh::Mesh( const MEDCoupling::MEDCouplingUMesh* mesh )
128 _spaceDim=mesh->getSpaceDimension();
129 _meshDim=mesh->getMeshDimension();
131 double* Box0=new double[2*_spaceDim];
132 mesh->getBoundingBox(Box0);
133 _name=mesh->getName();
134 _indexFacePeriodicSet=false;
149 _mesh=mesh->deepCopy();
154 //----------------------------------------------------------------------
155 Mesh::Mesh( const Mesh& mesh )
156 //----------------------------------------------------------------------
158 _spaceDim = mesh.getSpaceDimension() ;
159 _meshDim = mesh.getMeshDimension() ;
160 _name=mesh.getName();
161 _xMax=mesh.getXMax();
162 _yMin=mesh.getYMin();
163 _yMax=mesh.getYMax();
164 _zMin=mesh.getZMin();
165 _zMax=mesh.getZMax();
166 _isStructured=mesh.isStructured();
169 _nxyz = mesh.getCellGridStructure() ;
170 _dxyz=mesh.getDXYZ();
172 _numberOfNodes = mesh.getNumberOfNodes();
173 _numberOfFaces = mesh.getNumberOfFaces();
174 _numberOfCells = mesh.getNumberOfCells();
175 _numberOfEdges = mesh.getNumberOfEdges();
177 _faceGroupNames = mesh.getNameOfFaceGroups() ;
178 _faceGroups = mesh.getFaceGroups() ;
179 _nodeGroupNames = mesh.getNameOfNodeGroups() ;
180 _nodeGroups = mesh.getNodeGroups() ;
182 _nodes = new Node[_numberOfNodes] ;
183 _faces = new Face[_numberOfFaces] ;
184 _cells = new Cell[_numberOfCells] ;
186 //memcpy(_nodes,mesh.getNodes(),sizeof(*mesh.getNodes())) ;
187 //memcpy(_cells,mesh.getCells(),sizeof(*mesh.getCells())) ;
188 //memcpy(_faces,mesh.getFaces(),sizeof(*mesh.getFaces())) ;
190 for (int i=0;i<_numberOfNodes;i++)
191 _nodes[i]=mesh.getNode(i);
193 for (int i=0;i<_numberOfFaces;i++)
194 _faces[i]=mesh.getFace(i);
196 for (int i=0;i<_numberOfCells;i++)
197 _cells[i]=mesh.getCell(i);
199 _indexFacePeriodicSet= mesh.isIndexFacePeriodicSet();
200 if(_indexFacePeriodicSet)
201 _indexFacePeriodicMap=mesh.getIndexFacePeriodic();
203 _boundaryFaceIds=mesh.getBoundaryFaceIds();
204 _boundaryNodeIds=mesh.getBoundaryNodeIds();
206 _eltsTypes=mesh.getElementTypes();
207 _eltsTypesNames=mesh.getElementTypesNames();
209 MCAuto<MEDCouplingMesh> m1=mesh.getMEDCouplingMesh()->deepCopy();
213 //----------------------------------------------------------------------
214 Mesh::Mesh( const std::string filename, int meshLevel )
215 //----------------------------------------------------------------------
217 readMeshMed(filename, meshLevel);
220 //----------------------------------------------------------------------
222 Mesh::readMeshMed( const std::string filename, const int meshLevel)
223 //----------------------------------------------------------------------
225 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)
226 _mesh=m->getMeshAtLevel(meshLevel);
227 _mesh->checkConsistencyLight();
228 _mesh->setName(_mesh->getName());
229 _meshDim=_mesh->getMeshDimension();
230 _spaceDim=_mesh->getSpaceDimension();
231 _name=_mesh->getName();
232 _indexFacePeriodicSet=false;
233 MEDCoupling::MEDCouplingIMesh* structuredMesh = dynamic_cast<MEDCoupling::MEDCouplingIMesh*> (_mesh.retn());
237 _dxyz=structuredMesh->getDXYZ();
238 _nxyz=structuredMesh->getCellGridStructure();
239 double* Box0=new double[2*_spaceDim];
240 structuredMesh->getBoundingBox(Box0);
258 MEDCouplingUMesh* mu = setMesh();
261 cout<<endl<< "Loaded file "<< filename<<endl;
262 cout<<"Mesh name= "<<m->getName()<<", mesh dim="<< _meshDim<< ", space dim="<< _spaceDim<< ", nb cells= "<<getNumberOfCells()<< ", nb nodes= "<<getNumberOfNodes()<<endl;
269 Mesh::setGroupAtFaceByCoords(double x, double y, double z, double eps, std::string groupName)
271 int nbFace=getNumberOfFaces();
273 for (int iface=0;iface<nbFace;iface++)
275 double FX=_faces[iface].x();
276 double FY=_faces[iface].y();
277 double FZ=_faces[iface].z();
278 if (abs(FX-x)<eps && abs(FY-y)<eps && abs(FZ-z)<eps)
280 _faces[iface].setGroupName(groupName);
281 std::vector< int > nodesID= _faces[iface].getNodesId();
282 int nbNodes = _faces[iface].getNumberOfNodes();
283 for(int inode=0 ; inode<nbNodes ; inode++)
284 _nodes[nodesID[inode]].setGroupName(groupName);
291 _faceGroupNames.insert(_faceGroupNames.begin(),groupName);
292 _nodeGroupNames.insert(_nodeGroupNames.begin(),groupName);
293 //To do : update _faceGroups and _nodeGroups
298 Mesh::setGroupAtPlan(double value, int direction, double eps, std::string groupName)
300 int nbFace=getNumberOfFaces();
302 for (int iface=0;iface<nbFace;iface++)
304 double cord=_faces[iface].getBarryCenter()[direction];
305 if (abs(cord-value)<eps)
307 _faces[iface].setGroupName(groupName);
308 std::vector< int > nodesID= _faces[iface].getNodesId();
309 int nbNodes = _faces[iface].getNumberOfNodes();
310 for(int inode=0 ; inode<nbNodes ; inode++)
312 _nodes[nodesID[inode]].setGroupName(groupName);
320 _faceGroupNames.insert(_faceGroupNames.begin(),groupName);
321 _nodeGroupNames.insert(_nodeGroupNames.begin(),groupName);
322 //To do : update _faceGroups, _nodeGroups
327 Mesh::setBoundaryNodesFromFaces()
329 for (int iface=0;iface<_boundaryFaceIds.size();iface++)
331 std::vector< int > nodesID= _faces[_boundaryFaceIds[iface]].getNodesId();
332 int nbNodes = _faces[_boundaryFaceIds[iface]].getNumberOfNodes();
333 for(int inode=0 ; inode<nbNodes ; inode++)
335 std::vector<int>::const_iterator it = std::find(_boundaryNodeIds.begin(),_boundaryNodeIds.end(),nodesID[inode]);
336 if( it != _boundaryNodeIds.end() )
337 _boundaryNodeIds.push_back(nodesID[inode]);
343 Mesh::getIndexFacePeriodic( void ) const
345 return _indexFacePeriodicMap;
349 Mesh::setPeriodicFaces(bool check_groups, bool use_central_inversion)
351 if(_indexFacePeriodicSet)
356 for (int indexFace=0;indexFace<_boundaryFaceIds.size() ; indexFace++)
358 Face my_face=_faces[_boundaryFaceIds[indexFace]];
362 for (int iface=0;iface<_boundaryFaceIds.size() ; iface++)
365 iface_perio=_boundaryFaceIds[iface];
371 double x=my_face.x();
372 double y=my_face.y();
374 for (int iface=0;iface<_boundaryFaceIds.size() ; iface++)
376 Face face_i=_faces[_boundaryFaceIds[iface]];
377 double xi=face_i.x();
378 double yi=face_i.y();
379 if ( (abs(y-yi)<eps || abs(x-xi)<eps )// Case of a square geometry
380 && ( !check_groups || my_face.getGroupName()!=face_i.getGroupName()) //In case groups need to be checked
381 && ( !use_central_inversion || abs(y+yi) + abs(x+xi)<eps ) // Case of a central inversion
382 && fabs(my_face.getMeasure() - face_i.getMeasure())<eps
383 && fabs(my_face.getXN() + face_i.getXN())<eps
384 && fabs(my_face.getYN() + face_i.getYN())<eps
385 && fabs(my_face.getZN() + face_i.getZN())<eps )
387 iface_perio=_boundaryFaceIds[iface];
394 double x=my_face.x();
395 double y=my_face.y();
396 double z=my_face.z();
398 for (int iface=0;iface<_boundaryFaceIds.size() ; iface++)
400 Face face_i=_faces[_boundaryFaceIds[iface]];
401 double xi=face_i.x();
402 double yi=face_i.y();
403 double zi=face_i.z();
404 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
405 && ( !check_groups || my_face.getGroupName()!=face_i.getGroupName()) //In case groups need to be checked
406 && ( !use_central_inversion || abs(y+yi) + abs(x+xi) + abs(z+zi)<eps )// Case of a central inversion
407 && fabs(my_face.getMeasure() - face_i.getMeasure())<eps
408 && fabs(my_face.getXN() + face_i.getXN())<eps
409 && fabs(my_face.getYN() + face_i.getYN())<eps
410 && fabs(my_face.getZN() + face_i.getZN())<eps )
412 iface_perio=_boundaryFaceIds[iface];
418 throw CdmathException("Mesh::setPeriodicFaces: Mesh dimension should be 1, 2 or 3");
421 throw CdmathException("Mesh::setPeriodicFaces: periodic face not found, iface_perio==-1 " );
423 _indexFacePeriodicMap[_boundaryFaceIds[indexFace]]=iface_perio;
425 _indexFacePeriodicSet=true;
429 Mesh::getIndexFacePeriodic(int indexFace, bool check_groups, bool use_central_inversion)
431 if (!_faces[indexFace].isBorder())
433 cout<<"Pb with indexFace= "<<indexFace<<endl;
434 throw CdmathException("Mesh::getIndexFacePeriodic: not a border face" );
437 if(!_indexFacePeriodicSet)
438 setPeriodicFaces(check_groups, use_central_inversion);
440 std::map<int,int>::const_iterator it = _indexFacePeriodicMap.find(indexFace);
441 if( it != _indexFacePeriodicMap.end() )
445 cout<<"Pb with indexFace= "<<indexFace<<endl;
446 throw CdmathException("Mesh::getIndexFacePeriodic: not a periodic face" );
451 Mesh::isBorderNode(int nodeid) const
453 return _nodes[nodeid].isBorder();
457 Mesh::isBorderFace(int faceid) const
459 return _faces[faceid].isBorder();
463 Mesh::getBoundaryFaceIds() const
465 return _boundaryFaceIds;
469 Mesh::getBoundaryNodeIds() const
471 return _boundaryNodeIds;
475 Mesh::isTriangular() const
477 return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_TRI3;
480 Mesh::isTetrahedral() const
482 return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_TETRA4;
485 Mesh::isQuadrangular() const
487 return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_QUAD4;
490 Mesh::isHexahedral() const
492 return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_HEXA8;
495 Mesh::isStructured() const
497 return _isStructured;
500 std::vector< INTERP_KERNEL::NormalizedCellType >
501 Mesh::getElementTypes() const
506 std::vector< string >
507 Mesh::getElementTypesNames() const
509 std::vector< string > result(0);
510 for(int i=0; i< _eltsTypes.size(); i++)
512 if( _eltsTypes[i]==INTERP_KERNEL::NORM_POINT1)
513 result.push_back("Points");
514 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_SEG2)
515 result.push_back("Segments");
516 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_TRI3)
517 result.push_back("Triangles");
518 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_QUAD4)
519 result.push_back("Quadrangles");
520 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_POLYGON)
521 result.push_back("Polygons");
522 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_TETRA4)
523 result.push_back("Tetrahedra");
524 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_PYRA5)
525 result.push_back("Pyramids");
526 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_PENTA6)
527 result.push_back("Pentahedra");
528 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_HEXA8)
529 result.push_back("Hexahedra");
530 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_POLYHED)
531 result.push_back("Polyhedrons");
534 cout<< "Mesh " + _name + " contains an element of type " <<endl;
535 cout<< _eltsTypes[i]<<endl;
536 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");
543 Mesh::setGroups( const MEDFileUMesh* medmesh, MEDCouplingUMesh* mu)
545 //Searching for face groups
546 vector<string> faceGroups=medmesh->getGroupsNames() ;
548 for (unsigned int i=0;i<faceGroups.size();i++ )
550 string groupName=faceGroups[i];
551 vector<mcIdType> nonEmptyGrp(medmesh->getGrpNonEmptyLevels(groupName));
552 //We check if the group has a relative dimension equal to -1
553 //before call to the function getGroup(-1,groupName.c_str())
554 vector<mcIdType>::iterator it = find(nonEmptyGrp.begin(), nonEmptyGrp.end(), -1);
555 if (it != nonEmptyGrp.end())
557 cout<<"Boundary face group named "<< groupName << " found"<<endl;
558 MEDCouplingUMesh *m=medmesh->getGroup(-1,groupName.c_str());
560 _faceGroups.insert(_faceGroups.begin(),m);
561 _faceGroupNames.insert(_faceGroupNames.begin(),groupName);
562 DataArrayDouble *baryCell = m->computeIsoBarycenterOfNodesPerCell();//computeCellCenterOfMass() ;
563 const double *coorBary=baryCell->getConstPointer();
565 int nbCellsSubMesh=m->getNumberOfCells();
566 for (int ic(0), k(0); ic<nbCellsSubMesh; ic++, k+=_spaceDim)
568 vector<double> coorBaryXyz(3,0);
569 for (int d=0; d<_spaceDim; d++)
570 coorBaryXyz[d] = coorBary[k+d];
571 Point p1(coorBaryXyz[0],coorBaryXyz[1],coorBaryXyz[2]) ;
574 /* double min=INFINITY;
577 for (int iface=0;iface<_numberOfFaces;iface++ )
579 Point p2=_faces[iface].getBarryCenter();
580 /*if(groupName=="Wall" and p1.distance(p2)<min)
586 if(p1.distance(p2)<1.E-10)
588 _faces[iface].setGroupName(groupName);
593 /* if(groupName=="Wall" and min>1.E-10)
596 cout<<"Subcell number "<< ic <<" Min distance to Wall = "<<min <<" p= "<<p1[0]<<" , "<<p1[1]<<" , "<<p1[2]<<endl;
597 cout<<" attained at face "<< closestFaceNb << " p_min= "<<p3[0]<<" , "<<p3[1]<<" , "<<p3[2]<<endl;
600 throw CdmathException("No face belonging to group " + groupName + " found");
607 //Searching for node groups
608 vector<string> nodeGroups=medmesh->getGroupsOnSpecifiedLev(1) ;
610 for (unsigned int i=0;i<nodeGroups.size();i++ )
612 string groupName=nodeGroups[i];
613 DataArrayIdType * nodeGroup=medmesh->getNodeGroupArr( groupName );
614 const mcIdType *nodeids=nodeGroup->getConstPointer();
618 cout<<"Boundary node group named "<< groupName << " found"<<endl;
620 _nodeGroups.insert(_nodeGroups.begin(),nodeGroup);
621 _nodeGroupNames.insert(_nodeGroupNames.begin(),groupName);
623 int nbNodesSubMesh=nodeGroup->getNumberOfTuples();//nodeGroup->getNbOfElems();
625 DataArrayDouble *coo = mu->getCoords() ;
626 const double *cood=coo->getConstPointer();
628 for (int ic(0); ic<nbNodesSubMesh; ic++)
630 vector<double> coorP(3,0);
631 for (int d=0; d<_spaceDim; d++)
632 coorP[d] = cood[nodeids[ic]*_spaceDim+d];
633 Point p1(coorP[0],coorP[1],coorP[2]) ;
636 for (int inode=0;inode<_numberOfNodes;inode++ )
638 Point p2=_nodes[inode].getPoint();
639 if(p1.distance(p2)<1.E-10)
641 _nodes[inode].setGroupName(groupName);
647 throw CdmathException("No node belonging to group " + groupName + " found");
653 //----------------------------------------------------------------------
655 Mesh::setMesh( void )
656 //----------------------------------------------------------------------
658 DataArrayIdType *desc = DataArrayIdType::New();
659 DataArrayIdType *descI = DataArrayIdType::New();
660 DataArrayIdType *revDesc = DataArrayIdType::New();
661 DataArrayIdType *revDescI = DataArrayIdType::New();
662 MEDCouplingUMesh* mu = _mesh->buildUnstructured();
665 /* Save the boundary mesh for later use*/
666 _boundaryMesh = mu->computeSkin();
668 MEDCouplingUMesh* mu2=mu->buildDescendingConnectivity2(desc,descI,revDesc,revDescI);//mesh of dimension N-1 containing the cell interfaces
670 const mcIdType *tmp = desc->getConstPointer();
671 const mcIdType *tmpI=descI->getConstPointer();
673 const mcIdType *tmpA =revDesc->getConstPointer();
674 const mcIdType *tmpAI=revDescI->getConstPointer();
676 //const int *work=tmp+tmpI[id];//corresponds to buildDescendingConnectivity
678 //Test du type d'éléments contenus dans le maillage afin d'éviter les éléments contenant des points de gauss
679 _eltsTypes=mu->getAllGeoTypesSorted();
680 for(int i=0; i<_eltsTypes.size();i++)
683 _eltsTypes[i]!= INTERP_KERNEL::NORM_POINT1 && _eltsTypes[i]!= INTERP_KERNEL::NORM_SEG2
684 && _eltsTypes[i]!= INTERP_KERNEL::NORM_TRI3 && _eltsTypes[i]!= INTERP_KERNEL::NORM_QUAD4
685 && _eltsTypes[i]!= INTERP_KERNEL::NORM_TETRA4 && _eltsTypes[i]!= INTERP_KERNEL::NORM_PYRA5
686 && _eltsTypes[i]!= INTERP_KERNEL::NORM_PENTA6 && _eltsTypes[i]!= INTERP_KERNEL::NORM_HEXA8
687 && _eltsTypes[i]!= INTERP_KERNEL::NORM_POLYGON&& _eltsTypes[i]!= INTERP_KERNEL::NORM_POLYHED
690 cout<< "Mesh " + mu->getName() + " contains an element of type " <<endl;
691 cout<< _eltsTypes[i]<<endl;
692 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");
696 DataArrayDouble *baryCell = mu->computeCellCenterOfMass() ;
697 const double *coorBary=baryCell->getConstPointer();
699 MEDCouplingFieldDouble* fields=mu->getMeasureField(true);
700 DataArrayDouble *surface = fields->getArray();
701 const double *surf=surface->getConstPointer();
703 DataArrayDouble *coo = mu->getCoords() ;
704 const double *cood=coo->getConstPointer();
706 DataArrayIdType *revNode =DataArrayIdType::New();
707 DataArrayIdType *revNodeI=DataArrayIdType::New();
708 mu->getReverseNodalConnectivity(revNode,revNodeI) ;
709 const mcIdType *tmpN =revNode->getConstPointer();
710 const mcIdType *tmpNI=revNodeI->getConstPointer();
712 DataArrayIdType *revCell =DataArrayIdType::New();
713 DataArrayIdType *revCellI=DataArrayIdType::New();
714 mu2->getReverseNodalConnectivity(revCell,revCellI) ;
715 const mcIdType *tmpC =revCell->getConstPointer();
716 const mcIdType *tmpCI=revCellI->getConstPointer();
718 const DataArrayIdType *nodal = mu2->getNodalConnectivity() ;
719 const DataArrayIdType *nodalI = mu2->getNodalConnectivityIndex() ;
720 const mcIdType *tmpNE =nodal->getConstPointer();
721 const mcIdType *tmpNEI=nodalI->getConstPointer();
723 _numberOfCells = mu->getNumberOfCells() ;
724 _cells = new Cell[_numberOfCells] ;
726 _numberOfNodes = mu->getNumberOfNodes() ;
727 _nodes = new Node[_numberOfNodes] ;
729 _numberOfFaces = mu2->getNumberOfCells();
730 _faces = new Face[_numberOfFaces] ;
732 _indexFacePeriodicSet=false;
734 //Definition used if _meshDim =3 to determine the edges
735 DataArrayIdType *desc2 =DataArrayIdType::New();
736 DataArrayIdType *descI2=DataArrayIdType::New();
737 DataArrayIdType *revDesc2 =DataArrayIdType::New();
738 DataArrayIdType *revDescI2=DataArrayIdType::New();
739 DataArrayIdType *revNode2 =DataArrayIdType::New();
740 DataArrayIdType *revNodeI2=DataArrayIdType::New();
741 const mcIdType *tmpN2 ;
742 const mcIdType *tmpNI2;
743 MEDCouplingUMesh* mu3;
746 _numberOfEdges = mu->getNumberOfCells();
747 else if (_meshDim == 2)
748 _numberOfEdges = mu2->getNumberOfCells();
751 mu3=mu2->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2);//1D mesh of segments
752 _numberOfEdges = mu3->getNumberOfCells();
753 mu3->getReverseNodalConnectivity(revNode2,revNodeI2) ;
754 tmpN2 =revNode2->getConstPointer();
755 tmpNI2=revNodeI2->getConstPointer();
758 // _cells, _nodes and _faces initialization:
761 for( int id=0;id<_numberOfCells;id++ )
763 const mcIdType *work=tmp+tmpI[id];
764 int nbFaces=tmpI[id+1]-tmpI[id];
766 int nbVertices=mu->getNumberOfNodesInCell(id) ;
768 Cell ci( nbVertices, nbFaces, surf[id], Point(coorBary[id], 0.0, 0.0) ) ;
770 std::vector<mcIdType> nodeIdsOfCell ;
771 mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
772 for( int el=0;el<nbVertices;el++ )
774 ci.addNodeId(el,nodeIdsOfCell[el]) ;
775 ci.addFaceId(el,nodeIdsOfCell[el]) ;
778 //This assumes that in a 1D mesh the cell numbers form a consecutive sequence
779 //going either from left to right (xn=-1) or right to left (xn=1)
780 double xn = (cood[nodeIdsOfCell[nbVertices-1]] - cood[nodeIdsOfCell[0]] > 0.0) ? -1.0 : 1.0;
782 for( int el=0;el<nbFaces;el++ )
784 ci.addNormalVector(el,xn,0.0,0.0) ;
785 int indexFace=abs(work[el])-1;
786 ci.addFaceId(el,indexFace) ;
792 for( int id(0), k(0); id<_numberOfNodes; id++, k+=_spaceDim)
794 Point p(cood[k], 0.0, 0.0) ;
795 const mcIdType *workc=tmpN+tmpNI[id];
796 mcIdType nbCells=tmpNI[id+1]-tmpNI[id];
798 const mcIdType *workf=tmpC+tmpCI[id];
799 mcIdType nbFaces=tmpCI[id+1]-tmpCI[id];
800 const mcIdType *workn=tmpA+tmpAI[id];
801 mcIdType nbNeighbourNodes=tmpAI[id+1]-tmpAI[id];
802 Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
804 for( int el=0;el<nbCells;el++ )
805 vi.addCellId(el,workc[el]) ;
806 for( int el=0;el<nbFaces;el++ )
807 vi.addFaceId(el,workf[el]) ;
808 for( int el=0;el<nbNeighbourNodes;el++ )
809 vi.addNeighbourNodeId(el,workn[el]) ;
812 _boundaryNodeIds.push_back(0);
813 _boundaryNodeIds.push_back(_numberOfNodes-1);
815 for(int id(0), k(0); id<_numberOfFaces; id++, k+=_spaceDim)
817 Point p(cood[k], 0.0, 0.0) ;
818 const mcIdType *workc=tmpA+tmpAI[id];
819 mcIdType nbCells=tmpAI[id+1]-tmpAI[id];
821 const mcIdType *workv=tmpNE+tmpNEI[id]+1;
822 Face fi( 1, nbCells, 1.0, p, 1.0, 0.0, 0.0) ;
823 fi.addNodeId(0,workv[0]) ;
825 for(int idCell=0; idCell<nbCells; idCell++)
826 fi.addCellId(idCell,workc[idCell]) ;
830 _boundaryFaceIds.push_back(0);
831 _boundaryFaceIds.push_back(_numberOfFaces-1);
833 else if(_spaceDim==2 || _spaceDim==3)
835 DataArrayDouble *barySeg = mu2->computeIsoBarycenterOfNodesPerCell();//computeCellCenterOfMass() ;
836 const double *coorBarySeg=barySeg->getConstPointer();
838 MEDCouplingFieldDouble* fieldn;
839 DataArrayDouble *normal;
840 const double *tmpNormal;
842 if(_spaceDim==_meshDim)
843 fieldn = mu2->buildOrthogonalField();//Compute the normal to each cell interface
845 fieldn = mu->buildOrthogonalField();//compute the 3D normal vector to the 2D cell
847 normal = fieldn->getArray();
848 tmpNormal = normal->getConstPointer();
850 /*Building mesh cells */
851 for(int id(0), k(0); id<_numberOfCells; id++, k+=_spaceDim)
853 const mcIdType *work=tmp+tmpI[id];
854 mcIdType nbFaces=tmpI[id+1]-tmpI[id];
856 mcIdType nbVertices=mu->getNumberOfNodesInCell(id) ;
858 vector<double> coorBaryXyz(3,0);
859 for (int d=0; d<_spaceDim; d++)
860 coorBaryXyz[d] = coorBary[k+d];
862 Point p(coorBaryXyz[0],coorBaryXyz[1],coorBaryXyz[2]) ;
863 Cell ci( nbVertices, nbFaces, surf[id], p ) ;
865 /* Filling cell nodes */
866 std::vector<mcIdType> nodeIdsOfCell ;
867 mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
868 for( int el=0;el<nbVertices;el++ )
869 ci.addNodeId(el,nodeIdsOfCell[el]) ;
871 /* Filling cell faces */
872 if(_spaceDim==_meshDim)//use the normal field generated by buildOrthogonalField()
873 for( int el=0;el<nbFaces;el++ )
875 mcIdType faceIndex=(abs(work[el])-1);//=work[el] since Fortran type numbering was used, and negative sign means anticlockwise numbering
876 vector<double> xyzn(3,0);//Outer normal to the cell
878 for (int d=0; d<_spaceDim; d++)
879 xyzn[d] = -tmpNormal[_spaceDim*faceIndex+d];
881 for (int d=0; d<_spaceDim; d++)
882 xyzn[d] = +tmpNormal[_spaceDim*faceIndex+d];
883 ci.addNormalVector(el,xyzn[0],xyzn[1],xyzn[2]) ;
884 ci.addFaceId(el,faceIndex) ;
886 else//build normals associated to the couple (cell id, face el)
888 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
889 {//work[0]= first face global number, work[1]= second face global number
890 mcIdType indexFace0=abs(work[0])-1;//=work[0] since Fortran type numbering was used, and negative sign means anticlockwise numbering
891 mcIdType indexFace1=abs(work[1])-1;//=work[1] since Fortran type numbering was used, and negative sign means anticlockwise numbering
892 mcIdType idNodeA=(tmpNE+tmpNEI[indexFace0]+1)[0];//global number of the first face node work[0]=(abs(work[0])-1)
893 mcIdType idNodeB=(tmpNE+tmpNEI[indexFace1]+1)[0];//global number of the second face node work[1]=(abs(work[1])-1)
895 for(int i=0;i<_spaceDim;i++)
896 vecAB[i]=coo->getIJ(idNodeB,i) - coo->getIJ(idNodeA,i);
898 ci.addNormalVector(0,-vecAB[0],-vecAB[1],-vecAB[2]) ;
899 ci.addNormalVector(1,vecAB[0],vecAB[1],vecAB[2]) ;
900 ci.addFaceId(0,indexFace0) ;
901 ci.addFaceId(1,indexFace1) ;
903 else//_meshDim==2, number of faces around the cell id is variable, each face is composed of two nodes
906 for (int d=0; d<_spaceDim; d++)
907 xyzn[d] = tmpNormal[_spaceDim*id+d];
908 for( int el=0;el<nbFaces;el++ )
910 int faceIndex=(abs(work[el])-1);//=work[el] since Fortran type numbering was used, and negative sign means anticlockwise numbering
911 const mcIdType *workv=tmpNE+tmpNEI[faceIndex]+1;
912 mcIdType nbNodes= tmpNEI[faceIndex+1]-tmpNEI[faceIndex]-1;
913 if(nbNodes!=2)//We want to compute the normal to a straight line, not a curved interface composed of more thant 2 points
915 cout<<"Mesh name "<< mu->getName()<< " space dim= "<< _spaceDim <<" mesh dim= "<< _meshDim <<endl;
916 cout<<"For cell id "<<id<<" and local face number "<<el<<", the number of nodes is "<< nbNodes<< ", total number of faces is "<< nbFaces <<endl;
917 throw CdmathException("Mesh::setMesh number of nodes around a face should be 2");
920 mcIdType idNodeA=workv[0];
921 mcIdType idNodeB=workv[1];
922 vector<double> nodeA(_spaceDim), nodeB(_spaceDim), nodeP(_spaceDim);
923 for(int i=0;i<_spaceDim;i++)
925 nodeA[i]=coo->getIJ(idNodeA,i);
926 nodeB[i]=coo->getIJ(idNodeB,i);
927 nodeP[i]=coorBary[_spaceDim*id+i];
929 //Let P be the barycenter of the cell id
930 Vector vecAB(3), vecPA(3);
931 for(int i=0;i<_spaceDim;i++)
933 vecAB[i]=coo->getIJ(idNodeB,i) - coo->getIJ(idNodeA,i);
934 vecPA[i]=coo->getIJ(idNodeA,i) - coorBary[_spaceDim*id+i];
937 Vector normale = xyzn % vecAB;//Normal to the edge
938 normale/=normale.norm();
941 ci.addNormalVector(el,normale[0],normale[1],normale[2]) ;
943 ci.addNormalVector(el,-normale[0],-normale[1],-normale[2]) ;
944 ci.addFaceId(el,faceIndex) ;
951 MEDCouplingFieldDouble* fieldl=mu2->getMeasureField(true);
952 DataArrayDouble *longueur = fieldl->getArray();
953 const double *lon=longueur->getConstPointer();
955 if(_spaceDim!=_meshDim)
957 /* Since spaceDim!=meshDim, don't build normal to faces */
963 /*Building mesh faces */
964 for(int id(0), k(0); id<_numberOfFaces; id++, k+=_spaceDim)
966 vector<double> coorBarySegXyz(3,0);
967 for (int d=0; d<_spaceDim; d++)
968 coorBarySegXyz[d] = coorBarySeg[k+d];
969 Point p(coorBarySegXyz[0],coorBarySegXyz[1],coorBarySegXyz[2]) ;
970 const mcIdType *workc=tmpA+tmpAI[id];
971 mcIdType nbCells=tmpAI[id+1]-tmpAI[id];
973 if (nbCells>2 && _spaceDim==_meshDim)
975 cout<<"Warning : nbCells>2, numberOfFaces="<<_numberOfFaces<<endl;
976 cout<<"nbCells= "<<nbCells<<", _spaceDim="<<_spaceDim<<", _meshDim="<<_meshDim<<endl;
977 for(int icell=0; icell<nbCells; icell++)
978 cout<<workc[icell]<<", ";
980 throw CdmathException("Wrong mesh : nbCells>2 and spaceDim==meshDim");
983 _boundaryFaceIds.push_back(id);
985 const mcIdType *workv=tmpNE+tmpNEI[id]+1;
986 mcIdType nbNodes= tmpNEI[id+1]-tmpNEI[id]-1;
989 if(_spaceDim==_meshDim)//Euclidean flat mesh geometry
991 fi=Face( nbNodes, nbCells, lon[id], p, tmpNormal[k], tmpNormal[k+1], 0.0) ;
993 fi=Face( nbNodes, nbCells, lon[id], p, tmpNormal[k], tmpNormal[k+1], tmpNormal[k+2]) ;
994 else//Curved mesh geometry
995 fi=Face( nbNodes, nbCells, lon[id], p, 0.0, 0.0, 0.0) ;//Since spaceDim!=meshDim, normal to face is not defined
997 for(int node_id=0; node_id<nbNodes;node_id++)
998 fi.addNodeId(node_id,workv[node_id]) ;
1000 fi.addCellId(0,workc[0]) ;
1001 for(int cell_id=1; cell_id<nbCells;cell_id++)
1004 if (workc[cell_id]!=workc[cell_id-1])//For some meshes (bad ones) the same cell can appear several times
1006 fi.addCellId(cell_idx+1,workc[cell_id]) ;
1014 /*Building mesh nodes, should be done after building mesh faces in order to detect boundary nodes*/
1015 for(int id(0), k(0); id<_numberOfNodes; id++, k+=_spaceDim)
1017 vector<double> coorP(3,0);
1018 for (int d=0; d<_spaceDim; d++)
1019 coorP[d] = cood[k+d];
1020 Point p(coorP[0],coorP[1],coorP[2]) ;
1022 const mcIdType *workc=tmpN+tmpNI[id];
1023 mcIdType nbCells=tmpNI[id+1]-tmpNI[id];
1024 const mcIdType *workf=tmpC+tmpCI[id];
1025 mcIdType nbFaces=tmpCI[id+1]-tmpCI[id];
1026 const mcIdType *workn;
1027 mcIdType nbNeighbourNodes;
1030 workn=tmpA+tmpAI[id];
1031 nbNeighbourNodes=tmpAI[id+1]-tmpAI[id];
1033 else if (_meshDim == 2)
1035 workn=tmpC+tmpCI[id];
1036 nbNeighbourNodes=tmpCI[id+1]-tmpCI[id];
1040 workn=tmpN2+tmpNI2[id];
1041 nbNeighbourNodes=tmpNI2[id+1]-tmpNI2[id];
1043 Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
1045 for( int el=0;el<nbCells;el++ )
1046 vi.addCellId(el,workc[el]) ;
1047 for( int el=0;el<nbNeighbourNodes;el++ )
1048 vi.addNeighbourNodeId(el,workn[el]) ;
1049 //Detection of border nodes
1050 bool isBorder=false;
1051 for( int el=0;el<nbFaces;el++ )
1053 vi.addFaceId(el,workf[el],_faces[workf[el]].isBorder()) ;
1054 isBorder= isBorder || _faces[workf[el]].isBorder();
1057 _boundaryNodeIds.push_back(id);
1061 if(_spaceDim==_meshDim)
1067 throw CdmathException("Mesh::setMesh space dimension should be 1, 2 or 3");
1069 //definition of the bounding box for unstructured meshes
1070 if(!_isStructured)//Structured case is treated in function readMeshMed
1072 double Box0[2*_spaceDim];
1073 coo->getMinMaxPerComponent(Box0);
1089 //Set boundary groups
1090 _faceGroupNames.push_back("Boundary");
1091 _nodeGroupNames.push_back("Boundary");
1096 revDescI->decrRef();
1098 baryCell->decrRef();
1101 revNodeI->decrRef();
1103 revCellI->decrRef();
1107 revNode2->decrRef();
1108 revNodeI2->decrRef();
1111 revDesc2->decrRef();
1112 revDescI2->decrRef();
1119 //----------------------------------------------------------------------
1120 Mesh::Mesh( double xmin, double xmax, int nx, std::string meshName )
1121 //----------------------------------------------------------------------
1124 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx) : nx <= 0");
1126 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx) : xmin >= xmax");
1128 double dx = (xmax - xmin)/nx ;
1133 _isStructured = true;
1134 _indexFacePeriodicSet=false;
1143 _dxyz.resize(_spaceDim);
1145 _nxyz.resize(_spaceDim);
1148 double *originPtr = new double[_spaceDim];
1149 double *dxyzPtr = new double[_spaceDim];
1150 mcIdType *nodeStrctPtr = new mcIdType[_spaceDim];
1153 nodeStrctPtr[0]=nx+1;
1156 _mesh=MEDCouplingIMesh::New(meshName,
1159 nodeStrctPtr+_spaceDim,
1161 originPtr+_spaceDim,
1164 delete [] originPtr;
1166 delete [] nodeStrctPtr;
1168 DataArrayIdType *desc=DataArrayIdType::New();
1169 DataArrayIdType *descI=DataArrayIdType::New();
1170 DataArrayIdType *revDesc=DataArrayIdType::New();
1171 DataArrayIdType *revDescI=DataArrayIdType::New();
1172 MEDCouplingUMesh* mu=_mesh->buildUnstructured();
1173 MEDCouplingUMesh *mu2=mu->buildDescendingConnectivity(desc,descI,revDesc,revDescI);
1175 const mcIdType *tmp=desc->getConstPointer();
1176 const mcIdType *tmpI=descI->getConstPointer();
1178 const mcIdType *tmpA =revDesc->getConstPointer();
1179 const mcIdType *tmpAI=revDescI->getConstPointer();
1181 _eltsTypes=mu->getAllGeoTypesSorted();
1183 DataArrayDouble *baryCell = mu->computeCellCenterOfMass() ;
1184 const double *coorBary=baryCell->getConstPointer();
1186 _numberOfCells = _mesh->getNumberOfCells() ;
1187 _cells = new Cell[_numberOfCells] ;
1189 _numberOfNodes = mu->getNumberOfNodes() ;
1190 _nodes = new Node[_numberOfNodes] ;
1192 _numberOfFaces = _numberOfNodes;
1193 _faces = new Face[_numberOfFaces] ;
1195 _numberOfEdges = _numberOfCells;
1197 MEDCouplingFieldDouble* fieldl=mu->getMeasureField(true);
1198 DataArrayDouble *longueur = fieldl->getArray();
1199 const double *lon=longueur->getConstPointer();
1201 DataArrayDouble *coo = mu->getCoords() ;
1202 const double *cood=coo->getConstPointer();
1205 for( int id=0;id<_numberOfCells;id++ )
1207 int nbVertices=mu->getNumberOfNodesInCell(id) ;
1208 Point p(coorBary[id],0.0,0.0) ;
1209 Cell ci( nbVertices, nbVertices, lon[id], p ) ;
1211 std::vector<mcIdType> nodeIdsOfCell ;
1212 mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
1213 for( int el=0;el<nbVertices;el++ )
1215 ci.addNodeId(el,nodeIdsOfCell[el]) ;
1216 ci.addFaceId(el,nodeIdsOfCell[el]) ;
1219 double xn = (cood[nodeIdsOfCell[nbVertices-1]] - cood[nodeIdsOfCell[0]] > 0.0) ? -1.0 : 1.0;
1221 mcIdType nbFaces=tmpI[id+1]-tmpI[id];
1222 const mcIdType *work=tmp+tmpI[id];
1224 for( int el=0;el<nbFaces;el++ )
1226 ci.addNormalVector(el,xn,0.0,0.0) ;
1227 ci.addFaceId(el,work[el]) ;
1236 //Suppress the following since tmpN=tmpA
1237 DataArrayIdType *revNode=DataArrayIdType::New();
1238 DataArrayIdType *revNodeI=DataArrayIdType::New();
1239 mu->getReverseNodalConnectivity(revNode,revNodeI) ;
1240 const mcIdType *tmpN=revNode->getConstPointer();
1241 const mcIdType *tmpNI=revNodeI->getConstPointer();
1243 for( int id=0;id<_numberOfNodes;id++ )
1245 std::vector<double> coo ;
1246 mu->getCoordinatesOfNode(id,coo);
1247 Point p(coo[0],0.0,0.0) ;
1248 const mcIdType *workc=tmpN+tmpNI[id];
1249 mcIdType nbCells=tmpNI[id+1]-tmpNI[id];
1251 const mcIdType *workn=tmpA+tmpAI[id];
1252 mcIdType nbNeighbourNodes=tmpAI[id+1]-tmpAI[id];
1254 Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
1255 for( int el=0;el<nbCells;el++ )
1256 vi.addCellId(el,workc[el]) ;
1257 for( int el=0;el<nbFaces;el++ )
1258 vi.addFaceId(el,id) ;
1259 for( int el=0;el<nbNeighbourNodes;el++ )
1260 vi.addNeighbourNodeId(el,workn[el]) ;
1265 Face fi( nbVertices, nbCells, 1.0, p, 1., 0., 0. ) ;
1266 for( int el=0;el<nbVertices;el++ )
1267 fi.addNodeId(el,id) ;
1269 for( int el=0;el<nbCells;el++ )
1270 fi.addCellId(el,workc[el]) ;
1274 //Set boundary groups
1275 _faceGroupNames.push_back("Boundary");
1276 _nodeGroupNames.push_back("Boundary");
1277 _boundaryFaceIds.push_back(0);
1278 _boundaryFaceIds.push_back(_numberOfFaces-1);
1279 _boundaryNodeIds.push_back(0);
1280 _boundaryNodeIds.push_back(_numberOfNodes-1);
1283 baryCell->decrRef();
1287 revDescI->decrRef();
1289 revNodeI->decrRef();
1294 //----------------------------------------------------------------------
1295 Mesh::Mesh( std::vector<double> points, std::string meshName )
1296 //----------------------------------------------------------------------
1298 int nx=points.size();
1301 throw CdmathException("Mesh::Mesh( vector<double> points, string meshName) : nx < 2, vector should contain at least two values");
1303 while( i<nx-1 && points[i+1]>points[i] )
1307 //cout<< points << endl;
1308 throw CdmathException("Mesh::Mesh( vector<double> points, string meshName) : vector values should be sorted");
1321 _isStructured = false;
1322 _indexFacePeriodicSet=false;
1324 MEDCouplingUMesh * mesh1d = MEDCouplingUMesh::New(meshName, 1);
1325 mesh1d->allocateCells(nx - 1);
1326 double * coords = new double[nx];
1327 mcIdType * nodal_con = new mcIdType[2];
1328 coords[0]=points[0];
1329 for(int i=0; i<nx- 1 ; i++)
1333 mesh1d->insertNextCell(INTERP_KERNEL::NORM_SEG2, 2, nodal_con);
1334 coords[i+1]=points[i + 1];
1336 mesh1d->finishInsertingCells();
1338 DataArrayDouble * coords_arr = DataArrayDouble::New();
1339 coords_arr->alloc(nx,1);
1340 std::copy(coords,coords+nx,coords_arr->getPointer());
1341 mesh1d->setCoords(coords_arr);
1343 delete [] coords, nodal_con;
1344 coords_arr->decrRef();
1348 DataArrayIdType *desc=DataArrayIdType::New();
1349 DataArrayIdType *descI=DataArrayIdType::New();
1350 DataArrayIdType *revDesc=DataArrayIdType::New();
1351 DataArrayIdType *revDescI=DataArrayIdType::New();
1352 MEDCouplingUMesh* mu=_mesh->buildUnstructured();
1353 MEDCouplingUMesh *mu2=mu->buildDescendingConnectivity(desc,descI,revDesc,revDescI);
1355 DataArrayDouble *baryCell = mu->computeCellCenterOfMass() ;
1356 const double *coorBary=baryCell->getConstPointer();
1358 _numberOfCells = _mesh->getNumberOfCells() ;
1359 _cells = new Cell[_numberOfCells] ;
1361 _numberOfEdges = _numberOfCells;
1363 _eltsTypes=mu->getAllGeoTypesSorted();
1365 MEDCouplingFieldDouble* fieldl=mu->getMeasureField(true);
1366 DataArrayDouble *longueur = fieldl->getArray();
1367 const double *lon=longueur->getConstPointer();
1370 for( int id=0;id<_numberOfCells;id++ )
1372 int nbVertices=mu->getNumberOfNodesInCell(id) ;
1373 Point p(coorBary[id],0.0,0.0) ;
1374 Cell ci( nbVertices, nbVertices, lon[id], p ) ;
1376 std::vector<mcIdType> nodeIdsOfCell ;
1377 mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
1378 for( int el=0;el<nbVertices;el++ )
1380 ci.addNodeId(el,nodeIdsOfCell[el]) ;
1381 ci.addFaceId(el,nodeIdsOfCell[el]) ;
1388 //Suppress the following since tmpN=tmpA
1389 DataArrayIdType *revNode=DataArrayIdType::New();
1390 DataArrayIdType *revNodeI=DataArrayIdType::New();
1391 mu->getReverseNodalConnectivity(revNode,revNodeI) ;
1392 const mcIdType *tmpN=revNode->getConstPointer();
1393 const mcIdType *tmpNI=revNodeI->getConstPointer();
1395 _numberOfNodes = mu->getNumberOfNodes() ;
1396 _nodes = new Node[_numberOfNodes] ;
1397 _numberOfFaces = _numberOfNodes;
1398 _faces = new Face[_numberOfFaces] ;
1400 for( int id=0;id<_numberOfNodes;id++ )
1402 std::vector<double> coo ;
1403 mu->getCoordinatesOfNode(id,coo);
1404 Point p(coo[0],0.0,0.0) ;
1405 const mcIdType *workc=tmpN+tmpNI[id];
1406 mcIdType nbCells=tmpNI[id+1]-tmpNI[id];
1408 const mcIdType *workn=tmpN+tmpNI[id];
1409 mcIdType nbNeighbourNodes=tmpNI[id+1]-tmpNI[id];
1410 Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
1412 /* provisoire !!!!!!!!!!!!*/
1413 // Point pf(0.0,0.0,0.0) ;
1414 Face fi( nbVertices, nbCells, 0.0, p, 0., 0., 0. ) ;
1416 for( int el=0;el<nbCells;el++ )
1417 vi.addCellId(el,workc[el]) ;
1418 for( int el=0;el<nbFaces;el++ )
1419 vi.addFaceId(el,id) ;
1420 for( int el=0;el<nbNeighbourNodes;el++ )
1421 vi.addNeighbourNodeId(el,workn[el]) ;
1424 for( int el=0;el<nbVertices;el++ )
1425 fi.addNodeId(el,id) ;
1427 for( int el=0;el<nbCells;el++ )
1428 fi.addCellId(el,workc[el]) ;
1431 //Set boundary groups
1432 _faceGroupNames.push_back("Boundary");
1433 _nodeGroupNames.push_back("Boundary");
1434 _boundaryFaceIds.push_back(0);
1435 _boundaryFaceIds.push_back(_numberOfFaces-1);
1436 _boundaryNodeIds.push_back(0);
1437 _boundaryNodeIds.push_back(_numberOfNodes-1);
1440 baryCell->decrRef();
1444 revDescI->decrRef();
1446 revNodeI->decrRef();
1450 //----------------------------------------------------------------------
1451 Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, int split_to_triangles_policy, std::string meshName)
1452 //----------------------------------------------------------------------
1455 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : nx <= 0 or ny <= 0");
1457 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : xmin >= xmax");
1459 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : ymin >= ymax");
1469 double dx = (xmax - xmin)/nx ;
1470 double dy = (ymax - ymin)/ny ;
1475 _indexFacePeriodicSet=false;
1476 _isStructured = true;
1477 _nxyz.resize(_spaceDim);
1481 _dxyz.resize(_spaceDim);
1485 double *originPtr = new double[_spaceDim];
1486 double *dxyzPtr = new double[_spaceDim];
1487 mcIdType *nodeStrctPtr = new mcIdType[_spaceDim];
1491 nodeStrctPtr[0]=nx+1;
1492 nodeStrctPtr[1]=ny+1;
1496 _mesh=MEDCouplingIMesh::New(meshName,
1499 nodeStrctPtr+_spaceDim,
1501 originPtr+_spaceDim,
1505 if(split_to_triangles_policy==0 || split_to_triangles_policy==1)
1507 _mesh=_mesh->buildUnstructured();
1508 _mesh->simplexize(split_to_triangles_policy);
1510 else if (split_to_triangles_policy != -1)
1512 cout<< "split_to_triangles_policy = "<< split_to_triangles_policy << endl;
1513 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : Unknown splitting policy");
1516 delete [] originPtr;
1518 delete [] nodeStrctPtr;
1523 //----------------------------------------------------------------------
1524 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)
1525 //----------------------------------------------------------------------
1527 if(nx<=0 || ny<=0 || nz<=0)
1528 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");
1530 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : xmin >= xmax");
1532 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : ymin >= ymax");
1534 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : zmin >= zmax");
1539 _indexFacePeriodicSet=false;
1547 double dx = (xmax - xmin)/nx ;
1548 double dy = (ymax - ymin)/ny ;
1549 double dz = (zmax - zmin)/nz ;
1551 _isStructured = true;
1552 _dxyz.resize(_spaceDim);
1557 _nxyz.resize(_spaceDim);
1562 double *originPtr = new double[_spaceDim];
1563 double *dxyzPtr = new double[_spaceDim];
1564 mcIdType *nodeStrctPtr = new mcIdType[_spaceDim];
1569 nodeStrctPtr[0]=nx+1;
1570 nodeStrctPtr[1]=ny+1;
1571 nodeStrctPtr[2]=nz+1;
1576 _mesh=MEDCouplingIMesh::New(meshName,
1579 nodeStrctPtr+_spaceDim,
1581 originPtr+_spaceDim,
1585 if( split_to_tetrahedra_policy == 0 )
1587 _mesh=_mesh->buildUnstructured();
1588 _mesh->simplexize(INTERP_KERNEL::PLANAR_FACE_5);
1590 else if( split_to_tetrahedra_policy == 1 )
1592 _mesh=_mesh->buildUnstructured();
1593 _mesh->simplexize(INTERP_KERNEL::PLANAR_FACE_6);
1595 else if ( split_to_tetrahedra_policy != -1 )
1597 cout<< "split_to_tetrahedra_policy = "<< split_to_tetrahedra_policy << endl;
1598 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");
1601 delete [] originPtr;
1603 delete [] nodeStrctPtr;
1608 //----------------------------------------------------------------------
1610 Mesh::getSpaceDimension( void ) const
1611 //----------------------------------------------------------------------
1616 //----------------------------------------------------------------------
1618 Mesh::getMeshDimension( void ) const
1619 //----------------------------------------------------------------------
1625 Mesh::getDXYZ() const
1628 throw CdmathException("std::vector<double> Mesh::getDXYZ() : dx,dy and dz are defined only for structured meshes !");
1633 std::vector<mcIdType>
1634 Mesh::getCellGridStructure() const
1637 throw CdmathException("std::vector<int> Mesh::getCellGridStructure() : nx, ny and nz are defined only for structured meshes !");
1642 //----------------------------------------------------------------------
1644 Mesh::getNx( void ) const
1645 //----------------------------------------------------------------------
1648 throw CdmathException("int Mesh::getNx( void ) : Nx is defined only for structured meshes !");
1653 //----------------------------------------------------------------------
1655 Mesh::getNy( void ) const
1656 //----------------------------------------------------------------------
1659 throw CdmathException("int Mesh::getNy( void ) : Ny is defined only for structured meshes !");
1661 throw CdmathException("int double& Field::operator[ielem] : Ny is not defined in dimension < 2!");
1666 //----------------------------------------------------------------------
1668 Mesh::getNz( void ) const
1669 //----------------------------------------------------------------------
1672 throw CdmathException("int Mesh::getNz( void ) : Nz is defined only for structured meshes !");
1674 throw CdmathException("int Mesh::getNz( void ) : Nz is not defined in dimension < 3!");
1679 //----------------------------------------------------------------------
1681 Mesh::getXMin( void ) const
1682 //----------------------------------------------------------------------
1687 //----------------------------------------------------------------------
1689 Mesh::getXMax( void ) const
1690 //----------------------------------------------------------------------
1695 //----------------------------------------------------------------------
1697 Mesh::getYMin( void ) const
1698 //----------------------------------------------------------------------
1703 //----------------------------------------------------------------------
1705 Mesh::getYMax( void ) const
1706 //----------------------------------------------------------------------
1711 //----------------------------------------------------------------------
1713 Mesh::getZMin( void ) const
1714 //----------------------------------------------------------------------
1719 //----------------------------------------------------------------------
1721 Mesh::getZMax( void ) const
1722 //----------------------------------------------------------------------
1727 //----------------------------------------------------------------------
1728 MCAuto<MEDCouplingMesh>
1729 Mesh::getMEDCouplingMesh( void ) const
1730 //----------------------------------------------------------------------
1735 //----------------------------------------------------------------------
1737 Mesh::getNumberOfNodes ( void ) const
1738 //----------------------------------------------------------------------
1740 return _numberOfNodes ;
1743 //----------------------------------------------------------------------
1745 Mesh::getNumberOfCells ( void ) const
1746 //----------------------------------------------------------------------
1748 return _numberOfCells ;
1751 //----------------------------------------------------------------------
1753 Mesh::getNumberOfFaces ( void ) const
1754 //----------------------------------------------------------------------
1756 return _numberOfFaces ;
1759 //----------------------------------------------------------------------
1761 Mesh::getNumberOfEdges ( void ) const
1762 //----------------------------------------------------------------------
1764 return _numberOfEdges ;
1767 //----------------------------------------------------------------------
1769 Mesh::getFaces ( void ) const
1770 //----------------------------------------------------------------------
1775 //----------------------------------------------------------------------
1777 Mesh::getCells ( void ) const
1778 //----------------------------------------------------------------------
1783 //----------------------------------------------------------------------
1785 Mesh::getCell ( int i ) const
1786 //----------------------------------------------------------------------
1791 //----------------------------------------------------------------------
1793 Mesh::getFace ( int i ) const
1794 //----------------------------------------------------------------------
1799 //----------------------------------------------------------------------
1801 Mesh::getNode ( int i ) const
1802 //----------------------------------------------------------------------
1807 //----------------------------------------------------------------------
1809 Mesh::getNodes ( void ) const
1810 //----------------------------------------------------------------------
1816 Mesh::getNameOfFaceGroups( void ) const
1818 return _faceGroupNames;
1821 vector<MEDCoupling::MEDCouplingUMesh *>
1822 Mesh::getFaceGroups( void ) const
1828 Mesh::getNameOfNodeGroups( void ) const
1830 return _nodeGroupNames;
1833 vector<MEDCoupling::DataArrayIdType *>
1834 Mesh::getNodeGroups( void ) const
1839 //----------------------------------------------------------------------
1841 Mesh::operator= ( const Mesh& mesh )
1842 //----------------------------------------------------------------------
1844 _spaceDim = mesh.getSpaceDimension() ;
1845 _meshDim = mesh.getMeshDimension() ;
1846 _name = mesh.getName();
1847 _numberOfNodes = mesh.getNumberOfNodes();
1848 _numberOfFaces = mesh.getNumberOfFaces();
1849 _numberOfCells = mesh.getNumberOfCells();
1850 _numberOfEdges = mesh.getNumberOfEdges();
1852 _isStructured = mesh.isStructured();
1855 _nxyz = mesh.getCellGridStructure() ;
1856 _dxyz=mesh.getDXYZ();
1857 _xMin=mesh.getXMin();
1858 _xMax=mesh.getXMax();
1859 _yMin=mesh.getYMin();
1860 _yMax=mesh.getYMax();
1861 _zMin=mesh.getZMin();
1862 _zMax=mesh.getZMax();
1864 _faceGroupNames = mesh.getNameOfFaceGroups() ;
1865 _faceGroups = mesh.getFaceGroups() ;
1866 _nodeGroupNames = mesh.getNameOfNodeGroups() ;
1867 _nodeGroups = mesh.getNodeGroups() ;
1885 _nodes = new Node[_numberOfNodes] ;
1886 _faces = new Face[_numberOfFaces] ;
1887 _cells = new Cell[_numberOfCells] ;
1889 for (int i=0;i<_numberOfNodes;i++)
1890 _nodes[i]=mesh.getNode(i);
1892 for (int i=0;i<_numberOfFaces;i++)
1893 _faces[i]=mesh.getFace(i);
1895 for (int i=0;i<_numberOfCells;i++)
1896 _cells[i]=mesh.getCell(i);
1898 _indexFacePeriodicSet= mesh.isIndexFacePeriodicSet();
1899 if(_indexFacePeriodicSet)
1900 _indexFacePeriodicMap=mesh.getIndexFacePeriodic();
1902 _boundaryFaceIds=mesh.getBoundaryFaceIds();
1903 _boundaryNodeIds=mesh.getBoundaryNodeIds();
1905 _eltsTypes=mesh.getElementTypes();
1906 _eltsTypesNames=mesh.getElementTypesNames();
1908 MCAuto<MEDCouplingMesh> m1=mesh.getMEDCouplingMesh()->deepCopy();
1913 bool Mesh::isIndexFacePeriodicSet() const
1915 return _indexFacePeriodicSet;
1917 //----------------------------------------------------------------------
1919 Mesh::minRatioVolSurf()
1921 double dx_min = 1e30;
1922 for(int i=0; i<_numberOfCells; i++)
1924 Cell Ci = getCell(i);
1928 for(int k=0; k< Ci.getNumberOfFaces(); k++)
1930 int indexFace=Ci.getFacesId()[k];
1931 Face Fk = getFace(indexFace);
1932 perimeter+=Fk.getMeasure();
1934 dx_min = min(dx_min,Ci.getMeasure()/perimeter);
1937 dx_min = min(dx_min,Ci.getMeasure());
1944 Mesh::getBoundaryMesh ( void ) const
1946 return Mesh(_boundaryMesh);
1950 Mesh::getMaxNbNeighbours(EntityType type) const
1956 for(int i=0; i<_numberOfCells; i++)
1957 if(result < _cells[i].getNumberOfFaces())
1958 result=_cells[i].getNumberOfFaces();
1960 else if(type==NODES)
1962 for(int i=0; i<_numberOfNodes; i++)
1963 if(result < _nodes[i].getNumberOfEdges())
1964 result=_nodes[i].getNumberOfEdges();
1967 throw CdmathException("Mesh::getMaxNbNeighbours : entity type is not accepted. Should be CELLS or NODES");
1971 //----------------------------------------------------------------------
1973 Mesh::writeVTK ( const std::string fileName ) const
1974 //----------------------------------------------------------------------
1976 string fname=fileName+".vtu";
1977 _mesh->writeVTK(fname.c_str()) ;
1980 //----------------------------------------------------------------------
1982 Mesh::writeMED ( const std::string fileName ) const
1983 //----------------------------------------------------------------------
1985 string fname=fileName+".med";
1986 MEDCouplingUMesh* mu=_mesh->buildUnstructured();
1987 MEDCoupling::WriteUMesh(fname.c_str(),mu,true);
1989 //MEDFileUMesh meshMEDFile;
1990 //meshMEDFile.setMeshAtLevel(0,mu);
1991 //for(int i=0; i< _faceGroups.size(); i++)
1992 //meshMEDFile.setMeshAtLevel(-1,_faceGroups[i]);
1994 //MEDCoupling::meshMEDFile.write(fname.c_str(),2) ;
1996 //MEDCoupling::meshMEDFile.write(fname.c_str(),1) ;
2003 Mesh::getGroupFaceIds(std::string groupName) const
2005 if( std::find(_faceGroupNames.begin(),_faceGroupNames.end(),groupName) == _faceGroupNames.end() )
2007 cout<<"Mesh::getGroupFaceIds Error : face group " << groupName << " does not exist"<<endl;
2008 throw CdmathException("Required face group does not exist");
2012 std::vector< int > result(0);
2013 for(int i=0; i<_boundaryFaceIds.size(); i++)
2015 vector< std::string > groupNames = getFace(_boundaryFaceIds[i]).getGroupNames();
2016 if( std::find(groupNames.begin(), groupNames.end(),groupName) != groupNames.end() )
2017 result.push_back(_boundaryFaceIds[i]);
2024 Mesh::getGroupNodeIds(std::string groupName) const
2026 if( std::find(_nodeGroupNames.begin(),_nodeGroupNames.end(),groupName) == _nodeGroupNames.end() )
2028 cout<<"Mesh::getGroupNodeIds Error : node group " << groupName << " does not exist"<<endl;
2029 throw CdmathException("Required node group does not exist");
2033 std::vector< int > result(0);
2034 for(int i=0; i<_boundaryFaceIds.size(); i++)
2036 vector< std::string > groupNames = getNode(_boundaryFaceIds[i]).getGroupNames();
2037 if( std::find(groupNames.begin(), groupNames.end(),groupName) != groupNames.end() )
2038 result.push_back(_boundaryFaceIds[i]);
2045 Mesh::setFaceGroupByIds(std::vector< int > faceIds, std::string groupName)
2047 for(int i=0; i< faceIds.size(); i++)
2048 getFace(faceIds[i]).setGroupName(groupName);
2052 Mesh::setNodeGroupByIds(std::vector< int > nodeIds, std::string groupName)
2054 for(int i=0; i< nodeIds.size(); i++)
2055 getNode(nodeIds[i]).setGroupName(groupName);