1 // Copyright (C) 2007-2008 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
22 // NETGENPlugin : C++ implementation
23 // File : NETGENPlugin_Mesher.cxx
24 // Author : Michael Sazonov (OCN)
27 //=============================================================================
29 #include "NETGENPlugin_Mesher.hxx"
30 #include "NETGENPlugin_Hypothesis_2D.hxx"
31 #include "NETGENPlugin_SimpleHypothesis_3D.hxx"
33 #include <SMESH_Mesh.hxx>
34 #include <SMESH_Comment.hxx>
35 #include <SMESH_ComputeError.hxx>
36 #include <SMESH_subMesh.hxx>
37 #include <SMESH_MesherHelper.hxx>
38 #include <SMESHDS_Mesh.hxx>
39 #include <SMDS_MeshElement.hxx>
40 #include <SMDS_MeshNode.hxx>
41 #include <utilities.h>
46 #include <BRep_Tool.hxx>
47 #include <NCollection_Map.hxx>
48 #include <OSD_File.hxx>
49 #include <OSD_Path.hxx>
50 #include <Standard_ErrorHandler.hxx>
51 #include <Standard_ProgramError.hxx>
52 #include <TCollection_AsciiString.hxx>
54 #include <TopExp_Explorer.hxx>
55 #include <TopTools_DataMapOfShapeInteger.hxx>
56 #include <TopTools_DataMapIteratorOfDataMapOfShapeInteger.hxx>
57 #include <TopTools_ListIteratorOfListOfShape.hxx>
58 #include <TopTools_MapOfShape.hxx>
61 // Netgen include files
66 #include <occgeom.hpp>
67 #include <meshing.hpp>
68 //#include <ngexception.hpp>
70 extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, int, int, char*);
71 extern MeshingParameters mparam;
76 static void removeFile( const TCollection_AsciiString& fileName )
79 OSD_File( fileName ).Remove();
81 catch ( Standard_ProgramError ) {
82 MESSAGE("Can't remove file: " << fileName.ToCString() << " ; file does not exist or permission denied");
86 //=============================================================================
90 //=============================================================================
92 NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh* mesh,
93 const TopoDS_Shape& aShape,
104 //================================================================================
106 * \brief Initialize global NETGEN parameters with default values
108 //================================================================================
110 void NETGENPlugin_Mesher::defaultParameters()
113 netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters();
115 netgen::MeshingParameters& mparams = netgen::mparam;
117 // maximal mesh edge size
118 mparams.maxh = NETGENPlugin_Hypothesis::GetDefaultMaxSize();
119 // minimal number of segments per edge
120 mparams.segmentsperedge = NETGENPlugin_Hypothesis::GetDefaultNbSegPerEdge();
121 // rate of growth of size between elements
122 mparams.grading = NETGENPlugin_Hypothesis::GetDefaultGrowthRate();
123 // safety factor for curvatures (elements per radius)
124 mparams.curvaturesafety = NETGENPlugin_Hypothesis::GetDefaultNbSegPerRadius();
125 // create elements of second order
126 mparams.secondorder = NETGENPlugin_Hypothesis::GetDefaultSecondOrder() ? 1 : 0;
127 // quad-dominated surface meshing
131 mparams.quad = NETGENPlugin_Hypothesis_2D::GetDefaultQuadAllowed() ? 1 : 0;
134 //=============================================================================
136 * Pass parameters to NETGEN
138 //=============================================================================
139 void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp)
144 netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters();
146 netgen::MeshingParameters& mparams = netgen::mparam;
148 // Initialize global NETGEN parameters:
149 // maximal mesh segment size
150 mparams.maxh = hyp->GetMaxSize();
151 // minimal number of segments per edge
152 mparams.segmentsperedge = hyp->GetNbSegPerEdge();
153 // rate of growth of size between elements
154 mparams.grading = hyp->GetGrowthRate();
155 // safety factor for curvatures (elements per radius)
156 mparams.curvaturesafety = hyp->GetNbSegPerRadius();
157 // create elements of second order
158 mparams.secondorder = hyp->GetSecondOrder() ? 1 : 0;
159 // quad-dominated surface meshing
160 // only triangles are allowed for volumic mesh
162 mparams.quad = static_cast<const NETGENPlugin_Hypothesis_2D*>
163 (hyp)->GetQuadAllowed() ? 1 : 0;
164 _optimize = hyp->GetOptimize();
169 //=============================================================================
171 * Pass simple parameters to NETGEN
173 //=============================================================================
175 void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_SimpleHypothesis_2D* hyp)
182 //=============================================================================
184 * Link - a pair of integer numbers
186 //=============================================================================
190 Link(int _n1, int _n2) : n1(_n1), n2(_n2) {}
191 Link() : n1(0), n2(0) {}
194 int HashCode(const Link& aLink, int aLimit)
196 return HashCode(aLink.n1 + aLink.n2, aLimit);
199 Standard_Boolean IsEqual(const Link& aLink1, const Link& aLink2)
201 return (aLink1.n1 == aLink2.n1 && aLink1.n2 == aLink2.n2 ||
202 aLink1.n1 == aLink2.n2 && aLink1.n2 == aLink2.n1);
205 //================================================================================
207 * \brief Initialize netgen::OCCGeometry with OCCT shape
209 //================================================================================
211 void NETGENPlugin_Mesher::PrepareOCCgeometry(netgen::OCCGeometry& occgeo,
212 const TopoDS_Shape& shape,
214 list< SMESH_subMesh* > * meshedSM)
216 BRepTools::Clean (shape);
218 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
221 BRepMesh_IncrementalMesh::BRepMesh_IncrementalMesh (shape, 0.01, true);
222 } catch (Standard_Failure) {
225 BRepBndLib::Add (shape, bb);
226 double x1,y1,z1,x2,y2,z2;
227 bb.Get (x1,y1,z1,x2,y2,z2);
228 MESSAGE("shape bounding box:\n" <<
229 "(" << x1 << " " << y1 << " " << z1 << ") " <<
230 "(" << x2 << " " << y2 << " " << z2 << ")");
231 netgen::Point<3> p1 = netgen::Point<3> (x1,y1,z1);
232 netgen::Point<3> p2 = netgen::Point<3> (x2,y2,z2);
233 occgeo.boundingbox = netgen::Box<3> (p1,p2);
235 occgeo.shape = shape;
237 //occgeo.BuildFMap();
239 // fill maps of shapes of occgeo with not yet meshed subshapes
241 // get root submeshes
242 list< SMESH_subMesh* > rootSM;
243 if ( SMESH_subMesh* sm = mesh.GetSubMeshContaining( shape )) {
244 rootSM.push_back( sm );
247 for ( TopoDS_Iterator it( shape ); it.More(); it.Next() )
248 rootSM.push_back( mesh.GetSubMesh( it.Value() ));
251 // add subshapes of empty submeshes
252 list< SMESH_subMesh* >::iterator rootIt = rootSM.begin(), rootEnd = rootSM.end();
253 for ( ; rootIt != rootEnd; ++rootIt ) {
254 SMESH_subMesh * root = *rootIt;
255 SMESH_subMeshIteratorPtr smIt = root->getDependsOnIterator(/*includeSelf=*/true,
256 /*complexShapeFirst=*/true);
257 // to find a right orientation of subshapes (PAL20462)
258 TopTools_IndexedMapOfShape subShapes;
259 TopExp::MapShapes(root->GetSubShape(), subShapes);
260 while ( smIt->more() ) {
261 SMESH_subMesh* sm = smIt->next();
262 if ( !meshedSM || sm->IsEmpty() ) {
263 TopoDS_Shape shape = sm->GetSubShape();
264 if ( shape.ShapeType() != TopAbs_VERTEX )
265 shape = subShapes( subShapes.FindIndex( shape ));// - shape->index->oriented shape
266 switch ( shape.ShapeType() ) {
267 case TopAbs_FACE : occgeo.fmap.Add( shape ); break;
268 case TopAbs_EDGE : occgeo.emap.Add( shape ); break;
269 case TopAbs_VERTEX: occgeo.vmap.Add( shape ); break;
270 case TopAbs_SOLID :occgeo.somap.Add( shape ); break;
274 // collect submeshes of meshed shapes
276 meshedSM->push_back( sm );
280 occgeo.facemeshstatus.SetSize (occgeo.fmap.Extent());
281 occgeo.facemeshstatus = 0;
283 occgeo.face_maxh.SetSize(occgeo.fmap.Extent());
284 occgeo.face_maxh = netgen::mparam.maxh;
289 //================================================================================
291 * \brief return id of netgen point corresponding to SMDS node
293 //================================================================================
294 typedef map< const SMDS_MeshNode*, int > TNode2IdMap;
296 static int ngNodeId( const SMDS_MeshNode* node,
297 netgen::Mesh& ngMesh,
298 TNode2IdMap& nodeNgIdMap)
300 int newNgId = ngMesh.GetNP() + 1;
302 pair< TNode2IdMap::iterator, bool > it_isNew = nodeNgIdMap.insert( make_pair( node, newNgId ));
304 if ( it_isNew.second ) {
305 netgen::MeshPoint p( netgen::Point<3> (node->X(), node->Y(), node->Z()) );
306 ngMesh.AddPoint( p );
308 return it_isNew.first->second;
311 //================================================================================
313 * \brief fill ngMesh with nodes and elements of computed submeshes
315 //================================================================================
317 bool NETGENPlugin_Mesher::fillNgMesh(netgen::OCCGeometry& occgeom,
318 netgen::Mesh& ngMesh,
319 vector<SMDS_MeshNode*>& nodeVec,
320 const list< SMESH_subMesh* > & meshedSM)
322 TNode2IdMap nodeNgIdMap;
324 TopTools_MapOfShape visitedShapes;
326 SMESH_MesherHelper helper (*_mesh);
328 int faceID = occgeom.fmap.Extent();
329 _faceDescriptors.clear();
331 list< SMESH_subMesh* >::const_iterator smIt, smEnd = meshedSM.end();
332 for ( smIt = meshedSM.begin(); smIt != smEnd; ++smIt )
334 SMESH_subMesh* sm = *smIt;
335 if ( !visitedShapes.Add( sm->GetSubShape() ))
338 SMESHDS_SubMesh * smDS = sm->GetSubMeshDS();
340 switch ( sm->GetSubShape().ShapeType() )
342 case TopAbs_EDGE: { // EDGE
343 // ----------------------
344 const TopoDS_Edge& geomEdge = TopoDS::Edge( sm->GetSubShape() );
346 // Add ng segments for each not meshed face the edge bounds
347 TopTools_MapOfShape visitedAncestors;
348 const TopTools_ListOfShape& ancestors = _mesh->GetAncestors( geomEdge );
349 TopTools_ListIteratorOfListOfShape ancestorIt ( ancestors );
350 for ( ; ancestorIt.More(); ancestorIt.Next() )
352 const TopoDS_Shape & ans = ancestorIt.Value();
353 if ( ans.ShapeType() != TopAbs_FACE || !visitedAncestors.Add( ans ))
355 const TopoDS_Face& face = TopoDS::Face( ans );
357 int faceID = occgeom.fmap.FindIndex( face );
359 continue; // meshed face
361 // find out orientation of geomEdge within face
362 bool isForwad = false;
363 for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next() ) {
364 if ( geomEdge.IsSame( exp.Current() )) {
365 isForwad = ( exp.Current().Orientation() == geomEdge.Orientation() );
369 bool isQuad = smDS->GetElements()->next()->IsQuadratic();
371 // get all nodes from geomEdge
372 StdMeshers_FaceSide fSide( face, geomEdge, _mesh, isForwad, isQuad );
373 const vector<UVPtStruct>& points = fSide.GetUVPtStruct();
374 int i, nbSeg = fSide.NbSegments();
376 double otherSeamParam = 0;
377 helper.SetSubShape( face );
378 bool isSeam = helper.IsRealSeam( geomEdge );
381 helper.GetOtherParam( helper.GetPeriodicIndex() == 1 ? points[0].u : points[0].v );
385 int prevNgId = ngNodeId( points[0].node, ngMesh, nodeNgIdMap );
387 for ( i = 0; i < nbSeg; ++i )
389 const UVPtStruct& p1 = points[ i ];
390 const UVPtStruct& p2 = points[ i+1 ];
395 seg.pnums[0] = prevNgId;
396 seg.pnums[1] = prevNgId = ngNodeId( p2.node, ngMesh, nodeNgIdMap );
399 seg.p2 = prevNgId = ngNodeId( p2.node, ngMesh, nodeNgIdMap );
401 // node param on curve
402 seg.epgeominfo[ 0 ].dist = p1.param;
403 seg.epgeominfo[ 1 ].dist = p2.param;
405 seg.epgeominfo[ 0 ].u = p1.u;
406 seg.epgeominfo[ 0 ].v = p1.v;
407 seg.epgeominfo[ 1 ].u = p2.u;
408 seg.epgeominfo[ 1 ].v = p2.v;
410 //seg.epgeominfo[ iEnd ].edgenr = edgeID; // = geom.emap.FindIndex(edge);
411 seg.si = faceID; // = geom.fmap.FindIndex (face);
412 seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
413 ngMesh.AddSegment (seg);
417 if ( helper.GetPeriodicIndex() == 1 ) {
418 seg.epgeominfo[ 0 ].u = otherSeamParam;
419 seg.epgeominfo[ 1 ].u = otherSeamParam;
420 swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v);
422 seg.epgeominfo[ 0 ].v = otherSeamParam;
423 seg.epgeominfo[ 1 ].v = otherSeamParam;
424 swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u);
427 swap (seg.pnums[0], seg.pnums[1]);
429 swap (seg.p1, seg.p2);
431 swap (seg.epgeominfo[0].dist, seg.epgeominfo[1].dist);
432 seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
433 ngMesh.AddSegment (seg);
436 } // loop on geomEdge ancestors
439 } // case TopAbs_EDGE
441 case TopAbs_FACE: { // FACE
442 // ----------------------
443 const TopoDS_Face& geomFace = TopoDS::Face( sm->GetSubShape() );
444 helper.SetSubShape( geomFace );
446 // Find solids the geomFace bounds
447 int solidID1 = 0, solidID2 = 0;
448 const TopTools_ListOfShape& ancestors = _mesh->GetAncestors( geomFace );
449 TopTools_ListIteratorOfListOfShape ancestorIt ( ancestors );
450 for ( ; ancestorIt.More(); ancestorIt.Next() )
452 const TopoDS_Shape & solid = ancestorIt.Value();
453 if ( solid.ShapeType() == TopAbs_SOLID ) {
454 int id = occgeom.somap.FindIndex ( solid );
455 if ( solidID1 && id != solidID1 ) solidID2 = id;
460 _faceDescriptors[ faceID ].first = solidID1;
461 _faceDescriptors[ faceID ].second = solidID2;
463 // Orient the face correctly in solidID1 (issue 0020206)
464 bool reverse = false;
466 TopoDS_Shape solid = occgeom.somap( solidID1 );
467 for ( TopExp_Explorer f( solid, TopAbs_FACE ); f.More(); f.Next() ) {
468 if ( geomFace.IsSame( f.Current() )) {
469 reverse = SMESH_Algo::IsReversedSubMesh( TopoDS::Face( f.Current()), helper.GetMeshDS() );
475 // Add surface elements
476 SMDS_ElemIteratorPtr faces = smDS->GetElements();
477 while ( faces->more() ) {
479 const SMDS_MeshElement* f = faces->next();
480 if ( f->NbNodes() % 3 != 0 ) { // not triangle
481 for ( ancestorIt.Initialize(ancestors); ancestorIt.More(); ancestorIt.Next() )
482 if ( ancestorIt.Value().ShapeType() == TopAbs_SOLID ) {
483 sm = _mesh->GetSubMesh( ancestorIt.Value() );
486 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
487 smError.reset( new SMESH_ComputeError(COMPERR_BAD_INPUT_MESH,"Not triangle submesh"));
488 smError->myBadElements.push_back( f );
492 netgen::Element2d tri(3);
493 tri.SetIndex ( faceID );
495 for ( int i = 0; i < 3; ++i ) {
496 const SMDS_MeshNode* node = f->GetNode( i ), * inFaceNode=0;
497 if ( helper.IsSeamShape( node->GetPosition()->GetShapeId() ))
498 if ( helper.IsSeamShape( f->GetNodeWrap( i+1 )->GetPosition()->GetShapeId() ))
499 inFaceNode = f->GetNodeWrap( i-1 );
501 inFaceNode = f->GetNodeWrap( i+1 );
503 gp_XY uv = helper.GetNodeUV( geomFace, node, inFaceNode );
505 tri.GeomInfoPi(3-i).u = uv.X();
506 tri.GeomInfoPi(3-i).v = uv.Y();
507 tri.PNum (3-i) = ngNodeId( node, ngMesh, nodeNgIdMap );
509 tri.GeomInfoPi(i+1).u = uv.X();
510 tri.GeomInfoPi(i+1).v = uv.Y();
511 tri.PNum (i+1) = ngNodeId( node, ngMesh, nodeNgIdMap );
515 ngMesh.AddSurfaceElement (tri);
521 case TopAbs_VERTEX: { // VERTEX
522 // --------------------------
523 SMDS_NodeIteratorPtr nodeIt = smDS->GetNodes();
524 if ( nodeIt->more() )
525 ngNodeId( nodeIt->next(), ngMesh, nodeNgIdMap );
530 } // loop on submeshes
533 nodeVec.resize( ngMesh.GetNP() + 1 );
534 TNode2IdMap::iterator node_NgId, nodeNgIdEnd = nodeNgIdMap.end();
535 for ( node_NgId = nodeNgIdMap.begin(); node_NgId != nodeNgIdEnd; ++node_NgId)
536 nodeVec[ node_NgId->second ] = (SMDS_MeshNode*) node_NgId->first;
541 //=============================================================================
543 * Here we are going to use the NETGEN mesher
545 //=============================================================================
546 bool NETGENPlugin_Mesher::Compute()
549 netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters();
551 netgen::MeshingParameters& mparams = netgen::mparam;
553 MESSAGE("Compute with:\n"
554 " max size = " << mparams.maxh << "\n"
555 " segments per edge = " << mparams.segmentsperedge);
557 " growth rate = " << mparams.grading << "\n"
558 " elements per radius = " << mparams.curvaturesafety << "\n"
559 " second order = " << mparams.secondorder << "\n"
560 " quad allowed = " << mparams.quad);
562 SMESH_ComputeErrorPtr error = SMESH_ComputeError::New();
565 // -------------------------
566 // Prepare OCC geometry
567 // -------------------------
569 netgen::OCCGeometry occgeo;
570 list< SMESH_subMesh* > meshedSM;
571 PrepareOCCgeometry( occgeo, _shape, *_mesh, &meshedSM );
573 // -------------------------
575 // -------------------------
577 netgen::Mesh *ngMesh = NULL;
579 SMESH_Comment comment;
584 // vector of nodes in which node index == netgen ID
585 vector< SMDS_MeshNode* > nodeVec;
591 // pass 1D simple parameters to NETGEN
593 if ( int nbSeg = _simpleHyp->GetNumberOfSegments() ) {
595 mparams.segmentsperedge = nbSeg + 0.1;
596 mparams.maxh = occgeo.boundingbox.Diam();
597 mparams.grading = 0.01;
601 mparams.segmentsperedge = 1;
602 mparams.maxh = _simpleHyp->GetLocalLength();
605 // let netgen create ngMesh and calculate element size on not meshed shapes
607 int startWith = netgen::MESHCONST_ANALYSE;
608 int endWith = netgen::MESHCONST_ANALYSE;
609 err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
610 if (err) comment << "Error in netgen::OCCGenerateMesh() at MESHCONST_ANALYSE step";
612 // fill ngMesh with nodes and elements of computed submeshes
613 err = ! fillNgMesh(occgeo, *ngMesh, nodeVec, meshedSM);
614 nbInitNod = ngMesh->GetNP();
615 nbInitSeg = ngMesh->GetNSeg();
616 nbInitFac = ngMesh->GetNSE();
621 startWith = endWith = netgen::MESHCONST_MESHEDGES;
622 err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
623 if (err) comment << "Error in netgen::OCCGenerateMesh() at 1D mesh generation";
625 // ---------------------
626 // compute surface mesh
627 // ---------------------
630 // pass 2D simple parameters to NETGEN
632 if ( double area = _simpleHyp->GetMaxElementArea() ) {
634 mparams.maxh = sqrt(2. * area/sqrt(3.0));
635 mparams.grading = 0.4; // moderate size growth
640 TopTools_MapOfShape tmpMap;
641 for ( TopExp_Explorer exp( _shape, TopAbs_EDGE ); exp.More(); exp.Next() )
642 if( tmpMap.Add(exp.Current()) )
643 length += SMESH_Algo::EdgeLength( TopoDS::Edge( exp.Current() ));
645 if ( ngMesh->GetNSeg() ) {
646 // we have to multiply length by 2 since for each TopoDS_Edge there
647 // are double set of NETGEN edges or, in other words, we have to
648 // divide ngMesh->GetNSeg() on 2.
649 mparams.maxh = 2*length / ngMesh->GetNSeg();
653 mparams.grading = 0.2; // slow size growth
655 mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
656 ngMesh->SetGlobalH (mparams.maxh);
657 netgen::Box<3> bb = occgeo.GetBoundingBox();
658 bb.Increase (bb.Diam()/20);
659 ngMesh->SetLocalH (bb.PMin(), bb.PMax(), mparams.grading);
661 // let netgen compute 2D mesh
662 startWith = netgen::MESHCONST_MESHSURFACE;
663 endWith = _optimize ? netgen::MESHCONST_OPTSURFACE : netgen::MESHCONST_MESHSURFACE;
664 err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
665 if (err) comment << "Error in netgen::OCCGenerateMesh() at surface mesh generation";
667 // ---------------------
668 // generate volume mesh
669 // ---------------------
670 if (!err && _isVolume)
672 // add ng face descriptors of meshed faces
673 std::map< int, std::pair<int,int> >::iterator fId_soIds = _faceDescriptors.begin();
674 for ( ; fId_soIds != _faceDescriptors.end(); ++fId_soIds ) {
675 int faceID = fId_soIds->first;
676 int solidID1 = fId_soIds->second.first;
677 int solidID2 = fId_soIds->second.second;
678 ngMesh->AddFaceDescriptor (netgen::FaceDescriptor(faceID, solidID1, solidID2, 0));
680 // pass 3D simple parameters to NETGEN
681 const NETGENPlugin_SimpleHypothesis_3D* simple3d =
682 dynamic_cast< const NETGENPlugin_SimpleHypothesis_3D* > ( _simpleHyp );
684 if ( double vol = simple3d->GetMaxElementVolume() ) {
686 mparams.maxh = pow( 72, 1/6. ) * pow( vol, 1/3. );
687 mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
691 mparams.maxh = ngMesh->AverageH();
693 // netgen::ARRAY<double> maxhdom;
694 // maxhdom.SetSize (occgeo.NrSolids());
695 // maxhdom = mparams.maxh;
696 // ngMesh->SetMaxHDomain (maxhdom);
697 ngMesh->SetGlobalH (mparams.maxh);
698 mparams.grading = 0.4;
699 ngMesh->CalcLocalH();
701 // let netgen compute 3D mesh
702 startWith = netgen::MESHCONST_MESHVOLUME;
703 endWith = _optimize ? netgen::MESHCONST_OPTVOLUME : netgen::MESHCONST_MESHVOLUME;
704 err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
705 if (err) comment << "Error in netgen::OCCGenerateMesh()";
707 if (!err && mparams.secondorder > 0)
709 netgen::OCCRefinementSurfaces ref (occgeo);
710 ref.MakeSecondOrder (*ngMesh);
713 catch (netgen::NgException exc)
715 error->myName = err = COMPERR_ALGO_FAILED;
716 comment << exc.What();
719 int nbNod = ngMesh->GetNP();
720 int nbSeg = ngMesh->GetNSeg();
721 int nbFac = ngMesh->GetNSE();
722 int nbVol = ngMesh->GetNE();
724 MESSAGE((err ? "Mesh Generation failure" : "End of Mesh Generation") <<
725 ", nb nodes: " << nbNod <<
726 ", nb segments: " << nbSeg <<
727 ", nb faces: " << nbFac <<
728 ", nb volumes: " << nbVol);
730 // -----------------------------------------------------------
731 // Feed back the SMESHDS with the generated Nodes and Elements
732 // -----------------------------------------------------------
734 SMESHDS_Mesh* meshDS = _mesh->GetMeshDS();
735 bool isOK = ( !err && (_isVolume ? (nbVol > 0) : (nbFac > 0)) );
736 if ( true /*isOK*/ ) // get whatever built
738 // map of nodes assigned to submeshes
739 NCollection_Map<int> pindMap;
740 // create and insert nodes into nodeVec
741 nodeVec.resize( nbNod + 1 );
743 for (i = nbInitNod+1; i <= nbNod /*&& isOK*/; ++i )
745 const netgen::MeshPoint& ngPoint = ngMesh->Point(i);
746 SMDS_MeshNode* node = NULL;
747 bool newNodeOnVertex = false;
749 if (i-nbInitNod <= occgeo.vmap.Extent())
752 aVert = TopoDS::Vertex(occgeo.vmap(i-nbInitNod));
753 SMESHDS_SubMesh * submesh = meshDS->MeshElements(aVert);
756 SMDS_NodeIteratorPtr it = submesh->GetNodes();
759 node = const_cast<SMDS_MeshNode*> (it->next());
764 newNodeOnVertex = true;
768 node = meshDS->AddNode(ngPoint(0), ngPoint(1), ngPoint(2));
770 node = meshDS->AddNode(ngPoint.X(), ngPoint.Y(), ngPoint.Z());
774 MESSAGE("Cannot create a mesh node");
775 if ( !comment.size() ) comment << "Cannot create a mesh node";
776 nbSeg = nbFac = nbVol = isOK = 0;
779 nodeVec.at(i) = node;
783 meshDS->SetNodeOnVertex(node, aVert);
788 // create mesh segments along geometric edges
789 NCollection_Map<Link> linkMap;
790 for (i = nbInitSeg+1; i <= nbSeg/* && isOK*/; ++i )
792 const netgen::Segment& seg = ngMesh->LineSegment(i);
794 Link link(seg.pnums[0], seg.pnums[1]);
796 Link link(seg.p1, seg.p2);
798 if (linkMap.Contains(link))
803 int pinds[3] = { seg.pnums[0], seg.pnums[1], seg.pnums[2] };
805 int pinds[3] = { seg.p1, seg.p2, seg.pmid };
809 for (int j=0; j < 3; ++j)
812 if (pind <= 0) continue;
819 int aGeomEdgeInd = seg.epgeominfo[j].edgenr;
820 if (aGeomEdgeInd > 0 && aGeomEdgeInd <= occgeo.emap.Extent())
821 aEdge = TopoDS::Edge(occgeo.emap(aGeomEdgeInd));
823 param = seg.epgeominfo[j].dist;
827 param = param2 * 0.5;
828 if (pind <= nbInitNod || pindMap.Contains(pind))
832 meshDS->SetNodeOnEdge(nodeVec.at(pind), aEdge, param);
837 if (nbp < 3) // second order ?
838 edge = meshDS->AddEdge(nodeVec.at(pinds[0]), nodeVec.at(pinds[1]));
840 edge = meshDS->AddEdge(nodeVec.at(pinds[0]), nodeVec.at(pinds[1]),
841 nodeVec.at(pinds[2]));
844 if ( !comment.size() ) comment << "Cannot create a mesh edge";
845 MESSAGE("Cannot create a mesh edge");
846 nbSeg = nbFac = nbVol = isOK = 0;
850 meshDS->SetMeshElementOnShape(edge, aEdge);
853 // create mesh faces along geometric faces
854 for (i = nbInitFac+1; i <= nbFac/* && isOK*/; ++i )
856 const netgen::Element2d& elem = ngMesh->SurfaceElement(i);
857 int aGeomFaceInd = elem.GetIndex();
859 if (aGeomFaceInd > 0 && aGeomFaceInd <= occgeo.fmap.Extent())
860 aFace = TopoDS::Face(occgeo.fmap(aGeomFaceInd));
861 vector<SMDS_MeshNode*> nodes;
862 for (int j=1; j <= elem.GetNP(); ++j)
864 int pind = elem.PNum(j);
865 SMDS_MeshNode* node = nodeVec.at(pind);
866 nodes.push_back(node);
867 if (pind <= nbInitNod || pindMap.Contains(pind))
871 const netgen::PointGeomInfo& pgi = elem.GeomInfoPi(j);
872 meshDS->SetNodeOnFace(node, aFace, pgi.u, pgi.v);
876 SMDS_MeshFace* face = NULL;
877 switch (elem.GetType())
880 face = meshDS->AddFace(nodes[0],nodes[1],nodes[2]);
883 face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3]);
886 face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[5],nodes[3],nodes[4]);
889 face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3],
890 nodes[4],nodes[7],nodes[5],nodes[6]);
893 MESSAGE("NETGEN created a face of unexpected type, ignoring");
898 if ( !comment.size() ) comment << "Cannot create a mesh face";
899 MESSAGE("Cannot create a mesh face");
900 nbSeg = nbFac = nbVol = isOK = 0;
904 meshDS->SetMeshElementOnShape(face, aFace);
908 for (i = 1; i <= nbVol/* && isOK*/; ++i)
910 const netgen::Element& elem = ngMesh->VolumeElement(i);
911 int aSolidInd = elem.GetIndex();
913 if (aSolidInd > 0 && aSolidInd <= occgeo.somap.Extent())
914 aSolid = TopoDS::Solid(occgeo.somap(aSolidInd));
915 vector<SMDS_MeshNode*> nodes;
916 for (int j=1; j <= elem.GetNP(); ++j)
918 int pind = elem.PNum(j);
919 SMDS_MeshNode* node = nodeVec.at(pind);
920 nodes.push_back(node);
921 if (pind <= nbInitNod || pindMap.Contains(pind))
923 if (!aSolid.IsNull())
926 meshDS->SetNodeInVolume(node, aSolid);
930 SMDS_MeshVolume* vol = NULL;
931 switch (elem.GetType())
934 vol = meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3]);
937 vol = meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3],
938 nodes[4],nodes[7],nodes[5],nodes[6],nodes[8],nodes[9]);
941 MESSAGE("NETGEN created a volume of unexpected type, ignoring");
946 if ( !comment.size() ) comment << "Cannot create a mesh volume";
947 MESSAGE("Cannot create a mesh volume");
948 nbSeg = nbFac = nbVol = isOK = 0;
951 if (!aSolid.IsNull())
952 meshDS->SetMeshElementOnShape(vol, aSolid);
956 if ( error->IsOK() && ( !isOK || comment.size() > 0 ))
957 error->myName = COMPERR_ALGO_FAILED;
958 if ( !comment.empty() )
959 error->myComment = comment;
961 // set bad compute error to subshapes of all failed subshapes shapes
962 if ( !error->IsOK() && err )
964 for (int i = 1; i <= occgeo.fmap.Extent(); i++) {
965 int status = occgeo.facemeshstatus[i-1];
966 if (status == 1 ) continue;
967 if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( occgeo.fmap( i ))) {
968 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
969 if ( !smError || smError->IsOK() ) {
971 smError.reset( new SMESH_ComputeError( error->myName, error->myComment ));
973 smError.reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED, "Ignored" ));
979 nglib::Ng_DeleteMesh((nglib::Ng_Mesh*)ngMesh);
984 return error->IsOK();
987 //=============================================================================
991 //=============================================================================
992 bool NETGENPlugin_Mesher::Evaluate(MapShapeNbElems& aResMap)
995 netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters();
997 netgen::MeshingParameters& mparams = netgen::mparam;
1001 // -------------------------
1002 // Prepare OCC geometry
1003 // -------------------------
1004 netgen::OCCGeometry occgeo;
1005 PrepareOCCgeometry( occgeo, _shape, *_mesh );
1007 bool tooManyElems = false;
1008 const int hugeNb = std::numeric_limits<int>::max() / 100;
1013 // pass 1D simple parameters to NETGEN
1015 if ( int nbSeg = _simpleHyp->GetNumberOfSegments() ) {
1017 mparams.segmentsperedge = nbSeg + 0.1;
1018 mparams.maxh = occgeo.boundingbox.Diam();
1019 mparams.grading = 0.01;
1023 mparams.segmentsperedge = 1;
1024 mparams.maxh = _simpleHyp->GetLocalLength();
1027 // let netgen create ngMesh and calculate element size on not meshed shapes
1029 netgen::Mesh *ngMesh = NULL;
1031 int startWith = netgen::MESHCONST_ANALYSE;
1032 int endWith = netgen::MESHCONST_MESHEDGES;
1033 int err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
1035 if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( _shape ))
1036 sm->GetComputeError().reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED ));
1040 // calculate total nb of segments and length of edges
1041 double fullLen = 0.0;
1043 int entity = mparams.secondorder > 0 ? SMDSEntity_Quad_Edge : SMDSEntity_Edge;
1044 TopTools_DataMapOfShapeInteger Edge2NbSeg;
1045 for (TopExp_Explorer exp(_shape, TopAbs_EDGE); exp.More(); exp.Next())
1047 TopoDS_Edge E = TopoDS::Edge( exp.Current() );
1048 if( !Edge2NbSeg.Bind(E,0) )
1051 double aLen = SMESH_Algo::EdgeLength(E);
1054 vector<int>& aVec = aResMap[_mesh->GetSubMesh(E)];
1056 aVec.resize( SMDSEntity_Last, 0);
1058 fullNbSeg += aVec[ entity ];
1061 // store nb of segments computed by Netgen
1062 NCollection_Map<Link> linkMap;
1063 for (int i = 1; i <= ngMesh->GetNSeg(); ++i )
1065 const netgen::Segment& seg = ngMesh->LineSegment(i);
1067 Link link(seg.pnums[0], seg.pnums[1]);
1069 Link link(seg.p1, seg.p2);
1071 if ( !linkMap.Add( link )) continue;
1072 int aGeomEdgeInd = seg.epgeominfo[0].edgenr;
1073 if (aGeomEdgeInd > 0 && aGeomEdgeInd <= occgeo.emap.Extent())
1075 vector<int>& aVec = aResMap[_mesh->GetSubMesh(occgeo.emap(aGeomEdgeInd))];
1079 // store nb of nodes on edges computed by Netgen
1080 TopTools_DataMapIteratorOfDataMapOfShapeInteger Edge2NbSegIt(Edge2NbSeg);
1081 for (; Edge2NbSegIt.More(); Edge2NbSegIt.Next())
1083 vector<int>& aVec = aResMap[_mesh->GetSubMesh(Edge2NbSegIt.Key())];
1084 if ( aVec[ entity ] > 1 && aVec[ SMDSEntity_Node ] == 0 )
1085 aVec[SMDSEntity_Node] = mparams.secondorder > 0 ? 2*aVec[ entity ]-1 : aVec[ entity ]-1;
1087 fullNbSeg += aVec[ entity ];
1088 Edge2NbSeg( Edge2NbSegIt.Key() ) = aVec[ entity ];
1090 nglib::Ng_DeleteMesh((nglib::Ng_Mesh*)ngMesh);
1097 if ( double area = _simpleHyp->GetMaxElementArea() ) {
1099 mparams.maxh = sqrt(2. * area/sqrt(3.0));
1100 mparams.grading = 0.4; // moderate size growth
1103 // length from edges
1104 mparams.maxh = fullLen/fullNbSeg;
1105 mparams.grading = 0.2; // slow size growth
1108 mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
1109 mparams.maxh = min( mparams.maxh, fullLen/fullNbSeg * (1. + mparams.grading));
1111 for (TopExp_Explorer exp(_shape, TopAbs_FACE); exp.More(); exp.Next())
1113 TopoDS_Face F = TopoDS::Face( exp.Current() );
1114 SMESH_subMesh *sm = _mesh->GetSubMesh(F);
1116 BRepGProp::SurfaceProperties(F,G);
1117 double anArea = G.Mass();
1118 tooManyElems = tooManyElems || ( anArea/hugeNb > mparams.maxh*mparams.maxh );
1120 if ( !tooManyElems )
1122 TopTools_MapOfShape egdes;
1123 for (TopExp_Explorer exp1(F,TopAbs_EDGE); exp1.More(); exp1.Next())
1124 if ( egdes.Add( exp1.Current() ))
1125 nb1d += Edge2NbSeg.Find(exp1.Current());
1127 int nbFaces = tooManyElems ? hugeNb : int( 4*anArea / (mparams.maxh*mparams.maxh*sqrt(3.)));
1128 int nbNodes = tooManyElems ? hugeNb : (( nbFaces*3 - (nb1d-1)*2 ) / 6 + 1 );
1130 vector<int> aVec(SMDSEntity_Last, 0);
1131 if( mparams.secondorder > 0 ) {
1132 int nb1d_in = (nbFaces*3 - nb1d) / 2;
1133 aVec[SMDSEntity_Node] = nbNodes + nb1d_in;
1134 aVec[SMDSEntity_Quad_Triangle] = nbFaces;
1137 aVec[SMDSEntity_Node] = nbNodes;
1138 aVec[SMDSEntity_Triangle] = nbFaces;
1140 aResMap[sm].swap(aVec);
1147 // pass 3D simple parameters to NETGEN
1148 const NETGENPlugin_SimpleHypothesis_3D* simple3d =
1149 dynamic_cast< const NETGENPlugin_SimpleHypothesis_3D* > ( _simpleHyp );
1151 if ( double vol = simple3d->GetMaxElementVolume() ) {
1153 mparams.maxh = pow( 72, 1/6. ) * pow( vol, 1/3. );
1154 mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
1157 // using previous length from faces
1159 mparams.grading = 0.4;
1160 mparams.maxh = min( mparams.maxh, fullLen/fullNbSeg * (1. + mparams.grading));
1163 BRepGProp::VolumeProperties(_shape,G);
1164 double aVolume = G.Mass();
1165 double tetrVol = 0.1179*mparams.maxh*mparams.maxh*mparams.maxh;
1166 tooManyElems = tooManyElems || ( aVolume/hugeNb > tetrVol );
1167 int nbVols = tooManyElems ? hugeNb : int(aVolume/tetrVol);
1168 int nb1d_in = int(( nbVols*6 - fullNbSeg ) / 6 );
1169 vector<int> aVec(SMDSEntity_Last, 0 );
1170 if ( tooManyElems ) // avoid FPE
1172 aVec[SMDSEntity_Node] = hugeNb;
1173 aVec[ mparams.secondorder > 0 ? SMDSEntity_Quad_Tetra : SMDSEntity_Tetra] = hugeNb;
1177 if( mparams.secondorder > 0 ) {
1178 aVec[SMDSEntity_Node] = nb1d_in/3 + 1 + nb1d_in;
1179 aVec[SMDSEntity_Quad_Tetra] = nbVols;
1182 aVec[SMDSEntity_Node] = nb1d_in/3 + 1;
1183 aVec[SMDSEntity_Tetra] = nbVols;
1186 SMESH_subMesh *sm = _mesh->GetSubMesh(_shape);
1187 aResMap[sm].swap(aVec);
1193 //================================================================================
1195 * \brief Remove "test.out" and "problemfaces" files in current directory
1197 //================================================================================
1199 void NETGENPlugin_Mesher::RemoveTmpFiles()
1201 removeFile("test.out");
1202 removeFile("problemfaces");
1203 removeFile("occmesh.rep");