1 // Copyright (C) 2012-2015 ALNEOS
2 // Copyright (C) 2016-2023 EDF
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 // See http://www.alneos.com/ or email : contact@alneos.fr
19 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 #include "GMSHPlugin_Mesher.hxx"
22 #include "GMSHPlugin_Hypothesis_2D.hxx"
24 #include <SMDS_MeshElement.hxx>
25 #include <SMDS_MeshNode.hxx>
26 #include <SMESHDS_Mesh.hxx>
27 #include <SMESH_Comment.hxx>
28 #include <SMESH_ComputeError.hxx>
29 #include <SMESH_Gen_i.hxx>
30 #include <SMESH_Mesh.hxx>
31 #include <SMESH_MesherHelper.hxx>
32 #include <SMESH_subMesh.hxx>
33 #include <StdMeshers_FaceSide.hxx>
34 #include <utilities.h>
35 #include <StdMeshers_QuadToTriaAdaptor.hxx>
38 #include <BRep_Tool.hxx>
39 #include <GeomAPI_ProjectPointOnCurve.hxx>
47 #include <TopExp_Explorer.hxx>
51 #include <MTriangle.h>
52 #include <MQuadrangle.h>
53 #if GMSH_MAJOR_VERSION >=4
54 #include <GmshGlobal.h>
55 #include <gmsh/Context.h>
58 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
68 SBoundingBox3d _bounds;
72 //================================================================================
74 * \brief Retrieve ShapeBounds from a compound GEdge
76 //================================================================================
78 bool getBoundsOfShapes( GEdge* gEdge,
79 std::vector< ShapeBounds > & topoEdges )
82 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
83 for ( size_t i = 0; i < gEdge->compound.size(); ++i )
85 GEdge* gE = static_cast< GEdge* >( gEdge->compound[ i ]);
86 topoEdges.push_back( ShapeBounds{ gE->bounds(), *((TopoDS_Edge*)gE->getNativePtr()) });
89 for ( size_t i = 0; i < gEdge->_compound.size(); ++i )
91 GEdge* gE = static_cast< GEdge* >( gEdge->_compound[ i ]);
92 topoEdges.push_back( ShapeBounds{ gE->bounds(), *((TopoDS_Edge*)gE->getNativePtr()) });
95 return topoEdges.size();
98 //================================================================================
100 * \brief Retrieve ShapeBounds from a compound GFace
102 //================================================================================
104 bool getBoundsOfShapes( GFace* gFace,
105 std::vector< ShapeBounds > & topoFaces )
108 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
109 for ( size_t i = 0; i < gFace->compound.size(); ++i )
111 GFace* gF = static_cast< GFace* >( gFace->compound[ i ]);
112 topoFaces.push_back( ShapeBounds{ gF->bounds(), *((TopoDS_Face*)gF->getNativePtr()) });
115 for ( size_t i = 0; i < gFace->_compound.size(); ++i )
117 GFace* gF = static_cast< GFace* >( gFace->_compound[ i ]);
118 topoFaces.push_back( ShapeBounds{ gF->bounds(), *((TopoDS_Face*)gF->getNativePtr()) });
121 return topoFaces.size();
123 //================================================================================
125 * \brief Find a shape whose bounding box includes a given point
127 //================================================================================
129 TopoDS_Shape getShapeAtPoint( const SPoint3& point, const std::vector< ShapeBounds > & shapes )
132 float distmin = std::numeric_limits<float>::max();
133 for ( size_t i = 0; i < shapes.size(); ++i )
135 float dist = GMSHPlugin_Mesher::DistBoundingBox( shapes[i]._bounds, point );
138 shape = shapes[i]._shape;
147 double segmentSize( const UVPtStructVec& nodeParam, size_t i )
149 double l1 = SMESH_NodeXYZ( nodeParam[i].node ).Distance( nodeParam[i-1].node );
150 double l2 = SMESH_NodeXYZ( nodeParam[i].node ).Distance( nodeParam[i+1].node );
151 return 0.5 * ( l1 + l2 );
155 //=============================================================================
159 //=============================================================================
161 GMSHPlugin_Mesher::GMSHPlugin_Mesher (SMESH_Mesh* mesh,
162 const TopoDS_Shape& aShape,
170 // il faudra peut être mettre un truc par defaut si l'utilisateur ne rentre rien en para
171 //defaultParameters();
174 //void GMSHPlugin_Mesher::defaultParameters(){}
176 void GMSHPlugin_Mesher::SetParameters(const GMSHPlugin_Hypothesis* hyp)
180 _algo2d = hyp->Get2DAlgo();
181 _algo3d = hyp->Get3DAlgo();
182 _recomb2DAlgo = hyp->GetRecomb2DAlgo();
183 _recombineAll = hyp->GetRecombineAll();
184 _subdivAlgo = hyp->GetSubdivAlgo();
185 _remeshAlgo = hyp->GetRemeshAlgo();
186 _remeshPara = hyp->GetRemeshPara();
187 _smouthSteps = hyp->GetSmouthSteps();
188 _sizeFactor = hyp->GetSizeFactor();
189 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=10
190 _meshCurvatureSize = hyp->GetMeshCurvatureSize();
192 _minSize = hyp->GetMinSize();
193 _maxSize = hyp->GetMaxSize();
194 _secondOrder = hyp->GetSecondOrder();
195 _useIncomplElem = hyp->GetUseIncomplElem();
196 _verbLvl = hyp->GetVerbosityLevel();
197 _compounds = hyp->GetCompoundOnEntries();
198 // 6 in the enum corresponds to 99 in gmsh
207 _recombineAll = false;
213 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=10
214 _meshCurvatureSize = 0;
218 _secondOrder = false;
219 _useIncomplElem = true;
225 //================================================================================
227 * \brief Set maximum number of threads to be used by Gmsh
229 //================================================================================
231 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
232 void GMSHPlugin_Mesher::SetMaxThreadsGmsh()
234 MESSAGE("GMSHPlugin_Mesher::SetMaxThreadsGmsh");
235 // compound meshing (_compounds.size() > 0) and quad meshing (_algo2d >= 5) will
236 // not be multi-threaded
237 if (_compounds.size() > 0 || _algo2d >= 5){
241 _maxThreads = omp_get_max_threads();
246 //================================================================================
248 * \brief Check that the passed nodes are all IN the face
249 * \param element the element
250 * \param F the geom Face
251 * \param uvValues a vector of size elem->NbCornerNodes() to save the uv coordinate points on the face
252 * \return true if all the nodes are IN the face
254 //================================================================================
255 bool GMSHPlugin_Mesher::IsAllNodesInSameFace( const SMDS_MeshElement* element, const TopoDS_Face& F,
256 std::vector<gp_XY>& uvValues )
258 Handle(ShapeAnalysis_Surface) sprojector = new ShapeAnalysis_Surface( BRep_Tool::Surface( F ));
259 double tol = BRep_Tool::MaxTolerance( F, TopAbs_FACE );
260 int nbN = element->NbCornerNodes();
261 gp_Pnt surfPnt(0,0,0);
262 for ( int i = 0; i < nbN; ++i )
264 SMESH_NodeXYZ nXYZ( element->GetNode( i ) );
265 gp_XY uv = sprojector->ValueOfUV( nXYZ, tol ).XY();
266 surfPnt = sprojector->Value( uv );
267 double dist = surfPnt.Distance( nXYZ );
276 //================================================================================
278 * \brief Associate mesh elements to geometrical faces.
279 * \param the list of elements
280 * \return a map between faces (incremental order) and mesh elements found to be placed on the face
282 //================================================================================
283 std::map<int,std::vector<std::tuple<smIdType,bool,std::vector<gp_XY>>>> GMSHPlugin_Mesher::AssociateElementsToFaces( std::map<const SMDS_MeshElement*, bool, TIDCompare>& listElements )
285 // Map faces to elementId and uv of nodes.
287 // Index vector with element smIdType, [ gp_XY_0, gp_XY_1, gp_XY_2 ]
288 std::map<int,std::vector<std::tuple<smIdType,bool,std::vector<gp_XY>>>> elementToFaceMap;
289 SMESHDS_Mesh* meshDS = _mesh->GetMeshDS();
290 SMDS_ElemIteratorPtr iteratorElem = meshDS->elementsIterator(SMDSAbs_Face);
291 for ( auto const& [elem, IsReverse] : listElements ) // loop on elements on a geom face
293 int nbN = elem->NbCornerNodes();
294 std::vector<gp_XY> uvValues(nbN);
295 if ( nbN > 4 /*this restriction might be eliminated. Have to adapt FillGeomMapMeshUsing2DMeshIterator function too */)
296 throw std::string("Polygon sub-meshes not supported");
299 for( GModel::fiter it = _gModel->firstFace(); it != _gModel->lastFace(); ++it )
302 TopoDS_Face topoFace = *((TopoDS_Face*)gFace->getNativePtr());
303 if ( IsAllNodesInSameFace( elem, topoFace, uvValues ) )
305 elementToFaceMap[ faceId ].push_back( std::make_tuple( elem->GetID(), IsReverse, uvValues ) );
311 return elementToFaceMap;
314 //================================================================================
316 * \brief Add the elements found associated to the face as gmsh elements
318 //================================================================================
319 void GMSHPlugin_Mesher::Set2DMeshes( std::vector< const SMDS_MeshNode* >& nodeVec, std::map<const SMDS_MeshElement*, bool, TIDCompare>& listElements )
321 std::map< const SMDS_MeshNode* , const MVertex * > nodes2mvertMap;
322 SMESHDS_Mesh* meshDS = _mesh->GetMeshDS();
323 auto elementToFaceMap = AssociateElementsToFaces( listElements );
325 std::vector<MVertex *> mVertices;
327 for(GModel::fiter it = _gModel->firstFace(); it != _gModel->lastFace(); ++it)
332 auto element2uv = elementToFaceMap.find( faceId )->second;
334 const int numberOfEntries = element2uv.size();
336 for (int el = 0; el < numberOfEntries; el++)
338 const smIdType elementId = std::get<0>( element2uv[ el ] ); // smesh element id
339 bool isReverse = std::get<1>( element2uv[ el ] );
340 const SMDS_MeshElement* elem = meshDS->FindElement( elementId );
342 int nbN = elem->NbCornerNodes();
343 mVertices.resize( nbN );
345 for ( int i = 0; i < nbN; ++i )
347 const SMDS_MeshNode* n = elem->GetNode( i );
348 MVertex * mv = nullptr;
349 auto n2v = nodes2mvertMap.find( n );
350 if ( n2v != nodes2mvertMap.end() )
352 mv = const_cast< MVertex*>( n2v->second );
356 if ( n->GetPosition()->GetDim() < 2 )
357 throw std::string("Wrong mapping of edge nodes to GMSH nodes");
358 SMESH_NodeXYZ xyz = n;
359 gp_XY uv = std::get<2>(element2uv[ el ])[ i ];
360 mv = new MFaceVertex( xyz.X(), xyz.Y(), xyz.Z(), gFace, uv.X(), uv.Y() );
361 gFace->mesh_vertices.push_back( mv );
362 nodes2mvertMap.insert({ n, mv });
363 _nodeMap.insert ({ mv, n });
364 _premeshednodeMap.insert({ mv, n });
368 // create GMSH mesh faces
372 gFace->triangles.push_back (new MTriangle(mVertices[0], mVertices[2], mVertices[1]));
374 gFace->triangles.push_back (new MTriangle(mVertices[0], mVertices[1], mVertices[2]));
378 gFace->quadrangles.push_back (new MQuadrangle(mVertices[0], mVertices[3],
379 mVertices[2], mVertices[1]));
381 gFace->quadrangles.push_back (new MQuadrangle(mVertices[0], mVertices[1],
382 mVertices[2], mVertices[3]));
387 faceId++; // face counter
388 } // iterator in the face
391 nodeVec.resize( nodes2mvertMap.size() + 1, 0 );
393 for (auto k : nodes2mvertMap )
395 nodeVec[ count ] = k.first; // Index the node id to the smesh node itself
401 //================================================================================
403 * \brief Initialize GMSH model with mesh elements as geometry objects.
404 * Nodes are vertexes and element connections are geom lines
406 //================================================================================
407 void GMSHPlugin_Mesher::FillGeomMapMeshUsing2DMeshIterator( std::map<const SMDS_MeshElement*, bool, TIDCompare>& listElements )
410 gmsh::model::add("mesh");
412 typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
413 typedef TNodeToIDMap::value_type TN2ID;
414 typedef map<std::pair<int, int>, int> TLineToIDMap;
415 TNodeToIDMap aNodeToID;
416 TLineToIDMap aLineToID;
421 const int invalid_ID = -1;
422 std::vector<int> aTrinagle( 3, 0 );
424 // Playing around with SMESHDS_Mesh structure
425 SMESHDS_Mesh* meshDS = _mesh->GetMeshDS();
427 for ( auto const& [elem, IsReverse] : listElements ) // loop on elements on a geom face
429 if ( elem->NbCornerNodes() != 3 )
432 for (int iN = 0; iN < 3; ++iN)
434 const SMDS_MeshNode* aNode = elem->GetNode(iN);
436 int& ngID = aNodeToID.insert(TN2ID(aNode, invalid_ID)).first->second;
437 if (ngID == invalid_ID)
440 gmsh::model::occ::addPoint(aNode->X(), aNode->Y(), aNode->Z(), 1.e-2, ngID);
443 aTrinagle[ IsReverse ? 2 - iN : iN ] = ngID;
446 if ((aTrinagle[0] == aTrinagle[1] ||
447 aTrinagle[0] == aTrinagle[2] ||
448 aTrinagle[2] == aTrinagle[1]))
451 std::vector<int> LinesID(3, 0);
452 for (int anIndex = 0; anIndex < 3; ++anIndex)
454 int aNextIndex = (anIndex + 1) % 3;
455 if (aLineToID.find({ aTrinagle[anIndex], aTrinagle[aNextIndex] }) == aLineToID.end()
456 && aLineToID.find({ aTrinagle[aNextIndex], aTrinagle[anIndex] }) == aLineToID.end())
458 LinesID[anIndex] = aLineToID.insert({ { aTrinagle[aNextIndex], aTrinagle[anIndex] }, ++aNbOfLines }).first->second;
459 gmsh::model::occ::addLine(aTrinagle[anIndex], aTrinagle[aNextIndex], LinesID[anIndex]);
463 LinesID[anIndex] = aLineToID.find({ aTrinagle[anIndex], aTrinagle[aNextIndex] })->second;
464 if (LinesID[anIndex] == 0)
465 LinesID[anIndex] = aLineToID.find({ aTrinagle[aNextIndex], aTrinagle[anIndex] })->second;
469 // if (!aProxyMesh->IsTemporary(ls.first))
470 // swap(aTrinagle[1], aTrinagle[2]);
471 gmsh::model::occ::addCurveLoop(LinesID);
475 //================================================================================
477 * \brief Initialize GMSH model
479 //================================================================================
480 void GMSHPlugin_Mesher::FillGMSHMesh()
483 gmsh::model::add("mesh");
485 SMESHDS_Mesh* meshDS = _mesh->GetMeshDS();
489 std::vector<int> aTrinagle(3, 0);
491 const int invalid_ID = -1;
494 typedef map< const SMDS_MeshNode*, int, TIDCompare > TNodeToIDMap;
495 typedef TNodeToIDMap::value_type TN2ID;
496 typedef map<std::pair<int, int>, int> TLineToIDMap;
497 TNodeToIDMap aNodeToID;
498 TLineToIDMap aLineToID;
500 TopAbs_ShapeEnum aMainType = _mesh->GetShapeToMesh().ShapeType();
501 bool aCheckReverse = (aMainType == TopAbs_COMPOUND || aMainType == TopAbs_COMPSOLID);
503 SMESH_MesherHelper aHelper(*_mesh);
504 SMESH_ProxyMesh::Ptr aProxyMesh(new SMESH_ProxyMesh(*_mesh));
505 if (_mesh->NbQuadrangles() > 0)
507 StdMeshers_QuadToTriaAdaptor* Adaptor = new StdMeshers_QuadToTriaAdaptor;
508 Adaptor->Compute(*_mesh, _shape, aProxyMesh.get());
509 aProxyMesh.reset(Adaptor);
512 std::map<const SMDS_MeshElement*, bool, TIDCompare> listElements;
513 for (TopExp_Explorer exFa(_shape, TopAbs_FACE); exFa.More(); exFa.Next())
515 const TopoDS_Shape& aShapeFace = exFa.Current();
516 int faceID = meshDS->ShapeToIndex(aShapeFace);
518 if (aCheckReverse && aHelper.NbAncestors(aShapeFace, *_mesh, _shape.ShapeType()) > 1)
519 // IsReversedSubMesh() can work wrong on strongly curved faces,
520 // so we use it as less as possible
521 isRev = aHelper.IsReversedSubMesh(TopoDS::Face(aShapeFace));
523 const SMESHDS_SubMesh* aSubMeshDSFace = aProxyMesh->GetSubMesh(aShapeFace);
526 SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
527 if (aHelper.IsQuadraticSubMesh(_shape) &&
528 dynamic_cast<const SMESH_ProxyMesh::SubMesh*>(aSubMeshDSFace))
530 // add medium nodes of proxy triangles to helper
531 while (iteratorElem->more())
532 aHelper.AddTLinks(static_cast<const SMDS_MeshFace*>(iteratorElem->next()));
534 iteratorElem = aSubMeshDSFace->GetElements();
536 while (iteratorElem->more()) // loop on elements on a geom face
539 const SMDS_MeshElement* elem = iteratorElem->next();
542 if (elem->NbCornerNodes() != 3)
544 listElements[elem] = isRev;
548 for (auto const& ls : listElements)
550 // Add nodes of triangles and triangles them-selves to netgen mesh
551 // add three nodes of
552 bool hasDegen = false;
553 for (int iN = 0; iN < 3; ++iN)
555 const SMDS_MeshNode* aNode = ls.first->GetNode(iN);
556 const int shapeID = aNode->getshapeId();
557 if (aNode->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE &&
558 aHelper.IsDegenShape(shapeID))
560 // ignore all nodes on degeneraged edge and use node on its vertex instead
561 TopoDS_Shape vertex = TopoDS_Iterator(meshDS->IndexToShape(shapeID)).Value();
562 aNode = SMESH_Algo::VertexNode(TopoDS::Vertex(vertex), meshDS);
565 int& ngID = aNodeToID.insert(TN2ID(aNode, invalid_ID)).first->second;
566 if (ngID == invalid_ID)
569 gmsh::model::occ::addPoint(aNode->X(), aNode->Y(), aNode->Z(), 1.e-2, ngID);
571 aTrinagle[ls.second ? 2 - iN : iN] = ngID;
574 if (hasDegen && (aTrinagle[0] == aTrinagle[1] ||
575 aTrinagle[0] == aTrinagle[2] ||
576 aTrinagle[2] == aTrinagle[1]))
580 std::vector<int> LinesID(3, 0);
581 for (int anIndex = 0; anIndex < 3; ++anIndex)
583 int aNextIndex = (anIndex + 1) % 3;
584 if (aLineToID.find({ aTrinagle[anIndex], aTrinagle[aNextIndex] }) == aLineToID.end()
585 && aLineToID.find({ aTrinagle[aNextIndex], aTrinagle[anIndex] }) == aLineToID.end())
587 LinesID[anIndex] = aLineToID.insert({ { aTrinagle[aNextIndex], aTrinagle[anIndex] }, ++aNbOfLines }).first->second;
588 gmsh::model::occ::addLine(aTrinagle[anIndex], aTrinagle[aNextIndex], LinesID[anIndex]);
592 LinesID[anIndex] = aLineToID.find({ aTrinagle[anIndex], aTrinagle[aNextIndex] })->second;
593 if (LinesID[anIndex] == 0)
594 LinesID[anIndex] = aLineToID.find({ aTrinagle[aNextIndex], aTrinagle[anIndex] })->second;
598 if (!aProxyMesh->IsTemporary(ls.first))
599 swap(aTrinagle[1], aTrinagle[2]);
601 gmsh::model::occ::addCurveLoop(LinesID);
604 // Generate 1D and 2D mesh
605 _gModel->mesh( /*dim=*/ 1);
606 Set1DSubMeshes(_gModel);
607 _gModel->mesh( /*dim=*/ 2);
610 //================================================================================
612 * \brief Set Gmsh Options
614 //================================================================================
615 void GMSHPlugin_Mesher::SetGmshOptions()
617 MESSAGE("GMSHPlugin_Mesher::SetGmshOptions");
619 printf("We chose _algo2d %d \n", _algo2d );
620 printf("We chose _algo3d %d \n", _algo3d );
621 printf("We chose _recomb2DAlgo %d \n", _recomb2DAlgo );
622 printf("We chose _recombineAll %d \n", (_recombineAll)?1:0);
623 printf("We chose _subdivAlgo %d \n", _subdivAlgo );
624 printf("We chose _remeshAlgo %d \n", _remeshAlgo );
625 printf("We chose _remeshPara %d \n", _remeshPara );
626 printf("We chose _smouthSteps %e \n", _smouthSteps );
627 printf("We chose _sizeFactor %e \n", _sizeFactor );
628 printf("We chose _minSize %e \n", _minSize );
629 printf("We chose _maxSize %e \n", _maxSize );
630 printf("We chose _secondOrder %d \n", (_secondOrder)?1:0);
631 printf("We chose _useIncomplElem %d \n", (_useIncomplElem)?1:0);
632 printf("We are in dimension %d \n", (_is2d)?2:3);
635 std::map <int,double> mapAlgo2d;
636 mapAlgo2d[0]=2; // Automatic
637 mapAlgo2d[1]=1; // MeshAdapt
638 mapAlgo2d[2]=5; // Delaunay
639 mapAlgo2d[3]=6; // Frontal-Delaunay
640 mapAlgo2d[4]=8; // DelQuad (Frontal-Delaunay for Quads)
641 mapAlgo2d[5]=9; // Packing of parallelograms
642 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=10
643 mapAlgo2d[6]=11;// Quasistructured quads with cross-fields
646 std::map <int,double> mapAlgo3d;
647 mapAlgo3d[0]=1; // Delaunay
648 mapAlgo3d[1]=4; // Frontal
649 mapAlgo3d[2]=7; // MMG3D
650 mapAlgo3d[3]=9; // R-tree
651 mapAlgo3d[4]=10;// HXT
654 ok = GmshSetOption("Mesh", "Algorithm" , mapAlgo2d[_algo2d]) ;
658 ok = GmshSetOption("Mesh", "Algorithm3D" , mapAlgo3d[_algo3d]) ;
661 ok = GmshSetOption("Mesh", "RecombinationAlgorithm" , (double)_recomb2DAlgo) ;
663 ok = GmshSetOption("Mesh", "RecombineAll" , (_recombineAll)?1.:0.) ;
665 ok = GmshSetOption("Mesh", "SubdivisionAlgorithm" , (double)_subdivAlgo) ;
667 ok = GmshSetOption("Mesh", "RemeshAlgorithm" , (double)_remeshAlgo) ;
669 ok = GmshSetOption("Mesh", "RemeshParametrization" , (double)_remeshPara) ;
671 ok = GmshSetOption("Mesh", "Smoothing" , (double)_smouthSteps) ;
673 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=10
674 ok = GmshSetOption("Mesh", "MeshSizeFromCurvature" , _meshCurvatureSize) ;
677 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
678 ok = GmshSetOption("Mesh", "MeshSizeFactor" , _sizeFactor) ;
680 ok = GmshSetOption("Mesh", "MeshSizeMin" , _minSize) ;
682 ok = GmshSetOption("Mesh", "MeshSizeMax" , _maxSize) ;
685 ok = GmshSetOption("Mesh", "CharacteristicLengthFactor" , _sizeFactor) ;
687 ok = GmshSetOption("Mesh", "CharacteristicLengthMin" , _minSize) ;
689 ok = GmshSetOption("Mesh", "CharacteristicLengthMax" , _maxSize) ;
692 ok = GmshSetOption("Mesh", "ElementOrder" , (_secondOrder)?2.:1.) ;
696 ok = GmshSetOption("Mesh", "SecondOrderIncomplete" ,(_useIncomplElem)?1.:0.);
700 ok = GmshSetOption("General", "Verbosity" , (double) _verbLvl ) ; // Verbosity level
703 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
704 /*ok = GmshSetOption("Mesh", "MaxNumThreads1D" , 0. ) ; // Coarse-grain algo threads
706 ok = GmshSetOption("Mesh", "MaxNumThreads2D" , 0. ) ; // Coarse-grain algo threads
708 ok = GmshSetOption("Mesh", "MaxNumThreads3D" , 0. ) ; // Fine-grain algo threads HXT
711 ok = GmshSetOption("General", "NumThreads" , _maxThreads ) ; // system default i.e. OMP_NUM_THREADS
714 if ( GMSHPlugin_Hypothesis::Algo3D::hxt == _algo3d ){
715 MESSAGE("GMSHPlugin_Mesher::SetGmshOptions: HXT algorithm is being used. Setting number of threads to 1.");
716 ok = GmshSetOption("Mesh", "MaxNumThreads3D" , 1. );
723 //================================================================================
725 * \brief Create and add Compounds into GModel _gModel.
727 //================================================================================
729 void GMSHPlugin_Mesher::CreateGmshCompounds()
731 MESSAGE("GMSHPlugin_Mesher::CreateGmshCompounds");
733 SMESH_Gen_i* smeshGen_i = SMESH_Gen_i::GetSMESHGen();
735 OCC_Internals* occgeo = _gModel->getOCCInternals();
736 bool toSynchronize = false;
738 for(std::set<std::string>::const_iterator its = _compounds.begin();its != _compounds.end(); ++its )
740 GEOM::GEOM_Object_var aGeomObj;
741 TopoDS_Shape geomShape = TopoDS_Shape();
742 SALOMEDS::SObject_var aSObj = SMESH_Gen_i::GetSMESHGen()->getStudyServant()->FindObjectID( (*its).c_str() );
743 SALOMEDS::GenericAttribute_var anAttr;
744 if (!aSObj->_is_nil() && aSObj->FindAttribute(anAttr, "AttributeIOR"))
746 SALOMEDS::AttributeIOR_var anIOR = SALOMEDS::AttributeIOR::_narrow(anAttr);
747 CORBA::String_var aVal = anIOR->Value();
748 CORBA::Object_var obj = SMESH_Gen_i::GetSMESHGen()->getStudyServant()->ConvertIORToObject(aVal);
749 aGeomObj = GEOM::GEOM_Object::_narrow(obj);
751 geomShape = smeshGen_i->GeomObjectToShape( aGeomObj.in() );
752 if ( geomShape.IsNull() )
755 TopAbs_ShapeEnum geomType = geomShape.ShapeType();
756 if ( geomType == TopAbs_COMPOUND)// voir s'il ne faut pas mettre une erreur dans le cas contraire
758 MESSAGE("shapeType == TopAbs_COMPOUND");
759 TopoDS_Iterator it(geomShape);
762 TopAbs_ShapeEnum shapeType = it.Value().ShapeType();
763 std::vector< std::pair< int, int > > dimTags;
764 for ( ; it.More(); it.Next())
766 const TopoDS_Shape& topoShape = it.Value();
767 ASSERT(topoShape.ShapeType() == shapeType);
768 if ( _mesh->GetMeshDS()->ShapeToIndex( topoShape ) > 0 )
769 occgeo->importShapes( &topoShape, false, dimTags );
772 TopoDS_Shape face = TopExp_Explorer( _shape, shapeType ).Current();
773 SMESH_subMesh* sm = _mesh->GetSubMesh( face );
774 sm->GetComputeError() =
775 SMESH_ComputeError::New
776 ( COMPERR_WARNING, "Compound shape does not belong to the main geometry. Ignored");
779 std::vector<int> tags;
780 int dim = ( shapeType == TopAbs_EDGE ) ? 1 : 2;
781 for ( size_t i = 0; i < dimTags.size(); ++i )
783 if ( dimTags[i].first == dim )
784 tags.push_back( dimTags[i].second );
788 _gModel->getGEOInternals()->setCompoundMesh( dim, tags );
789 toSynchronize = true;
792 _gModel->getGEOInternals()->synchronize(_gModel);
797 //================================================================================
799 * \brief For a compound mesh set the mesh components to be transmitted to SMESH
801 //================================================================================
803 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
804 void GMSHPlugin_Mesher::SetCompoundMeshVisibility()
807 // Loop over all faces, if the face belongs to a compound entry then
808 // for all (boundary) edges whithin the face visibility is set to 0,
809 // if the face doesn't belong to a compound entry then visibility is
810 // set to 1 for all its (boundary) edges. Later, in FillSMesh() func
811 // getVisibility() (returns either 1 or 0) is used to decide weather
812 // the mesh of edge should be transmitted to SMESH or not.
814 for ( GModel::fiter itF = _gModel->firstFace(); itF != _gModel->lastFace(); ++itF )
816 std::vector< GEdge *> faceEdges = (*itF)->edges();
818 for ( auto itE = faceEdges.begin(); itE != faceEdges.end(); ++itE )
820 if ( ((*itF)->compound.size()) )
821 (*itE)->setVisibility(0);
823 (*itE)->setVisibility(1);
828 // Loop over all edges, if the edge belongs to a compound entry then
829 // for all (boundary) vertices whithin the edge visibility is set to
830 // 0, if the edge doesn't belong to a compound entry then visibility
831 // is set to 1 for all its (boundary) vertices. Later, in FillSMesh()
832 // func getVisibility() (returns either 1 or 0) is used to decide we-
833 // ather the mesh of vertices should be transmitted to SMESH or not.
835 for ( GModel::eiter itE = _gModel->firstEdge(); itE != _gModel->lastEdge(); ++itE )
837 std::vector<GVertex *> bndVerticies = (*itE)->vertices();
839 for( auto itV = bndVerticies.begin(); itV != bndVerticies.end(); ++itV )
841 if(((*itE)->compound.size()))
842 (*itV)->setVisibility(0);
844 (*itV)->setVisibility(1);
851 //================================================================================
853 * \brief Get a node by a GMSH mesh vertex
855 //================================================================================
857 const SMDS_MeshNode* GMSHPlugin_Mesher::Node( const MVertex* v )
859 std::map< const MVertex *, const SMDS_MeshNode* >::iterator v2n = _nodeMap.find( v );
860 if ( v2n != _nodeMap.end() )
866 //================================================================================
868 * \brief Get a node by a GMSH mesh vertex
870 //================================================================================
872 const SMDS_MeshNode* GMSHPlugin_Mesher::PremeshedNode( const MVertex* v )
874 std::map< const MVertex *, const SMDS_MeshNode* >::iterator v2n = _premeshednodeMap.find( v );
875 if ( v2n != _premeshednodeMap.end() )
882 //================================================================================
884 * \brief Return a corresponding sub-mesh if a shape is meshed
886 //================================================================================
888 SMESHDS_SubMesh* GMSHPlugin_Mesher::HasSubMesh( const TopoDS_Shape& s )
890 if ( SMESHDS_SubMesh* sm = _mesh->GetMeshDS()->MeshElements( s ))
892 if ( s.ShapeType() == TopAbs_VERTEX )
893 return ( sm->NbNodes() > 0 ) ? sm : nullptr;
895 return ( sm->NbElements() > 0 ) ? sm : nullptr;
900 //================================================================================
902 * \brief Write mesh from GModel instance to SMESH instance
904 //================================================================================
906 void GMSHPlugin_Mesher::FillSMesh()
908 SMESHDS_Mesh* meshDS = _mesh->GetMeshDS();
911 for ( GModel::viter it = _gModel->firstVertex(); it != _gModel->lastVertex(); ++it)
913 GVertex *gVertex = *it;
915 // GET topoVertex CORRESPONDING TO gVertex
916 TopoDS_Vertex topoVertex = *((TopoDS_Vertex*)gVertex->getNativePtr());
918 if (gVertex->getVisibility() == 0) // belongs to a compound
920 SMESH_subMesh* sm = _mesh->GetSubMesh(topoVertex);
921 sm->SetIsAlwaysComputed(true); // prevent from displaying errors
925 // FILL SMESH FOR topoVertex
927 for( size_t i = 0; i < gVertex->mesh_vertices.size(); i++)
929 MVertex *v = gVertex->mesh_vertices[i];
930 if(v->getIndex() >= 0)
932 if ( SMESHDS_SubMesh* sm = HasSubMesh( topoVertex ))
934 const SMDS_MeshNode *node = sm->GetNodes()->next();
935 _nodeMap.insert({ v, node });
939 SMDS_MeshNode *node = meshDS->AddNode( v->x(),v->y(),v->z() );
940 meshDS->SetNodeOnVertex( node, topoVertex );
941 _nodeMap.insert({ v, node });
945 // WE DON'T ADD 0D ELEMENTS because it does not follow the salome meshers philosophy
947 // for(unsigned int i = 0; i < gVertex->getNumMeshElements(); i++)
949 // MElement *e = gVertex->getMeshElement(i);
950 // std::vector<MVertex*> verts;
951 // e->getVertices(verts);
952 // ASSERT(verts.size()==1);
953 // SMDS_Mesh0DElement* zeroDElement = 0;
954 // zeroDElement = meshDS->Add0DElementWithID(verts[0]->getNum(),e->getNum());
955 // meshDS->SetMeshElementOnShape(zeroDElement, topoVertex);
960 for(GModel::eiter it = _gModel->firstEdge(); it != _gModel->lastEdge(); ++it)
964 // GET topoEdge CORRESPONDING TO gEdge
965 TopoDS_Edge topoEdge;
966 std::vector< ShapeBounds > topoEdges;
967 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
968 if(gEdge->haveParametrization())
970 if ( gEdge->geomType() != GEntity::CompoundCurve )
973 topoEdge = *((TopoDS_Edge*)gEdge->getNativePtr());
974 if (gEdge->getVisibility() == 0) // belongs to a compound
976 SMESH_subMesh* sm = _mesh->GetSubMesh(topoEdge);
977 sm->SetIsAlwaysComputed(true); // prevent from displaying errors
980 if ( HasSubMesh( topoEdge ))
981 continue; // a meshed sub-mesh
983 bool isCompound = getBoundsOfShapes( gEdge, topoEdges );
985 // FILL SMESH FOR topoEdge
987 for ( size_t i = 0; i < gEdge->mesh_vertices.size(); i++ )
989 MVertex *v = gEdge->mesh_vertices[i];
990 if ( v->getIndex() >= 0 )
992 SMDS_MeshNode *node = meshDS->AddNode( v->x(),v->y(),v->z() );
995 topoEdge = TopoDS::Edge( getShapeAtPoint( v->point(), topoEdges ));
997 // Based on BLSURFPlugin_BLSURF
998 gp_Pnt point3D( v->x(),v->y(),v->z() );
999 Standard_Real p0 = 0.0;
1000 Standard_Real p1 = 1.0;
1001 TopLoc_Location loc;
1002 Handle(Geom_Curve) curve = BRep_Tool::Curve(topoEdge, loc, p0, p1);
1004 if ( !curve.IsNull() )
1006 if ( !loc.IsIdentity() )
1007 point3D.Transform( loc.Transformation().Inverted() );
1009 GeomAPI_ProjectPointOnCurve proj(point3D, curve, p0, p1);
1012 if ( proj.NbPoints() > 0 )
1013 pa = (double)proj.LowerDistanceParameter();
1015 meshDS->SetNodeOnEdge( node, topoEdge, pa );
1019 meshDS->SetNodeOnEdge( node, topoEdge );
1021 //END on BLSURFPlugin_BLSURF
1024 _nodeMap.insert({ v, node });
1029 for ( GModel::eiter it = _gModel->firstEdge(); it != _gModel->lastEdge(); ++it )
1032 if ( gEdge->getVisibility() == 0) // belongs to a compound
1035 TopoDS_Edge topoEdge;
1036 std::vector< ShapeBounds > topoEdges;
1037 bool isCompound = getBoundsOfShapes( gEdge, topoEdges );
1039 topoEdge = *((TopoDS_Edge*)gEdge->getNativePtr());
1041 if ( HasSubMesh( topoEdge ))
1042 continue; // a meshed sub-mesh
1045 std::vector<MVertex*> verts(3);
1046 for ( size_t i = 0; i < gEdge->getNumMeshElements(); i++ )
1048 MElement *e = gEdge->getMeshElement(i);
1050 e->getVertices(verts);
1052 // if a node wasn't set, it is assigned here
1053 for ( size_t j = 0; j < verts.size(); j++ )
1055 if ( verts[j]->onWhat()->getVisibility() == 0 )
1057 SMDS_MeshNode *node = meshDS->AddNode(verts[j]->x(),verts[j]->y(),verts[j]->z() );
1059 gp_Pnt point3D( verts[j]->x(),verts[j]->y(),verts[j]->z() );
1060 Standard_Real p0 = 0.0;
1061 Standard_Real p1 = 1.0;
1062 TopLoc_Location loc;
1063 Handle(Geom_Curve) curve = BRep_Tool::Curve(topoEdge, loc, p0, p1);
1065 if ( !curve.IsNull() )
1067 if ( !loc.IsIdentity() )
1068 point3D.Transform( loc.Transformation().Inverted() );
1070 GeomAPI_ProjectPointOnCurve proj(point3D, curve, p0, p1);
1073 if ( proj.NbPoints() > 0 )
1074 pa = (double)proj.LowerDistanceParameter();
1076 meshDS->SetNodeOnEdge( node, topoEdge, pa );
1080 meshDS->SetNodeOnEdge( node, topoEdge );
1083 verts[j]->setEntity(gEdge);
1084 _nodeMap.insert({ verts[j], node });
1088 SMDS_MeshEdge* edge = 0;
1089 switch (verts.size())
1092 edge = meshDS->AddEdge(Node( verts[0]),
1096 edge = meshDS->AddEdge(Node( verts[0]),
1105 topoEdge = TopoDS::Edge( getShapeAtPoint( e->barycenter(), topoEdges ));
1107 meshDS->SetMeshElementOnShape( edge, topoEdge );
1112 for ( GModel::fiter it = _gModel->firstFace(); it != _gModel->lastFace(); ++it)
1116 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
1117 // Gmsh since version 4.3 is now producing extra surface and mesh when
1118 // compounds are involved. Since in Gmsh meshing procedure needs acess
1119 // to each of the original topology and the meshed topology. Hence we
1120 // bypass the additional mesh in case of compounds. Note, similar cri-
1121 // teria also occurs in the following 'for' loop.
1122 if ( _compounds.size() && gFace->geomType() == GEntity::DiscreteSurface )
1126 // GET topoFace CORRESPONDING TO gFace
1127 TopoDS_Face topoFace;
1128 std::vector< ShapeBounds > topoFaces;
1130 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
1131 if(gFace->haveParametrization())
1133 if ( gFace->geomType() != GEntity::CompoundSurface )
1136 topoFace = *((TopoDS_Face*)gFace->getNativePtr());
1137 if (gFace->getVisibility() == 0) // belongs to a compound
1139 SMESH_subMesh* sm = _mesh->GetSubMesh(topoFace);
1140 sm->SetIsAlwaysComputed(true); // prevent from displaying errors
1143 if ( HasSubMesh( topoFace ))
1144 continue; // a meshed sub-mesh
1146 bool isCompound = getBoundsOfShapes( gFace, topoFaces );
1148 // FILL SMESH FOR topoFace
1150 for ( size_t i = 0; i < gFace->mesh_vertices.size(); i++ )
1152 MVertex *v = gFace->mesh_vertices[i];
1153 if ( v->getIndex() >= 0 )
1155 SMDS_MeshNode *node = meshDS->AddNode( v->x(),v->y(),v->z() );
1158 topoFace = TopoDS::Face( getShapeAtPoint( v->point(), topoFaces ));
1160 meshDS->SetNodeOnFace( node, topoFace );
1161 _nodeMap.insert({ v, node });
1166 for ( GModel::fiter it = _gModel->firstFace(); it != _gModel->lastFace(); ++it)
1170 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
1171 if ( _compounds.size() && gFace->geomType() == GEntity::DiscreteSurface )
1174 bool isCompound = (!gFace->haveParametrization());
1176 bool isCompound = ( gFace->geomType() == GEntity::CompoundSurface );
1179 if ( !isCompound && gFace->getVisibility() == 0 )
1180 continue; // belongs to a compound
1182 TopoDS_Face topoFace;
1183 std::vector< ShapeBounds > topoFaces;
1185 getBoundsOfShapes( gFace, topoFaces );
1187 topoFace = *((TopoDS_Face*)gFace->getNativePtr());
1189 if ( HasSubMesh( topoFace ))
1190 continue; // a meshed sub-mesh
1193 std::vector<MVertex*> verts;
1194 for ( size_t i = 0; i < gFace->getNumMeshElements(); i++ )
1196 MElement *e = gFace->getMeshElement(i);
1198 e->getVertices(verts);
1199 SMDS_MeshFace* face = 0;
1201 // if a node wasn't set, it is assigned here
1202 for ( size_t j = 0; j < verts.size(); j++)
1204 if(verts[j]->onWhat()->getVisibility() == 0)
1206 SMDS_MeshNode *node = meshDS->AddNode(verts[j]->x(),verts[j]->y(),verts[j]->z());
1207 meshDS->SetNodeOnFace( node, topoFace );
1208 _nodeMap.insert({ verts[j], node });
1209 verts[j]->setEntity(gFace);
1212 switch (verts.size())
1215 face = meshDS->AddFace(Node( verts[0]),
1220 face = meshDS->AddFace(Node( verts[0]),
1226 face = meshDS->AddFace(Node( verts[0]),
1234 face = meshDS->AddFace(Node( verts[0]),
1244 face = meshDS->AddFace(Node( verts[0]),
1260 topoFace = TopoDS::Face( getShapeAtPoint( e->barycenter(), topoFaces ));
1262 meshDS->SetMeshElementOnShape(face, topoFace);
1267 for ( GModel::riter it = _gModel->firstRegion(); it != _gModel->lastRegion(); ++it)
1269 GRegion *gRegion = *it;
1270 if (gRegion->getVisibility() == 0)
1273 // GET topoSolid CORRESPONDING TO gRegion
1274 TopoDS_Solid topoSolid = *((TopoDS_Solid*)gRegion->getNativePtr());
1276 // FILL SMESH FOR topoSolid
1279 for( size_t i = 0; i < gRegion->mesh_vertices.size(); i++)
1281 MVertex *v = gRegion->mesh_vertices[i];
1282 if(v->getIndex() >= 0)
1284 SMDS_MeshNode *node = meshDS->AddNode( v->x(),v->y(),v->z() );
1285 meshDS->SetNodeInVolume( node, topoSolid );
1286 _nodeMap.insert({ v, node });
1291 std::vector<MVertex*> verts;
1292 for( size_t i = 0; i < gRegion->getNumMeshElements(); i++)
1294 MElement *e = gRegion->getMeshElement(i);
1296 e->getVertices(verts);
1297 SMDS_MeshVolume* volume = 0;
1298 switch (verts.size()){
1300 volume = meshDS->AddVolume(Node( verts[0]),
1306 volume = meshDS->AddVolume(Node( verts[0]),
1313 volume = meshDS->AddVolume(Node( verts[0]),
1321 volume = meshDS->AddVolume(Node( verts[0]),
1331 volume = meshDS->AddVolume(Node( verts[0]),
1343 volume = meshDS->AddVolume(Node( verts[0] ),
1357 case 14: // same as case 13, because no pyra14 in smesh
1358 volume = meshDS->AddVolume(Node( verts[0] ),
1373 volume = meshDS->AddVolume(Node( verts[0] ),
1389 case 18: // same as case 15, because no penta18 in smesh
1390 volume = meshDS->AddVolume(Node( verts[0] ),
1407 volume = meshDS->AddVolume(Node( verts[0] ),
1429 volume = meshDS->AddVolume(Node( verts[0] ),
1461 meshDS->SetMeshElementOnShape(volume, topoSolid);
1468 //================================================================================
1470 * \brief Find if SPoint point is in SBoundingBox3d bounds
1472 //================================================================================
1474 float GMSHPlugin_Mesher::DistBoundingBox(const SBoundingBox3d& bounds, const SPoint3& point)
1476 SPoint3 min = bounds.min();
1477 SPoint3 max = bounds.max();
1481 if (point.x() < min.x())
1482 x = min.x()-point.x();
1483 else if (point.x() > max.x())
1484 x = point.x()-max.x();
1488 if (point.y() < min.y())
1489 y = min.y()-point.y();
1490 else if (point.y() > max.y())
1491 y = point.y()-max.y();
1495 if (point.z() < min.z())
1496 z = min.z()-point.z();
1497 else if (point.z() > max.z())
1498 z = point.z()-max.z();
1504 //================================================================================
1506 * \brief Reimplemented GmshMessage call. Actions done if errors occurs
1507 * during gmsh meshing. We define here what to display and what to do.
1509 //================================================================================
1510 void GMSHPlugin_Mesher::mymsg::operator()(std::string level, std::string msg)
1512 //MESSAGE("level="<< level.c_str() << ", msg=" << msg.c_str()<< "\n");
1513 printf("level=%s msg=%s\n", level.c_str(), msg.c_str());
1515 if(level == "Fatal" || level == "Error")
1517 std::ostringstream oss;
1518 if (level == "Fatal")
1519 oss << "Fatal error during Generation of Gmsh Mesh\n";
1521 oss << "Error during Generation of Gmsh Mesh\n";
1522 oss << " " << msg.c_str() << "\n";
1523 GEntity *e = _gModel->getCurrentMeshEntity();
1526 oss << " error occurred while meshing entity:\n" <<
1527 " tag=" << e->tag() << "\n" <<
1528 " dimension=" << e->dim() << "\n" <<
1529 " native pointer=" << e->getNativePtr();
1530 //if(e->geomType() != GEntity::CompoundCurve and e->geomType() != GEntity::CompoundSurface)
1532 //SMESH_subMesh *sm = _mesh->GetSubMesh(*((TopoDS_Shape*)e->getNativePtr()));
1533 //SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1534 //SMESH_Comment comment;
1535 //comment << SMESH_Comment(oss.str);
1536 //std::string str = oss.str();
1537 //smError.reset( new SMESH_ComputeError( str ));
1539 // plutot que de faire de la merde ici, on pourait simplement
1540 // remplir une liste pour dire sur quelles entités gmsh se plante
1541 // (puis faire le fillsmesh)
1542 // puis faire une nouvelle routine qui réécrit les messages d'erreur
1543 // probleme : gmsh peut planté en Fatal, dans ce cas pas de fillsmesh
1546 if (level == "Fatal")
1548 CTX::instance()->lock = 0;
1552 printf("%s\n", oss.str().c_str());
1556 bool GMSHPlugin_Mesher::Compute3D( std::vector< const SMDS_MeshNode* >& nodeVec, std::map<const SMDS_MeshElement*, bool, TIDCompare>& listElements, bool addElements )
1558 MESSAGE("GMSHPlugin_Mesher::Compute3D");
1561 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=3
1565 char* argv[] = {"-noenv"};
1566 GmshInitialize(1,argv);
1568 _gModel = new GModel();
1570 GmshSetMessageHandler(&msg);
1572 _gModel->importOCCShape((void*)&_shape);
1575 HideComputedEntities( _gModel, true );
1576 // fill geometry with elements as geom objects
1577 FillGeomMapMeshUsing2DMeshIterator( listElements );
1578 Set2DMeshes( nodeVec, listElements );
1579 _gModel->mesh( /*dim=*/ 3);
1581 catch (std::string& str)
1589 MESSAGE("Unrecoverable error during Generation of Gmsh Mesh");
1594 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
1595 if (_compounds.size() > 0)
1596 SetCompoundMeshVisibility();
1602 MESSAGE("GMSHPlugin_Mesher::Compute3D:End");
1606 void GMSHPlugin_Mesher::finalizeGModel()
1610 GmshSetMessageHandler(nullptr);
1616 //=============================================================================
1618 * Here we are going to use the GMSH mesher
1620 //=============================================================================
1622 bool GMSHPlugin_Mesher::Compute()
1624 MESSAGE("GMSHPlugin_Mesher::Compute");
1628 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
1629 SetMaxThreadsGmsh();
1631 //RNV: to avoid modification of PATH and PYTHONPATH
1632 char* argv[] = {"-noenv"};
1633 GmshInitialize(1,argv);
1635 _gModel = new GModel();
1637 GmshSetMessageHandler(&msg);
1638 _gModel->importOCCShape((void*)&_shape);
1639 if (_compounds.size() > 0) CreateGmshCompounds();
1643 HideComputedEntities(_gModel);
1647 Set2DSubMeshes(_gModel);
1648 _gModel->mesh( /*dim=*/ 3);
1652 //CTX::instance()->mesh.maxNumThreads1D=1;
1653 _gModel->mesh( /*dim=*/ 1);
1655 Set1DSubMeshes(_gModel);
1657 //_gModel->writeUNV("/tmp/1D.unv", 1,0,0);
1658 //CTX::instance()->mesh.maxNumThreads2D=1;
1660 _gModel->mesh( /*dim=*/ 2);
1664 Set2DSubMeshes(_gModel);
1666 //CTX::instance()->mesh.maxNumThreads3D=1;
1668 _gModel->mesh( /*dim=*/ 3);
1670 RestoreVisibility(_gModel);
1672 #ifdef WITH_SMESH_CANCEL_COMPUTE
1676 catch (std::string& str)
1679 std::cerr << "GMSH: exception caught: " << str << std::endl;
1685 std::cerr << "GMSH: Unknown exception caught: " << std::endl;
1686 MESSAGE("Unrecoverable error during Generation of Gmsh Mesh");
1691 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
1692 if (_compounds.size() > 0)
1693 SetCompoundMeshVisibility();
1697 GmshSetMessageHandler(nullptr);
1700 MESSAGE("GMSHPlugin_Mesher::Compute:End");
1704 //================================================================================
1706 * \brief Set 1D sub-meshes to GModel. GMSH 1D mesh is made by now.
1707 * \param [inout] _gModel - GMSH model
1709 //================================================================================
1711 void GMSHPlugin_Mesher::Set1DSubMeshes( GModel* gModel )
1713 SMESHDS_Mesh* meshDS = _mesh->GetMeshDS();
1715 for(GModel::eiter it = gModel->firstEdge(); it != gModel->lastEdge(); ++it)
1719 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
1720 if ( !gEdge->haveParametrization())
1722 if ( gEdge->geomType() == GEntity::CompoundCurve )
1726 TopoDS_Edge topoEdge = *((TopoDS_Edge*)gEdge->getNativePtr());
1727 if ( !HasSubMesh( topoEdge ))
1728 continue; // empty sub-mesh
1730 gEdge->deleteMesh();
1732 // get node parameters on topoEdge
1733 StdMeshers_FaceSide side( TopoDS_Face(), topoEdge, _mesh, /*fwd=*/true, /*skpMedium=*/true);
1734 const UVPtStructVec& nodeParam = side.GetUVPtStruct();
1735 if ( nodeParam.empty() )
1736 throw std::string("Pb with StdMeshers_FaceSide::GetUVPtStruct()");
1738 // get GMSH mesh vertices on VERTEX'es
1739 std::vector<MVertex *> mVertices( nodeParam.size(), nullptr );
1740 GVertex * gV0 = gEdge->getBeginVertex(), *gV1 = gEdge->getEndVertex();
1741 mVertices[0] = gV0->mesh_vertices[ 0 ];
1742 mVertices.back() = gV1->mesh_vertices[ 0 ];
1743 TopoDS_Vertex v01 = *((TopoDS_Vertex*) gV0->getNativePtr());
1744 TopoDS_Shape v02 = SMESH_MesherHelper::GetSubShapeByNode( nodeParam[0].node, meshDS );
1745 bool reverse = !v01.IsSame( v02 );
1746 if ( mVertices[0] == mVertices.back() )
1747 reverse = ( nodeParam[0].param > nodeParam.back().param );
1748 const SMDS_MeshNode* n0 = reverse ? nodeParam.back().node : nodeParam[0].node;
1749 const SMDS_MeshNode* n1 = reverse ? nodeParam[0].node : nodeParam.back().node;
1750 _nodeMap.insert({ mVertices[ 0 ], n0 });
1751 _nodeMap.insert({ mVertices.back(), n1 });
1753 // create GMSH mesh vertices on gEdge
1754 for ( size_t i = 1; i < nodeParam.size() - 1; ++i )
1756 size_t iN = reverse ? ( nodeParam.size() - 1 - i ) : i;
1757 SMESH_NodeXYZ xyz = nodeParam[ iN ].node;
1758 double lc = segmentSize( nodeParam, iN );
1759 // SVector3 der = gEdge->firstDer(nodeParam[ iN ].param);
1760 // double lc = norm(der) / segmentSize( nodeParam, i );
1762 mVertices[ i ] = new MEdgeVertex( xyz.X(), xyz.Y(), xyz.Z(),
1763 gEdge, nodeParam[ iN ].param, 0, lc);
1764 gEdge->mesh_vertices.push_back( mVertices[ i ]);
1765 _nodeMap.insert({ mVertices[ i ], nodeParam[ iN ].node });
1767 // create GMSH mesh edges
1768 for ( size_t i = 1; i < mVertices.size(); ++i )
1770 gEdge->lines.push_back( new MLine( mVertices[ i - 1 ],
1774 cout << endl << "EDGE " << gEdge->tag() <<
1775 ( topoEdge.Orientation() == TopAbs_FORWARD ? " F" : " R") << endl;
1776 MVertex* mv = gV0->mesh_vertices[ 0 ];
1777 cout << "V0: " << mv->x() << ", " << mv->y() << ", " << mv->z() << ", "<<endl;
1778 for ( size_t i = 0; i < gEdge->mesh_vertices.size(); ++i )
1780 MEdgeVertex* mv = (MEdgeVertex*) gEdge->mesh_vertices[i];
1781 cout << i << ": " << mv->x() << ", " << mv->y() << ", " << mv->z() << ", ";
1783 mv->getParameter(0, t );
1784 cout << ":\t t = "<< t << " lc = " << mv->getLc() << endl;
1786 mv = gV1->mesh_vertices[ 0 ];
1787 cout << "V1: " << mv->x() << ", " << mv->y() << ", " << mv->z() << ", "<<endl;
1793 //================================================================================
1795 * \brief Set 2D sub-meshes to GModel. GMSH 2D mesh is made by now.
1796 * \param [inout] _gModel - GMSH model
1798 //================================================================================
1800 void GMSHPlugin_Mesher::Set2DSubMeshes( GModel* gModel )
1802 if ( _nodeMap.empty() )
1803 return; // no sub-meshes
1805 SMESH_MesherHelper helper( *_mesh );
1807 std::map< const SMDS_MeshNode* , const MVertex * > nodes2mvertMap;
1808 for ( auto & v2n : _nodeMap )
1809 nodes2mvertMap.insert({ v2n.second, v2n.first });
1811 std::vector<MVertex *> mVertices;
1813 for(GModel::fiter it = gModel->firstFace(); it != gModel->lastFace(); ++it)
1817 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
1818 if ( !gFace->haveParametrization())
1820 if ( gFace->geomType() == GEntity::CompoundSurface )
1824 TopoDS_Face topoFace = *((TopoDS_Face*)gFace->getNativePtr());
1825 SMESHDS_SubMesh* sm = HasSubMesh( topoFace );
1828 //_gModel->writeUNV("/tmp/befDEl.unv", 1,0,0);
1830 gFace->deleteMesh();
1832 bool reverse = false;
1833 if ( gFace->getRegion(0) )
1835 //GRegion * gRegion = gFace->getRegion(0);
1836 TopoDS_Shape topoSolid = *((TopoDS_Shape*)gFace->getNativePtr());
1837 TopAbs_Orientation faceOriInSolid = helper.GetSubShapeOri( topoSolid, topoFace );
1838 if ( faceOriInSolid >= 0 )
1840 helper.IsReversedSubMesh( TopoDS::Face( topoFace.Oriented( faceOriInSolid )));
1843 for ( SMDS_ElemIteratorPtr fIt = sm->GetElements(); fIt->more(); )
1845 const SMDS_MeshElement* f = fIt->next();
1847 int nbN = f->NbCornerNodes();
1849 throw std::string("Polygon sub-meshes not supported");
1851 mVertices.resize( nbN );
1852 for ( int i = 0; i < nbN; ++i )
1854 const SMDS_MeshNode* n = f->GetNode( i );
1855 MVertex * mv = nullptr;
1856 auto n2v = nodes2mvertMap.find( n );
1857 if ( n2v != nodes2mvertMap.end() )
1859 mv = const_cast< MVertex*>( n2v->second );
1863 if ( n->GetPosition()->GetDim() < 2 )
1864 throw std::string("Wrong mapping of edge nodes to GMSH nodes");
1865 SMESH_NodeXYZ xyz = n;
1867 gp_XY uv = helper.GetNodeUV( topoFace, n, nullptr, &ok );
1868 mv = new MFaceVertex( xyz.X(), xyz.Y(), xyz.Z(), gFace, uv.X(), uv.Y() );
1869 gFace->mesh_vertices.push_back( mv );
1870 nodes2mvertMap.insert({ n, mv });
1871 _nodeMap.insert ({ mv, n });
1873 mVertices[ i ] = mv;
1875 // create GMSH mesh faces
1879 gFace->triangles.push_back (new MTriangle(mVertices[0], mVertices[2], mVertices[1]));
1881 gFace->triangles.push_back (new MTriangle(mVertices[0], mVertices[1], mVertices[2]));
1885 gFace->quadrangles.push_back (new MQuadrangle(mVertices[0], mVertices[3],
1886 mVertices[2], mVertices[1]));
1888 gFace->quadrangles.push_back (new MQuadrangle(mVertices[0], mVertices[1],
1889 mVertices[2], mVertices[3]));
1894 } // loop on GMSH faces
1899 //================================================================================
1901 * \brief Set visibility 0 to already computed geom entities
1902 * to prevent their meshing
1904 //================================================================================
1906 void GMSHPlugin_Mesher::HideComputedEntities( GModel* gModel, bool hideAnyway )
1908 CTX::instance()->mesh.meshOnlyVisible = true;
1910 for(GModel::eiter it = gModel->firstEdge(); it != gModel->lastEdge(); ++it)
1914 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
1915 if ( !gEdge->haveParametrization())
1917 if ( gEdge->geomType() == GEntity::CompoundCurve )
1921 TopoDS_Edge topoEdge = *((TopoDS_Edge*)gEdge->getNativePtr());
1922 if ( HasSubMesh( topoEdge ) || hideAnyway )
1923 gEdge->setVisibility(0);
1927 for(GModel::fiter it = gModel->firstFace(); it != gModel->lastFace(); ++it)
1931 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
1932 if ( !gFace->haveParametrization())
1934 if ( gFace->geomType() == GEntity::CompoundSurface )
1938 TopoDS_Face topoFace = *((TopoDS_Face*)gFace->getNativePtr());
1939 if ( HasSubMesh( topoFace ) || hideAnyway )
1940 gFace->setVisibility(0);
1944 //================================================================================
1946 * \brief Restore visibility of all geom entities
1948 //================================================================================
1950 void GMSHPlugin_Mesher::RestoreVisibility( GModel* gModel )
1952 for(GModel::eiter it = gModel->firstEdge(); it != gModel->lastEdge(); ++it)
1955 gEdge->setVisibility(1);
1957 for(GModel::fiter it = gModel->firstFace(); it != gModel->lastFace(); ++it)
1960 gFace->setVisibility(1);
1965 void GMSHPlugin_Mesher::toPython( GModel* )
1967 const char* pyFile = "/tmp/gMesh.py";
1968 ofstream outfile( pyFile, ios::out );
1969 if ( !outfile ) return;
1971 outfile << "import salome, SMESH" << std::endl
1972 << "from salome.smesh import smeshBuilder" << std::endl
1973 << "smesh = smeshBuilder.New()" << std::endl
1974 << "mesh = smesh.Mesh()" << std::endl << std::endl;
1976 outfile << "## VERTICES" << endl;
1977 for ( GModel::viter it = _gModel->firstVertex(); it != _gModel->lastVertex(); ++it)
1979 GVertex *gVertex = *it;
1981 for(unsigned int i = 0; i < gVertex->mesh_vertices.size(); i++)
1983 MVertex *v = gVertex->mesh_vertices[i];
1984 if ( v->getIndex() >= 0)
1986 outfile << "n" << v->getNum() << " = mesh.AddNode("
1987 << v->x() << ", " << v->y() << ", " << v->z()<< ")"
1988 << " ## tag = " << gVertex->tag() << endl;
1993 for(GModel::eiter it = _gModel->firstEdge(); it != _gModel->lastEdge(); ++it)
1996 outfile << "## GEdge " << gEdge->tag() << endl;
1998 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
1999 if(gEdge->haveParametrization())
2001 if ( gEdge->geomType() != GEntity::CompoundCurve )
2003 for ( size_t i = 0; i < gEdge->mesh_vertices.size(); i++ )
2005 MVertex *v = gEdge->mesh_vertices[i];
2006 if ( v->getIndex() >= 0 )
2008 outfile << "n" << v->getNum() << " = mesh.AddNode("
2009 << v->x() << ", " << v->y() << ", " << v->z()<< ")"
2010 << " ## tag = " << gEdge->tag() << endl;
2015 for ( GModel::fiter it = _gModel->firstFace(); it != _gModel->lastFace(); ++it)
2018 if ( _compounds.size() && gFace->geomType() == GEntity::DiscreteSurface )
2020 outfile << "## GFace " << gFace->tag() << endl;
2022 for ( size_t i = 0; i < gFace->mesh_vertices.size(); i++ )
2024 MVertex *v = gFace->mesh_vertices[i];
2025 if ( v->getIndex() >= 0 )
2027 outfile << "n" << v->getNum() << " = mesh.AddNode("
2028 << v->x() << ", " << v->y() << ", " << v->z()<< ")"
2029 << " ## tag = " << gFace->tag() << endl;
2034 std::vector<MVertex*> verts(3);
2035 for(GModel::eiter it = _gModel->firstEdge(); it != _gModel->lastEdge(); ++it)
2038 outfile << "## GEdge " << gEdge->tag() << endl;
2040 #if GMSH_MAJOR_VERSION >=4 && GMSH_MINOR_VERSION >=8
2041 if(gEdge->haveParametrization())
2043 if ( gEdge->geomType() != GEntity::CompoundCurve )
2045 for ( size_t i = 0; i < gEdge->getNumMeshElements(); i++ )
2047 MElement *e = gEdge->getMeshElement(i);
2049 e->getVertices(verts);
2051 outfile << "e" << e->getNum() << " = mesh.AddEdge(["
2052 << "n" << verts[0]->getNum() << ","
2053 << "n" << verts[1]->getNum();
2054 if ( verts.size() == 3 )
2055 outfile << "n" << verts[2]->getNum();
2056 outfile << "])"<< endl;
2060 for ( GModel::fiter it = _gModel->firstFace(); it != _gModel->lastFace(); ++it)
2063 if ( _compounds.size() && gFace->geomType() == GEntity::DiscreteSurface )
2065 outfile << "## GFace " << gFace->tag() << endl;
2067 for ( size_t i = 0; i < gFace->getNumMeshElements(); i++ )
2069 MElement *e = gFace->getMeshElement(i);
2071 e->getVertices(verts);
2073 outfile << "f" << e->getNum() << " = mesh.AddFace([";
2074 for ( size_t j = 0; j < verts.size(); j++)
2076 outfile << "n" << verts[j]->getNum();
2077 if ( j < verts.size()-1 )
2080 outfile << "])" << endl;
2083 std::cout << "Write " << pyFile << std::endl;