4 * Created on: 22 janv. 2012
13 #include "CdmathException.hxx"
14 #include "MEDFileMesh.hxx"
15 #include "MEDLoader.hxx"
25 using namespace MEDCoupling;
28 //----------------------------------------------------------------------
30 //----------------------------------------------------------------------
33 _meshNotDeleted=false;
53 _boundaryFaceIds.resize(0);
54 _boundaryNodeIds.resize(0);
55 _faceGroupNames.resize(0);
56 _faceGroups.resize(0);
57 _faceGroupsIds.resize(0);
58 _nodeGroupNames.resize(0);
59 _nodeGroups.resize(0);
60 _nodeGroupsIds.resize(0);
61 _indexFacePeriodicSet=false;
66 //----------------------------------------------------------------------
68 //----------------------------------------------------------------------
74 //for(int i=0; i< _faceGroups.size(); i++)
75 // _faceGroups[i]->decrRef();
76 // _nodeGroups[i]->decrRef();
78 (_mesh.retn())->decrRef();
82 Mesh::getName( void ) const
87 Mesh::Mesh( const MEDCoupling::MEDCouplingIMesh* mesh )
89 _spaceDim=mesh->getSpaceDimension();
90 _meshDim=mesh->getMeshDimension();
91 _dxyz=mesh->getDXYZ();
92 _nxyz=mesh->getCellGridStructure();
93 double* Box0=new double[2*_spaceDim];
94 mesh->getBoundingBox(Box0);
95 _name=mesh->getName();
97 _indexFacePeriodicSet=false;
112 double *originPtr = new double[_spaceDim];
113 double *dxyzPtr = new double[_spaceDim];
114 mcIdType *nodeStrctPtr = new mcIdType[_spaceDim];
116 for(int i=0;i<_spaceDim;i++)
118 originPtr[i]=Box0[2*i];
119 nodeStrctPtr[i]=_nxyz[i]+1;
122 _mesh=MEDCouplingIMesh::New(_name,
125 nodeStrctPtr+_spaceDim,
130 _meshNotDeleted=true;
135 delete [] nodeStrctPtr;
141 Mesh::Mesh( const MEDCoupling::MEDCouplingUMesh* mesh )
143 _spaceDim=mesh->getSpaceDimension();
144 _meshDim=mesh->getMeshDimension();
145 double* Box0=new double[2*_spaceDim];
146 mesh->getBoundingBox(Box0);
147 _name=mesh->getName();
149 _indexFacePeriodicSet=false;
164 _mesh=mesh->deepCopy();
165 _mesh=mesh->buildUnstructured();
166 _meshNotDeleted=true;
173 //----------------------------------------------------------------------
174 Mesh::Mesh( const std::string filename, int meshLevel )
175 //----------------------------------------------------------------------
177 readMeshMed(filename, meshLevel);
180 //----------------------------------------------------------------------
181 Mesh::Mesh( const Mesh& mesh )
182 //----------------------------------------------------------------------
184 _spaceDim = mesh.getSpaceDimension() ;
185 _meshDim = mesh.getMeshDimension() ;
186 _name=mesh.getName();
187 _epsilon=mesh.getComparisonEpsilon();
188 _xMax=mesh.getXMax();
189 _yMin=mesh.getYMin();
190 _yMax=mesh.getYMax();
191 _zMin=mesh.getZMin();
192 _zMax=mesh.getZMax();
193 _isStructured=mesh.isStructured();
196 _nxyz = mesh.getCellGridStructure() ;
197 _dxyz=mesh.getDXYZ();
199 _numberOfNodes = mesh.getNumberOfNodes();
200 _numberOfFaces = mesh.getNumberOfFaces();
201 _numberOfCells = mesh.getNumberOfCells();
202 _numberOfEdges = mesh.getNumberOfEdges();
204 _faceGroupNames = mesh.getNameOfFaceGroups() ;
205 _faceGroups = mesh.getFaceGroups() ;
206 _nodeGroupNames = mesh.getNameOfNodeGroups() ;
207 _nodeGroups = mesh.getNodeGroups() ;
209 _nodes = new Node[_numberOfNodes] ;
210 _faces = new Face[_numberOfFaces] ;
211 _cells = new Cell[_numberOfCells] ;
213 //memcpy(_nodes,mesh.getNodes(),sizeof(*mesh.getNodes())) ;
214 //memcpy(_cells,mesh.getCells(),sizeof(*mesh.getCells())) ;
215 //memcpy(_faces,mesh.getFaces(),sizeof(*mesh.getFaces())) ;
217 for (int i=0;i<_numberOfNodes;i++)
218 _nodes[i]=mesh.getNode(i);
220 for (int i=0;i<_numberOfFaces;i++)
221 _faces[i]=mesh.getFace(i);
223 for (int i=0;i<_numberOfCells;i++)
224 _cells[i]=mesh.getCell(i);
226 _indexFacePeriodicSet= mesh.isIndexFacePeriodicSet();
227 if(_indexFacePeriodicSet)
228 _indexFacePeriodicMap=mesh.getIndexFacePeriodic();
230 _boundaryFaceIds=mesh.getBoundaryFaceIds();
231 _boundaryNodeIds=mesh.getBoundaryNodeIds();
233 _eltsTypes=mesh.getElementTypes();
234 _eltsTypesNames=mesh.getElementTypesNames();
236 MCAuto<MEDCouplingMesh> m1=mesh.getMEDCouplingMesh()->deepCopy();
238 _meshNotDeleted=mesh.meshNotDeleted();
241 //----------------------------------------------------------------------
243 Mesh::readMeshMed( const std::string filename, const int meshLevel)
244 //----------------------------------------------------------------------
246 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)
247 _mesh=m->getMeshAtLevel(meshLevel);
248 _mesh->checkConsistencyLight();
249 _mesh->setName(_mesh->getName());
250 _meshDim=_mesh->getMeshDimension();
251 _spaceDim=_mesh->getSpaceDimension();
252 _name=_mesh->getName();
254 _indexFacePeriodicSet=false;
255 MEDCoupling::MEDCouplingIMesh* structuredMesh = dynamic_cast<MEDCoupling::MEDCouplingIMesh*> (_mesh.retn());
259 _meshNotDeleted=false;
260 _dxyz=structuredMesh->getDXYZ();
261 _nxyz=structuredMesh->getCellGridStructure();
262 double* Box0=new double[2*_spaceDim];
263 structuredMesh->getBoundingBox(Box0);
281 _meshNotDeleted=true;
284 MEDCouplingUMesh* mu = setMesh();
287 cout<<endl<< "Loaded file "<< filename<<endl;
288 cout<<"Mesh name = "<<m->getName()<<", mesh dim = "<< _meshDim<< ", space dim = "<< _spaceDim<< ", nb cells= "<<getNumberOfCells()<< ", nb nodes= "<<getNumberOfNodes()<<endl;
294 Mesh::setGroupAtFaceByCoords(double x, double y, double z, double eps, std::string groupName, bool isBoundaryGroup)
296 std::vector< int > faceIds(0);
300 /* Construction of the face group */
302 for(int i=0; i<_boundaryFaceIds.size(); i++)
304 Fi=_faces[_boundaryFaceIds[i]];
308 if (abs(FX-x)<eps && abs(FY-y)<eps && abs(FZ-z)<eps)
310 faceIds.insert(faceIds.end(),_boundaryFaceIds[i]);
311 _faces[_boundaryFaceIds[i]].setGroupName(groupName);
315 for (int iface=0;iface<_numberOfFaces;iface++)
321 if (abs(FX-x)<eps && abs(FY-y)<eps && abs(FZ-z)<eps)
323 faceIds.insert(faceIds.end(),iface);
324 _faces[iface].setGroupName(groupName);
328 if (faceIds.size()>0)
330 std::vector< std::string >::iterator it = std::find(_faceGroupNames.begin(), _faceGroupNames.end(), groupName);
331 if(it == _faceGroupNames.end())//No group named groupName
333 _faceGroupNames.insert(_faceGroupNames.end(),groupName);
334 _faceGroupsIds.insert( _faceGroupsIds.end(),faceIds);
335 _faceGroups.insert( _faceGroups.end(), NULL);//No mesh created. Create one ?
339 std::vector< int > faceGroupIds = _faceGroupsIds[it-_faceGroupNames.begin()];
340 faceGroupIds.insert( faceGroupIds.end(), faceIds.begin(), faceIds.end());
341 /* Detect and erase duplicates face ids */
342 sort( faceGroupIds.begin(), faceGroupIds.end() );
343 faceGroupIds.erase( unique( faceGroupIds.begin(), faceGroupIds.end() ), faceGroupIds.end() );
344 _faceGroupsIds[it-_faceGroupNames.begin()] = faceGroupIds;
350 Mesh::setGroupAtNodeByCoords(double x, double y, double z, double eps, std::string groupName, bool isBoundaryGroup)
352 std::vector< int > nodeIds(0);
356 /* Construction of the node group */
358 for(int i=0; i<_boundaryNodeIds.size(); i++)
360 Ni=_nodes[_boundaryNodeIds[i]];
364 if (abs(NX-x)<eps && abs(NY-y)<eps && abs(NZ-z)<eps)
366 nodeIds.insert(nodeIds.end(),_boundaryNodeIds[i]);
367 _nodes[_boundaryNodeIds[i]].setGroupName(groupName);
371 for (int inode=0;inode<_numberOfNodes;inode++)
373 NX=_nodes[inode].x();
374 NY=_nodes[inode].y();
375 NZ=_nodes[inode].z();
376 if (abs(NX-x)<eps && abs(NY-y)<eps && abs(NZ-z)<eps)
378 nodeIds.insert(nodeIds.end(),inode);
379 _nodes[inode].setGroupName(groupName);
383 if (nodeIds.size()>0)
385 std::vector< std::string >::iterator it = std::find(_nodeGroupNames.begin(), _nodeGroupNames.end(), groupName);
386 if(it == _nodeGroupNames.end())//No group named groupName
388 _nodeGroupNames.insert(_nodeGroupNames.end(),groupName);
389 _nodeGroupsIds.insert( _nodeGroupsIds.end(),nodeIds);
390 _nodeGroups.insert( _nodeGroups.end(), NULL);//No mesh created. Create one ?
394 std::vector< int > nodeGroupIds = _nodeGroupsIds[it-_nodeGroupNames.begin()];
395 nodeGroupIds.insert( nodeGroupIds.end(), nodeIds.begin(), nodeIds.end());
396 /* Detect and erase duplicates node ids */
397 sort( nodeGroupIds.begin(), nodeGroupIds.end() );
398 nodeGroupIds.erase( unique( nodeGroupIds.begin(), nodeGroupIds.end() ), nodeGroupIds.end() );
399 _nodeGroupsIds[it-_nodeGroupNames.begin()] = nodeGroupIds;
405 Mesh::setGroupAtPlan(double value, int direction, double eps, std::string groupName, bool isBoundaryGroup)
407 std::vector< int > faceIds(0), nodeIds(0);
410 /* Construction of the face group */
412 for(int i=0; i<_boundaryFaceIds.size(); i++)
414 cord=_faces[_boundaryFaceIds[i]].getBarryCenter()[direction];
415 if (abs(cord-value)<eps)
417 faceIds.insert(faceIds.end(),_boundaryFaceIds[i]);
418 _faces[_boundaryFaceIds[i]].setGroupName(groupName);
422 for (int iface=0;iface<_numberOfFaces;iface++)
424 cord=_faces[iface].getBarryCenter()[direction];
425 if (abs(cord-value)<eps)
427 faceIds.insert(faceIds.end(),iface);
428 _faces[iface].setGroupName(groupName);
432 /* Construction of the node group */
434 for(int i=0; i<_boundaryNodeIds.size(); i++)
436 cord=_nodes[_boundaryNodeIds[i]].getPoint()[direction];
437 if (abs(cord-value)<eps)
439 nodeIds.insert(nodeIds.end(),_boundaryNodeIds[i]);
440 _nodes[_boundaryNodeIds[i]].setGroupName(groupName);
444 for (int inode=0;inode<_numberOfNodes;inode++)
446 cord=_nodes[inode].getPoint()[direction];
447 if (abs(cord-value)<eps)
449 nodeIds.insert(nodeIds.end(),inode);
450 _nodes[inode].setGroupName(groupName);
454 if (faceIds.size()>0)
456 std::vector< std::string >::iterator it = std::find(_faceGroupNames.begin(), _faceGroupNames.end(), groupName);
457 if(it == _faceGroupNames.end())//No group named groupName
459 _faceGroupNames.insert(_faceGroupNames.end(),groupName);
460 _faceGroupsIds.insert( _faceGroupsIds.end(),faceIds);
461 _faceGroups.insert( _faceGroups.end(), NULL);//No mesh created. Create one ?
465 std::vector< int > faceGroupIds = _faceGroupsIds[it-_faceGroupNames.begin()];
466 faceGroupIds.insert( faceGroupIds.end(), faceIds.begin(), faceIds.end());
467 /* Detect and erase duplicates face ids */
468 sort( faceGroupIds.begin(), faceGroupIds.end() );
469 faceGroupIds.erase( unique( faceGroupIds.begin(), faceGroupIds.end() ), faceGroupIds.end() );
470 _faceGroupsIds[it-_faceGroupNames.begin()] = faceGroupIds;
473 if (nodeIds.size()>0)
475 std::vector< std::string >::iterator it = std::find(_nodeGroupNames.begin(), _nodeGroupNames.end(), groupName);
476 if(it == _nodeGroupNames.end())//No group named groupName
478 _nodeGroupNames.insert(_nodeGroupNames.end(),groupName);
479 _nodeGroupsIds.insert( _nodeGroupsIds.end(),nodeIds);
480 _nodeGroups.insert( _nodeGroups.end(), NULL);//No mesh created. Create one ?
484 std::vector< int > nodeGroupIds = _nodeGroupsIds[it-_nodeGroupNames.begin()];
485 nodeGroupIds.insert( nodeGroupIds.end(), nodeIds.begin(), nodeIds.end());
486 /* Detect and erase duplicates node ids */
487 sort( nodeGroupIds.begin(), nodeGroupIds.end() );
488 nodeGroupIds.erase( unique( nodeGroupIds.begin(), nodeGroupIds.end() ), nodeGroupIds.end() );
489 _nodeGroupsIds[it-_nodeGroupNames.begin()] = nodeGroupIds;
495 Mesh::setBoundaryNodesFromFaces()
497 for (int iface=0;iface<_boundaryFaceIds.size();iface++)
499 std::vector< int > nodesID= _faces[_boundaryFaceIds[iface]].getNodesId();
500 int nbNodes = _faces[_boundaryFaceIds[iface]].getNumberOfNodes();
501 for(int inode=0 ; inode<nbNodes ; inode++)
503 std::vector<int>::const_iterator it = std::find(_boundaryNodeIds.begin(),_boundaryNodeIds.end(),nodesID[inode]);
504 if( it != _boundaryNodeIds.end() )
505 _boundaryNodeIds.push_back(nodesID[inode]);
511 Mesh::getIndexFacePeriodic( void ) const
513 return _indexFacePeriodicMap;
517 Mesh::setPeriodicFaces(bool check_groups, bool use_central_inversion)
519 if(_indexFacePeriodicSet)
522 for (int indexFace=0;indexFace<_boundaryFaceIds.size() ; indexFace++)
524 Face my_face=_faces[_boundaryFaceIds[indexFace]];
528 for (int iface=0;iface<_boundaryFaceIds.size() ; iface++)
531 iface_perio=_boundaryFaceIds[iface];
537 double x=my_face.x();
538 double y=my_face.y();
540 for (int iface=0;iface<_boundaryFaceIds.size() ; iface++)
542 Face face_i=_faces[_boundaryFaceIds[iface]];
543 double xi=face_i.x();
544 double yi=face_i.y();
545 if ( (abs(y-yi)<_epsilon || abs(x-xi)<_epsilon )// Case of a square geometry
546 && ( !check_groups || my_face.getGroupName()!=face_i.getGroupName()) //In case groups need to be checked
547 && ( !use_central_inversion || abs(y+yi) + abs(x+xi)<_epsilon ) // Case of a central inversion
548 && fabs(my_face.getMeasure() - face_i.getMeasure())<_epsilon
549 && fabs(my_face.getXN() + face_i.getXN())<_epsilon
550 && fabs(my_face.getYN() + face_i.getYN())<_epsilon
551 && fabs(my_face.getZN() + face_i.getZN())<_epsilon )
553 iface_perio=_boundaryFaceIds[iface];
560 double x=my_face.x();
561 double y=my_face.y();
562 double z=my_face.z();
564 for (int iface=0;iface<_boundaryFaceIds.size() ; iface++)
566 Face face_i=_faces[_boundaryFaceIds[iface]];
567 double xi=face_i.x();
568 double yi=face_i.y();
569 double zi=face_i.z();
570 if ( ((abs(y-yi)<_epsilon && abs(x-xi)<_epsilon) || (abs(x-xi)<_epsilon && abs(z-zi)<_epsilon) || (abs(y-yi)<_epsilon && abs(z-zi)<_epsilon))// Case of a cube geometry
571 && ( !check_groups || my_face.getGroupName()!=face_i.getGroupName()) //In case groups need to be checked
572 && ( !use_central_inversion || abs(y+yi) + abs(x+xi) + abs(z+zi)<_epsilon )// Case of a central inversion
573 && fabs(my_face.getMeasure() - face_i.getMeasure())<_epsilon
574 && fabs(my_face.getXN() + face_i.getXN())<_epsilon
575 && fabs(my_face.getYN() + face_i.getYN())<_epsilon
576 && fabs(my_face.getZN() + face_i.getZN())<_epsilon )
578 iface_perio=_boundaryFaceIds[iface];
584 throw CdmathException("Mesh::setPeriodicFaces: Mesh dimension should be 1, 2 or 3");
587 throw CdmathException("Mesh::setPeriodicFaces: periodic face not found, iface_perio==-1 " );
589 _indexFacePeriodicMap[_boundaryFaceIds[indexFace]]=iface_perio;
591 _indexFacePeriodicSet=true;
595 Mesh::getIndexFacePeriodic(int indexFace, bool check_groups, bool use_central_inversion)
597 if (!_faces[indexFace].isBorder())
599 cout<<"Pb with indexFace= "<<indexFace<<endl;
600 throw CdmathException("Mesh::getIndexFacePeriodic: not a border face" );
603 if(!_indexFacePeriodicSet)
604 setPeriodicFaces(check_groups, use_central_inversion);
606 std::map<int,int>::const_iterator it = _indexFacePeriodicMap.find(indexFace);
607 if( it != _indexFacePeriodicMap.end() )
611 cout<<"Pb with indexFace= "<<indexFace<<endl;
612 throw CdmathException("Mesh::getIndexFacePeriodic: not a periodic face" );
617 Mesh::isBorderNode(int nodeid) const
619 return _nodes[nodeid].isBorder();
623 Mesh::isBorderFace(int faceid) const
625 return _faces[faceid].isBorder();
629 Mesh::getBoundaryFaceIds() const
631 return _boundaryFaceIds;
635 Mesh::getBoundaryNodeIds() const
637 return _boundaryNodeIds;
641 Mesh::isTriangular() const
643 return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_TRI3;
646 Mesh::isTetrahedral() const
648 return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_TETRA4;
651 Mesh::isQuadrangular() const
653 return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_QUAD4;
656 Mesh::isHexahedral() const
658 return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_HEXA8;
661 Mesh::isStructured() const
663 return _isStructured;
666 std::vector< INTERP_KERNEL::NormalizedCellType >
667 Mesh::getElementTypes() const
672 std::vector< string >
673 Mesh::getElementTypesNames() const
675 std::vector< string > result(0);
676 for(int i=0; i< _eltsTypes.size(); i++)
678 if( _eltsTypes[i]==INTERP_KERNEL::NORM_POINT1)
679 result.push_back("Points");
680 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_SEG2)
681 result.push_back("Segments");
682 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_TRI3)
683 result.push_back("Triangles");
684 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_QUAD4)
685 result.push_back("Quadrangles");
686 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_POLYGON)
687 result.push_back("Polygons");
688 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_TETRA4)
689 result.push_back("Tetrahedra");
690 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_PYRA5)
691 result.push_back("Pyramids");
692 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_PENTA6)
693 result.push_back("Pentahedra");
694 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_HEXA8)
695 result.push_back("Hexahedra");
696 else if( _eltsTypes[i]==INTERP_KERNEL::NORM_POLYHED)
697 result.push_back("Polyhedrons");
700 cout<< "Mesh " + _name + " contains an element of type " <<endl;
701 cout<< _eltsTypes[i]<<endl;
702 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");
709 Mesh::setGroups( const MEDFileUMesh* medmesh, MEDCouplingUMesh* mu)
711 int nbCellsSubMesh, nbNodesSubMesh;
712 bool foundFace, foundNode;
714 /* Searching for face groups */
715 vector<string> faceGroups=medmesh->getGroupsNames() ;
717 for (unsigned int i=0;i<faceGroups.size();i++ )
719 string groupName=faceGroups[i];
720 vector<mcIdType> nonEmptyGrp(medmesh->getGrpNonEmptyLevels(groupName));
721 //We check if the group has a relative dimension equal to -1
722 //before call to the function getGroup(-1,groupName.c_str())
723 vector<mcIdType>::iterator it = find(nonEmptyGrp.begin(), nonEmptyGrp.end(), -1);
724 if (it != nonEmptyGrp.end())
726 MEDCouplingUMesh *m=medmesh->getGroup(-1,groupName.c_str());
728 m->sortCellsInMEDFileFrmt( );
729 nbCellsSubMesh=m->getNumberOfCells();
731 _faceGroups.insert(_faceGroups.end(),m);//Vector of group meshes
732 _faceGroupNames.insert(_faceGroupNames.end(),groupName);//Vector of group names
733 _faceGroupsIds.insert(_faceGroupsIds.end(),std::vector<int>(nbCellsSubMesh));//Vector of group face Ids. The filling of the face ids takes place below.
735 DataArrayDouble *baryCell = m->computeIsoBarycenterOfNodesPerCell();//computeCellCenterOfMass() ;
736 const double *coorBary=baryCell->getConstPointer();
737 /* Face identification */
738 for (int ic(0), k(0); ic<nbCellsSubMesh; ic++, k+=_spaceDim)
740 vector<double> coorBaryXyz(3,0);
741 for (int d=0; d<_spaceDim; d++)
742 coorBaryXyz[d] = coorBary[k+d];
743 Point p1(coorBaryXyz[0],coorBaryXyz[1],coorBaryXyz[2]) ;
746 for (int iface=0;iface<_numberOfFaces;iface++ )
748 Point p2=_faces[iface].getBarryCenter();
749 if(p1.distance(p2)<_epsilon)
751 _faces[iface].setGroupName(groupName);
752 _faceGroupsIds[_faceGroupsIds.size()-1][ic]=iface;
758 throw CdmathException("No face found for group " + groupName );
765 /* Searching for node groups */
766 vector<string> nodeGroups=medmesh->getGroupsOnSpecifiedLev(1) ;
768 for (unsigned int i=0;i<nodeGroups.size();i++ )
770 string groupName=nodeGroups[i];
771 DataArrayIdType * nodeGroup=medmesh->getNodeGroupArr( groupName );
772 const mcIdType *nodeids=nodeGroup->getConstPointer();
776 _nodeGroups.insert(_nodeGroups.end(),nodeGroup);
777 _nodeGroupNames.insert(_nodeGroupNames.end(),groupName);
779 nbNodesSubMesh=nodeGroup->getNumberOfTuples();//nodeGroup->getNbOfElems();
781 DataArrayDouble *coo = mu->getCoords() ;
782 const double *cood=coo->getConstPointer();
784 _nodeGroupsIds.insert(_nodeGroupsIds.end(),std::vector<int>(nbNodesSubMesh));//Vector of boundary faces
785 /* Node identification */
786 for (int ic(0); ic<nbNodesSubMesh; ic++)
788 vector<double> coorP(3,0);
789 for (int d=0; d<_spaceDim; d++)
790 coorP[d] = cood[nodeids[ic]*_spaceDim+d];
791 Point p1(coorP[0],coorP[1],coorP[2]) ;
794 for (int inode=0;inode<_numberOfNodes;inode++ )
796 Point p2=_nodes[inode].getPoint();
797 if(p1.distance(p2)<_epsilon)
799 _nodes[inode].setGroupName(groupName);
800 _nodeGroupsIds[_nodeGroupsIds.size()-1][ic]=inode;
806 throw CdmathException("No node found for group " + groupName );
812 //----------------------------------------------------------------------
814 Mesh::setMesh( void )
815 //----------------------------------------------------------------------
817 /* This is the main function translating medcouplingumesh info into Mesh class to be used when designing numerical methods
818 * We need the level 0 mesh to extract the cell-node connectvity
819 * We need the level -1 mesh to extract the cell-face and face-node connectivities (use o build descending connectivity)
820 * Be careful : the nodes in the medcoupling mesh are not necessarily all conected to a cell/face.
821 * Mesh class discard isolated nodes, hence the number of nodes in Mesh class can be lower than the number of nodes in medcouplingumesh.
824 DataArrayIdType *desc = DataArrayIdType::New();
825 DataArrayIdType *descI = DataArrayIdType::New();
826 DataArrayIdType *revDesc = DataArrayIdType::New();
827 DataArrayIdType *revDescI = DataArrayIdType::New();
828 MEDCouplingUMesh* mu = _mesh->buildUnstructured();
829 MEDCouplingUMesh* mu2;//mesh of dimension N-1 containing the cell interfaces->cell/face connectivity
832 mu->sortCellsInMEDFileFrmt( );
835 mu2=mu->buildDescendingConnectivity(desc,descI,revDesc,revDescI);
837 mu2=mu->buildDescendingConnectivity2(desc,descI,revDesc,revDescI);
839 const mcIdType *tmp = desc->getConstPointer();//Lists the faces surrounding each cell
840 const mcIdType *tmpI=descI->getConstPointer();
842 const mcIdType *tmpA =revDesc->getConstPointer();//Lists the cells surrounding each face
843 const mcIdType *tmpAI=revDescI->getConstPointer();
845 //Test du type d'éléments contenus dans le maillage afin d'éviter les éléments contenant des points de gauss
846 _eltsTypes=mu->getAllGeoTypesSorted();
847 for(int i=0; i<_eltsTypes.size();i++)
850 _eltsTypes[i]!= INTERP_KERNEL::NORM_POINT1 && _eltsTypes[i]!= INTERP_KERNEL::NORM_SEG2
851 && _eltsTypes[i]!= INTERP_KERNEL::NORM_TRI3 && _eltsTypes[i]!= INTERP_KERNEL::NORM_QUAD4
852 && _eltsTypes[i]!= INTERP_KERNEL::NORM_TETRA4 && _eltsTypes[i]!= INTERP_KERNEL::NORM_PYRA5
853 && _eltsTypes[i]!= INTERP_KERNEL::NORM_PENTA6 && _eltsTypes[i]!= INTERP_KERNEL::NORM_HEXA8
854 && _eltsTypes[i]!= INTERP_KERNEL::NORM_POLYGON&& _eltsTypes[i]!= INTERP_KERNEL::NORM_POLYHED
857 cout<< "Mesh " + mu->getName() + " contains an element of type " <<endl;
858 cout<< _eltsTypes[i]<<endl;
859 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");
863 DataArrayDouble *baryCell = mu->computeCellCenterOfMass() ;
864 const double *coorBary=baryCell->getConstPointer();//Used for cell center coordinates
866 MEDCouplingFieldDouble* fields=mu->getMeasureField(true);
867 DataArrayDouble *surface = fields->getArray();
868 const double *surf=surface->getConstPointer();//Used for cell lenght/surface/volume
870 DataArrayDouble *coo = mu->getCoords() ;
871 const double *cood=coo->getConstPointer();//Used for nodes coordinates
873 DataArrayIdType *revNode =DataArrayIdType::New();
874 DataArrayIdType *revNodeI=DataArrayIdType::New();
875 mu->getReverseNodalConnectivity(revNode,revNodeI) ;
876 const mcIdType *tmpN =revNode->getConstPointer();//Used to know which cells surround a given node
877 const mcIdType *tmpNI=revNodeI->getConstPointer();
879 DataArrayIdType *revCell =DataArrayIdType::New();
880 DataArrayIdType *revCellI=DataArrayIdType::New();
881 mu2->getReverseNodalConnectivity(revCell,revCellI);
882 const mcIdType *tmpC =revCell->getConstPointer();//Used to know which faces surround a given node
883 const mcIdType *tmpCI=revCellI->getConstPointer();
885 const DataArrayIdType *nodal = mu2->getNodalConnectivity() ;
886 const DataArrayIdType *nodalI = mu2->getNodalConnectivityIndex() ;
887 const mcIdType *tmpNE =nodal->getConstPointer();//Used to know which nodes surround a given face
888 const mcIdType *tmpNEI=nodalI->getConstPointer();
890 _numberOfCells = mu->getNumberOfCells() ;
891 _cells = new Cell[_numberOfCells] ;
893 _numberOfNodes = mu->getNumberOfNodes() ;//This number may include isolated nodes that will not be loaded. The number will be updated during nodes constructions
894 _nodes = new Node[_numberOfNodes] ;//This array may be resized if isolated nodes are found
896 _numberOfFaces = mu2->getNumberOfCells();
897 _faces = new Face[_numberOfFaces] ;
899 _indexFacePeriodicSet=false;
901 //Definition used if _meshDim =3 to determine the edges
902 DataArrayIdType *desc2 =DataArrayIdType::New();
903 DataArrayIdType *descI2=DataArrayIdType::New();
904 DataArrayIdType *revDesc2 =DataArrayIdType::New();
905 DataArrayIdType *revDescI2=DataArrayIdType::New();
906 DataArrayIdType *revNode2 =DataArrayIdType::New();
907 DataArrayIdType *revNodeI2=DataArrayIdType::New();
908 const mcIdType *tmpN2 ;
909 const mcIdType *tmpNI2;
910 MEDCouplingUMesh* mu3;
913 _numberOfEdges = mu->getNumberOfCells();
914 else if (_meshDim == 2)
915 _numberOfEdges = mu2->getNumberOfCells();
918 mu3=mu2->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2);//1D mesh of segments
919 _numberOfEdges = mu3->getNumberOfCells();
920 mu3->getReverseNodalConnectivity(revNode2,revNodeI2) ;
921 tmpN2 =revNode2->getConstPointer();
922 tmpNI2=revNodeI2->getConstPointer();
925 // _cells, _nodes and _faces initialization:
928 double xn, yn=0., zn=0.;//Components of the normal vector at a cell interface
930 for( int id=0;id<_numberOfCells;id++ )
932 Point p(0.0,0.0,0.0) ;
933 for(int idim=0; idim<_spaceDim; idim++)
934 p[idim]=coorBary[id*_spaceDim+idim];
936 mcIdType nbVertices=mu->getNumberOfNodesInCell(id) ;//should be equal to 2
937 assert( nbVertices==2);
938 std::vector<mcIdType> nodeIdsOfCell ;
939 mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
941 mcIdType nbFaces=tmpI[id+1]-tmpI[id];//should be equal to 2
943 const mcIdType *work=tmp+tmpI[id];
945 /* compute the normal to the face */
946 xn = cood[nodeIdsOfCell[0]*_spaceDim ] - cood[nodeIdsOfCell[nbFaces-1]*_spaceDim ];
948 yn = cood[nodeIdsOfCell[0]*_spaceDim+1] - cood[nodeIdsOfCell[nbFaces-1]*_spaceDim+1];
950 zn = cood[nodeIdsOfCell[0]*_spaceDim+2] - cood[nodeIdsOfCell[nbFaces-1]*_spaceDim+2];
951 norm = sqrt(xn*xn+yn*yn+zn*zn);
953 throw CdmathException("!!! Mesh::setMesh Normal vector has norm 0 !!!");
961 Cell ci( nbVertices, nbFaces, surf[id], p ) ;//nbCells=nbFaces=2
962 for( int el=0;el<nbFaces;el++ )
964 ci.addNodeId(el,nodeIdsOfCell[el]) ;//global node number
965 ci.addNormalVector(el,xn,yn,zn) ;
966 ci.addFaceId(el,work[el]) ;
967 xn = - xn; yn=-yn; zn=-zn;
972 for( int id(0); id<_numberOfFaces; id++ )
974 const mcIdType *workv=tmpNE+tmpNEI[id]+1;
975 mcIdType nbNodes= tmpNEI[id+1]-tmpNEI[id]-1;//Normally equal to 1.
978 std::vector<double> coo(0) ;
979 mu2->getCoordinatesOfNode(workv[0],coo);
981 for(int idim=0; idim<_spaceDim; idim++)
984 const mcIdType *workc=tmpA+tmpAI[id];
985 mcIdType nbCells=tmpAI[id+1]-tmpAI[id];
986 assert( nbCells>0);//To make sure our face is not located on an isolated node
988 Face fi( nbNodes, nbCells, 1.0, p, 1., 0., 0. ) ;
989 for(int node_id=0; node_id<nbNodes;node_id++)//This loop could b deleted since nbNodes=1. Trying to merge with setMesh
990 fi.addNodeId(node_id,workv[node_id]) ;//global node number
992 fi.addCellId(0,workc[0]) ;
993 for(int cell_id=1; cell_id<nbCells;cell_id++)
996 if (workc[cell_id]!=workc[cell_id-1])//For some meshes (bad ones) the same cell can appear several times
998 fi.addCellId(cell_idx+1,workc[cell_id]) ;
1003 _boundaryFaceIds.push_back(id);
1007 int correctNbNodes=0;
1008 for( int id=0;id<_numberOfNodes;id++ )
1010 const mcIdType *workc=tmpN+tmpNI[id];
1011 mcIdType nbCells=tmpNI[id+1]-tmpNI[id];
1013 if( nbCells>0)//To make sure this is not an isolated node
1016 std::vector<double> coo(0) ;
1017 mu->getCoordinatesOfNode(id,coo);
1018 Point p(0,0.0,0.0) ;
1019 for(int idim=0; idim<_spaceDim; idim++)
1022 const mcIdType *workf=tmpC+tmpCI[id];
1023 mcIdType nbFaces=tmpCI[id+1]-tmpCI[id];
1024 assert( nbFaces==1);
1026 const mcIdType *workn=tmpN+tmpNI[id];
1027 mcIdType nbNeighbourNodes=tmpNI[id+1]-tmpNI[id];
1029 Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
1030 for( int el=0;el<nbCells;el++ )
1031 vi.addCellId(el,workc[el]) ;
1032 for( int el=0;el<nbNeighbourNodes;el++ )
1033 vi.addNeighbourNodeId(el,workn[el]) ;//global node number
1034 for( int el=0;el<nbFaces;el++ )
1035 vi.addFaceId(el,workf[el],_faces[workf[el]].isBorder()) ;
1037 _boundaryNodeIds.push_back(id);
1041 if( _numberOfNodes!=correctNbNodes)
1042 cout<<"Found isolated nodes : correctNbNodes= "<<correctNbNodes<<", _numberOfNodes= "<<_numberOfNodes<<endl;
1044 else if(_meshDim==2 || _meshDim==3)
1046 DataArrayDouble *barySeg = mu2->computeIsoBarycenterOfNodesPerCell();//computeCellCenterOfMass() ;//Used as face center
1047 const double *coorBarySeg=barySeg->getConstPointer();
1049 MEDCouplingFieldDouble* fieldl=mu2->getMeasureField(true);
1050 DataArrayDouble *longueur = fieldl->getArray();
1051 const double *lon=longueur->getConstPointer();//The lenght/area of each face
1053 MEDCouplingFieldDouble* fieldn;//The normal to each face
1054 DataArrayDouble *normal;
1055 const double *tmpNormal;
1057 if(_spaceDim==_meshDim)
1058 fieldn = mu2->buildOrthogonalField();//Compute the normal to each cell interface
1060 fieldn = mu->buildOrthogonalField();//compute the 3D normal vector to the 2D cell
1062 normal = fieldn->getArray();
1063 tmpNormal = normal->getConstPointer();
1065 /*Building mesh cells */
1066 for(int id(0), k(0); id<_numberOfCells; id++, k+=_spaceDim)
1068 const mcIdType *work=tmp+tmpI[id];
1069 mcIdType nbFaces=tmpI[id+1]-tmpI[id];
1071 mcIdType nbVertices=mu->getNumberOfNodesInCell(id) ;
1073 vector<double> coorBaryXyz(3,0);
1074 for (int d=0; d<_spaceDim; d++)
1075 coorBaryXyz[d] = coorBary[k+d];
1077 Point p(coorBaryXyz[0],coorBaryXyz[1],coorBaryXyz[2]) ;
1078 Cell ci( nbVertices, nbFaces, surf[id], p ) ;
1080 /* Filling cell nodes */
1081 std::vector<mcIdType> nodeIdsOfCell ;
1082 mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
1083 for( int el=0;el<nbVertices;el++ )
1084 ci.addNodeId(el,nodeIdsOfCell[el]) ;
1086 /* Filling cell faces */
1087 if(_spaceDim==_meshDim)//use the normal field generated by buildOrthogonalField()
1088 for( int el=0;el<nbFaces;el++ )
1090 mcIdType faceIndex=(abs(work[el])-1);//=work[el] since Fortran type numbering was used, and negative sign means anticlockwise numbering
1091 vector<double> xyzn(3,0);//Outer normal to the cell
1093 for (int d=0; d<_spaceDim; d++)
1094 xyzn[d] = -tmpNormal[_spaceDim*faceIndex+d];
1096 for (int d=0; d<_spaceDim; d++)
1097 xyzn[d] = +tmpNormal[_spaceDim*faceIndex+d];
1098 ci.addNormalVector(el,xyzn[0],xyzn[1],xyzn[2]) ;
1099 ci.addFaceId(el,faceIndex) ;
1101 else//build normals associated to the couple (cell id, face el)
1102 {//Case _meshDim=1 should be moved upper since we are in the 2D/3D branch
1103 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
1104 {//work[0]= first face global number, work[1]= second face global number
1105 mcIdType indexFace0=abs(work[0])-1;//=work[0] since Fortran type numbering was used, and negative sign means anticlockwise numbering
1106 mcIdType indexFace1=abs(work[1])-1;//=work[1] since Fortran type numbering was used, and negative sign means anticlockwise numbering
1107 mcIdType idNodeA=(tmpNE+tmpNEI[indexFace0]+1)[0];//global number of the first face node work[0]=(abs(work[0])-1)
1108 mcIdType idNodeB=(tmpNE+tmpNEI[indexFace1]+1)[0];//global number of the second face node work[1]=(abs(work[1])-1)
1110 for(int i=0;i<_spaceDim;i++)
1111 vecAB[i]=coo->getIJ(idNodeB,i) - coo->getIJ(idNodeA,i);
1112 vecAB/=vecAB.norm();
1113 ci.addNormalVector(0,-vecAB[0],-vecAB[1],-vecAB[2]) ;
1114 ci.addNormalVector(1,vecAB[0],vecAB[1],vecAB[2]) ;
1115 ci.addFaceId(0,indexFace0) ;
1116 ci.addFaceId(1,indexFace1) ;
1118 else//_meshDim==2, number of faces around the cell id is variable, each face is composed of two nodes
1121 for (int d=0; d<_spaceDim; d++)
1122 xyzn[d] = tmpNormal[_spaceDim*id+d];
1123 for( int el=0;el<nbFaces;el++ )
1125 int faceIndex=(abs(work[el])-1);//=work[el] since Fortran type numbering was used, and negative sign means anticlockwise numbering
1126 const mcIdType *workv=tmpNE+tmpNEI[faceIndex]+1;
1127 mcIdType nbNodes= tmpNEI[faceIndex+1]-tmpNEI[faceIndex]-1;
1128 if(nbNodes!=2)//We want to compute the normal to a straight line, not a curved interface composed of more thant 2 points
1130 cout<<"Mesh name "<< mu->getName()<< " space dim= "<< _spaceDim <<" mesh dim= "<< _meshDim <<endl;
1131 cout<<"For cell id "<<id<<" and local face number "<<el<<", the number of nodes is "<< nbNodes<< ", total number of faces is "<< nbFaces <<endl;
1132 throw CdmathException("Mesh::setMesh number of nodes around a face should be 2");
1135 mcIdType idNodeA=workv[0];
1136 mcIdType idNodeB=workv[1];
1137 vector<double> nodeA(_spaceDim), nodeB(_spaceDim), nodeP(_spaceDim);
1138 for(int i=0;i<_spaceDim;i++)
1140 nodeA[i]=coo->getIJ(idNodeA,i);
1141 nodeB[i]=coo->getIJ(idNodeB,i);
1142 nodeP[i]=coorBary[_spaceDim*id+i];
1144 //Let P be the barycenter of the cell id
1145 Vector vecAB(3), vecPA(3);
1146 for(int i=0;i<_spaceDim;i++)
1148 vecAB[i]=coo->getIJ(idNodeB,i) - coo->getIJ(idNodeA,i);
1149 vecPA[i]=coo->getIJ(idNodeA,i) - coorBary[_spaceDim*id+i];
1152 Vector normale = xyzn % vecAB;//Normal to the edge
1153 normale/=normale.norm();
1156 ci.addNormalVector(el,normale[0],normale[1],normale[2]) ;
1158 ci.addNormalVector(el,-normale[0],-normale[1],-normale[2]) ;
1159 ci.addFaceId(el,faceIndex) ;
1166 if(_spaceDim!=_meshDim)
1168 /* Since spaceDim!=meshDim, don't build normal to faces */
1174 /*Building mesh faces */
1175 for(int id(0), k(0); id<_numberOfFaces; id++, k+=_spaceDim)
1177 vector<double> coorBarySegXyz(3);
1178 for (int d=0; d<_spaceDim; d++)
1179 coorBarySegXyz[d] = coorBarySeg[k+d];
1180 Point p(coorBarySegXyz[0],coorBarySegXyz[1],coorBarySegXyz[2]) ;
1181 const mcIdType *workc=tmpA+tmpAI[id];
1182 mcIdType nbCells=tmpAI[id+1]-tmpAI[id];
1184 if (nbCells>2 && _spaceDim==_meshDim)
1186 cout<<"Warning : nbCells>2, numberOfFaces="<<_numberOfFaces<<endl;
1187 cout<<"nbCells= "<<nbCells<<", _spaceDim="<<_spaceDim<<", _meshDim="<<_meshDim<<endl;
1188 for(int icell=0; icell<nbCells; icell++)
1189 cout<<workc[icell]<<", ";
1191 throw CdmathException("Wrong mesh : nbCells>2 and spaceDim==meshDim");
1194 _boundaryFaceIds.push_back(id);
1196 const mcIdType *workv=tmpNE+tmpNEI[id]+1;
1197 mcIdType nbNodes= tmpNEI[id+1]-tmpNEI[id]-1;
1200 if(_spaceDim==_meshDim)//Euclidean flat mesh geometry
1202 fi=Face( nbNodes, nbCells, lon[id], p, tmpNormal[k], tmpNormal[k+1], 0.0) ;
1204 fi=Face( nbNodes, nbCells, lon[id], p, tmpNormal[k], tmpNormal[k+1], tmpNormal[k+2]) ;
1205 else//Curved mesh geometry
1206 fi=Face( nbNodes, nbCells, lon[id], p, 0.0, 0.0, 0.0) ;//Since spaceDim!=meshDim, normal to face is not defined
1208 for(int node_id=0; node_id<nbNodes;node_id++)
1209 fi.addNodeId(node_id,workv[node_id]) ;
1211 fi.addCellId(0,workc[0]) ;
1212 for(int cell_id=1; cell_id<nbCells;cell_id++)
1215 if (workc[cell_id]!=workc[cell_id-1])//For some meshes (bad ones) the same cell can appear several times
1217 fi.addCellId(cell_idx+1,workc[cell_id]) ;
1225 /*Building mesh nodes, should be done after building mesh faces in order to detect boundary nodes*/
1226 int correctNbNodes=0;
1227 for(int id(0), k(0); id<_numberOfNodes; id++, k+=_spaceDim)
1229 const mcIdType *workc=tmpN+tmpNI[id];
1230 mcIdType nbCells=tmpNI[id+1]-tmpNI[id];
1232 if( nbCells>0)//To make sure this is not an isolated node
1235 vector<double> coorP(3);
1236 for (int d=0; d<_spaceDim; d++)
1237 coorP[d] = cood[k+d];
1238 Point p(coorP[0],coorP[1],coorP[2]) ;
1240 const mcIdType *workf=tmpC+tmpCI[id];
1241 mcIdType nbFaces=tmpCI[id+1]-tmpCI[id];
1242 const mcIdType *workn;
1243 mcIdType nbNeighbourNodes;
1246 workn=tmpA+tmpAI[id];
1247 nbNeighbourNodes=tmpAI[id+1]-tmpAI[id];
1249 else if (_meshDim == 2)
1251 workn=tmpC+tmpCI[id];
1252 nbNeighbourNodes=tmpCI[id+1]-tmpCI[id];
1256 workn=tmpN2+tmpNI2[id];
1257 nbNeighbourNodes=tmpNI2[id+1]-tmpNI2[id];
1259 Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
1261 for( int el=0;el<nbCells;el++ )
1262 vi.addCellId(el,workc[el]) ;
1263 for( int el=0;el<nbNeighbourNodes;el++ )
1264 vi.addNeighbourNodeId(el,workn[el]) ;
1265 //Detection of border nodes
1266 for( int el=0;el<nbFaces;el++ )
1267 vi.addFaceId(el,workf[el],_faces[workf[el]].isBorder()) ;
1269 _boundaryNodeIds.push_back(id);
1273 if( _numberOfNodes!=correctNbNodes)
1274 cout<<"Found isolated nodes : correctNbNodes= "<<correctNbNodes<<", _numberOfNodes= "<<_numberOfNodes<<endl;
1276 if(_spaceDim==_meshDim)
1282 throw CdmathException("Mesh::setMesh space dimension should be 1, 2 or 3");
1284 //definition of the bounding box for unstructured meshes
1285 if(!_isStructured)//Structured case is treated in function readMeshMed
1287 double Box0[2*_spaceDim];
1288 coo->getMinMaxPerComponent(Box0);
1304 //Set boundary groups
1305 _faceGroupNames.push_back("Boundary");
1306 _nodeGroupNames.push_back("Boundary");
1307 _faceGroupsIds.push_back(_boundaryFaceIds);
1308 _nodeGroupsIds.push_back(_boundaryNodeIds);
1310 { //Set face boundary group
1311 _boundaryMesh = mu->computeSkin();
1312 _faceGroups.push_back(_boundaryMesh);
1315 _faceGroups.push_back(NULL);
1316 _nodeGroups.push_back(NULL);
1321 revDescI->decrRef();
1323 baryCell->decrRef();
1326 revNodeI->decrRef();
1328 revCellI->decrRef();
1332 revNode2->decrRef();
1333 revNodeI2->decrRef();
1336 revDesc2->decrRef();
1337 revDescI2->decrRef();
1344 //----------------------------------------------------------------------
1345 Mesh::Mesh( std::vector<double> points, std::string meshName )
1346 //----------------------------------------------------------------------
1348 int nx=points.size();
1351 throw CdmathException("Mesh::Mesh( vector<double> points, string meshName) : nx < 2, vector should contain at least two values");
1353 while( i<nx-1 && points[i+1]>points[i] )
1357 //cout<< points << endl;
1358 throw CdmathException("Mesh::Mesh( vector<double> points, string meshName) : vector values should be sorted");
1372 _isStructured = false;
1373 _indexFacePeriodicSet=false;
1375 MEDCouplingUMesh * mesh1d = MEDCouplingUMesh::New(meshName, 1);
1376 mesh1d->allocateCells(nx - 1);
1377 double * coords = new double[nx];
1378 mcIdType * nodal_con = new mcIdType[2];
1379 coords[0]=points[0];
1380 for(int i=0; i<nx- 1 ; i++)
1384 mesh1d->insertNextCell(INTERP_KERNEL::NORM_SEG2, 2, nodal_con);
1385 coords[i+1]=points[i + 1];
1387 mesh1d->finishInsertingCells();
1389 DataArrayDouble * coords_arr = DataArrayDouble::New();
1390 coords_arr->alloc(nx,1);
1391 std::copy(coords,coords+nx,coords_arr->getPointer());
1392 mesh1d->setCoords(coords_arr);
1394 delete [] coords, nodal_con;
1395 coords_arr->decrRef();
1397 _mesh=mesh1d->buildUnstructured();//To enable writeMED. Because we declared the mesh as unstructured, we decide to build the unstructured data (not mandatory)
1398 _meshNotDeleted=true;
1403 //----------------------------------------------------------------------
1404 Mesh::Mesh( double xmin, double xmax, int nx, std::string meshName )
1405 //----------------------------------------------------------------------
1408 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx) : nx <= 0");
1410 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx) : xmin >= xmax");
1412 double dx = (xmax - xmin)/nx ;
1418 _isStructured = true;
1419 _indexFacePeriodicSet=false;
1428 _dxyz.resize(_spaceDim);
1430 _nxyz.resize(_spaceDim);
1433 double *originPtr = new double[_spaceDim];
1434 double *dxyzPtr = new double[_spaceDim];
1435 mcIdType *nodeStrctPtr = new mcIdType[_spaceDim];
1438 nodeStrctPtr[0]=nx+1;
1441 _mesh=MEDCouplingIMesh::New(meshName,
1444 nodeStrctPtr+_spaceDim,
1446 originPtr+_spaceDim,
1449 _meshNotDeleted=true;//Because the mesh is structured cartesian : no data in memory. No nodes and cell coordinates stored
1451 delete [] originPtr;
1453 delete [] nodeStrctPtr;
1458 //----------------------------------------------------------------------
1459 Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, int split_to_triangles_policy, std::string meshName)
1460 //----------------------------------------------------------------------
1463 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : nx <= 0 or ny <= 0");
1465 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : xmin >= xmax");
1467 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : ymin >= ymax");
1477 double dx = (xmax - xmin)/nx ;
1478 double dy = (ymax - ymin)/ny ;
1484 _indexFacePeriodicSet=false;
1485 _nxyz.resize(_spaceDim);
1489 _dxyz.resize(_spaceDim);
1493 double *originPtr = new double[_spaceDim];
1494 double *dxyzPtr = new double[_spaceDim];
1495 mcIdType *nodeStrctPtr = new mcIdType[_spaceDim];
1499 nodeStrctPtr[0]=nx+1;
1500 nodeStrctPtr[1]=ny+1;
1504 _mesh=MEDCouplingIMesh::New(meshName,
1507 nodeStrctPtr+_spaceDim,
1509 originPtr+_spaceDim,
1512 _meshNotDeleted=true;//Because the mesh is structured cartesian : no data in memory. No nodes and cell coordinates stored
1513 _isStructured = true;
1515 if(split_to_triangles_policy==0 || split_to_triangles_policy==1)
1517 _mesh=_mesh->buildUnstructured();
1518 _mesh->simplexize(split_to_triangles_policy);
1519 _meshNotDeleted=true;//Now the mesh is unstructured and stored with nodes and cell coordinates
1520 _isStructured = false;
1522 else if (split_to_triangles_policy != -1)
1524 cout<< "split_to_triangles_policy = "<< split_to_triangles_policy << endl;
1525 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : Unknown splitting policy");
1528 delete [] originPtr;
1530 delete [] nodeStrctPtr;
1535 //----------------------------------------------------------------------
1536 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)
1537 //----------------------------------------------------------------------
1539 if(nx<=0 || ny<=0 || nz<=0)
1540 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");
1542 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : xmin >= xmax");
1544 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : ymin >= ymax");
1546 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : zmin >= zmax");
1559 double dx = (xmax - xmin)/nx ;
1560 double dy = (ymax - ymin)/ny ;
1561 double dz = (zmax - zmin)/nz ;
1563 _dxyz.resize(_spaceDim);
1568 _nxyz.resize(_spaceDim);
1573 double *originPtr = new double[_spaceDim];
1574 double *dxyzPtr = new double[_spaceDim];
1575 mcIdType *nodeStrctPtr = new mcIdType[_spaceDim];
1580 nodeStrctPtr[0]=nx+1;
1581 nodeStrctPtr[1]=ny+1;
1582 nodeStrctPtr[2]=nz+1;
1587 _mesh=MEDCouplingIMesh::New(meshName,
1590 nodeStrctPtr+_spaceDim,
1592 originPtr+_spaceDim,
1595 _meshNotDeleted=true;//Because the mesh is structured cartesian : no data in memory. Nno nodes and cell coordinates stored
1596 _isStructured = true;
1598 if( split_to_tetrahedra_policy == 0 )
1600 _mesh=_mesh->buildUnstructured();
1601 _mesh->simplexize(INTERP_KERNEL::PLANAR_FACE_5);
1602 _meshNotDeleted=true;//Now the mesh is unstructured and stored with nodes and cell coordinates
1603 _isStructured = false;
1605 else if( split_to_tetrahedra_policy == 1 )
1607 _mesh=_mesh->buildUnstructured();
1608 _mesh->simplexize(INTERP_KERNEL::PLANAR_FACE_6);
1609 _meshNotDeleted=true;//Now the mesh is unstructured and stored with nodes and cell coordinates
1610 _isStructured = false;
1612 else if ( split_to_tetrahedra_policy != -1 )
1614 cout<< "split_to_tetrahedra_policy = "<< split_to_tetrahedra_policy << endl;
1615 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");
1618 delete [] originPtr;
1620 delete [] nodeStrctPtr;
1625 //----------------------------------------------------------------------
1627 Mesh::getSpaceDimension( void ) const
1628 //----------------------------------------------------------------------
1633 //----------------------------------------------------------------------
1635 Mesh::getMeshDimension( void ) const
1636 //----------------------------------------------------------------------
1642 Mesh::getDXYZ() const
1645 throw CdmathException("std::vector<double> Mesh::getDXYZ() : dx,dy and dz are defined only for structured meshes !");
1650 std::vector<mcIdType>
1651 Mesh::getCellGridStructure() const
1654 throw CdmathException("std::vector<int> Mesh::getCellGridStructure() : nx, ny and nz are defined only for structured meshes !");
1659 //----------------------------------------------------------------------
1661 Mesh::getNx( void ) const
1662 //----------------------------------------------------------------------
1665 throw CdmathException("int Mesh::getNx( void ) : Nx is defined only for structured meshes !");
1670 //----------------------------------------------------------------------
1672 Mesh::getNy( void ) const
1673 //----------------------------------------------------------------------
1676 throw CdmathException("int Mesh::getNy( void ) : Ny is defined only for structured meshes !");
1678 throw CdmathException("int double& Field::operator[ielem] : Ny is not defined in dimension < 2!");
1683 //----------------------------------------------------------------------
1685 Mesh::getNz( void ) const
1686 //----------------------------------------------------------------------
1689 throw CdmathException("int Mesh::getNz( void ) : Nz is defined only for structured meshes !");
1691 throw CdmathException("int Mesh::getNz( void ) : Nz is not defined in dimension < 3!");
1696 //----------------------------------------------------------------------
1698 Mesh::getXMin( void ) const
1699 //----------------------------------------------------------------------
1704 //----------------------------------------------------------------------
1706 Mesh::getXMax( void ) const
1707 //----------------------------------------------------------------------
1712 //----------------------------------------------------------------------
1714 Mesh::getYMin( void ) const
1715 //----------------------------------------------------------------------
1720 //----------------------------------------------------------------------
1722 Mesh::getYMax( void ) const
1723 //----------------------------------------------------------------------
1728 //----------------------------------------------------------------------
1730 Mesh::getZMin( void ) const
1731 //----------------------------------------------------------------------
1736 //----------------------------------------------------------------------
1738 Mesh::getZMax( void ) const
1739 //----------------------------------------------------------------------
1744 //----------------------------------------------------------------------
1745 MCAuto<MEDCouplingMesh>
1746 Mesh::getMEDCouplingMesh( void ) const
1747 //----------------------------------------------------------------------
1752 //----------------------------------------------------------------------
1754 Mesh::getNumberOfNodes ( void ) const
1755 //----------------------------------------------------------------------
1757 return _numberOfNodes ;
1760 //----------------------------------------------------------------------
1762 Mesh::getNumberOfCells ( void ) const
1763 //----------------------------------------------------------------------
1765 return _numberOfCells ;
1768 //----------------------------------------------------------------------
1770 Mesh::getNumberOfFaces ( void ) const
1771 //----------------------------------------------------------------------
1773 return _numberOfFaces ;
1776 //----------------------------------------------------------------------
1778 Mesh::getNumberOfEdges ( void ) const
1779 //----------------------------------------------------------------------
1781 return _numberOfEdges ;
1784 //----------------------------------------------------------------------
1786 Mesh::getFaces ( void ) const
1787 //----------------------------------------------------------------------
1792 //----------------------------------------------------------------------
1794 Mesh::getCells ( void ) const
1795 //----------------------------------------------------------------------
1800 //----------------------------------------------------------------------
1802 Mesh::getCell ( int i ) const
1803 //----------------------------------------------------------------------
1808 //----------------------------------------------------------------------
1810 Mesh::getFace ( int i ) const
1811 //----------------------------------------------------------------------
1816 //----------------------------------------------------------------------
1818 Mesh::getNode ( int i ) const
1819 //----------------------------------------------------------------------
1824 //----------------------------------------------------------------------
1826 Mesh::getNodes ( void ) const
1827 //----------------------------------------------------------------------
1833 Mesh::getNameOfFaceGroups( void ) const
1835 return _faceGroupNames;
1838 vector<MEDCoupling::MEDCouplingUMesh *>
1839 Mesh::getFaceGroups( void ) const
1845 Mesh::getNameOfNodeGroups( void ) const
1847 return _nodeGroupNames;
1850 vector<MEDCoupling::DataArrayIdType *>
1851 Mesh::getNodeGroups( void ) const
1856 //----------------------------------------------------------------------
1858 Mesh::operator= ( const Mesh& mesh )
1859 //----------------------------------------------------------------------
1861 _spaceDim = mesh.getSpaceDimension() ;
1862 _meshDim = mesh.getMeshDimension() ;
1863 _name = mesh.getName();
1864 _epsilon=mesh.getComparisonEpsilon();
1865 _numberOfNodes = mesh.getNumberOfNodes();
1866 _numberOfFaces = mesh.getNumberOfFaces();
1867 _numberOfCells = mesh.getNumberOfCells();
1868 _numberOfEdges = mesh.getNumberOfEdges();
1870 _isStructured = mesh.isStructured();
1871 _meshNotDeleted = mesh.meshNotDeleted();
1875 _nxyz = mesh.getCellGridStructure() ;
1876 _dxyz=mesh.getDXYZ();
1877 _xMin=mesh.getXMin();
1878 _xMax=mesh.getXMax();
1879 _yMin=mesh.getYMin();
1880 _yMax=mesh.getYMax();
1881 _zMin=mesh.getZMin();
1882 _zMax=mesh.getZMax();
1884 _faceGroupNames = mesh.getNameOfFaceGroups() ;
1885 _faceGroups = mesh.getFaceGroups() ;
1886 _nodeGroupNames = mesh.getNameOfNodeGroups() ;
1887 _nodeGroups = mesh.getNodeGroups() ;
1905 _nodes = new Node[_numberOfNodes] ;
1906 _faces = new Face[_numberOfFaces] ;
1907 _cells = new Cell[_numberOfCells] ;
1909 for (int i=0;i<_numberOfNodes;i++)
1910 _nodes[i]=mesh.getNode(i);
1912 for (int i=0;i<_numberOfFaces;i++)
1913 _faces[i]=mesh.getFace(i);
1915 for (int i=0;i<_numberOfCells;i++)
1916 _cells[i]=mesh.getCell(i);
1918 _indexFacePeriodicSet= mesh.isIndexFacePeriodicSet();
1919 if(_indexFacePeriodicSet)
1920 _indexFacePeriodicMap=mesh.getIndexFacePeriodic();
1922 _boundaryFaceIds=mesh.getBoundaryFaceIds();
1923 _boundaryNodeIds=mesh.getBoundaryNodeIds();
1925 _eltsTypes=mesh.getElementTypes();
1926 _eltsTypesNames=mesh.getElementTypesNames();
1928 MCAuto<MEDCouplingMesh> m1=mesh.getMEDCouplingMesh()->deepCopy();
1933 bool Mesh::isIndexFacePeriodicSet() const
1935 return _indexFacePeriodicSet;
1937 //----------------------------------------------------------------------
1939 Mesh::minRatioVolSurf() const
1941 double dx_min = 1e30;
1942 for(int i=0; i<_numberOfCells; i++)
1944 Cell Ci = getCell(i);
1948 for(int k=0; k< Ci.getNumberOfFaces(); k++)
1950 int indexFace=Ci.getFacesId()[k];
1951 Face Fk = getFace(indexFace);
1952 perimeter+=Fk.getMeasure();
1954 dx_min = min(dx_min,Ci.getMeasure()/perimeter);
1957 dx_min = min(dx_min,Ci.getMeasure());
1964 Mesh::getBoundaryMesh ( void ) const
1966 return Mesh(_boundaryMesh);
1970 Mesh::getBoundaryGroupMesh ( std::string groupName, int nth_occurence ) const
1972 //count occurences of groupName in known group name list
1973 int count_occurences = std::count (_faceGroupNames.begin(),_faceGroupNames.end(),groupName);
1975 if( count_occurences ==0 )//No group found
1977 cout<<"Mesh::getBoundaryGroupMesh Error : face group " << groupName << " does not exist"<<endl;
1978 cout<<"Known face group names are " ;
1979 for(int i=0; i<_faceGroupNames.size(); i++)
1980 cout<< _faceGroupNames[i]<<" ";
1982 throw CdmathException("Required face group does not exist");
1984 else if ( count_occurences <= nth_occurence)//Too many groups found
1986 cout<<"Mesh::getBoundaryGroupMesh Error : "<<count_occurences<<" groups have name " << groupName<<", but you asked fo occurencer number "<<nth_occurence<<"which is too large"<<endl;
1987 cout<<"Call function getBoundaryGroupMesh ( string groupName, int nth_group_match ) with nth_group_match between 0 and "<<count_occurences-1<<" to discriminate them "<<endl ;
1988 throw CdmathException("Several face groups have the same name but you asked for an occurence that does not exsit");
1990 else if( count_occurences >1 )//Wrning several groups found
1991 cout<<"Warning : "<<count_occurences<<" groups have name " << groupName<<". Searching occurence number "<<nth_occurence<<endl;
1993 //search occurence of group name in known group name list
1994 std::vector<std::string>::const_iterator it = _faceGroupNames.begin();
1995 for (int i = 0; i<nth_occurence+1; i++)
1996 it = std::find(it,_faceGroupNames.end(),groupName);
1998 return Mesh(_faceGroups[it - _faceGroupNames.begin()]);
2002 Mesh::getMaxNbNeighbours(EntityType type) const
2008 int nbNeib;//local number of neighbours
2009 for(int i=0; i<_numberOfCells; i++)
2011 Cell Ci = _cells[i];
2012 //Careful with mesh with junctions : do not just take Ci.getNumberOfFaces()
2014 for(int j=0; j<Ci.getNumberOfFaces(); j++)
2015 nbNeib+=_faces[Ci.getFacesId()[j]].getNumberOfCells()-1;//Without junction this would be +=1
2021 else if(type==NODES)
2023 for(int i=0; i<_numberOfNodes; i++)
2024 if(result < _nodes[i].getNumberOfEdges())
2025 result=_nodes[i].getNumberOfEdges();
2028 throw CdmathException("Mesh::getMaxNbNeighbours : entity type is not accepted. Should be CELLS or NODES");
2032 //----------------------------------------------------------------------
2034 Mesh::writeVTK ( const std::string fileName ) const
2035 //----------------------------------------------------------------------
2037 if( !_isStructured && !_meshNotDeleted )
2038 throw CdmathException("Mesh::writeVTK : Cannot save mesh : no MEDCouplingUMesh loaded");
2040 string fname=fileName+".vtu";
2041 _mesh->writeVTK(fname.c_str()) ;
2044 //----------------------------------------------------------------------
2046 Mesh::writeMED ( const std::string fileName ) const
2047 //----------------------------------------------------------------------
2049 if( !_meshNotDeleted )
2050 throw CdmathException("Mesh::writeMED : Cannot save mesh : no MEDCouplingUMesh loaded");
2052 string fname=fileName+".med";
2053 MEDCouplingUMesh* mu=_mesh->buildUnstructured();
2054 MEDCoupling::WriteUMesh(fname.c_str(),mu,true);
2056 /* Try to save mesh groups */
2057 //MEDFileUMesh meshMEDFile;
2058 //meshMEDFile.setMeshAtLevel(0,mu);
2059 //for(int i=0; i< _faceGroups.size(); i++)
2060 //meshMEDFile.setMeshAtLevel(-1,_faceGroups[i]);
2062 //MEDCoupling::meshMEDFile.write(fname.c_str(),2) ;
2064 //MEDCoupling::meshMEDFile.write(fname.c_str(),1) ;
2068 Mesh::getFaceGroupIds(std::string groupName, bool isBoundaryGroup) const
2070 std::vector<std::string>::const_iterator it = std::find(_faceGroupNames.begin(),_faceGroupNames.end(),groupName);
2071 if( it == _faceGroupNames.end() )
2073 cout<<"Mesh::getGroupFaceIds Error : face group " << groupName << " does not exist"<<endl;
2074 throw CdmathException("Required face group does not exist");
2078 return _faceGroupsIds[it-_faceGroupNames.begin()];
2083 Mesh::getNodeGroupIds(std::string groupName, bool isBoundaryGroup) const
2085 std::vector<std::string>::const_iterator it = std::find(_nodeGroupNames.begin(),_nodeGroupNames.end(),groupName);
2086 if( it == _nodeGroupNames.end() )
2088 cout<<"Mesh::getGroupNodeIds Error : node group " << groupName << " does not exist"<<endl;
2089 throw CdmathException("Required node group does not exist");
2092 return _nodeGroupsIds[it-_nodeGroupNames.begin()];
2096 Mesh::setFaceGroupByIds(std::vector< int > faceIds, std::string groupName)
2098 for(int i=0; i< faceIds.size(); i++)
2099 getFace(faceIds[i]).setGroupName(groupName);
2101 _faceGroupsIds.insert( _faceGroupsIds.end(),faceIds);
2102 _faceGroupNames.insert(_faceGroupNames.end(), groupName);
2103 _faceGroups.insert( _faceGroups.end(), NULL);
2107 Mesh::setNodeGroupByIds(std::vector< int > nodeIds, std::string groupName)
2109 for(int i=0; i< nodeIds.size(); i++)
2110 getNode(nodeIds[i]).setGroupName(groupName);
2113 void Mesh::deleteMEDCouplingUMesh()
2117 (_mesh.retn())->decrRef();
2118 _meshNotDeleted=false;
2121 throw CdmathException("Mesh::deleteMEDCouplingMesh() : mesh is not loaded");