1 // Copyright (C) 2007-2024 CEA, EDF, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License, or (at your option) any later version.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
22 // File : NETGENPlugin_Remesher_2D.cxx
23 // Created : Thu Sep 21 16:48:46 2017
24 // Author : Edward AGAPOV (eap)
27 #include "NETGENPlugin_Remesher_2D.hxx"
29 #include "NETGENPlugin_Mesher.hxx"
30 #include "NETGENPlugin_Hypothesis_2D.hxx"
32 #include <SMDS_SetIterator.hxx>
33 #include <SMESHDS_Group.hxx>
34 #include <SMESHDS_Mesh.hxx>
35 #include <SMESH_ControlsDef.hxx>
36 #include <SMESH_Gen.hxx>
37 #include <SMESH_MeshAlgos.hxx>
38 #include <SMESH_MesherHelper.hxx>
39 #include <SMESH_Group.hxx>
40 #include <SMESH_MeshEditor.hxx>
41 #include <SMESH_subMesh.hxx>
43 #include <Bnd_B3d.hxx>
44 #include <Precision.hxx>
46 #include <occgeom.hpp>
47 #include <meshing.hpp>
48 #include <stlgeom.hpp>
49 //#include <stltool.hpp>
51 #include <boost/container/flat_set.hpp>
53 using namespace nglib;
57 NETGENPLUGIN_DLL_HEADER
58 extern MeshingParameters mparam;
60 NETGENPLUGIN_DLL_HEADER
61 extern STLParameters stlparam;
63 NETGENPLUGIN_DLL_HEADER
64 extern STLDoctorParams stldoctor;
68 NETGENPLUGIN_DLL_HEADER
70 extern netgen::NgArray<netgen::Point<3> > readedges;
72 extern netgen::Array<netgen::Point<3> > readedges;
78 //=============================================================================
80 * \brief Fill holes in the mesh, since netgen can remesh only a closed shell mesh.
81 * At destruction, remove triangles filling the holes
86 HoleFiller( SMESH_Mesh& meshDS );
88 void AddHoleBordersAndEdges( Ng_STL_Geometry * ngStlGeo, bool toAddEdges );
89 void KeepHole() { myHole.clear(); myCapElems.clear(); }
90 void ClearCapElements() { myCapElems.clear(); }
93 SMESHDS_Mesh* myMeshDS;
94 std::vector< std::vector< gp_XYZ > > myHole; // initial border nodes
95 std::vector< gp_XYZ > myInHolePos; // position inside each hole
96 std::vector< const SMDS_MeshElement* > myCapElems; // elements closing holes
99 //================================================================================
101 * \brief Fill holes in the mesh
103 //================================================================================
105 HoleFiller::HoleFiller( SMESH_Mesh& theMesh ):
106 myMeshDS( theMesh.GetMeshDS() )
108 SMESH_MeshEditor editor( &theMesh );
111 // const double tol = Max( 0.1 * netgen::mparam.minh, Precision::Confusion() );
112 // TIDSortedNodeSet allNodes;
113 // SMESH_MeshEditor::TListOfListOfNodes equalNodes;
114 // editor.FindCoincidentNodes( allNodes, tol, equalNodes, true );
115 // editor.MergeNodes( equalNodes, /*noHoles=*/false );
119 SMESH_MeshAlgos::TFreeBorderVec holes;
120 bool isManifold = true, isGoodOri = true;
121 SMESH_MeshAlgos::FindFreeBorders( *myMeshDS, holes, /*closedOnly=*/true,
122 &isManifold, &isGoodOri );
126 // set bad faces into a compute error
127 const char* text = "Non-manifold mesh. Only manifold mesh can be re-meshed";
128 SMESH_BadInputElements* error =
129 new SMESH_BadInputElements( myMeshDS, COMPERR_BAD_INPUT_MESH, text );
130 SMESH::Controls::MultiConnection2D fun;
131 fun.SetMesh( myMeshDS );
132 SMDS_ElemIteratorPtr fIt = myMeshDS->elementsIterator( SMDSAbs_Face );
133 while ( fIt->more() )
135 const SMDS_MeshElement* f = fIt->next();
136 if ( fun.GetValue( f->GetID() ) > 2 )
137 error->myBadElements.push_back( f );
139 theMesh.GetSubMesh( theMesh.GetShapeToMesh() )->GetComputeError().reset( error );
140 throw SALOME_Exception( text );
143 // don't want to sew coincident borders
144 if ( !holes.empty() )
147 double tol, len, sumLen = 0, minLen = 1e100;
149 for ( size_t i = 0; i < holes.size(); ++i )
151 nbSeg += holes[i].size();
152 SMESH_NodeXYZ p1 = holes[i][0];
153 for ( size_t iP = 1; iP < holes[i].size(); ++iP )
155 SMESH_NodeXYZ p2 = holes[i][iP];
156 len = ( p1 - p2 ).Modulus();
158 minLen = Min( minLen, len );
162 double avgLen = sumLen / double( nbSeg );
163 if ( minLen > 1e-5 * avgLen )
164 tol = 0.1 * minLen; // minLen is not degenerate
168 SMESH_MeshAlgos::CoincidentFreeBorders freeBords;
169 SMESH_MeshAlgos::FindCoincidentFreeBorders( *myMeshDS, tol, freeBords );
170 if ( !freeBords._coincidentGroups.empty() )
172 const char* text = "Can't re-meshed a mesh with coincident free edges";
173 SMESH_BadInputElements* error =
174 new SMESH_BadInputElements( myMeshDS, COMPERR_BAD_INPUT_MESH, text );
175 for ( size_t i = 0; i < freeBords._borders.size(); ++i )
176 error->myBadElements.insert( error->myBadElements.end(),
177 freeBords._borders[i].begin(),
178 freeBords._borders[i].end() );
179 theMesh.GetSubMesh( theMesh.GetShapeToMesh() )->GetComputeError().reset( error );
180 throw SALOME_Exception( text );
185 myHole.resize( holes.size() );
186 myInHolePos.resize( holes.size() );
187 std::vector<const SMDS_MeshElement*> newFaces;
188 for ( size_t i = 0; i < holes.size(); ++i )
191 SMESH_MeshAlgos::FillHole( holes[i], *myMeshDS, newFaces );
193 // keep data to be able to remove hole filling faces after remeshing
194 if ( !newFaces.empty() )
196 myHole[i].resize( holes[i].size() );
197 for ( size_t iP = 0; iP < holes[i].size(); ++iP )
198 myHole[i][iP] = SMESH_NodeXYZ( holes[i][iP] );
200 myInHolePos[i] = ( SMESH_NodeXYZ( newFaces[0]->GetNode(0)) +
201 SMESH_NodeXYZ( newFaces[0]->GetNode(1)) +
202 SMESH_NodeXYZ( newFaces[0]->GetNode(2)) ) / 3.;
203 myCapElems.insert( myCapElems.end(), newFaces.begin(), newFaces.end() );
204 // unmark to be able to remove them if meshing is canceled
205 // for ( size_t iF = 0; iF < newFaces.size(); ++iF )
206 // newFaces[iF]->setIsMarked( false );
212 SMDS_ElemIteratorPtr fIt = myMeshDS->elementsIterator( SMDSAbs_Face );
213 while ( fIt->more() )
215 const SMDS_MeshElement* f = fIt->next();
217 if ( SMESH_MeshAlgos::FaceNormal( f, normal ))
219 TIDSortedElemSet allFaces, refFaces = { f };
220 editor.Reorient2D( allFaces, normal, refFaces, /*allowNonManifold=*/true );
226 //================================================================================
228 * \brief Add hole borders to be kept in a new mesh
230 //================================================================================
232 void HoleFiller::AddHoleBordersAndEdges( Ng_STL_Geometry * ngStlGeo, bool toAddEdges )
234 nglib::readedges.SetSize(0);
236 for ( size_t i = 0; i < myHole.size(); ++i )
237 for ( size_t iP = 1; iP < myHole[i].size(); ++iP )
239 Ng_STL_AddEdge( ngStlGeo,
240 myHole[i][iP-1].ChangeData(),
241 myHole[i][iP-0].ChangeData() );
246 std::vector<const SMDS_MeshNode *> nodes(2);
247 std::vector<const SMDS_MeshElement *> faces(2);
248 SMDS_EdgeIteratorPtr eIt = myMeshDS->edgesIterator();
249 while ( eIt->more() )
251 const SMDS_MeshElement* edge = eIt->next();
252 nodes[0] = edge->GetNode(0);
253 nodes[1] = edge->GetNode(1);
254 // check that an edge is a face border
255 if ( myMeshDS->GetElementsByNodes( nodes, faces, SMDSAbs_Face ))
257 Ng_STL_AddEdge( ngStlGeo,
258 SMESH_NodeXYZ( nodes[0] ).ChangeData(),
259 SMESH_NodeXYZ( nodes[1] ).ChangeData() );
265 //================================================================================
267 * \brief Remove triangles filling the holes
269 //================================================================================
271 HoleFiller::~HoleFiller()
273 if ( myMeshDS->NbNodes() < 3 )
276 if ( !myCapElems.empty() ) // old mesh not removed; simply remove myCapElems
278 for ( size_t i = 0; i < myCapElems.size(); ++i )
279 myMeshDS->RemoveFreeElement( myCapElems[i], /*sm=*/0 );
283 bool hasOrphanNodes = true;
285 const double tol = Max( 1e-3 * netgen::mparam.minh, Precision::Confusion() );
287 for ( size_t i = 0; i < myHole.size(); ++i )
289 std::vector< gp_XYZ >& borderPnt = myHole[i];
290 const gp_XYZ& inHolePos = myInHolePos[i];
291 if ( borderPnt.empty() ) continue;
292 borderPnt.pop_back(); // first point repeated at end
294 // mark all nodes located on the hole border
296 // new nodeSearcher for each hole, otherwise it contains removed nodes for i > 0
297 SMESHUtils::Deleter< SMESH_NodeSearcher > nodeSearcher;
298 if ( hasOrphanNodes )
300 std::vector< const SMDS_MeshNode* > sharedNodes;
301 sharedNodes.reserve( myMeshDS->NbNodes() );
302 SMDS_NodeIteratorPtr nIt = myMeshDS->nodesIterator();
303 while ( nIt->more() )
305 const SMDS_MeshNode* n = nIt->next();
306 if ( n->NbInverseElements() )
307 sharedNodes.push_back( n );
309 hasOrphanNodes = ((int) sharedNodes.size() < myMeshDS->NbNodes() );
310 SMDS_ElemIteratorPtr elemIt( new SMDS_NodeVectorElemIterator( sharedNodes.begin(),
311 sharedNodes.end() ));
312 nodeSearcher._obj = SMESH_MeshAlgos::GetNodeSearcher( elemIt );
316 nodeSearcher._obj = SMESH_MeshAlgos::GetNodeSearcher( *myMeshDS );
319 std::vector< const SMDS_MeshElement* > edgesToRemove;
320 edgesToRemove.reserve( borderPnt.size() );
322 // look for a border point coincident with a node
324 SMESH_NodeXYZ bordNode1;
325 for ( ; iP < borderPnt.size(); ++iP )
327 bordNode1 = nodeSearcher->FindClosestTo( borderPnt[iP] );
328 if (( bordNode1 - borderPnt[iP] ).SquareModulus() < tol * tol )
332 bordNode1._node->setIsMarked( true );
334 // find the rest nodes located on the hole border
335 boost::container::flat_set< const SMDS_MeshNode* > checkedNodes;
336 gp_XYZ p1 = bordNode1;
337 for ( size_t j = 0; j < borderPnt.size()+1; ++j, iP = ( iP+1 ) % borderPnt.size() )
339 // among nodes surrounding bordNode1 find one most close to vec12
340 gp_XYZ vec12 = borderPnt[iP] - p1;
341 bool pntReached = false; // last found node is at iP
342 while ( !pntReached )
344 const SMDS_MeshNode* bordNode = bordNode1._node;
345 SMDS_ElemIteratorPtr fIt = bordNode->GetInverseElementIterator( SMDSAbs_Face );
346 double minArea = 1e100;
347 checkedNodes.clear();
348 checkedNodes.insert( bordNode );
349 while ( fIt->more() )
351 const SMDS_MeshElement* f = fIt->next();
352 for ( int iN = 0, nbN = f->NbNodes(); iN < nbN; ++iN )
354 const SMDS_MeshNode* n = f->GetNode( iN );
355 if ( !checkedNodes.insert( n ).second )
357 SMESH_NodeXYZ pn = n;
358 gp_XYZ vecPN = pn - bordNode1;
359 if ( vecPN * vec12 <= 0 )
361 gp_XYZ vec1N = pn - p1;
362 double a = vec12.CrossSquareMagnitude( vec1N );
369 if ( minArea < std::numeric_limits<double>::min() )
372 if ( bordNode == bordNode1._node )
373 return; // bug in the loop above
375 SMESH_NodeXYZ bordNode2 = bordNode;
376 gp_XYZ vec1N = bordNode2 - p1;
377 double u = ( vec12 * vec1N ) / vec12.SquareModulus(); // param [0,1] of bordNode on vec12
380 bordNode->setIsMarked( true );
381 //cout << bordNode->GetID() << " ";
383 if ( const SMDS_MeshElement* edge = myMeshDS->FindEdge( bordNode1._node, bordNode ))
384 edgesToRemove.push_back( edge );
386 edgesToRemove.push_back( myMeshDS->AddEdge( bordNode1._node, bordNode ));
387 edgesToRemove.back()->setIsMarked( true );
389 if ( minArea > std::numeric_limits<double>::min() &&
390 minArea / vec12.SquareModulus() > tol * tol )
392 // node is far from the border, move it
393 gp_XYZ p = p1 + u * vec12;
394 myMeshDS->MoveNode( bordNode, p.X(), p.Y(), p.Z() );
396 bordNode1 = bordNode2;
398 //else -- there must be another border point between bordNode1 and bordNode
399 pntReached = ( u > 1 - tol );
404 //cout << endl << endl;
406 // remove all faces starting from inHolePos
408 // get a starting face
409 std::vector< const SMDS_MeshNode* > nodesToRemove;
410 std::vector< const SMDS_MeshElement* > facesToRemove;
411 const SMDS_MeshNode* inHoleNode = nodeSearcher->FindClosestTo( inHolePos );
412 if ( inHoleNode && ! inHoleNode->isMarked() )
414 SMDS_ElemIteratorPtr fIt = inHoleNode->GetInverseElementIterator( SMDSAbs_Face );
415 while ( fIt->more() )
416 facesToRemove.push_back( fIt->next() );
420 SMESHUtils::Deleter< SMESH_ElementSearcher > faceSearcher
421 ( SMESH_MeshAlgos::GetElementSearcher( *myMeshDS ));
422 if ( const SMDS_MeshElement* f = faceSearcher->FindClosestTo( inHolePos, SMDSAbs_Face ))
423 facesToRemove.push_back( f );
427 for ( size_t iF = 0; iF < facesToRemove.size(); ++iF )
428 facesToRemove[iF]->setIsMarked( true );
430 // remove faces and nodes
431 TIDSortedElemSet elemSet, avoidSet;
432 const SMDS_MeshElement* e;
433 while ( !facesToRemove.empty() )
435 const SMDS_MeshElement* inHoleFace = facesToRemove.back();
436 facesToRemove.pop_back();
438 // add adjacent faces into facesToRemove
439 for ( int iN = 0, nbN = inHoleFace->NbNodes(); iN < nbN; ++iN )
441 const SMDS_MeshNode* n1 = inHoleFace->GetNode( iN );
442 if ( !n1->isMarked() )
444 SMDS_ElemIteratorPtr eIt = n1->GetInverseElementIterator();
445 while ( eIt->more() )
448 if ( e->GetType() == SMDSAbs_Face )
450 if ( !e->isMarked() )
451 facesToRemove.push_back( e );
452 e->setIsMarked( true );
454 else if ( e->GetType() == SMDSAbs_Edge )
456 myMeshDS->RemoveFreeElement( e, 0, /*fromGroups=*/false );
459 if ( n1->NbInverseElements() == 1 )
460 nodesToRemove.push_back( n1 );
464 const SMDS_MeshNode* n2 = inHoleFace->GetNodeWrap( iN+1 );
465 if (( n2->isMarked() ) &&
466 ( !(e = myMeshDS->FindEdge( n1, n2 )) || !e->isMarked() )) // n1-n2 not hole border
468 if ( e ) // remove edge
469 myMeshDS->RemoveFreeElement( e, 0, /*fromGroups=*/false );
471 avoidSet.insert( inHoleFace );
472 if (( e = SMESH_MeshAlgos::FindFaceInSet( n1, n2, elemSet, avoidSet )))
474 if ( !e->isMarked() )
475 facesToRemove.push_back( e );
476 e->setIsMarked( true );
481 myMeshDS->RemoveFreeElement( inHoleFace, 0, /*fromGroups=*/false );
483 for ( size_t iN = 0; iN < nodesToRemove.size(); ++iN )
484 myMeshDS->RemoveFreeNode( nodesToRemove[iN], 0, /*fromGroups=*/false );
485 nodesToRemove.clear();
488 // remove edges from the hole border
489 // for ( size_t iE = 0; iE < edgesToRemove.size(); ++iE )
490 // myMeshDS->RemoveFreeElement( edgesToRemove[iE], 0, /*fromGroups=*/false );
498 //================================================================================
500 * \brief Fix nodes of a group
502 //================================================================================
504 void fixNodes( SMESHDS_GroupBase* group, netgen::STLGeometry* stlGeom )
506 SMESH_MeshAlgos::MarkElemNodes( group->GetElements(), false ); // un-mark nodes
508 for ( SMDS_ElemIteratorPtr eIt = group->GetElements(); eIt->more(); )
510 const SMDS_MeshElement* e = eIt->next();
511 for ( SMDS_NodeIteratorPtr nIt = e->nodeIterator(); nIt->more(); )
513 const SMDS_MeshNode* n = nIt->next();
516 n->setIsMarked( true );
518 SMESH_NodeXYZ p( n );
519 int id = stlGeom->GetPointNum( netgen::Point<3>( p.X(),p.Y(),p.Z() ));
521 stlGeom->SetLineEndPoint( id );
528 //=============================================================================
532 //=============================================================================
534 NETGENPlugin_Remesher_2D::NETGENPlugin_Remesher_2D(int hypId, SMESH_Gen* gen)
535 : SMESH_2D_Algo(hypId, gen)
537 _name = "NETGEN_Remesher_2D";
538 _shapeType = (1 << TopAbs_FACE); // 1 bit /shape type
539 _compatibleHypothesis.push_back("NETGEN_RemesherParameters_2D");
540 _requireShape = false;
545 //=============================================================================
547 * Check assigned hypotheses
549 //=============================================================================
551 bool NETGENPlugin_Remesher_2D::CheckHypothesis (SMESH_Mesh& theMesh,
552 const TopoDS_Shape& theShape,
553 Hypothesis_Status& theStatus)
557 // can work with no hypothesis
558 theStatus = SMESH_Hypothesis::HYP_OK;
560 const list<const SMESHDS_Hypothesis*>& hyps =
561 GetUsedHypothesis( theMesh, theShape, /*skipAux=*/true );
563 switch ( hyps.size() ) {
567 _hypothesis = hyps.front();
570 theStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
573 return theStatus == SMESH_Hypothesis::HYP_OK;
576 //=============================================================================
578 * Compute mesh on an input mesh
580 //=============================================================================
582 bool NETGENPlugin_Remesher_2D::Compute(SMESH_Mesh& theMesh,
583 SMESH_MesherHelper* theHelper)
585 if ( theMesh.NbFaces() == 0 )
586 return !error( COMPERR_WARNING, "No faces in input mesh");
588 NETGENPlugin_Mesher mesher( &theMesh, theMesh.GetShapeToMesh(), /*isVol=*/false);
589 NETGENPlugin_NetgenLibWrapper ngLib;
590 netgen::Mesh * ngMesh = (netgen::Mesh*) ngLib._ngMesh;
591 Ng_STL_Geometry * ngStlGeo = Ng_STL_NewGeometry();
592 netgen::STLTopology* stlTopo = (netgen::STLTopology*) ngStlGeo;
593 netgen::multithread.terminate = 0;
595 const NETGENPlugin_RemesherHypothesis_2D* hyp =
596 dynamic_cast<const NETGENPlugin_RemesherHypothesis_2D*>( _hypothesis );
597 mesher.SetParameters( hyp );// for holeFiller
599 SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
600 HoleFiller holeFiller( theMesh );
601 //theHelper->SetIsQuadratic( theMesh.NbFaces( ORDER_QUADRATIC ));
603 // fill ngStlGeo with triangles
604 SMDS_ElemIteratorPtr fIt = meshDS->elementsIterator( SMDSAbs_Face );
605 while ( fIt->more() )
607 const SMDS_MeshElement* f = fIt->next();
608 SMESH_NodeXYZ n1 = f->GetNode( 0 );
609 SMESH_NodeXYZ n2 = f->GetNode( 1 );
610 SMESH_NodeXYZ n3 = f->GetNode( 2 );
611 Ng_STL_AddTriangle( ngStlGeo,
615 if ( f->NbNodes() > 3 )
617 n2.Set( f->GetNode( 3 ));
618 Ng_STL_AddTriangle( ngStlGeo,
625 bool toAddExistingEdges = ( hyp && hyp->GetKeepExistingEdges() );
626 holeFiller.AddHoleBordersAndEdges( ngStlGeo, toAddExistingEdges );
629 //netgen::stldoctor.geom_tol_fact = 1e-12; // pointtol=boundingbox.Diam()*stldoctor.geom_tol_fact
630 Ng_Result ng_res = Ng_STL_InitSTLGeometry( ngStlGeo );
631 if ( ng_res != NG_OK )
634 holeFiller.KeepHole();
636 std::string txt = "Error Initialising the STL Geometry";
637 if ( !stlTopo->GetStatusText().empty() )
638 txt += ". " + stlTopo->GetStatusText();
639 return error( COMPERR_BAD_INPUT_MESH, txt );
643 Ng_Meshing_Parameters ngParams;
646 ngParams.maxh = hyp->GetMaxSize();
647 ngParams.minh = hyp->GetMinSize();
648 ngParams.meshsize_filename = (char*) hyp->GetMeshSizeFile().c_str();
649 ngParams.quad_dominated = hyp->GetQuadAllowed();
651 netgen::stlparam.yangle = hyp->GetRidgeAngle();
652 netgen::stlparam.edgecornerangle = hyp->GetEdgeCornerAngle();
653 netgen::stlparam.chartangle = hyp->GetChartAngle();
654 netgen::stlparam.outerchartangle = hyp->GetOuterChartAngle();
655 netgen::stlparam.resthchartdistfac = hyp->GetRestHChartDistFactor();
656 netgen::stlparam.resthchartdistenable = hyp->GetRestHChartDistEnable();
657 netgen::stlparam.resthlinelengthfac = hyp->GetRestHLineLengthFactor();
658 netgen::stlparam.resthlinelengthenable = hyp->GetRestHLineLengthEnable();
660 netgen::stlparam.resthcloseedgefac = hyp->GetRestHCloseEdgeFactor();
661 netgen::stlparam.resthcloseedgeenable = hyp->GetRestHCloseEdgeEnable();
663 netgen::stlparam.resthsurfcurvfac = hyp->GetRestHSurfCurvFactor();
664 netgen::stlparam.resthsurfcurvenable = hyp->GetRestHSurfCurvEnable();
665 netgen::stlparam.resthedgeanglefac = hyp->GetRestHEdgeAngleFactor();
666 netgen::stlparam.resthedgeangleenable = hyp->GetRestHEdgeAngleEnable();
667 netgen::stlparam.resthsurfmeshcurvfac = hyp->GetRestHSurfMeshCurvFactor();
668 netgen::stlparam.resthsurfmeshcurvenable = hyp->GetRestHSurfMeshCurvEnable();
670 mesher.SetParameters( hyp );
674 double diagSize = Dist( stlTopo->GetBoundingBox().PMin(), stlTopo->GetBoundingBox().PMax());
675 netgen::mparam.maxh = diagSize / GetGen()->GetBoundaryBoxSegmentation();
676 netgen::mparam.minh = netgen::mparam.maxh;
679 // save netgen::mparam as Ng_STL_MakeEdges() modify it by Ng_Meshing_Parameters
680 netgen::MeshingParameters savedParams = netgen::mparam;
682 if ( SMESH_Group* fixedEdges = ( hyp ? hyp->GetFixedEdgeGroup( theMesh ) : 0 ))
684 netgen::STLGeometry* stlGeom = (netgen::STLGeometry*)ngStlGeo;
686 // the following code is taken from STLMeshing() method
689 stlGeom->BuildEdges( netgen::stlparam );
690 stlGeom->MakeAtlas( *ngMesh, netgen::mparam, netgen::stlparam );
691 stlGeom->CalcFaceNums();
692 stlGeom->AddFaceEdges();
693 fixNodes( fixedEdges->GetGroupDS(), stlGeom );
694 stlGeom->LinkEdges( netgen::stlparam );
697 stlGeom->BuildEdges();
698 stlGeom->MakeAtlas( *ngMesh );
699 stlGeom->CalcFaceNums();
700 stlGeom->AddFaceEdges();
701 fixNodes( fixedEdges->GetGroupDS(), stlGeom );
702 stlGeom->LinkEdges();
704 ngMesh->ClearFaceDescriptors();
705 for (int i = 1; i <= stlGeom->GetNOFaces(); i++)
706 ngMesh->AddFaceDescriptor (netgen::FaceDescriptor (i, 1, 0, 0));
708 stlGeom->edgesfound = 1;
712 Ng_STL_MakeEdges( ngStlGeo, ngLib.ngMesh(), &ngParams );
715 netgen::mparam = savedParams;
717 double h = netgen::mparam.maxh;
718 ngMesh->SetGlobalH( h );
719 ngMesh->SetMinimalH( netgen::mparam.minh );
720 ngMesh->SetLocalH( stlTopo->GetBoundingBox().PMin() - netgen::Vec3d(h, h, h),
721 stlTopo->GetBoundingBox().PMax() + netgen::Vec3d(h, h, h),
722 netgen::mparam.grading );
723 ngMesh->LoadLocalMeshSize( ngParams.meshsize_filename );
725 netgen::OCCGeometry occgeo;
726 mesher.SetLocalSize( occgeo, *ngMesh );
731 ng_res = Ng_STL_GenerateSurfaceMesh( ngStlGeo, ngLib.ngMesh(), &ngParams );
733 catch (netgen::NgException & ex)
735 if ( netgen::multithread.terminate )
736 if ( !hyp || !hyp->GetLoadMeshOnCancel() )
739 if ( ng_res != NG_OK )
740 return error( "Error in Surface Meshing" );
742 int nbN = ngMesh->GetNP();
743 int nbE = ngMesh->GetNSeg();
744 int nbF = ngMesh->GetNSE();
746 return error( "Error in Surface Meshing" );
748 // remove existing mesh
749 holeFiller.ClearCapElements();
750 SMDS_ElemIteratorPtr eIt = meshDS->elementsIterator();
751 while ( eIt->more() )
752 meshDS->RemoveFreeElement( eIt->next(), /*sm=*/0 );
753 SMDS_NodeIteratorPtr nIt = meshDS->nodesIterator();
754 while ( nIt->more() )
755 meshDS->RemoveFreeNode( nIt->next(), /*sm=*/0 );
760 std::vector< const SMDS_MeshNode* > newNodes( nbN+1 );
761 for ( int i = 1; i <= nbN; ++i )
763 const netgen::MeshPoint& p = ngMesh->Point(i);
764 newNodes[i] = meshDS->AddNode( p(0),p(1),p(2) );
768 std::vector<const SMDS_MeshNode*> nodes(4);
769 for ( int i = 1; i <= nbE; ++i )
771 const netgen::Segment& seg = ngMesh->LineSegment(i);
773 for ( int j = 0; j < 2; ++j )
775 size_t pind = seg.pnums[j];
776 if ( pind > 0 && pind < newNodes.size() )
777 nodes.push_back( newNodes[ pind ]);
781 if ( nodes.size() == 2 && !meshDS->FindEdge( nodes[0], nodes[1] ))
782 meshDS->AddEdge( nodes[0], nodes[1] );
785 // find existing groups
786 const char* theNamePrefix = "Surface_";
787 const size_t theNamePrefixLen = strlen( theNamePrefix );
788 std::vector< SMESHDS_Group* > groups;
789 if ( hyp && hyp->GetMakeGroupsOfSurfaces() )
791 SMESH_Mesh::GroupIteratorPtr grIt = theMesh.GetGroups();
792 while ( grIt->more() )
794 SMESH_Group* group = grIt->next();
795 SMESHDS_Group* groupDS;
796 if (( group->GetGroupDS()->GetType() == SMDSAbs_Face ) &&
797 ( strncmp( group->GetName(), theNamePrefix, theNamePrefixLen ) == 0 ) &&
798 ( groupDS = dynamic_cast<SMESHDS_Group*>( group->GetGroupDS() )))
799 groups.push_back( groupDS );
804 for ( int i = 1; i <= nbF; ++i )
806 const netgen::Element2d& elem = ngMesh->SurfaceElement(i);
808 for ( int j = 1; j <= elem.GetNP(); ++j )
810 size_t pind = elem.PNum(j);
811 if ( pind > 0 && pind < newNodes.size() )
812 nodes.push_back( newNodes[ pind ]);
816 const SMDS_MeshElement* newFace = 0;
817 switch( nodes.size() )
819 case 3: newFace = meshDS->AddFace( nodes[0], nodes[1], nodes[2] ); break;
820 case 4: newFace = meshDS->AddFace( nodes[0], nodes[1], nodes[2], nodes[3] ); break;
823 // add newFace to a group
824 if ( newFace && hyp && hyp->GetMakeGroupsOfSurfaces() )
826 if ((size_t) elem.GetIndex()-1 >= groups.size() )
827 groups.resize( elem.GetIndex(), 0 );
829 SMESHDS_Group* & group = groups[ elem.GetIndex()-1 ];
832 SMESH_Group* gr = theMesh.AddGroup( SMDSAbs_Face, "");
833 group = static_cast<SMESHDS_Group*>( gr->GetGroupDS() );
835 group->SMDSGroup().Add( newFace );
841 for ( size_t i = 0; i < groups.size(); ++i )
845 if ( groups[i]->IsEmpty() )
847 theMesh.RemoveGroup( groups[i]->GetID() );
849 else if ( SMESH_Group* g = theMesh.GetGroup( groups[i]->GetID() ))
851 g->SetName( SMESH_Comment( theNamePrefix ) << groupIndex++ );
855 // as we don't assign the new triangles to a shape (the pseudo-shape),
856 // to avoid their removal at hypothesis modification,
857 // we mark the shape as always computed to avoid the error messages
858 // that no elements assigned to the shape
859 theMesh.GetSubMesh( theHelper->GetSubShape() )->SetIsAlwaysComputed( true );
864 //=============================================================================
866 * Do not compute mesh on geometry
868 //=============================================================================
870 bool NETGENPlugin_Remesher_2D::Compute(SMESH_Mesh& /*theMesh*/,
871 const TopoDS_Shape& /*theShape*/)
876 //=============================================================================
878 * Terminate Compute()
880 //=============================================================================
882 void NETGENPlugin_Remesher_2D::CancelCompute()
884 SMESH_Algo::CancelCompute();
885 netgen::multithread.terminate = 1;
888 //================================================================================
890 * \brief Return progress of Compute() [0.,1]
892 //================================================================================
894 double NETGENPlugin_Remesher_2D::GetProgress() const
896 return netgen::multithread.percent / 100.;
899 //=============================================================================
903 //=============================================================================
905 bool NETGENPlugin_Remesher_2D::Evaluate(SMESH_Mesh& /*aMesh*/,
906 const TopoDS_Shape& /*aShape*/,
907 MapShapeNbElems& /*aResMap*/)