1 // Copyright (C) 2007-2010 CEA/DEN, EDF R&D, 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.
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
23 // NETGENPlugin : C++ implementation
24 // File : NETGENPlugin_Mesher.cxx
25 // Author : Michael Sazonov (OCN)
28 //=============================================================================
30 #include "NETGENPlugin_Mesher.hxx"
31 #include "NETGENPlugin_Hypothesis_2D.hxx"
32 #include "NETGENPlugin_SimpleHypothesis_3D.hxx"
34 #include <SMDS_FaceOfNodes.hxx>
35 #include <SMDS_MeshElement.hxx>
36 #include <SMDS_MeshNode.hxx>
37 #include <SMESHDS_Mesh.hxx>
38 #include <SMESH_Comment.hxx>
39 #include <SMESH_ComputeError.hxx>
40 #include <SMESH_File.hxx>
41 #include <SMESH_Mesh.hxx>
42 #include <SMESH_MesherHelper.hxx>
43 #include <SMESH_subMesh.hxx>
44 #include <SMESH_Gen_i.hxx>
45 #include <utilities.h>
50 #include <BRep_Tool.hxx>
51 #include <NCollection_Map.hxx>
52 #include <OSD_File.hxx>
53 #include <OSD_Path.hxx>
54 #include <Standard_ErrorHandler.hxx>
55 #include <Standard_ProgramError.hxx>
56 #include <TCollection_AsciiString.hxx>
58 #include <TopExp_Explorer.hxx>
59 #include <TopTools_DataMapIteratorOfDataMapOfShapeInteger.hxx>
60 #include <TopTools_DataMapIteratorOfDataMapOfShapeShape.hxx>
61 #include <TopTools_DataMapOfShapeInteger.hxx>
62 #include <TopTools_DataMapOfShapeShape.hxx>
63 #include <TopTools_ListIteratorOfListOfShape.hxx>
64 #include <TopTools_MapOfShape.hxx>
66 #include <GeomAdaptor_Curve.hxx>
67 #include <GCPnts_AbscissaPoint.hxx>
69 // Netgen include files
73 #include <occgeom.hpp>
74 #include <meshing.hpp>
75 //#include <ngexception.hpp>
77 extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, int, int, char*);
78 extern MeshingParameters mparam;
81 using namespace nglib;
85 #define nodeVec_ACCESS(index) (SMDS_MeshNode*) nodeVec.at((index))
87 #define nodeVec_ACCESS(index) (SMDS_MeshNode*) nodeVec[index]
91 #define NGPOINT_COORDS(p) p(0),p(1),p(2)
93 #define NGPOINT_COORDS(p) p.X(),p.Y(),p.Z()
96 // dump elements added to ng mesh
97 //#define DUMP_SEGMENTS
98 //#define DUMP_TRIANGLES
99 //#define DUMP_TRIANGLES_SCRIPT "/tmp/trias.py" //!< debug addIntVerticesInSolids()
101 TopTools_IndexedMapOfShape ShapesWithLocalSize;
102 std::map<int,double> VertexId2LocalSize;
103 std::map<int,double> EdgeId2LocalSize;
104 std::map<int,double> FaceId2LocalSize;
106 //=============================================================================
110 //=============================================================================
112 NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh* mesh,
113 const TopoDS_Shape& aShape,
122 ShapesWithLocalSize.Clear();
123 VertexId2LocalSize.clear();
124 EdgeId2LocalSize.clear();
125 FaceId2LocalSize.clear();
128 //================================================================================
130 * \brief Initialize global NETGEN parameters with default values
132 //================================================================================
134 void NETGENPlugin_Mesher::defaultParameters()
136 netgen::MeshingParameters& mparams = netgen::mparam;
137 // maximal mesh edge size
138 mparams.maxh = NETGENPlugin_Hypothesis::GetDefaultMaxSize();
139 // minimal number of segments per edge
140 mparams.segmentsperedge = NETGENPlugin_Hypothesis::GetDefaultNbSegPerEdge();
141 // rate of growth of size between elements
142 mparams.grading = NETGENPlugin_Hypothesis::GetDefaultGrowthRate();
143 // safety factor for curvatures (elements per radius)
144 mparams.curvaturesafety = NETGENPlugin_Hypothesis::GetDefaultNbSegPerRadius();
145 // create elements of second order
146 mparams.secondorder = NETGENPlugin_Hypothesis::GetDefaultSecondOrder() ? 1 : 0;
147 // quad-dominated surface meshing
151 mparams.quad = NETGENPlugin_Hypothesis_2D::GetDefaultQuadAllowed() ? 1 : 0;
154 //=============================================================================
158 //=============================================================================
159 void SetLocalSize(TopoDS_Shape GeomShape, double LocalSize)
161 TopAbs_ShapeEnum GeomType = GeomShape.ShapeType();
162 if (GeomType == TopAbs_COMPOUND) {
163 for (TopoDS_Iterator it (GeomShape); it.More(); it.Next()) {
164 SetLocalSize(it.Value(), LocalSize);
169 if (! ShapesWithLocalSize.Contains(GeomShape))
170 key = ShapesWithLocalSize.Add(GeomShape);
172 key = ShapesWithLocalSize.FindIndex(GeomShape);
173 if (GeomType == TopAbs_VERTEX) {
174 VertexId2LocalSize[key] = LocalSize;
175 } else if (GeomType == TopAbs_EDGE) {
176 EdgeId2LocalSize[key] = LocalSize;
177 } else if (GeomType == TopAbs_FACE) {
178 FaceId2LocalSize[key] = LocalSize;
182 //=============================================================================
184 * Pass parameters to NETGEN
186 //=============================================================================
187 void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp)
191 netgen::MeshingParameters& mparams = netgen::mparam;
192 // Initialize global NETGEN parameters:
193 // maximal mesh segment size
194 mparams.maxh = hyp->GetMaxSize();
195 // minimal number of segments per edge
196 mparams.segmentsperedge = hyp->GetNbSegPerEdge();
197 // rate of growth of size between elements
198 mparams.grading = hyp->GetGrowthRate();
199 // safety factor for curvatures (elements per radius)
200 mparams.curvaturesafety = hyp->GetNbSegPerRadius();
201 // create elements of second order
202 mparams.secondorder = hyp->GetSecondOrder() ? 1 : 0;
203 // quad-dominated surface meshing
204 // only triangles are allowed for volumic mesh
206 mparams.quad = static_cast<const NETGENPlugin_Hypothesis_2D*>
207 (hyp)->GetQuadAllowed() ? 1 : 0;
208 _optimize = hyp->GetOptimize();
211 SMESH_Gen_i* smeshGen_i = SMESH_Gen_i::GetSMESHGen();
212 CORBA::Object_var anObject = smeshGen_i->GetNS()->Resolve("/myStudyManager");
213 SALOMEDS::StudyManager_var aStudyMgr = SALOMEDS::StudyManager::_narrow(anObject);
214 SALOMEDS::Study_var myStudy = aStudyMgr->GetStudyByID(hyp->GetStudyId());
216 const NETGENPlugin_Hypothesis::TLocalSize localSizes = hyp->GetLocalSizesAndEntries();
217 NETGENPlugin_Hypothesis::TLocalSize::const_iterator it = localSizes.begin();
218 for (it ; it != localSizes.end() ; it++)
220 std::string entry = (*it).first;
221 double val = (*it).second;
223 GEOM::GEOM_Object_var aGeomObj;
224 TopoDS_Shape S = TopoDS_Shape();
225 SALOMEDS::SObject_var aSObj = myStudy->FindObjectID( entry.c_str() );
226 SALOMEDS::GenericAttribute_var anAttr;
227 if (!aSObj->_is_nil() && aSObj->FindAttribute(anAttr, "AttributeIOR")) {
228 SALOMEDS::AttributeIOR_var anIOR = SALOMEDS::AttributeIOR::_narrow(anAttr);
229 CORBA::String_var aVal = anIOR->Value();
230 CORBA::Object_var obj = myStudy->ConvertIORToObject(aVal);
231 aGeomObj = GEOM::GEOM_Object::_narrow(obj);
233 if ( !aGeomObj->_is_nil() )
234 S = smeshGen_i->GeomObjectToShape( aGeomObj.in() );
236 SetLocalSize(S, val);
241 //=============================================================================
243 * Pass simple parameters to NETGEN
245 //=============================================================================
247 void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_SimpleHypothesis_2D* hyp)
254 //=============================================================================
256 * Link - a pair of integer numbers
258 //=============================================================================
262 Link(int _n1, int _n2) : n1(_n1), n2(_n2) {}
263 Link() : n1(0), n2(0) {}
266 int HashCode(const Link& aLink, int aLimit)
268 return HashCode(aLink.n1 + aLink.n2, aLimit);
271 Standard_Boolean IsEqual(const Link& aLink1, const Link& aLink2)
273 return (aLink1.n1 == aLink2.n1 && aLink1.n2 == aLink2.n2 ||
274 aLink1.n1 == aLink2.n2 && aLink1.n2 == aLink2.n1);
277 //================================================================================
279 * \brief Initialize netgen::OCCGeometry with OCCT shape
281 //================================================================================
283 void NETGENPlugin_Mesher::PrepareOCCgeometry(netgen::OCCGeometry& occgeo,
284 const TopoDS_Shape& shape,
286 list< SMESH_subMesh* > * meshedSM,
287 NETGENPlugin_Internals* intern)
289 BRepTools::Clean (shape);
291 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
294 BRepMesh_IncrementalMesh::BRepMesh_IncrementalMesh (shape, 0.01, true);
295 } catch (Standard_Failure) {
298 BRepBndLib::Add (shape, bb);
299 double x1,y1,z1,x2,y2,z2;
300 bb.Get (x1,y1,z1,x2,y2,z2);
301 MESSAGE("shape bounding box:\n" <<
302 "(" << x1 << " " << y1 << " " << z1 << ") " <<
303 "(" << x2 << " " << y2 << " " << z2 << ")");
304 netgen::Point<3> p1 = netgen::Point<3> (x1,y1,z1);
305 netgen::Point<3> p2 = netgen::Point<3> (x2,y2,z2);
306 occgeo.boundingbox = netgen::Box<3> (p1,p2);
308 occgeo.shape = shape;
311 // fill maps of shapes of occgeo with not yet meshed subshapes
313 // get root submeshes
314 list< SMESH_subMesh* > rootSM;
315 if ( SMESH_subMesh* sm = mesh.GetSubMeshContaining( shape )) {
316 rootSM.push_back( sm );
319 for ( TopoDS_Iterator it( shape ); it.More(); it.Next() )
320 rootSM.push_back( mesh.GetSubMesh( it.Value() ));
323 // add subshapes of empty submeshes
324 list< SMESH_subMesh* >::iterator rootIt = rootSM.begin(), rootEnd = rootSM.end();
325 for ( ; rootIt != rootEnd; ++rootIt ) {
326 SMESH_subMesh * root = *rootIt;
327 SMESH_subMeshIteratorPtr smIt = root->getDependsOnIterator(/*includeSelf=*/true,
328 /*complexShapeFirst=*/true);
329 // to find a right orientation of subshapes (PAL20462)
330 TopTools_IndexedMapOfShape subShapes;
331 TopExp::MapShapes(root->GetSubShape(), subShapes);
332 while ( smIt->more() )
334 SMESH_subMesh* sm = smIt->next();
335 TopoDS_Shape shape = sm->GetSubShape();
336 if ( intern && intern->isShapeToPrecompute( shape ))
338 if ( !meshedSM || sm->IsEmpty() )
340 if ( shape.ShapeType() != TopAbs_VERTEX )
341 shape = subShapes( subShapes.FindIndex( shape ));// shape -> index -> oriented shape
342 if ( shape.Orientation() >= TopAbs_INTERNAL )
343 shape.Orientation( TopAbs_FORWARD ); // isuue 0020676
344 switch ( shape.ShapeType() ) {
345 case TopAbs_FACE : occgeo.fmap.Add( shape ); break;
346 case TopAbs_EDGE : occgeo.emap.Add( shape ); break;
347 case TopAbs_VERTEX: occgeo.vmap.Add( shape ); break;
348 case TopAbs_SOLID :occgeo.somap.Add( shape ); break;
352 // collect submeshes of meshed shapes
355 meshedSM->push_back( sm );
359 occgeo.facemeshstatus.SetSize (occgeo.fmap.Extent());
360 occgeo.facemeshstatus = 0;
362 occgeo.face_maxh.SetSize(occgeo.fmap.Extent());
363 occgeo.face_maxh = netgen::mparam.maxh;
364 occgeo.face_maxh_modified.SetSize(occgeo.fmap.Extent());
365 occgeo.face_maxh_modified = 0;
372 //================================================================================
374 * \brief return id of netgen point corresponding to SMDS node
376 //================================================================================
377 typedef map< const SMDS_MeshNode*, int > TNode2IdMap;
379 int ngNodeId( const SMDS_MeshNode* node,
380 netgen::Mesh& ngMesh,
381 TNode2IdMap& nodeNgIdMap)
383 int newNgId = ngMesh.GetNP() + 1;
385 TNode2IdMap::iterator node_id = nodeNgIdMap.insert( make_pair( node, newNgId )).first;
387 if ( node_id->second == newNgId)
389 #if defined(DUMP_SEGMENTS) || defined(DUMP_TRIANGLES)
390 cout << "Ng " << newNgId << " - " << node;
392 netgen::MeshPoint p( netgen::Point<3> (node->X(), node->Y(), node->Z()) );
393 ngMesh.AddPoint( p );
395 return node_id->second;
399 //================================================================================
401 * \brief fill ngMesh with nodes and elements of computed submeshes
403 //================================================================================
405 bool NETGENPlugin_Mesher::fillNgMesh(const netgen::OCCGeometry& occgeom,
406 netgen::Mesh& ngMesh,
407 vector<const SMDS_MeshNode*>& nodeVec,
408 const list< SMESH_subMesh* > & meshedSM)
410 TNode2IdMap nodeNgIdMap;
411 if ( !nodeVec.empty() )
412 for ( int i = 1; i < nodeVec.size(); ++i )
413 nodeNgIdMap.insert( make_pair( nodeVec[i], i ));
415 TopTools_MapOfShape visitedShapes;
417 SMESH_MesherHelper helper (*_mesh);
419 int faceID = occgeom.fmap.Extent();
421 list< SMESH_subMesh* >::const_iterator smIt, smEnd = meshedSM.end();
422 for ( smIt = meshedSM.begin(); smIt != smEnd; ++smIt )
424 SMESH_subMesh* sm = *smIt;
425 if ( !visitedShapes.Add( sm->GetSubShape() ))
428 SMESHDS_SubMesh * smDS = sm->GetSubMeshDS();
429 if ( !smDS ) continue;
431 switch ( sm->GetSubShape().ShapeType() )
433 case TopAbs_EDGE: { // EDGE
434 // ----------------------
435 TopoDS_Edge geomEdge = TopoDS::Edge( sm->GetSubShape() );
436 if ( geomEdge.Orientation() >= TopAbs_INTERNAL )
437 geomEdge.Orientation( TopAbs_FORWARD ); // isuue 0020676
439 // Add ng segments for each not meshed face the edge bounds
440 TopTools_MapOfShape visitedAncestors;
441 PShapeIteratorPtr fIt = helper.GetAncestors( geomEdge, *sm->GetFather(), TopAbs_FACE );
442 while ( const TopoDS_Shape * anc = fIt->next() )
444 if ( !visitedAncestors.Add( *anc )) continue;
445 TopoDS_Face face = TopoDS::Face( *anc );
446 if ( face.Orientation() >= TopAbs_INTERNAL )
447 face.Orientation( TopAbs_FORWARD ); // isuue 0020676
449 int faceID = occgeom.fmap.FindIndex( face );
451 continue; // meshed face
453 // find out orientation of geomEdge within face
454 TopAbs_Orientation fOri = helper.GetSubShapeOri( face, geomEdge );
456 // get all nodes from geomEdge
457 bool isForwad = ( fOri == geomEdge.Orientation() );
458 bool isQuad = smDS->NbElements() ? smDS->GetElements()->next()->IsQuadratic() : false;
459 StdMeshers_FaceSide fSide( face, geomEdge, _mesh, isForwad, isQuad );
460 const vector<UVPtStruct>& points = fSide.GetUVPtStruct();
461 int i, nbSeg = fSide.NbSegments();
463 double otherSeamParam = 0;
464 helper.SetSubShape( face );
465 bool isSeam = helper.IsRealSeam( geomEdge );
468 helper.GetOtherParam( helper.GetPeriodicIndex() == 1 ? points[0].u : points[0].v );
472 int prevNgId = ngNodeId( points[0].node, ngMesh, nodeNgIdMap );
474 for ( i = 0; i < nbSeg; ++i )
476 const UVPtStruct& p1 = points[ i ];
477 const UVPtStruct& p2 = points[ i+1 ];
482 seg[1] = prevNgId = ngNodeId( p2.node, ngMesh, nodeNgIdMap );
483 // node param on curve
484 seg.epgeominfo[ 0 ].dist = p1.param;
485 seg.epgeominfo[ 1 ].dist = p2.param;
487 seg.epgeominfo[ 0 ].u = p1.u;
488 seg.epgeominfo[ 0 ].v = p1.v;
489 seg.epgeominfo[ 1 ].u = p2.u;
490 seg.epgeominfo[ 1 ].v = p2.v;
492 //seg.epgeominfo[ iEnd ].edgenr = edgeID; // = geom.emap.FindIndex(edge);
493 seg.si = faceID; // = geom.fmap.FindIndex (face);
494 seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
495 ngMesh.AddSegment (seg);
497 cout << "Segment: " << seg.edgenr << " on SMESH face " << helper.GetMeshDS()->ShapeToIndex( face ) << endl
498 << "\tface index: " << seg.si << endl
499 << "\tp1: " << seg.p1 << endl
500 << "\tp2: " << seg.p2 << endl
501 << "\tp0 param: " << seg.epgeominfo[ 0 ].dist << endl
502 << "\tp0 uv: " << seg.epgeominfo[ 0 ].u <<", "<< seg.epgeominfo[ 0 ].v << endl
503 << "\tp0 edge: " << seg.epgeominfo[ 0 ].edgenr << endl
504 << "\tp1 param: " << seg.epgeominfo[ 1 ].dist << endl
505 << "\tp1 uv: " << seg.epgeominfo[ 1 ].u <<", "<< seg.epgeominfo[ 1 ].v << endl
506 << "\tp1 edge: " << seg.epgeominfo[ 1 ].edgenr << endl;
510 if ( helper.GetPeriodicIndex() == 1 ) {
511 seg.epgeominfo[ 0 ].u = otherSeamParam;
512 seg.epgeominfo[ 1 ].u = otherSeamParam;
513 swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v);
515 seg.epgeominfo[ 0 ].v = otherSeamParam;
516 seg.epgeominfo[ 1 ].v = otherSeamParam;
517 swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u);
519 swap (seg[0], seg[1]);
520 swap (seg.epgeominfo[0].dist, seg.epgeominfo[1].dist);
521 seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
522 ngMesh.AddSegment (seg);
524 else if ( fOri == TopAbs_INTERNAL )
526 swap (seg[0], seg[1]);
527 swap( seg.epgeominfo[0], seg.epgeominfo[1] );
528 seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
529 ngMesh.AddSegment (seg);
531 cout << "Segment: " << seg.edgenr << endl << "\t is REVERSE of the previous" << endl;
535 } // loop on geomEdge ancestors
538 } // case TopAbs_EDGE
540 case TopAbs_FACE: { // FACE
541 // ----------------------
542 const TopoDS_Face& geomFace = TopoDS::Face( sm->GetSubShape() );
543 helper.SetSubShape( geomFace );
544 bool isInternalFace = ( geomFace.Orientation() == TopAbs_INTERNAL );
546 // Find solids the geomFace bounds
547 int solidID1 = 0, solidID2 = 0;
548 PShapeIteratorPtr solidIt = helper.GetAncestors( geomFace, *sm->GetFather(), TopAbs_SOLID);
549 while ( const TopoDS_Shape * solid = solidIt->next() )
551 int id = occgeom.somap.FindIndex ( *solid );
552 if ( solidID1 && id != solidID1 ) solidID2 = id;
556 _faceDescriptors[ faceID ].first = solidID1;
557 _faceDescriptors[ faceID ].second = solidID2;
559 // Orient the face correctly in solidID1 (issue 0020206)
560 bool reverse = false;
562 TopoDS_Shape solid = occgeom.somap( solidID1 );
563 TopAbs_Orientation faceOriInSolid = helper.GetSubShapeOri( solid, geomFace );
564 if ( faceOriInSolid >= 0 )
565 reverse = SMESH_Algo::IsReversedSubMesh
566 ( TopoDS::Face( geomFace.Oriented( faceOriInSolid )), helper.GetMeshDS() );
569 // Add surface elements
571 netgen::Element2d tri(3);
572 tri.SetIndex ( faceID );
575 #ifdef DUMP_TRIANGLES
576 cout << "SMESH face " << helper.GetMeshDS()->ShapeToIndex( geomFace )
577 << " internal="<<isInternalFace<< " border="<<isBorderFace << endl;
579 SMDS_ElemIteratorPtr faces = smDS->GetElements();
580 while ( faces->more() )
582 const SMDS_MeshElement* f = faces->next();
583 if ( f->NbNodes() % 3 != 0 ) // not triangle
585 PShapeIteratorPtr solidIt=helper.GetAncestors(geomFace,*sm->GetFather(),TopAbs_SOLID);
586 if ( const TopoDS_Shape * solid = solidIt->next() )
587 sm = _mesh->GetSubMesh( *solid );
588 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
589 smError.reset( new SMESH_ComputeError(COMPERR_BAD_INPUT_MESH,"Not triangle submesh"));
590 smError->myBadElements.push_back( f );
594 for ( int i = 0; i < 3; ++i )
596 const SMDS_MeshNode* node = f->GetNode( i ), * inFaceNode=0;
598 // get node UV on face
599 int shapeID = node->GetPosition()->GetShapeId();
600 if ( helper.IsSeamShape( shapeID ))
601 if ( helper.IsSeamShape( f->GetNodeWrap( i+1 )->GetPosition()->GetShapeId() ))
602 inFaceNode = f->GetNodeWrap( i-1 );
604 inFaceNode = f->GetNodeWrap( i+1 );
605 gp_XY uv = helper.GetNodeUV( geomFace, node, inFaceNode );
607 int ind = reverse ? 3-i : i+1;
608 tri.GeomInfoPi(ind).u = uv.X();
609 tri.GeomInfoPi(ind).v = uv.Y();
610 tri.PNum (ind) = ngNodeId( node, ngMesh, nodeNgIdMap );
613 ngMesh.AddSurfaceElement (tri);
614 #ifdef DUMP_TRIANGLES
618 if ( isInternalFace )
620 swap( tri[1], tri[2] );
621 ngMesh.AddSurfaceElement (tri);
622 #ifdef DUMP_TRIANGLES
628 } // case TopAbs_FACE
630 case TopAbs_VERTEX: { // VERTEX
631 // --------------------------
632 SMDS_NodeIteratorPtr nodeIt = smDS->GetNodes();
633 if ( nodeIt->more() )
634 ngNodeId( nodeIt->next(), ngMesh, nodeNgIdMap );
639 } // loop on submeshes
642 nodeVec.resize( ngMesh.GetNP() + 1 );
643 TNode2IdMap::iterator node_NgId, nodeNgIdEnd = nodeNgIdMap.end();
644 for ( node_NgId = nodeNgIdMap.begin(); node_NgId != nodeNgIdEnd; ++node_NgId)
645 nodeVec[ node_NgId->second ] = (SMDS_MeshNode*) node_NgId->first;
650 //================================================================================
652 * \brief Duplicate mesh faces on internal geom faces
654 //================================================================================
656 void NETGENPlugin_Mesher::fixIntFaces(const netgen::OCCGeometry& occgeom,
657 netgen::Mesh& ngMesh,
658 NETGENPlugin_Internals& internalShapes)
660 SMESHDS_Mesh* meshDS = internalShapes.getMesh().GetMeshDS();
662 // find ng indices of internal faces
664 for ( int ngFaceID = 1; ngFaceID <= occgeom.fmap.Extent(); ++ngFaceID )
666 int smeshID = meshDS->ShapeToIndex( occgeom.fmap( ngFaceID ));
667 if ( internalShapes.isInternalShape( smeshID ))
668 ngFaceIds.insert( ngFaceID );
670 if ( !ngFaceIds.empty() )
673 int i, nbFaces = ngMesh.GetNSE();
674 for (int i = 1; i <= nbFaces; ++i)
676 netgen::Element2d elem = ngMesh.SurfaceElement(i);
677 if ( ngFaceIds.count( elem.GetIndex() ))
679 swap( elem[1], elem[2] );
680 ngMesh.AddSurfaceElement (elem);
688 //================================================================================
689 // define gp_XY_Subtracted pointer to function calling gp_XY::Subtracted(gp_XY)
690 gp_XY_FunPtr(Subtracted);
691 //gp_XY_FunPtr(Added);
693 //================================================================================
695 * \brief Evaluate distance between two 2d points along the surface
697 //================================================================================
699 double evalDist( const gp_XY& uv1,
701 const Handle(Geom_Surface)& surf,
702 const int stopHandler=-1)
704 if ( stopHandler > 0 ) // continue recursion
706 gp_XY mid = SMESH_MesherHelper::GetMiddleUV( surf, uv1, uv2 );
707 return evalDist( uv1,mid, surf, stopHandler-1 ) + evalDist( mid,uv2, surf, stopHandler-1 );
709 double dist3D = surf->Value( uv1.X(), uv1.Y() ).Distance( surf->Value( uv2.X(), uv2.Y() ));
710 if ( stopHandler == 0 ) // stop recursion
713 // start recursion if necessary
714 double dist2D = SMESH_MesherHelper::applyIn2D(surf, uv1, uv2, gp_XY_Subtracted, 0).Modulus();
715 if ( fabs( dist3D - dist2D ) < dist2D * 1e-10 )
716 return dist3D; // equal parametrization of a planar surface
718 return evalDist( uv1, uv2, surf, 3 ); // start recursion
721 //================================================================================
723 * \brief Data of vertex internal in geom face
725 //================================================================================
729 gp_XY uv; //!< UV in face parametric space
730 int ngId; //!< ng id of corrsponding node
731 gp_XY uvClose; //!< UV of closest boundary node
732 int ngIdClose; //!< ng id of closest boundary node
735 //================================================================================
737 * \brief Data of vertex internal in solid
739 //================================================================================
743 int ngId; //!< ng id of corrsponding node
744 int ngIdClose; //!< ng id of closest 2d mesh element
745 int ngIdCloseN; //!< ng id of closest node of the closest 2d mesh element
748 inline double dist2(const netgen::MeshPoint& p1, const netgen::MeshPoint& p2)
750 return gp_Pnt( NGPOINT_COORDS(p1)).SquareDistance( gp_Pnt( NGPOINT_COORDS(p2)));
754 //================================================================================
756 * \brief Make netgen take internal vertices in faces into account by adding
757 * segments including internal vertices
759 * This function works in supposition that 1D mesh is already computed in ngMesh
761 //================================================================================
763 void NETGENPlugin_Mesher::addIntVerticesInFaces(const netgen::OCCGeometry& occgeom,
764 netgen::Mesh& ngMesh,
765 vector<const SMDS_MeshNode*>& nodeVec,
766 NETGENPlugin_Internals& internalShapes)
768 if ( nodeVec.size() < ngMesh.GetNP() )
769 nodeVec.resize( ngMesh.GetNP(), 0 );
771 SMESHDS_Mesh* meshDS = internalShapes.getMesh().GetMeshDS();
772 SMESH_MesherHelper helper( internalShapes.getMesh() );
774 const map<int,list<int> >& face2Vert = internalShapes.getFacesWithVertices();
775 map<int,list<int> >::const_iterator f2v = face2Vert.begin();
776 for ( ; f2v != face2Vert.end(); ++f2v )
778 const TopoDS_Face& face = TopoDS::Face( meshDS->IndexToShape( f2v->first ));
779 if ( face.IsNull() ) continue;
780 int faceNgID = occgeom.fmap.FindIndex (face);
781 if ( faceNgID < 0 ) continue;
784 Handle(Geom_Surface) surf = BRep_Tool::Surface(face,loc);
786 helper.SetSubShape( face );
787 helper.SetElementsOnShape( true );
789 // Get data of internal vertices and add them to ngMesh
791 multimap< double, TIntVData > dist2VData; // sort vertices by distance from boundary nodes
793 int i, nbSegInit = ngMesh.GetNSeg();
795 // boundary characteristics
796 double totSegLen2D = 0;
799 const list<int>& iVertices = f2v->second;
800 list<int>::const_iterator iv = iVertices.begin();
801 for ( int nbV = 0; iv != iVertices.end(); ++iv, nbV++ )
804 // get node on vertex
805 const TopoDS_Vertex V = TopoDS::Vertex( meshDS->IndexToShape( *iv ));
806 const SMDS_MeshNode * nV = SMESH_Algo::VertexNode( V, meshDS );
809 SMESH_subMesh* sm = helper.GetMesh()->GetSubMesh( V );
810 sm->ComputeStateEngine( SMESH_subMesh::COMPUTE );
811 nV = SMESH_Algo::VertexNode( V, meshDS );
815 netgen::MeshPoint mp( netgen::Point<3> (nV->X(), nV->Y(), nV->Z()) );
816 ngMesh.AddPoint ( mp, 1, netgen::EDGEPOINT );
817 vData.ngId = ngMesh.GetNP();
818 nodeVec.push_back( nV );
822 vData.uv = helper.GetNodeUV( face, nV, 0, &uvOK );
823 if ( !uvOK ) helper.CheckNodeUV( face, nV, vData.uv, BRep_Tool::Tolerance(V),/*force=*/1);
825 // loop on all segments of the face to find the node closest to vertex and to count
826 // average segment 2d length
827 double closeDist2 = numeric_limits<double>::max(), dist2;
829 for (i = 1; i <= ngMesh.GetNSeg(); ++i)
831 netgen::Segment & seg = ngMesh.LineSegment(i);
832 if ( seg.si != faceNgID ) continue;
834 for ( int iEnd = 0; iEnd < 2; ++iEnd)
836 uv[iEnd].SetCoord( seg.epgeominfo[iEnd].u, seg.epgeominfo[iEnd].v );
837 if ( ngIdLast == seg[ iEnd ] ) continue;
838 dist2 = helper.applyIn2D(surf, uv[iEnd], vData.uv, gp_XY_Subtracted,0).SquareModulus();
839 if ( dist2 < closeDist2 )
840 vData.ngIdClose = seg[ iEnd ], vData.uvClose = uv[iEnd], closeDist2 = dist2;
841 ngIdLast = seg[ iEnd ];
845 totSegLen2D += helper.applyIn2D(surf, uv[0], uv[1], gp_XY_Subtracted, false).Modulus();
849 dist2VData.insert( make_pair( closeDist2, vData ));
852 if ( totNbSeg == 0 ) break;
853 double avgSegLen2d = totSegLen2D / totNbSeg;
855 // Loop on vertices to add segments
857 multimap< double, TIntVData >::iterator dist_vData = dist2VData.begin();
858 for ( ; dist_vData != dist2VData.end(); ++dist_vData )
860 double closeDist2 = dist_vData->first, dist2;
861 TIntVData & vData = dist_vData->second;
863 // try to find more close node among segments added for internal vertices
864 for (i = nbSegInit+1; i <= ngMesh.GetNSeg(); ++i)
866 netgen::Segment & seg = ngMesh.LineSegment(i);
867 if ( seg.si != faceNgID ) continue;
869 for ( int iEnd = 0; iEnd < 2; ++iEnd)
871 uv[iEnd].SetCoord( seg.epgeominfo[iEnd].u, seg.epgeominfo[iEnd].v );
872 dist2 = helper.applyIn2D(surf, uv[iEnd], vData.uv, gp_XY_Subtracted,0).SquareModulus();
873 if ( dist2 < closeDist2 )
874 vData.ngIdClose = seg[ iEnd ], vData.uvClose = uv[iEnd], closeDist2 = dist2;
877 // decide whether to use the closest node as the second end of segment or to
878 // create a new point
879 int segEnd1 = vData.ngId;
880 int segEnd2 = vData.ngIdClose; // to use closest node
881 gp_XY uvV = vData.uv, uvP = vData.uvClose;
882 double segLenHint = ngMesh.GetH( ngMesh.Point( vData.ngId ));
883 double nodeDist2D = sqrt( closeDist2 );
884 double nodeDist3D = evalDist( vData.uv, vData.uvClose, surf );
885 bool avgLenOK = ( avgSegLen2d < 0.75 * nodeDist2D );
886 bool hintLenOK = ( segLenHint < 0.75 * nodeDist3D );
887 //cout << "uvV " << uvV.X() <<","<<uvV.Y() << " ";
888 if ( hintLenOK || avgLenOK )
890 // create a point between the closest node and V
893 double r = min( 0.5, ( hintLenOK ? segLenHint/nodeDist3D : avgSegLen2d/nodeDist2D ));
894 // direction from V to closet node in 2D
895 gp_Dir2d v2n( helper.applyIn2D(surf, uvP, uvV, gp_XY_Subtracted, false ));
897 uvP = vData.uv + r * nodeDist2D * v2n.XY();
898 gp_Pnt P = surf->Value( uvP.X(), uvP.Y() ).Transformed( loc );
900 netgen::MeshPoint mp( netgen::Point<3> (P.X(), P.Y(), P.Z()));
901 ngMesh.AddPoint ( mp, 1, netgen::EDGEPOINT );
902 segEnd2 = ngMesh.GetNP();
903 //cout << "Middle " << r << " uv " << uvP.X() << "," << uvP.Y() << "( " << ngMesh.Point(segEnd2).X()<<","<<ngMesh.Point(segEnd2).Y()<<","<<ngMesh.Point(segEnd2).Z()<<" )"<< endl;
904 SMDS_MeshNode * nP = helper.AddNode(P.X(), P.Y(), P.Z());
905 nodeVec.push_back( nP );
907 //else cout << "at Node " << " uv " << uvP.X() << "," << uvP.Y() << endl;
912 if ( segEnd1 > segEnd2 ) swap( segEnd1, segEnd2 ), swap( uvV, uvP );
913 seg[0] = segEnd1; // ng node id
914 seg[1] = segEnd2; // ng node id
915 seg.edgenr = ngMesh.GetNSeg() + 1;// segment id
918 seg.epgeominfo[ 0 ].dist = 0; // param on curve
919 seg.epgeominfo[ 0 ].u = uvV.X();
920 seg.epgeominfo[ 0 ].v = uvV.Y();
921 seg.epgeominfo[ 1 ].dist = 1; // param on curve
922 seg.epgeominfo[ 1 ].u = uvP.X();
923 seg.epgeominfo[ 1 ].v = uvP.Y();
925 // seg.epgeominfo[ 0 ].edgenr = 10; // = geom.emap.FindIndex(edge);
926 // seg.epgeominfo[ 1 ].edgenr = 10; // = geom.emap.FindIndex(edge);
928 ngMesh.AddSegment (seg);
930 // add reverse segment
931 swap (seg[0], seg[1]);
932 swap( seg.epgeominfo[0], seg.epgeominfo[1] );
933 seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
934 ngMesh.AddSegment (seg);
940 //================================================================================
942 * \brief Make netgen take internal vertices in solids into account by adding
943 * faces including internal vertices
945 * This function works in supposition that 2D mesh is already computed in ngMesh
947 //================================================================================
949 void NETGENPlugin_Mesher::addIntVerticesInSolids(const netgen::OCCGeometry& occgeom,
950 netgen::Mesh& ngMesh,
951 vector<const SMDS_MeshNode*>& nodeVec,
952 NETGENPlugin_Internals& internalShapes)
954 #ifdef DUMP_TRIANGLES_SCRIPT
955 // create a python script making a mesh containing triangles added for internal vertices
956 ofstream py(DUMP_TRIANGLES_SCRIPT);
957 py << "from smesh import * "<< endl
958 << "m = Mesh(name='triangles')" << endl;
960 if ( nodeVec.size() < ngMesh.GetNP() )
961 nodeVec.resize( ngMesh.GetNP(), 0 );
963 SMESHDS_Mesh* meshDS = internalShapes.getMesh().GetMeshDS();
964 SMESH_MesherHelper helper( internalShapes.getMesh() );
966 const map<int,list<int> >& so2Vert = internalShapes.getSolidsWithVertices();
967 map<int,list<int> >::const_iterator s2v = so2Vert.begin();
968 for ( ; s2v != so2Vert.end(); ++s2v )
970 const TopoDS_Shape& solid = meshDS->IndexToShape( s2v->first );
971 if ( solid.IsNull() ) continue;
972 int solidNgID = occgeom.somap.FindIndex (solid);
973 if ( solidNgID < 0 && !occgeom.somap.IsEmpty() ) continue;
975 helper.SetSubShape( solid );
976 helper.SetElementsOnShape( true );
978 // find ng indices of faces within the solid
980 for (TopExp_Explorer fExp(solid, TopAbs_FACE); fExp.More(); fExp.Next() )
981 ngFaceIds.insert( occgeom.fmap.FindIndex( fExp.Current() ));
982 if ( ngFaceIds.size() == 1 && *ngFaceIds.begin() == 0 )
983 ngFaceIds.insert( 1 );
985 // Get data of internal vertices and add them to ngMesh
987 multimap< double, TIntVSoData > dist2VData; // sort vertices by distance from ng faces
989 int i, nbFaceInit = ngMesh.GetNSE();
991 // boundary characteristics
992 double totSegLen = 0;
995 const list<int>& iVertices = s2v->second;
996 list<int>::const_iterator iv = iVertices.begin();
997 for ( int nbV = 0; iv != iVertices.end(); ++iv, nbV++ )
1000 const TopoDS_Vertex V = TopoDS::Vertex( meshDS->IndexToShape( *iv ));
1002 // get node on vertex
1003 const SMDS_MeshNode * nV = SMESH_Algo::VertexNode( V, meshDS );
1006 SMESH_subMesh* sm = helper.GetMesh()->GetSubMesh( V );
1007 sm->ComputeStateEngine( SMESH_subMesh::COMPUTE );
1008 nV = SMESH_Algo::VertexNode( V, meshDS );
1009 if ( !nV ) continue;
1012 netgen::MeshPoint mpV( netgen::Point<3> (nV->X(), nV->Y(), nV->Z()) );
1013 ngMesh.AddPoint ( mpV, 1, netgen::FIXEDPOINT );
1014 vData.ngId = ngMesh.GetNP();
1015 nodeVec.push_back( nV );
1017 // loop on all 2d elements to find the one closest to vertex and to count
1018 // average segment length
1019 double closeDist2 = numeric_limits<double>::max(), avgDist2;
1020 for (i = 1; i <= ngMesh.GetNSE(); ++i)
1022 const netgen::Element2d& elem = ngMesh.SurfaceElement(i);
1023 if ( !ngFaceIds.count( elem.GetIndex() )) continue;
1025 multimap< double, int> dist2nID; // sort nodes of element by distance from V
1026 for ( int j = 0; j < elem.GetNP(); ++j)
1028 netgen::MeshPoint mp = ngMesh.Point( elem[j] );
1029 double d2 = dist2( mpV, mp );
1030 dist2nID.insert( make_pair( d2, elem[j] ));
1031 avgDist2 += d2 / elem.GetNP();
1033 totNbSeg++, totSegLen+= sqrt( dist2( mp, ngMesh.Point( elem[(j+1)%elem.GetNP()])));
1035 double dist = dist2nID.begin()->first; //avgDist2;
1036 if ( dist < closeDist2 )
1037 vData.ngIdClose= i, vData.ngIdCloseN= dist2nID.begin()->second, closeDist2= dist;
1039 dist2VData.insert( make_pair( closeDist2, vData ));
1042 if ( totNbSeg == 0 ) break;
1043 double avgSegLen = totSegLen / totNbSeg;
1045 // Loop on vertices to add triangles
1047 multimap< double, TIntVSoData >::iterator dist_vData = dist2VData.begin();
1048 for ( ; dist_vData != dist2VData.end(); ++dist_vData )
1050 double closeDist2 = dist_vData->first;
1051 TIntVSoData & vData = dist_vData->second;
1053 const netgen::MeshPoint& mpV = ngMesh.Point( vData.ngId );
1055 // try to find more close face among ones added for internal vertices
1056 for (i = nbFaceInit+1; i <= ngMesh.GetNSE(); ++i)
1058 double avgDist2 = 0;
1059 multimap< double, int> dist2nID;
1060 const netgen::Element2d& elem = ngMesh.SurfaceElement(i);
1061 for ( int j = 0; j < elem.GetNP(); ++j)
1063 double d = dist2( mpV, ngMesh.Point( elem[j] ));
1064 dist2nID.insert( make_pair( d, elem[j] ));
1065 avgDist2 += d / elem.GetNP();
1066 if ( avgDist2 < closeDist2 )
1067 vData.ngIdClose= i, vData.ngIdCloseN= dist2nID.begin()->second, closeDist2= avgDist2;
1070 // sort nodes of the closest face by angle with vector from V to the closest node
1071 const double tol = numeric_limits<double>::min();
1072 map< double, int > angle2ID;
1073 const netgen::Element2d& closeFace = ngMesh.SurfaceElement( vData.ngIdClose );
1074 netgen::MeshPoint mp[2];
1075 mp[0] = ngMesh.Point( vData.ngIdCloseN );
1076 gp_XYZ p1( NGPOINT_COORDS( mp[0] ));
1077 gp_XYZ pV( NGPOINT_COORDS( mpV ));
1078 gp_Vec v2p1( pV, p1 );
1079 double distN1 = v2p1.Magnitude();
1080 if ( distN1 <= tol ) continue;
1082 for ( int j = 0; j < closeFace.GetNP(); ++j)
1084 mp[1] = ngMesh.Point( closeFace[j] );
1085 gp_Vec v2p( pV, gp_Pnt( NGPOINT_COORDS( mp[1] )) );
1086 angle2ID.insert( make_pair( v2p1.Angle( v2p ), closeFace[j]));
1088 // get node with angle of 60 degrees or greater
1089 map< double, int >::iterator angle_id = angle2ID.lower_bound( 60*PI180 );
1090 if ( angle_id == angle2ID.end() ) angle_id = --angle2ID.end();
1091 const double minAngle = 30 * PI180;
1092 const double angle = angle_id->first;
1093 bool angleOK = ( angle > minAngle );
1095 // find points to create a triangle
1096 netgen::Element2d tri(3);
1098 tri[0] = vData.ngId;
1099 tri[1] = vData.ngIdCloseN; // to use the closest nodes
1100 tri[2] = angle_id->second; // to use the node with best angle
1102 // decide whether to use the closest node and the node with best angle or to create new ones
1103 for ( int isBestAngleN = 0; isBestAngleN < 2; ++isBestAngleN )
1105 bool createNew = !angleOK, distOK = true;
1107 int triInd = isBestAngleN ? 2 : 1;
1108 mp[isBestAngleN] = ngMesh.Point( tri[triInd] );
1113 double distN2 = sqrt( dist2( mpV, mp[isBestAngleN]));
1114 createNew = ( fabs( distN2 - distN1 ) > 0.25 * distN1 );
1116 else if ( angle < tol )
1118 v2p1.SetX( v2p1.X() + 1e-3 );
1124 double segLenHint = ngMesh.GetH( ngMesh.Point( vData.ngId ));
1125 bool avgLenOK = ( avgSegLen < 0.75 * distN1 );
1126 bool hintLenOK = ( segLenHint < 0.75 * distN1 );
1127 createNew = (createNew || avgLenOK || hintLenOK );
1128 // we create a new node not closer than 0.5 to the closest face
1129 // in order not to clash with other close face
1130 double r = min( 0.5, ( hintLenOK ? segLenHint : avgSegLen ) / distN1 );
1131 distFromV = r * distN1;
1135 // create a new point, between the node and the vertex if angleOK
1136 gp_XYZ p( NGPOINT_COORDS( mp[isBestAngleN] ));
1137 gp_Vec v2p( pV, p ); v2p.Normalize();
1138 if ( isBestAngleN && !angleOK )
1139 p = p1 + gp_Dir( v2p.XYZ() - v2p1.XYZ()).XYZ() * distN1 * 0.95;
1141 p = pV + v2p.XYZ() * distFromV;
1143 if ( !isBestAngleN ) p1 = p, distN1 = distFromV;
1145 mp[isBestAngleN].SetPoint( netgen::Point<3> (p.X(), p.Y(), p.Z()));
1146 ngMesh.AddPoint ( mp[isBestAngleN], 1, netgen::SURFACEPOINT );
1147 tri[triInd] = ngMesh.GetNP();
1148 nodeVec.push_back( helper.AddNode( p.X(), p.Y(), p.Z()) );
1151 ngMesh.AddSurfaceElement (tri);
1152 swap( tri[1], tri[2] );
1153 ngMesh.AddSurfaceElement (tri);
1155 #ifdef DUMP_TRIANGLES_SCRIPT
1156 py << "n1 = m.AddNode( "<< mpV.X()<<", "<< mpV.Y()<<", "<< mpV.Z()<<") "<< endl
1157 << "n2 = m.AddNode( "<< mp[0].X()<<", "<< mp[0].Y()<<", "<< mp[0].Z()<<") "<< endl
1158 << "n3 = m.AddNode( "<< mp[1].X()<<", "<< mp[1].Y()<<", "<< mp[1].Z()<<" )" << endl
1159 << "m.AddFace([n1,n2,n3])" << endl;
1161 } // loop on internal vertices of a solid
1163 } // loop on solids with internal vertices
1166 //================================================================================
1168 * \brief Fill SMESH mesh according to contents of netgen mesh
1169 * \param occgeo - container of OCCT geometry to mesh
1170 * \param ngMesh - netgen mesh
1171 * \param initState - bn of entities in netgen mesh before computing
1172 * \param sMesh - SMESH mesh to fill in
1173 * \param nodeVec - vector of nodes in which node index == netgen ID
1174 * \retval int - error
1176 //================================================================================
1178 int NETGENPlugin_Mesher::FillSMesh(const netgen::OCCGeometry& occgeo,
1179 const netgen::Mesh& ngMesh,
1180 const NETGENPlugin_ngMeshInfo& initState,
1182 std::vector<const SMDS_MeshNode*>& nodeVec,
1183 SMESH_Comment& comment)
1185 int nbNod = ngMesh.GetNP();
1186 int nbSeg = ngMesh.GetNSeg();
1187 int nbFac = ngMesh.GetNSE();
1188 int nbVol = ngMesh.GetNE();
1190 SMESHDS_Mesh* meshDS = sMesh.GetMeshDS();
1192 // map of nodes assigned to submeshes
1193 NCollection_Map<int> pindMap;
1194 // create and insert nodes into nodeVec
1195 nodeVec.resize( nbNod + 1 );
1196 int i, nbInitNod = initState._nbNodes;
1197 for (i = nbInitNod+1; i <= nbNod; ++i )
1199 const netgen::MeshPoint& ngPoint = ngMesh.Point(i);
1200 SMDS_MeshNode* node = NULL;
1201 TopoDS_Vertex aVert;
1202 // First, netgen creates nodes on vertices in occgeo.vmap,
1203 // so node index corresponds to vertex index
1204 // but (isuue 0020776) netgen does not create nodes with equal coordinates
1205 if ( i-nbInitNod <= occgeo.vmap.Extent() )
1207 gp_Pnt p ( NGPOINT_COORDS(ngPoint) );
1208 for (int iV = i-nbInitNod; aVert.IsNull() && iV <= occgeo.vmap.Extent(); ++iV)
1210 aVert = TopoDS::Vertex( occgeo.vmap( iV ) );
1211 gp_Pnt pV = BRep_Tool::Pnt( aVert );
1212 if ( p.SquareDistance( pV ) > 1e-20 )
1215 node = const_cast<SMDS_MeshNode*>( SMESH_Algo::VertexNode( aVert, meshDS ));
1218 if (node) // node found on vertex
1222 node = meshDS->AddNode( NGPOINT_COORDS( ngPoint ));
1223 if (!aVert.IsNull())
1226 meshDS->SetNodeOnVertex(node, aVert);
1233 // create mesh segments along geometric edges
1234 NCollection_Map<Link> linkMap;
1235 int nbInitSeg = initState._nbSegments;
1236 for (i = nbInitSeg+1; i <= nbSeg; ++i )
1238 const netgen::Segment& seg = ngMesh.LineSegment(i);
1239 Link link(seg[0], seg[1]);
1240 if (!linkMap.Add(link))
1244 int pinds[3] = { seg.pnums[0], seg.pnums[1], seg.pnums[2] };
1246 int pinds[3] = { seg.p1, seg.p2, seg.pmid };
1250 for (int j=0; j < 3; ++j)
1252 int pind = pinds[j];
1253 if (pind <= 0) continue;
1260 int aGeomEdgeInd = seg.epgeominfo[j].edgenr;
1261 if (aGeomEdgeInd > 0 && aGeomEdgeInd <= occgeo.emap.Extent())
1262 aEdge = TopoDS::Edge(occgeo.emap(aGeomEdgeInd));
1264 param = seg.epgeominfo[j].dist;
1268 param = param2 * 0.5;
1269 if (pind <= nbInitNod || pindMap.Contains(pind))
1271 if (!aEdge.IsNull())
1273 meshDS->SetNodeOnEdge(nodeVec_ACCESS(pind), aEdge, param);
1277 SMDS_MeshEdge* edge;
1278 if (nbp < 3) // second order ?
1279 edge = meshDS->AddEdge(nodeVec_ACCESS(pinds[0]), nodeVec_ACCESS(pinds[1]));
1281 edge = meshDS->AddEdge(nodeVec_ACCESS(pinds[0]), nodeVec_ACCESS(pinds[1]),
1282 nodeVec_ACCESS(pinds[2]));
1285 if ( !comment.size() ) comment << "Cannot create a mesh edge";
1286 MESSAGE("Cannot create a mesh edge");
1287 nbSeg = nbFac = nbVol = 0;
1290 if (!aEdge.IsNull())
1291 meshDS->SetMeshElementOnShape(edge, aEdge);
1294 // create mesh faces along geometric faces
1295 int nbInitFac = initState._nbFaces;
1296 for (i = nbInitFac+1; i <= nbFac; ++i )
1298 const netgen::Element2d& elem = ngMesh.SurfaceElement(i);
1299 int aGeomFaceInd = elem.GetIndex();
1301 if (aGeomFaceInd > 0 && aGeomFaceInd <= occgeo.fmap.Extent())
1302 aFace = TopoDS::Face(occgeo.fmap(aGeomFaceInd));
1303 vector<SMDS_MeshNode*> nodes;
1304 for (int j=1; j <= elem.GetNP(); ++j)
1306 int pind = elem.PNum(j);
1307 SMDS_MeshNode* node = nodeVec_ACCESS(pind);
1308 nodes.push_back(node);
1309 if (pind <= nbInitNod || pindMap.Contains(pind))
1311 if (!aFace.IsNull())
1313 const netgen::PointGeomInfo& pgi = elem.GeomInfoPi(j);
1314 meshDS->SetNodeOnFace(node, aFace, pgi.u, pgi.v);
1318 SMDS_MeshFace* face = NULL;
1319 switch (elem.GetType())
1322 face = meshDS->AddFace(nodes[0],nodes[1],nodes[2]);
1325 face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3]);
1328 face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[5],nodes[3],nodes[4]);
1331 face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3],
1332 nodes[4],nodes[7],nodes[5],nodes[6]);
1335 MESSAGE("NETGEN created a face of unexpected type, ignoring");
1340 if ( !comment.size() ) comment << "Cannot create a mesh face";
1341 MESSAGE("Cannot create a mesh face");
1342 nbSeg = nbFac = nbVol = 0;
1345 if (!aFace.IsNull())
1346 meshDS->SetMeshElementOnShape(face, aFace);
1349 // create tetrahedra
1350 for (i = 1; i <= nbVol/* && isOK*/; ++i)
1352 const netgen::Element& elem = ngMesh.VolumeElement(i);
1353 int aSolidInd = elem.GetIndex();
1354 TopoDS_Solid aSolid;
1355 if (aSolidInd > 0 && aSolidInd <= occgeo.somap.Extent())
1356 aSolid = TopoDS::Solid(occgeo.somap(aSolidInd));
1357 vector<SMDS_MeshNode*> nodes;
1358 for (int j=1; j <= elem.GetNP(); ++j)
1360 int pind = elem.PNum(j);
1361 SMDS_MeshNode* node = nodeVec_ACCESS(pind);
1362 nodes.push_back(node);
1363 if (pind <= nbInitNod || pindMap.Contains(pind))
1365 if (!aSolid.IsNull())
1368 meshDS->SetNodeInVolume(node, aSolid);
1372 SMDS_MeshVolume* vol = NULL;
1373 switch (elem.GetType())
1376 vol = meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3]);
1379 vol = meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3],
1380 nodes[4],nodes[7],nodes[5],nodes[6],nodes[8],nodes[9]);
1383 MESSAGE("NETGEN created a volume of unexpected type, ignoring");
1388 if ( !comment.size() ) comment << "Cannot create a mesh volume";
1389 MESSAGE("Cannot create a mesh volume");
1390 nbSeg = nbFac = nbVol = 0;
1393 if (!aSolid.IsNull())
1394 meshDS->SetMeshElementOnShape(vol, aSolid);
1396 return comment.empty() ? 0 : 1;
1399 //=============================================================================
1401 * Here we are going to use the NETGEN mesher
1403 //=============================================================================
1405 bool NETGENPlugin_Mesher::Compute()
1407 NETGENPlugin_NetgenLibWrapper ngLib;
1409 netgen::MeshingParameters& mparams = netgen::mparam;
1410 MESSAGE("Compute with:\n"
1411 " max size = " << mparams.maxh << "\n"
1412 " segments per edge = " << mparams.segmentsperedge);
1414 " growth rate = " << mparams.grading << "\n"
1415 " elements per radius = " << mparams.curvaturesafety << "\n"
1416 " second order = " << mparams.secondorder << "\n"
1417 " quad allowed = " << mparams.quad);
1419 SMESH_ComputeErrorPtr error = SMESH_ComputeError::New();
1422 // -------------------------
1423 // Prepare OCC geometry
1424 // -------------------------
1426 netgen::OCCGeometry occgeo;
1427 list< SMESH_subMesh* > meshedSM;
1428 NETGENPlugin_Internals internals( *_mesh, _shape, _isVolume );
1429 PrepareOCCgeometry( occgeo, _shape, *_mesh, &meshedSM, &internals );
1431 // -------------------------
1432 // Local size on faces
1433 // -------------------------
1438 for(std::map<int,double>::const_iterator it=FaceId2LocalSize.begin(); it!=FaceId2LocalSize.end(); it++)
1440 int key = (*it).first;
1441 double val = (*it).second;
1442 const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
1443 int faceID = occgeo.fmap.FindIndex(shape);
1444 occgeo.SetFaceMaxH(faceID, val);
1449 // -------------------------
1450 // Generate the mesh
1451 // -------------------------
1453 netgen::Mesh *ngMesh = NULL;
1454 NETGENPlugin_ngMeshInfo initState;
1456 SMESH_Comment comment;
1459 // vector of nodes in which node index == netgen ID
1460 vector< const SMDS_MeshNode* > nodeVec;
1466 // Pass 1D simple parameters to NETGEN
1468 if ( int nbSeg = _simpleHyp->GetNumberOfSegments() ) {
1470 mparams.segmentsperedge = nbSeg + 0.1;
1471 mparams.maxh = occgeo.boundingbox.Diam();
1472 mparams.grading = 0.01;
1476 mparams.segmentsperedge = 1;
1477 mparams.maxh = _simpleHyp->GetLocalLength();
1480 // Let netgen create ngMesh and calculate element size on not meshed shapes
1482 int startWith = netgen::MESHCONST_ANALYSE;
1483 int endWith = netgen::MESHCONST_ANALYSE;
1484 err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
1485 if (err) comment << "Error in netgen::OCCGenerateMesh() at MESHCONST_ANALYSE step";
1486 ngLib.setMesh(( Ng_Mesh*) ngMesh );
1488 // --------------------------------
1489 // Local size on vertices and edges
1490 // --------------------------------
1494 for(std::map<int,double>::const_iterator it=EdgeId2LocalSize.begin(); it!=EdgeId2LocalSize.end(); it++)
1496 int key = (*it).first;
1497 double hi = (*it).second;
1498 const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
1499 const TopoDS_Edge& e = TopoDS::Edge(shape);
1500 Standard_Real u1, u2;
1501 Handle(Geom_Curve) curve = BRep_Tool::Curve(e, u1, u2);
1502 GeomAdaptor_Curve AdaptCurve(curve);
1503 double length = GCPnts_AbscissaPoint::Length(AdaptCurve, u1, u2);
1506 Standard_Real delta = (u2-u1)/nb;
1507 for(int i=0; i<nb; i++)
1509 Standard_Real u = u1 + delta*i;
1510 gp_Pnt p = curve->Value(u);
1511 netgen::Point3d pi(p.X(), p.Y(), p.Z());
1512 ngMesh->RestrictLocalH(pi, hi);
1515 for(std::map<int,double>::const_iterator it=VertexId2LocalSize.begin(); it!=VertexId2LocalSize.end(); it++)
1517 int key = (*it).first;
1518 double hi = (*it).second;
1519 const TopoDS_Shape& shape = ShapesWithLocalSize.FindKey(key);
1520 const TopoDS_Vertex& v = TopoDS::Vertex(shape);
1521 gp_Pnt p = BRep_Tool::Pnt(v);
1522 netgen::Point3d pi(p.X(), p.Y(), p.Z());
1523 ngMesh->RestrictLocalH(pi, hi);
1527 // Precompute internal edges (issue 0020676) in order to
1528 // add mesh on them correctly (twice) to netgen mesh
1529 if ( !err && internals.hasInternalEdges() )
1531 // load internal shapes into OCCGeometry
1532 netgen::OCCGeometry intOccgeo;
1533 internals.getInternalEdges( intOccgeo.fmap, intOccgeo.emap, intOccgeo.vmap, meshedSM);
1534 intOccgeo.boundingbox = occgeo.boundingbox;
1535 intOccgeo.shape = occgeo.shape;
1537 // let netgen compute element size by the main geometry in temporary mesh
1538 netgen::Mesh *tmpNgMesh = NULL;
1539 netgen::OCCGenerateMesh(occgeo, tmpNgMesh, startWith, endWith, optstr);
1540 // compute mesh on internal edges
1541 endWith = netgen::MESHCONST_MESHEDGES;
1542 err = netgen::OCCGenerateMesh(intOccgeo, tmpNgMesh, startWith, endWith, optstr);
1543 if (err) comment << "Error in netgen::OCCGenerateMesh() at meshing internal edges";
1545 // fill SMESH by netgen mesh
1546 vector< const SMDS_MeshNode* > tmpNodeVec;
1547 FillSMesh( intOccgeo, *tmpNgMesh, initState, *_mesh, tmpNodeVec, comment );
1548 err = ( !comment.empty() );
1550 nglib::Ng_DeleteMesh((nglib::Ng_Mesh*)tmpNgMesh);
1553 // Fill ngMesh with nodes and elements of computed submeshes
1556 _faceDescriptors.clear();
1557 err = ! fillNgMesh(occgeo, *ngMesh, nodeVec, meshedSM);
1559 initState = NETGENPlugin_ngMeshInfo(ngMesh);
1564 startWith = endWith = netgen::MESHCONST_MESHEDGES;
1565 err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
1566 if (err) comment << "Error in netgen::OCCGenerateMesh() at 1D mesh generation";
1568 // ---------------------
1569 // compute surface mesh
1570 // ---------------------
1573 // Pass 2D simple parameters to NETGEN
1575 if ( double area = _simpleHyp->GetMaxElementArea() ) {
1577 mparams.maxh = sqrt(2. * area/sqrt(3.0));
1578 mparams.grading = 0.4; // moderate size growth
1581 // length from edges
1582 if ( ngMesh->GetNSeg() ) {
1583 double edgeLength = 0;
1584 TopTools_MapOfShape visitedEdges;
1585 for ( TopExp_Explorer exp( _shape, TopAbs_EDGE ); exp.More(); exp.Next() )
1586 if( visitedEdges.Add(exp.Current()) )
1587 edgeLength += SMESH_Algo::EdgeLength( TopoDS::Edge( exp.Current() ));
1588 // we have to multiply length by 2 since for each TopoDS_Edge there
1589 // are double set of NETGEN edges or, in other words, we have to
1590 // divide ngMesh->GetNSeg() by 2.
1591 mparams.maxh = 2*edgeLength / ngMesh->GetNSeg();
1594 mparams.maxh = 1000;
1596 mparams.grading = 0.2; // slow size growth
1598 mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
1599 ngMesh->SetGlobalH (mparams.maxh);
1600 netgen::Box<3> bb = occgeo.GetBoundingBox();
1601 bb.Increase (bb.Diam()/20);
1602 ngMesh->SetLocalH (bb.PMin(), bb.PMax(), mparams.grading);
1605 // Care of vertices internal in faces (issue 0020676)
1606 if ( internals.hasInternalVertexInFace() )
1608 // store computed segments in SMESH in order not to create SMESH
1609 // edges for ng segments added by addIntVerticesInFaces()
1610 FillSMesh( occgeo, *ngMesh, initState, *_mesh, nodeVec, comment );
1611 // add segments to faces with internal vertices
1612 addIntVerticesInFaces( occgeo, *ngMesh, nodeVec, internals );
1613 initState = NETGENPlugin_ngMeshInfo(ngMesh);
1616 // Precompute internal faces (issue 0020676) in order to
1617 // add mesh on them correctly (twice to emulate the crack) to netgen mesh
1618 //if ( internals.hasInternalFaces() )
1620 // // fill SMESH with generated segments
1621 // FillSMesh( occgeo, *ngMesh, initState, *_mesh, nodeVec, comment );
1623 // // load internal shapes into a separate OCCGeometry
1624 // netgen::OCCGeometry intOccgeo;
1625 // list< SMESH_subMesh* > boundarySM;
1626 // internals.getInternalFaces( intOccgeo.fmap, intOccgeo.emap, meshedSM, boundarySM);
1627 // intOccgeo.boundingbox = occgeo.boundingbox;
1628 // intOccgeo.shape = occgeo.shape;
1629 // intOccgeo.facemeshstatus.SetSize (intOccgeo.fmap.Extent());
1630 // intOccgeo.facemeshstatus = 0;
1632 // // let netgen compute element size by the main geometry in temporary mesh
1633 // int start = netgen::MESHCONST_ANALYSE, end = netgen::MESHCONST_ANALYSE;
1634 // netgen::Mesh *tmpNgMesh = NULL;
1635 // netgen::OCCGenerateMesh(occgeo, tmpNgMesh, start, end, optstr);
1637 // // add already computed elements from submeshes of internal faces to tmpNgMesh
1638 // vector< const SMDS_MeshNode* > tmpNodeVec;
1639 // fillNgMesh(intOccgeo, *tmpNgMesh, tmpNodeVec, boundarySM);
1640 // addIntVerticesInFaces( intOccgeo, *tmpNgMesh, tmpNodeVec, internals );
1642 // // compute mesh on internal faces
1643 // NETGENPlugin_ngMeshInfo prevState(tmpNgMesh);
1644 // start = netgen::MESHCONST_MESHEDGES;
1645 // end = netgen::MESHCONST_MESHSURFACE;
1646 // err = netgen::OCCGenerateMesh(intOccgeo, tmpNgMesh, start, end, optstr);
1647 // if (err) comment << "Error in netgen::OCCGenerateMesh() at meshing internal faces";
1649 // // fill SMESH with computed elements
1650 // FillSMesh( intOccgeo, *tmpNgMesh, prevState, *_mesh, tmpNodeVec, comment );
1651 // err = ( !comment.empty() );
1653 // // finally, correctly add elements on internal faces to netgen mesh
1654 // err = ! fillNgMesh(occgeo, *ngMesh, nodeVec, meshedSM);
1655 // initState = NETGENPlugin_ngMeshInfo(ngMesh);
1657 // nglib::Ng_DeleteMesh((nglib::Ng_Mesh*)tmpNgMesh);
1660 // Let netgen compute 2D mesh
1661 startWith = netgen::MESHCONST_MESHSURFACE;
1662 endWith = _optimize ? netgen::MESHCONST_OPTSURFACE : netgen::MESHCONST_MESHSURFACE;
1663 err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
1664 if (err) comment << "Error in netgen::OCCGenerateMesh() at surface mesh generation";
1666 // ---------------------
1667 // generate volume mesh
1668 // ---------------------
1669 if (!err && _isVolume)
1671 // Add ng face descriptors of meshed faces
1672 map< int, pair<int,int> >::iterator fId_soIds = _faceDescriptors.begin();
1673 for ( ; fId_soIds != _faceDescriptors.end(); ++fId_soIds ) {
1674 int faceID = fId_soIds->first;
1675 int solidID1 = fId_soIds->second.first;
1676 int solidID2 = fId_soIds->second.second;
1677 ngMesh->AddFaceDescriptor (netgen::FaceDescriptor(faceID, solidID1, solidID2, 0));
1679 // Pass 3D simple parameters to NETGEN
1680 const NETGENPlugin_SimpleHypothesis_3D* simple3d =
1681 dynamic_cast< const NETGENPlugin_SimpleHypothesis_3D* > ( _simpleHyp );
1683 if ( double vol = simple3d->GetMaxElementVolume() ) {
1685 mparams.maxh = pow( 72, 1/6. ) * pow( vol, 1/3. );
1686 mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
1689 // length from faces
1690 mparams.maxh = ngMesh->AverageH();
1692 // netgen::ARRAY<double> maxhdom;
1693 // maxhdom.SetSize (occgeo.NrSolids());
1694 // maxhdom = mparams.maxh;
1695 // ngMesh->SetMaxHDomain (maxhdom);
1696 ngMesh->SetGlobalH (mparams.maxh);
1697 mparams.grading = 0.4;
1698 ngMesh->CalcLocalH();
1700 // Care of vertices internal in solids and internal faces (issue 0020676)
1701 if ( internals.hasInternalVertexInSolid() || internals.hasInternalFaces() )
1703 // store computed faces in SMESH in order not to create SMESH
1704 // faces for ng faces added here
1705 FillSMesh( occgeo, *ngMesh, initState, *_mesh, nodeVec, comment );
1706 // add ng faces to solids with internal vertices
1707 addIntVerticesInSolids( occgeo, *ngMesh, nodeVec, internals );
1708 // duplicate mesh faces on internal faces
1709 fixIntFaces( occgeo, *ngMesh, internals );
1710 initState = NETGENPlugin_ngMeshInfo(ngMesh);
1712 // Let netgen compute 3D mesh
1713 startWith = netgen::MESHCONST_MESHVOLUME;
1714 endWith = _optimize ? netgen::MESHCONST_OPTVOLUME : netgen::MESHCONST_MESHVOLUME;
1715 err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
1716 if (err) comment << "Error in netgen::OCCGenerateMesh()";
1718 if (!err && mparams.secondorder > 0)
1720 netgen::OCCRefinementSurfaces ref (occgeo);
1721 ref.MakeSecondOrder (*ngMesh);
1724 catch (netgen::NgException exc)
1726 error->myName = err = COMPERR_ALGO_FAILED;
1727 comment << exc.What();
1730 int nbNod = ngMesh->GetNP();
1731 int nbSeg = ngMesh->GetNSeg();
1732 int nbFac = ngMesh->GetNSE();
1733 int nbVol = ngMesh->GetNE();
1734 bool isOK = ( !err && (_isVolume ? (nbVol > 0) : (nbFac > 0)) );
1736 MESSAGE((err ? "Mesh Generation failure" : "End of Mesh Generation") <<
1737 ", nb nodes: " << nbNod <<
1738 ", nb segments: " << nbSeg <<
1739 ", nb faces: " << nbFac <<
1740 ", nb volumes: " << nbVol);
1742 // ------------------------------------------------------------
1743 // Feed back the SMESHDS with the generated Nodes and Elements
1744 // ------------------------------------------------------------
1746 if ( true /*isOK*/ ) // get whatever built
1747 FillSMesh( occgeo, *ngMesh, initState, *_mesh, nodeVec, comment ); //!<
1749 SMESH_ComputeErrorPtr readErr = readErrors(nodeVec);
1750 if ( readErr && !readErr->myBadElements.empty() )
1753 if ( error->IsOK() && ( !isOK || comment.size() > 0 ))
1754 error->myName = COMPERR_ALGO_FAILED;
1755 if ( !comment.empty() )
1756 error->myComment = comment;
1758 // set bad compute error to subshapes of all failed subshapes shapes
1759 if ( !error->IsOK() && err )
1762 for (int i = 1; i <= occgeo.fmap.Extent(); i++) {
1763 int status = occgeo.facemeshstatus[i-1];
1764 if (status == 1 ) continue;
1766 if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( occgeo.fmap( i ))) {
1767 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1768 if ( !smError || smError->IsOK() ) {
1770 smError.reset( new SMESH_ComputeError( *error ));
1772 smError.reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED, "Ignored" ));
1776 if ( !pb2D ) // all faces are OK
1777 for (int i = 1; i <= occgeo.somap.Extent(); i++)
1778 if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( occgeo.somap( i )))
1780 bool smComputed = !sm->IsEmpty();
1781 if ( smComputed && internals.hasInternalVertexInSolid( sm->GetId() ))
1783 int nbIntV = internals.getSolidsWithVertices().find( sm->GetId() )->second.size();
1784 SMESHDS_SubMesh* smDS = sm->GetSubMeshDS();
1785 smComputed = ( smDS->NbElements() > 0 || smDS->NbNodes() > nbIntV );
1787 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
1788 if ( !smComputed && ( !smError || smError->IsOK() ))
1789 smError.reset( new SMESH_ComputeError( *error ));
1793 return error->IsOK();
1796 //=============================================================================
1800 //=============================================================================
1801 bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
1803 netgen::MeshingParameters& mparams = netgen::mparam;
1806 // -------------------------
1807 // Prepare OCC geometry
1808 // -------------------------
1809 netgen::OCCGeometry occgeo;
1810 PrepareOCCgeometry( occgeo, _shape, *_mesh );
1812 bool tooManyElems = false;
1813 const int hugeNb = std::numeric_limits<int>::max() / 100;
1818 // pass 1D simple parameters to NETGEN
1820 if ( int nbSeg = _simpleHyp->GetNumberOfSegments() ) {
1822 mparams.segmentsperedge = nbSeg + 0.1;
1823 mparams.maxh = occgeo.boundingbox.Diam();
1824 mparams.grading = 0.01;
1828 mparams.segmentsperedge = 1;
1829 mparams.maxh = _simpleHyp->GetLocalLength();
1832 // let netgen create ngMesh and calculate element size on not meshed shapes
1833 NETGENPlugin_NetgenLibWrapper ngLib;
1834 netgen::Mesh *ngMesh = NULL;
1836 int startWith = netgen::MESHCONST_ANALYSE;
1837 int endWith = netgen::MESHCONST_MESHEDGES;
1838 int err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
1839 ngLib.setMesh(( Ng_Mesh*) ngMesh );
1841 if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( _shape ))
1842 sm->GetComputeError().reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED ));
1846 // calculate total nb of segments and length of edges
1847 double fullLen = 0.0;
1849 int entity = mparams.secondorder > 0 ? SMDSEntity_Quad_Edge : SMDSEntity_Edge;
1850 TopTools_DataMapOfShapeInteger Edge2NbSeg;
1851 for (TopExp_Explorer exp(_shape, TopAbs_EDGE); exp.More(); exp.Next())
1853 TopoDS_Edge E = TopoDS::Edge( exp.Current() );
1854 if( !Edge2NbSeg.Bind(E,0) )
1857 double aLen = SMESH_Algo::EdgeLength(E);
1860 vector<int>& aVec = aResMap[_mesh->GetSubMesh(E)];
1862 aVec.resize( SMDSEntity_Last, 0);
1864 fullNbSeg += aVec[ entity ];
1867 // store nb of segments computed by Netgen
1868 NCollection_Map<Link> linkMap;
1869 for (int i = 1; i <= ngMesh->GetNSeg(); ++i )
1871 const netgen::Segment& seg = ngMesh->LineSegment(i);
1872 Link link(seg[0], seg[1]);
1873 if ( !linkMap.Add( link )) continue;
1874 int aGeomEdgeInd = seg.epgeominfo[0].edgenr;
1875 if (aGeomEdgeInd > 0 && aGeomEdgeInd <= occgeo.emap.Extent())
1877 vector<int>& aVec = aResMap[_mesh->GetSubMesh(occgeo.emap(aGeomEdgeInd))];
1881 // store nb of nodes on edges computed by Netgen
1882 TopTools_DataMapIteratorOfDataMapOfShapeInteger Edge2NbSegIt(Edge2NbSeg);
1883 for (; Edge2NbSegIt.More(); Edge2NbSegIt.Next())
1885 vector<int>& aVec = aResMap[_mesh->GetSubMesh(Edge2NbSegIt.Key())];
1886 if ( aVec[ entity ] > 1 && aVec[ SMDSEntity_Node ] == 0 )
1887 aVec[SMDSEntity_Node] = mparams.secondorder > 0 ? 2*aVec[ entity ]-1 : aVec[ entity ]-1;
1889 fullNbSeg += aVec[ entity ];
1890 Edge2NbSeg( Edge2NbSegIt.Key() ) = aVec[ entity ];
1897 if ( double area = _simpleHyp->GetMaxElementArea() ) {
1899 mparams.maxh = sqrt(2. * area/sqrt(3.0));
1900 mparams.grading = 0.4; // moderate size growth
1903 // length from edges
1904 mparams.maxh = fullLen/fullNbSeg;
1905 mparams.grading = 0.2; // slow size growth
1908 mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
1909 mparams.maxh = min( mparams.maxh, fullLen/fullNbSeg * (1. + mparams.grading));
1911 for (TopExp_Explorer exp(_shape, TopAbs_FACE); exp.More(); exp.Next())
1913 TopoDS_Face F = TopoDS::Face( exp.Current() );
1914 SMESH_subMesh *sm = _mesh->GetSubMesh(F);
1916 BRepGProp::SurfaceProperties(F,G);
1917 double anArea = G.Mass();
1918 tooManyElems = tooManyElems || ( anArea/hugeNb > mparams.maxh*mparams.maxh );
1920 if ( !tooManyElems )
1922 TopTools_MapOfShape egdes;
1923 for (TopExp_Explorer exp1(F,TopAbs_EDGE); exp1.More(); exp1.Next())
1924 if ( egdes.Add( exp1.Current() ))
1925 nb1d += Edge2NbSeg.Find(exp1.Current());
1927 int nbFaces = tooManyElems ? hugeNb : int( 4*anArea / (mparams.maxh*mparams.maxh*sqrt(3.)));
1928 int nbNodes = tooManyElems ? hugeNb : (( nbFaces*3 - (nb1d-1)*2 ) / 6 + 1 );
1930 vector<int> aVec(SMDSEntity_Last, 0);
1931 if( mparams.secondorder > 0 ) {
1932 int nb1d_in = (nbFaces*3 - nb1d) / 2;
1933 aVec[SMDSEntity_Node] = nbNodes + nb1d_in;
1934 aVec[SMDSEntity_Quad_Triangle] = nbFaces;
1937 aVec[SMDSEntity_Node] = nbNodes;
1938 aVec[SMDSEntity_Triangle] = nbFaces;
1940 aResMap[sm].swap(aVec);
1947 // pass 3D simple parameters to NETGEN
1948 const NETGENPlugin_SimpleHypothesis_3D* simple3d =
1949 dynamic_cast< const NETGENPlugin_SimpleHypothesis_3D* > ( _simpleHyp );
1951 if ( double vol = simple3d->GetMaxElementVolume() ) {
1953 mparams.maxh = pow( 72, 1/6. ) * pow( vol, 1/3. );
1954 mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
1957 // using previous length from faces
1959 mparams.grading = 0.4;
1960 mparams.maxh = min( mparams.maxh, fullLen/fullNbSeg * (1. + mparams.grading));
1963 BRepGProp::VolumeProperties(_shape,G);
1964 double aVolume = G.Mass();
1965 double tetrVol = 0.1179*mparams.maxh*mparams.maxh*mparams.maxh;
1966 tooManyElems = tooManyElems || ( aVolume/hugeNb > tetrVol );
1967 int nbVols = tooManyElems ? hugeNb : int(aVolume/tetrVol);
1968 int nb1d_in = int(( nbVols*6 - fullNbSeg ) / 6 );
1969 vector<int> aVec(SMDSEntity_Last, 0 );
1970 if ( tooManyElems ) // avoid FPE
1972 aVec[SMDSEntity_Node] = hugeNb;
1973 aVec[ mparams.secondorder > 0 ? SMDSEntity_Quad_Tetra : SMDSEntity_Tetra] = hugeNb;
1977 if( mparams.secondorder > 0 ) {
1978 aVec[SMDSEntity_Node] = nb1d_in/3 + 1 + nb1d_in;
1979 aVec[SMDSEntity_Quad_Tetra] = nbVols;
1982 aVec[SMDSEntity_Node] = nb1d_in/3 + 1;
1983 aVec[SMDSEntity_Tetra] = nbVols;
1986 SMESH_subMesh *sm = _mesh->GetSubMesh(_shape);
1987 aResMap[sm].swap(aVec);
1993 //================================================================================
1995 * \brief Remove "test.out" and "problemfaces" files in current directory
1997 //================================================================================
1999 void NETGENPlugin_Mesher::RemoveTmpFiles()
2001 SMESH_File("test.out").remove();
2002 SMESH_File("problemfaces").remove();
2003 SMESH_File("occmesh.rep").remove();
2006 //================================================================================
2008 * \brief Read mesh entities preventing successful computation from "test.out" file
2010 //================================================================================
2012 SMESH_ComputeErrorPtr
2013 NETGENPlugin_Mesher::readErrors(const vector<const SMDS_MeshNode* >& nodeVec)
2015 SMESH_ComputeErrorPtr err = SMESH_ComputeError::New
2016 (COMPERR_BAD_INPUT_MESH, "Some edges multiple times in surface mesh");
2017 SMESH_File file("test.out");
2019 const char* badEdgeStr = " multiple times in surface mesh";
2020 const int badEdgeStrLen = strlen( badEdgeStr );
2021 while( !file.eof() )
2023 if ( strncmp( file, "Edge ", 5 ) == 0 &&
2024 file.getInts( two ) &&
2025 strncmp( file, badEdgeStr, badEdgeStrLen ) == 0 &&
2026 two[0] < nodeVec.size() && two[1] < nodeVec.size())
2028 err->myBadElements.push_back( new SMDS_MeshEdge( nodeVec[ two[0]], nodeVec[ two[1]] ));
2029 file += badEdgeStrLen;
2031 else if ( strncmp( file, "Intersecting: ", 14 ) == 0 )
2034 // openelement 18 with open element 126
2037 vector<int> three1(3), three2(3);
2039 const char* pos = file;
2040 bool ok = ( strncmp( file, "openelement ", 12 ) == 0 );
2041 ok = ok && file.getInts( two );
2042 ok = ok && file.getInts( three1 );
2043 ok = ok && file.getInts( three2 );
2044 for ( int i = 0; ok && i < 3; ++i )
2045 ok = ( three1[i] < nodeVec.size() && nodeVec[ three1[i]]);
2046 for ( int i = 0; ok && i < 3; ++i )
2047 ok = ( three2[i] < nodeVec.size() && nodeVec[ three2[i]]);
2050 err->myBadElements.push_back( new SMDS_FaceOfNodes( nodeVec[ three1[0]],
2051 nodeVec[ three1[1]],
2052 nodeVec[ three1[2]]));
2053 err->myBadElements.push_back( new SMDS_FaceOfNodes( nodeVec[ three2[0]],
2054 nodeVec[ three2[1]],
2055 nodeVec[ three2[2]]));
2056 err->myComment = "Intersecting triangles";
2071 //================================================================================
2073 * \brief Constructor of NETGENPlugin_ngMeshInfo
2075 //================================================================================
2077 NETGENPlugin_ngMeshInfo::NETGENPlugin_ngMeshInfo( netgen::Mesh* ngMesh)
2081 _nbNodes = ngMesh->GetNP();
2082 _nbSegments = ngMesh->GetNSeg();
2083 _nbFaces = ngMesh->GetNSE();
2084 _nbVolumes = ngMesh->GetNE();
2088 _nbNodes = _nbSegments = _nbFaces = _nbVolumes = 0;
2092 //================================================================================
2094 * \brief Find "internal" sub-shapes
2096 //================================================================================
2098 NETGENPlugin_Internals::NETGENPlugin_Internals( SMESH_Mesh& mesh,
2099 const TopoDS_Shape& shape,
2101 : _mesh( mesh ), _is3D( is3D )
2103 SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
2105 TopExp_Explorer f,e;
2106 for ( f.Init( shape, TopAbs_FACE ); f.More(); f.Next() )
2108 int faceID = meshDS->ShapeToIndex( f.Current() );
2110 // find not computed internal edges
2112 for ( e.Init( f.Current().Oriented(TopAbs_FORWARD), TopAbs_EDGE ); e.More(); e.Next() )
2113 if ( e.Current().Orientation() == TopAbs_INTERNAL )
2115 SMESH_subMesh* eSM = mesh.GetSubMesh( e.Current() );
2116 if ( eSM->IsEmpty() )
2118 _e2face.insert( make_pair( eSM->GetId(), faceID ));
2119 for ( TopoDS_Iterator v(e.Current()); v.More(); v.Next() )
2120 _e2face.insert( make_pair( meshDS->ShapeToIndex( v.Value() ), faceID ));
2124 // find internal vertices in a face
2125 set<int> intVV; // issue 0020850 where same vertex is twice in a face
2126 for ( TopoDS_Iterator fSub( f.Current() ); fSub.More(); fSub.Next())
2127 if ( fSub.Value().ShapeType() == TopAbs_VERTEX )
2129 int vID = meshDS->ShapeToIndex( fSub.Value() );
2130 if ( intVV.insert( vID ).second )
2131 _f2v[ faceID ].push_back( vID );
2136 // find internal faces and their subshapes where nodes are to be doubled
2137 // to make a crack with non-sewed borders
2139 if ( f.Current().Orientation() == TopAbs_INTERNAL )
2141 _intShapes.insert( meshDS->ShapeToIndex( f.Current() ));
2144 list< TopoDS_Shape > edges;
2145 for ( e.Init( f.Current(), TopAbs_EDGE ); e.More(); e.Next())
2146 if ( SMESH_MesherHelper::NbAncestors( e.Current(), mesh, TopAbs_FACE ) > 1 )
2148 _intShapes.insert( meshDS->ShapeToIndex( e.Current() ));
2149 edges.push_back( e.Current() );
2150 // find border faces
2151 PShapeIteratorPtr fIt =
2152 SMESH_MesherHelper::GetAncestors( edges.back(),mesh,TopAbs_FACE );
2153 while ( const TopoDS_Shape* pFace = fIt->next() )
2154 if ( !pFace->IsSame( f.Current() ))
2155 _borderFaces.insert( meshDS->ShapeToIndex( *pFace ));
2158 // we consider vertex internal if it is shared by more than one internal edge
2159 list< TopoDS_Shape >::iterator edge = edges.begin();
2160 for ( ; edge != edges.end(); ++edge )
2161 for ( TopoDS_Iterator v( *edge ); v.More(); v.Next() )
2163 set<int> internalEdges;
2164 PShapeIteratorPtr eIt =
2165 SMESH_MesherHelper::GetAncestors( v.Value(),mesh,TopAbs_EDGE );
2166 while ( const TopoDS_Shape* pEdge = eIt->next() )
2168 int edgeID = meshDS->ShapeToIndex( *pEdge );
2169 if ( isInternalShape( edgeID ))
2170 internalEdges.insert( edgeID );
2172 if ( internalEdges.size() > 1 )
2173 _intShapes.insert( meshDS->ShapeToIndex( v.Value() ));
2177 } // loop on geom faces
2179 // find vertices internal in solids
2182 for ( TopExp_Explorer so(shape, TopAbs_SOLID); so.More(); so.Next())
2184 int soID = meshDS->ShapeToIndex( so.Current() );
2185 for ( TopoDS_Iterator soSub( so.Current() ); soSub.More(); soSub.Next())
2186 if ( soSub.Value().ShapeType() == TopAbs_VERTEX )
2187 _s2v[ soID ].push_back( meshDS->ShapeToIndex( soSub.Value() ));
2192 //================================================================================
2194 * \brief Find mesh faces on non-internal geom faces sharing internal edge
2195 * some nodes of which are to be doubled to make the second border of the "crack"
2197 //================================================================================
2199 void NETGENPlugin_Internals::findBorderElements( TIDSortedElemSet & borderElems )
2201 if ( _intShapes.empty() ) return;
2203 SMESH_Mesh& mesh = const_cast<SMESH_Mesh&>(_mesh);
2204 SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
2206 // loop on internal geom edges
2207 set<int>::const_iterator intShapeId = _intShapes.begin();
2208 for ( ; intShapeId != _intShapes.end(); ++intShapeId )
2210 const TopoDS_Shape& s = meshDS->IndexToShape( *intShapeId );
2211 if ( s.ShapeType() != TopAbs_EDGE ) continue;
2213 // get internal and non-internal geom faces sharing the internal edge <s>
2215 set<int>::iterator bordFace = _borderFaces.end();
2216 PShapeIteratorPtr faces = SMESH_MesherHelper::GetAncestors( s, _mesh, TopAbs_FACE );
2217 while ( const TopoDS_Shape* pFace = faces->next() )
2219 int faceID = meshDS->ShapeToIndex( *pFace );
2220 if ( isInternalShape( faceID ))
2223 bordFace = _borderFaces.insert( faceID ).first;
2225 if ( bordFace == _borderFaces.end() || !intFace ) continue;
2227 // get all links of mesh faces on internal geom face sharing nodes on edge <s>
2228 set< SMESH_OrientedLink > links; //!< links of faces on internal geom face
2229 list<const SMDS_MeshElement*> suspectFaces[2]; //!< mesh faces on border geom faces
2230 int nbSuspectFaces = 0;
2231 SMESHDS_SubMesh* intFaceSM = meshDS->MeshElements( intFace );
2232 if ( !intFaceSM || intFaceSM->NbElements() == 0 ) continue;
2233 SMESH_subMeshIteratorPtr smIt = mesh.GetSubMesh( s )->getDependsOnIterator(true,true);
2234 while ( smIt->more() )
2236 SMESHDS_SubMesh* sm = smIt->next()->GetSubMeshDS();
2237 if ( !sm ) continue;
2238 SMDS_NodeIteratorPtr nIt = sm->GetNodes();
2239 while ( nIt->more() )
2241 const SMDS_MeshNode* nOnEdge = nIt->next();
2242 SMDS_ElemIteratorPtr fIt = nOnEdge->GetInverseElementIterator(SMDSAbs_Face);
2243 while ( fIt->more() )
2245 const SMDS_MeshElement* f = fIt->next();
2246 int nbNodes = f->NbNodes() / ( f->IsQuadratic() ? 2 : 1 );
2247 if ( intFaceSM->Contains( f ))
2249 for ( int i = 0; i < nbNodes; ++i )
2250 links.insert( SMESH_OrientedLink( f->GetNode(i), f->GetNode((i+1)%nbNodes)));
2255 for ( int i = 0; i < nbNodes; ++i )
2256 nbDblNodes += isInternalShape( f->GetNode(i)->GetPosition()->GetShapeId() );
2258 suspectFaces[ nbDblNodes < 2 ].push_back( f );
2264 // suspectFaces[0] having link with same orientation as mesh faces on
2265 // the internal geom face are <borderElems>. suspectFaces[1] have
2266 // only one node on edge <s>, we decide on them later (at the 2nd loop)
2267 // by links of <borderElems> found at the 1st and 2nd loops
2268 set< SMESH_OrientedLink > borderLinks;
2269 for ( int isPostponed = 0; isPostponed < 2; ++isPostponed )
2271 list<const SMDS_MeshElement*>::iterator fIt = suspectFaces[isPostponed].begin();
2272 for ( int nbF = 0; fIt != suspectFaces[isPostponed].end(); ++fIt, ++nbF )
2274 const SMDS_MeshElement* f = *fIt;
2275 bool isBorder = false, linkFound = false, borderLinkFound = false;
2276 list< SMESH_OrientedLink > faceLinks;
2277 int nbNodes = f->NbNodes() / ( f->IsQuadratic() ? 2 : 1 );
2278 for ( int i = 0; i < nbNodes; ++i )
2280 SMESH_OrientedLink link( f->GetNode(i), f->GetNode((i+1)%nbNodes));
2281 faceLinks.push_back( link );
2284 set< SMESH_OrientedLink >::iterator foundLink = links.find( link );
2285 if ( foundLink != links.end() )
2288 isBorder = ( foundLink->_reversed == link._reversed );
2289 if ( !isBorder && !isPostponed ) break;
2290 faceLinks.pop_back();
2292 else if ( isPostponed && !borderLinkFound )
2294 foundLink = borderLinks.find( link );
2295 if ( foundLink != borderLinks.end() )
2297 borderLinkFound = true;
2298 isBorder = ( foundLink->_reversed != link._reversed );
2305 borderElems.insert( f );
2306 borderLinks.insert( faceLinks.begin(), faceLinks.end() );
2308 else if ( !linkFound && !borderLinkFound )
2310 suspectFaces[1].push_back( f );
2311 if ( nbF > 2 * nbSuspectFaces )
2312 break; // dead loop protection
2319 //================================================================================
2321 * \brief put internal shapes in maps and fill in submeshes to precompute
2323 //================================================================================
2325 void NETGENPlugin_Internals::getInternalEdges( TopTools_IndexedMapOfShape& fmap,
2326 TopTools_IndexedMapOfShape& emap,
2327 TopTools_IndexedMapOfShape& vmap,
2328 list< SMESH_subMesh* >& smToPrecompute)
2330 if ( !hasInternalEdges() ) return;
2331 map<int,int>::const_iterator ev_face = _e2face.begin();
2332 for ( ; ev_face != _e2face.end(); ++ev_face )
2334 const TopoDS_Shape& ev = _mesh.GetMeshDS()->IndexToShape( ev_face->first );
2335 const TopoDS_Shape& face = _mesh.GetMeshDS()->IndexToShape( ev_face->second );
2337 ( ev.ShapeType() == TopAbs_EDGE ? emap : vmap ).Add( ev );
2339 //cout<<"INTERNAL EDGE or VERTEX "<<ev_face->first<<" on face "<<ev_face->second<<endl;
2341 smToPrecompute.push_back( _mesh.GetSubMeshContaining( ev_face->first ));
2345 //================================================================================
2347 * \brief return shapes and submeshes to be meshed and already meshed boundary submeshes
2349 //================================================================================
2351 void NETGENPlugin_Internals::getInternalFaces( TopTools_IndexedMapOfShape& fmap,
2352 TopTools_IndexedMapOfShape& emap,
2353 list< SMESH_subMesh* >& intFaceSM,
2354 list< SMESH_subMesh* >& boundarySM)
2356 if ( !hasInternalFaces() ) return;
2358 // <fmap> and <emap> are for not yet meshed shapes
2359 // <intFaceSM> is for submeshes of faces
2360 // <boundarySM> is for meshed edges and vertices
2365 set<int> shapeIDs ( _intShapes );
2366 if ( !_borderFaces.empty() )
2367 shapeIDs.insert( _borderFaces.begin(), _borderFaces.end() );
2369 set<int>::const_iterator intS = shapeIDs.begin();
2370 for ( ; intS != shapeIDs.end(); ++intS )
2372 SMESH_subMesh* sm = _mesh.GetSubMeshContaining( *intS );
2374 if ( sm->GetSubShape().ShapeType() != TopAbs_FACE ) continue;
2376 intFaceSM.push_back( sm );
2378 // add submeshes of not computed internal faces
2379 if ( !sm->IsEmpty() ) continue;
2381 SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(true,true);
2382 while ( smIt->more() )
2385 const TopoDS_Shape& s = sm->GetSubShape();
2387 if ( sm->IsEmpty() )
2390 switch ( s.ShapeType() ) {
2391 case TopAbs_FACE: fmap.Add ( s ); break;
2392 case TopAbs_EDGE: emap.Add ( s ); break;
2398 if ( s.ShapeType() != TopAbs_FACE )
2399 boundarySM.push_back( sm );
2405 //================================================================================
2407 * \brief Return true if given shape is to be precomputed in order to be correctly
2408 * added to netgen mesh
2410 //================================================================================
2412 bool NETGENPlugin_Internals::isShapeToPrecompute(const TopoDS_Shape& s)
2414 int shapeID = _mesh.GetMeshDS()->ShapeToIndex( s );
2415 switch ( s.ShapeType() ) {
2416 case TopAbs_FACE : break; //return isInternalShape( shapeID ) || isBorderFace( shapeID );
2417 case TopAbs_EDGE : return isInternalEdge( shapeID );
2418 case TopAbs_VERTEX: break;
2424 //================================================================================
2426 * \brief Return SMESH
2428 //================================================================================
2430 SMESH_Mesh& NETGENPlugin_Internals::getMesh() const
2432 return const_cast<SMESH_Mesh&>( _mesh );
2435 //================================================================================
2437 * \brief Initialize netgen library
2439 //================================================================================
2441 NETGENPlugin_NetgenLibWrapper::NETGENPlugin_NetgenLibWrapper()
2444 _ngMesh = Ng_NewMesh();
2447 //================================================================================
2449 * \brief Finish using netgen library
2451 //================================================================================
2453 NETGENPlugin_NetgenLibWrapper::~NETGENPlugin_NetgenLibWrapper()
2455 Ng_DeleteMesh( _ngMesh );
2457 NETGENPlugin_Mesher::RemoveTmpFiles();
2460 //================================================================================
2462 * \brief Set netgen mesh to delete at destruction
2464 //================================================================================
2466 void NETGENPlugin_NetgenLibWrapper::setMesh( Ng_Mesh* mesh )
2469 Ng_DeleteMesh( _ngMesh );