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)
28 //=============================================================================
30 #include "NETGENPlugin_Mesher.hxx"
31 #include "NETGENPlugin_Hypothesis_2D.hxx"
32 #include "NETGENPlugin_SimpleHypothesis_3D.hxx"
34 #include <SMESH_Mesh.hxx>
35 #include <SMESH_Comment.hxx>
36 #include <SMESH_ComputeError.hxx>
37 #include <SMESH_subMesh.hxx>
38 #include <SMESH_MesherHelper.hxx>
39 #include <SMESHDS_Mesh.hxx>
40 #include <SMDS_MeshElement.hxx>
41 #include <SMDS_MeshNode.hxx>
42 #include <utilities.h>
46 #include <BRep_Tool.hxx>
48 #include <TopExp_Explorer.hxx>
50 #include <NCollection_Map.hxx>
51 #include <OSD_Path.hxx>
52 #include <OSD_File.hxx>
53 #include <TCollection_AsciiString.hxx>
54 #include <TopTools_ListIteratorOfListOfShape.hxx>
55 #include <Standard_ErrorHandler.hxx>
57 // Netgen include files
62 #include <occgeom.hpp>
63 #include <meshing.hpp>
64 //#include <ngexception.hpp>
66 extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, int, int, char*);
67 extern MeshingParameters mparam;
72 //=============================================================================
76 //=============================================================================
78 NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh* mesh,
79 const TopoDS_Shape& aShape,
90 //================================================================================
92 * \brief Initialize global NETGEN parameters with default values
94 //================================================================================
96 void NETGENPlugin_Mesher::defaultParameters()
99 netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters();
101 netgen::MeshingParameters& mparams = netgen::mparam;
103 // maximal mesh edge size
104 mparams.maxh = NETGENPlugin_Hypothesis::GetDefaultMaxSize();
105 // minimal number of segments per edge
106 mparams.segmentsperedge = NETGENPlugin_Hypothesis::GetDefaultNbSegPerEdge();
107 // rate of growth of size between elements
108 mparams.grading = NETGENPlugin_Hypothesis::GetDefaultGrowthRate();
109 // safety factor for curvatures (elements per radius)
110 mparams.curvaturesafety = NETGENPlugin_Hypothesis::GetDefaultNbSegPerRadius();
111 // create elements of second order
112 mparams.secondorder = NETGENPlugin_Hypothesis::GetDefaultSecondOrder() ? 1 : 0;
113 // quad-dominated surface meshing
117 mparams.quad = NETGENPlugin_Hypothesis_2D::GetDefaultQuadAllowed() ? 1 : 0;
120 //=============================================================================
122 * Pass parameters to NETGEN
124 //=============================================================================
125 void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp)
130 netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters();
132 netgen::MeshingParameters& mparams = netgen::mparam;
134 // Initialize global NETGEN parameters:
135 // maximal mesh segment size
136 mparams.maxh = hyp->GetMaxSize();
137 // minimal number of segments per edge
138 mparams.segmentsperedge = hyp->GetNbSegPerEdge();
139 // rate of growth of size between elements
140 mparams.grading = hyp->GetGrowthRate();
141 // safety factor for curvatures (elements per radius)
142 mparams.curvaturesafety = hyp->GetNbSegPerRadius();
143 // create elements of second order
144 mparams.secondorder = hyp->GetSecondOrder() ? 1 : 0;
145 // quad-dominated surface meshing
146 // only triangles are allowed for volumic mesh
148 mparams.quad = static_cast<const NETGENPlugin_Hypothesis_2D*>
149 (hyp)->GetQuadAllowed() ? 1 : 0;
150 _optimize = hyp->GetOptimize();
155 //=============================================================================
157 * Pass simple parameters to NETGEN
159 //=============================================================================
161 void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_SimpleHypothesis_2D* hyp)
168 //=============================================================================
170 * Link - a pair of integer numbers
172 //=============================================================================
176 Link(int _n1, int _n2) : n1(_n1), n2(_n2) {}
177 Link() : n1(0), n2(0) {}
180 int HashCode(const Link& aLink, int aLimit)
182 return HashCode(aLink.n1 + aLink.n2, aLimit);
185 Standard_Boolean IsEqual(const Link& aLink1, const Link& aLink2)
187 return (aLink1.n1 == aLink2.n1 && aLink1.n2 == aLink2.n2 ||
188 aLink1.n1 == aLink2.n2 && aLink1.n2 == aLink2.n1);
191 //================================================================================
193 * \brief Initialize netgen::OCCGeometry with OCCT shape
195 //================================================================================
197 void NETGENPlugin_Mesher::PrepareOCCgeometry(netgen::OCCGeometry& occgeo,
198 const TopoDS_Shape& shape,
200 list< SMESH_subMesh* > * meshedSM)
202 BRepTools::Clean (shape);
204 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
207 BRepMesh_IncrementalMesh::BRepMesh_IncrementalMesh (shape, 0.01, true);
208 } catch (Standard_Failure) {
211 BRepBndLib::Add (shape, bb);
212 double x1,y1,z1,x2,y2,z2;
213 bb.Get (x1,y1,z1,x2,y2,z2);
214 MESSAGE("shape bounding box:\n" <<
215 "(" << x1 << " " << y1 << " " << z1 << ") " <<
216 "(" << x2 << " " << y2 << " " << z2 << ")");
217 netgen::Point<3> p1 = netgen::Point<3> (x1,y1,z1);
218 netgen::Point<3> p2 = netgen::Point<3> (x2,y2,z2);
219 occgeo.boundingbox = netgen::Box<3> (p1,p2);
221 occgeo.shape = shape;
223 //occgeo.BuildFMap();
225 // fill maps of shapes of occgeo with not yet meshed subshapes
227 // get root submeshes
228 list< SMESH_subMesh* > rootSM;
229 if ( SMESH_subMesh* sm = mesh.GetSubMeshContaining( shape )) {
230 rootSM.push_back( sm );
233 for ( TopoDS_Iterator it( shape ); it.More(); it.Next() )
234 rootSM.push_back( mesh.GetSubMesh( it.Value() ));
237 // add subshapes of empty submeshes
238 list< SMESH_subMesh* >::iterator rootIt = rootSM.begin(), rootEnd = rootSM.end();
239 for ( ; rootIt != rootEnd; ++rootIt ) {
240 SMESH_subMesh * root = *rootIt;
241 SMESH_subMeshIteratorPtr smIt = root->getDependsOnIterator(/*includeSelf=*/true,
242 /*complexShapeFirst=*/true);
243 // to find a right orientation of subshapes (PAL20462)
244 TopTools_IndexedMapOfShape subShapes;
245 TopExp::MapShapes(root->GetSubShape(), subShapes);
246 while ( smIt->more() ) {
247 SMESH_subMesh* sm = smIt->next();
248 if ( sm->IsEmpty() ) {
249 TopoDS_Shape shape = sm->GetSubShape();
250 if ( shape.ShapeType() != TopAbs_VERTEX )
251 shape = subShapes( subShapes.FindIndex( shape ));// - shape->index->oriented shape
252 switch ( shape.ShapeType() ) {
253 case TopAbs_FACE : occgeo.fmap.Add( shape ); break;
254 case TopAbs_EDGE : occgeo.emap.Add( shape ); break;
255 case TopAbs_VERTEX: occgeo.vmap.Add( shape ); break;
256 case TopAbs_SOLID :occgeo.somap.Add( shape ); break;
260 // collect submeshes of meshed shapes
262 meshedSM->push_back( sm );
266 occgeo.facemeshstatus.SetSize (occgeo.fmap.Extent());
267 occgeo.facemeshstatus = 0;
271 //================================================================================
273 * \brief return id of netgen point corresponding to SMDS node
275 //================================================================================
277 static int ngNodeId( const SMDS_MeshNode* node,
278 netgen::Mesh& ngMesh,
279 map< const SMDS_MeshNode*, int >& nodeNgIdMap)
281 int newNgId = ngMesh.GetNP() + 1;
283 pair< map< const SMDS_MeshNode*, int >::iterator, bool > it_isNew =
284 nodeNgIdMap.insert( make_pair( node, newNgId ));
286 if ( it_isNew.second ) {
287 netgen::MeshPoint p( netgen::Point<3> (node->X(), node->Y(), node->Z()) );
288 ngMesh.AddPoint( p );
290 return it_isNew.first->second;
293 //================================================================================
295 * \brief fill ngMesh with nodes and elements of computed submeshes
297 //================================================================================
299 bool NETGENPlugin_Mesher::fillNgMesh(netgen::OCCGeometry& occgeom,
300 netgen::Mesh& ngMesh,
301 vector<SMDS_MeshNode*>& nodeVec,
302 const list< SMESH_subMesh* > & meshedSM)
304 map< const SMDS_MeshNode*, int > nodeNgIdMap;
306 TopTools_MapOfShape visitedShapes;
308 SMESH_MesherHelper helper (*_mesh);
310 int faceID = occgeom.fmap.Extent();
311 _faceDescriptors.clear();
313 list< SMESH_subMesh* >::const_iterator smIt, smEnd = meshedSM.end();
314 for ( smIt = meshedSM.begin(); smIt != smEnd; ++smIt )
316 SMESH_subMesh* sm = *smIt;
317 if ( !visitedShapes.Add( sm->GetSubShape() ))
320 SMESHDS_SubMesh * smDS = sm->GetSubMeshDS();
322 switch ( sm->GetSubShape().ShapeType() )
324 case TopAbs_EDGE: { // EDGE
325 // ----------------------
326 const TopoDS_Edge& geomEdge = TopoDS::Edge( sm->GetSubShape() );
328 // Add ng segments for each not meshed face the edge bounds
329 TopTools_MapOfShape visitedAncestors;
330 const TopTools_ListOfShape& ancestors = _mesh->GetAncestors( geomEdge );
331 TopTools_ListIteratorOfListOfShape ancestorIt ( ancestors );
332 for ( ; ancestorIt.More(); ancestorIt.Next() )
334 const TopoDS_Shape & ans = ancestorIt.Value();
335 if ( ans.ShapeType() != TopAbs_FACE || !visitedAncestors.Add( ans ))
337 const TopoDS_Face& face = TopoDS::Face( ans );
339 int faceID = occgeom.fmap.FindIndex( face );
341 continue; // meshed face
343 // find out orientation of geomEdge within face
344 bool isForwad = false;
345 for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next() ) {
346 if ( geomEdge.IsSame( exp.Current() )) {
347 isForwad = ( exp.Current().Orientation() == geomEdge.Orientation() );
351 bool isQuad = smDS->GetElements()->next()->IsQuadratic();
353 // get all nodes from geomEdge
354 StdMeshers_FaceSide fSide( face, geomEdge, _mesh, isForwad, isQuad );
355 const vector<UVPtStruct>& points = fSide.GetUVPtStruct();
356 int i, nbSeg = fSide.NbSegments();
358 double otherSeamParam = 0;
359 helper.SetSubShape( face );
360 bool isSeam = helper.IsRealSeam( geomEdge );
363 helper.GetOtherParam( helper.GetPeriodicIndex() == 1 ? points[0].u : points[0].v );
367 int prevNgId = ngNodeId( points[0].node, ngMesh, nodeNgIdMap );
369 for ( i = 0; i < nbSeg; ++i )
371 const UVPtStruct& p1 = points[ i ];
372 const UVPtStruct& p2 = points[ i+1 ];
377 seg.p2 = prevNgId = ngNodeId( p2.node, ngMesh, nodeNgIdMap );
378 // node param on curve
379 seg.epgeominfo[ 0 ].dist = p1.param;
380 seg.epgeominfo[ 1 ].dist = p2.param;
382 seg.epgeominfo[ 0 ].u = p1.u;
383 seg.epgeominfo[ 0 ].v = p1.v;
384 seg.epgeominfo[ 1 ].u = p2.u;
385 seg.epgeominfo[ 1 ].v = p2.v;
387 //seg.epgeominfo[ iEnd ].edgenr = edgeID; // = geom.emap.FindIndex(edge);
388 seg.si = faceID; // = geom.fmap.FindIndex (face);
389 seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
390 ngMesh.AddSegment (seg);
394 if ( helper.GetPeriodicIndex() == 1 ) {
395 seg.epgeominfo[ 0 ].u = otherSeamParam;
396 seg.epgeominfo[ 1 ].u = otherSeamParam;
397 swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v);
399 seg.epgeominfo[ 0 ].v = otherSeamParam;
400 seg.epgeominfo[ 1 ].v = otherSeamParam;
401 swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u);
403 swap (seg.p1, seg.p2);
404 swap (seg.epgeominfo[0].dist, seg.epgeominfo[1].dist);
405 seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
406 ngMesh.AddSegment (seg);
409 } // loop on geomEdge ancestors
412 } // case TopAbs_EDGE
414 case TopAbs_FACE: { // FACE
415 // ----------------------
416 const TopoDS_Face& geomFace = TopoDS::Face( sm->GetSubShape() );
417 helper.SetSubShape( geomFace );
419 // find solids geomFace bounds
420 int solidID1 = 0, solidID2 = 0;
421 const TopTools_ListOfShape& ancestors = _mesh->GetAncestors( geomFace );
422 TopTools_ListIteratorOfListOfShape ancestorIt ( ancestors );
423 for ( ; ancestorIt.More(); ancestorIt.Next() )
425 const TopoDS_Shape & solid = ancestorIt.Value();
426 if ( solid.ShapeType() == TopAbs_SOLID ) {
427 int id = occgeom.somap.FindIndex ( solid );
428 if ( solidID1 && id != solidID1 ) solidID2 = id;
433 _faceDescriptors[ faceID ].first = solidID1;
434 _faceDescriptors[ faceID ].second = solidID2;
436 // add surface elements
437 SMDS_ElemIteratorPtr faces = smDS->GetElements();
438 while ( faces->more() ) {
440 const SMDS_MeshElement* f = faces->next();
441 if ( f->NbNodes() % 3 != 0 ) { // not triangle
442 for ( ancestorIt.Initialize(ancestors); ancestorIt.More(); ancestorIt.Next() )
443 if ( ancestorIt.Value().ShapeType() == TopAbs_SOLID ) {
444 sm = _mesh->GetSubMesh( ancestorIt.Value() );
447 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
448 smError.reset( new SMESH_ComputeError(COMPERR_BAD_INPUT_MESH,"Not triangle submesh"));
449 smError->myBadElements.push_back( f );
453 netgen::Element2d tri(3);
454 tri.SetIndex ( faceID );
456 for ( int i = 0; i < 3; ++i ) {
457 const SMDS_MeshNode* node = f->GetNode( i ), * inFaceNode=0;
458 if ( helper.IsSeamShape( node->GetPosition()->GetShapeId() ))
459 if ( helper.IsSeamShape( f->GetNode( i+1 )->GetPosition()->GetShapeId() ))
460 inFaceNode = f->GetNode( i-1 );
462 inFaceNode = f->GetNode( i+1 );
464 gp_XY uv = helper.GetNodeUV( geomFace, node, inFaceNode );
465 tri.GeomInfoPi(i+1).u = uv.X();
466 tri.GeomInfoPi(i+1).v = uv.Y();
467 tri.PNum(i+1) = ngNodeId( node, ngMesh, nodeNgIdMap );
470 ngMesh.AddSurfaceElement (tri);
476 case TopAbs_VERTEX: { // VERTEX
477 // --------------------------
478 SMDS_NodeIteratorPtr nodeIt = smDS->GetNodes();
479 if ( nodeIt->more() )
480 ngNodeId( nodeIt->next(), ngMesh, nodeNgIdMap );
485 } // loop on submeshes
488 nodeVec.resize( ngMesh.GetNP() + 1 );
489 map< const SMDS_MeshNode*, int >::iterator node_NgId, nodeNgIdEnd = nodeNgIdMap.end();
490 for ( node_NgId = nodeNgIdMap.begin(); node_NgId != nodeNgIdEnd; ++node_NgId)
491 nodeVec[ node_NgId->second ] = (SMDS_MeshNode*) node_NgId->first;
496 //=============================================================================
498 * Here we are going to use the NETGEN mesher
500 //=============================================================================
501 bool NETGENPlugin_Mesher::Compute()
504 netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters();
506 netgen::MeshingParameters& mparams = netgen::mparam;
508 MESSAGE("Compute with:\n"
509 " max size = " << mparams.maxh << "\n"
510 " segments per edge = " << mparams.segmentsperedge);
512 " growth rate = " << mparams.grading << "\n"
513 " elements per radius = " << mparams.curvaturesafety << "\n"
514 " second order = " << mparams.secondorder << "\n"
515 " quad allowed = " << mparams.quad);
517 SMESH_ComputeErrorPtr error = SMESH_ComputeError::New();
520 // -------------------------
521 // Prepare OCC geometry
522 // -------------------------
524 netgen::OCCGeometry occgeo;
525 list< SMESH_subMesh* > meshedSM;
526 PrepareOCCgeometry( occgeo, _shape, *_mesh, &meshedSM );
528 // -------------------------
530 // -------------------------
532 netgen::Mesh *ngMesh = NULL;
534 SMESH_Comment comment;
539 // vector of nodes in which node index == netgen ID
540 vector< SMDS_MeshNode* > nodeVec;
546 // pass 1D simple parameters to NETGEN
548 if ( int nbSeg = _simpleHyp->GetNumberOfSegments() ) {
550 mparams.segmentsperedge = nbSeg + 0.1;
551 mparams.maxh = occgeo.boundingbox.Diam();
556 mparams.segmentsperedge = 1;
557 mparams.maxh = _simpleHyp->GetLocalLength();
560 // let netgen create ngMesh and calculate element size on not meshed shapes
562 int startWith = netgen::MESHCONST_ANALYSE;
563 int endWith = netgen::MESHCONST_ANALYSE;
564 err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
565 if (err) comment << "Error in netgen::OCCGenerateMesh() at MESHCONST_ANALYSE step";
567 // fill ngMesh with nodes and elements of computed submeshes
568 err = ! fillNgMesh(occgeo, *ngMesh, nodeVec, meshedSM);
569 nbInitNod = ngMesh->GetNP();
570 nbInitSeg = ngMesh->GetNSeg();
571 nbInitFac = ngMesh->GetNSE();
576 startWith = endWith = netgen::MESHCONST_MESHEDGES;
577 err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
578 if (err) comment << "Error in netgen::OCCGenerateMesh() at 1D mesh generation";
580 // ---------------------
581 // compute surface mesh
582 // ---------------------
585 // pass 2D simple parameters to NETGEN
587 if ( double area = _simpleHyp->GetMaxElementArea() ) {
589 mparams.maxh = sqrt(2. * area/sqrt(3.0));
590 mparams.grading = 0.4; // moderate size growth
595 for ( TopExp_Explorer exp( _shape, TopAbs_EDGE ); exp.More(); exp.Next() )
596 length += SMESH_Algo::EdgeLength( TopoDS::Edge( exp.Current() ));
597 if ( ngMesh->GetNSeg() )
598 mparams.maxh = length / ngMesh->GetNSeg();
601 mparams.grading = 0.2; // slow size growth
603 mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
604 ngMesh->SetGlobalH (mparams.maxh);
605 netgen::Box<3> bb = occgeo.GetBoundingBox();
606 bb.Increase (bb.Diam()/20);
607 ngMesh->SetLocalH (bb.PMin(), bb.PMax(), mparams.grading);
609 // let netgen compute 2D mesh
610 startWith = netgen::MESHCONST_MESHSURFACE;
611 endWith = _optimize ? netgen::MESHCONST_OPTSURFACE : netgen::MESHCONST_MESHSURFACE;
612 err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
613 if (err) comment << "Error in netgen::OCCGenerateMesh() at surface mesh generation";
615 // ---------------------
616 // generate volume mesh
617 // ---------------------
618 if (!err && _isVolume)
620 // add ng face descriptors of meshed faces
621 std::map< int, std::pair<int,int> >::iterator fId_soIds = _faceDescriptors.begin();
622 for ( ; fId_soIds != _faceDescriptors.end(); ++fId_soIds ) {
623 int faceID = fId_soIds->first;
624 int solidID1 = fId_soIds->second.first;
625 int solidID2 = fId_soIds->second.second;
626 ngMesh->AddFaceDescriptor (netgen::FaceDescriptor(faceID, solidID1, solidID2, 0));
628 // pass 3D simple parameters to NETGEN
629 const NETGENPlugin_SimpleHypothesis_3D* simple3d =
630 dynamic_cast< const NETGENPlugin_SimpleHypothesis_3D* > ( _simpleHyp );
632 if ( double vol = simple3d->GetMaxElementVolume() ) {
634 mparams.maxh = pow( 72 * vol * vol, 1/3. );
635 mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
639 mparams.maxh = ngMesh->AverageH();
641 // netgen::ARRAY<double> maxhdom;
642 // maxhdom.SetSize (occgeo.NrSolids());
643 // maxhdom = mparams.maxh;
644 // ngMesh->SetMaxHDomain (maxhdom);
645 ngMesh->SetGlobalH (mparams.maxh);
646 mparams.grading = 0.4;
647 ngMesh->CalcLocalH();
649 // let netgen compute 3D mesh
650 startWith = netgen::MESHCONST_MESHVOLUME;
651 endWith = _optimize ? netgen::MESHCONST_OPTVOLUME : netgen::MESHCONST_MESHVOLUME;
652 err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
653 if (err) comment << "Error in netgen::OCCGenerateMesh()";
655 if (!err && mparams.secondorder > 0)
657 netgen::OCCRefinementSurfaces ref (occgeo);
658 ref.MakeSecondOrder (*ngMesh);
661 catch (netgen::NgException exc)
663 error->myName = err = COMPERR_ALGO_FAILED;
664 comment << exc.What();
667 int nbNod = ngMesh->GetNP();
668 int nbSeg = ngMesh->GetNSeg();
669 int nbFac = ngMesh->GetNSE();
670 int nbVol = ngMesh->GetNE();
672 MESSAGE((err ? "Mesh Generation failure" : "End of Mesh Generation") <<
673 ", nb nodes: " << nbNod <<
674 ", nb segments: " << nbSeg <<
675 ", nb faces: " << nbFac <<
676 ", nb volumes: " << nbVol);
678 // -----------------------------------------------------------
679 // Feed back the SMESHDS with the generated Nodes and Elements
680 // -----------------------------------------------------------
682 SMESHDS_Mesh* meshDS = _mesh->GetMeshDS();
683 bool isOK = ( !err && (_isVolume ? (nbVol > 0) : (nbFac > 0)) );
684 if ( true /*isOK*/ ) // get whatever built
686 // map of nodes assigned to submeshes
687 NCollection_Map<int> pindMap;
688 // create and insert nodes into nodeVec
689 nodeVec.resize( nbNod + 1 );
691 for (i = nbInitNod+1; i <= nbNod /*&& isOK*/; ++i )
693 const netgen::MeshPoint& ngPoint = ngMesh->Point(i);
694 SMDS_MeshNode* node = NULL;
695 bool newNodeOnVertex = false;
697 if (i-nbInitNod <= occgeo.vmap.Extent())
700 aVert = TopoDS::Vertex(occgeo.vmap(i-nbInitNod));
701 SMESHDS_SubMesh * submesh = meshDS->MeshElements(aVert);
704 SMDS_NodeIteratorPtr it = submesh->GetNodes();
707 node = const_cast<SMDS_MeshNode*> (it->next());
712 newNodeOnVertex = true;
715 node = meshDS->AddNode(ngPoint.X(), ngPoint.Y(), ngPoint.Z());
718 MESSAGE("Cannot create a mesh node");
719 if ( !comment.size() ) comment << "Cannot create a mesh node";
720 nbSeg = nbFac = nbVol = isOK = 0;
723 nodeVec.at(i) = node;
727 meshDS->SetNodeOnVertex(node, aVert);
732 // create mesh segments along geometric edges
733 NCollection_Map<Link> linkMap;
734 for (i = nbInitSeg+1; i <= nbSeg/* && isOK*/; ++i )
736 const netgen::Segment& seg = ngMesh->LineSegment(i);
737 Link link(seg.p1, seg.p2);
738 if (linkMap.Contains(link))
742 int pinds[3] = { seg.p1, seg.p2, seg.pmid };
745 for (int j=0; j < 3; ++j)
748 if (pind <= 0) continue;
755 int aGeomEdgeInd = seg.epgeominfo[j].edgenr;
756 if (aGeomEdgeInd > 0 && aGeomEdgeInd <= occgeo.emap.Extent())
757 aEdge = TopoDS::Edge(occgeo.emap(aGeomEdgeInd));
759 param = seg.epgeominfo[j].dist;
763 param = param2 * 0.5;
764 if (pind <= nbInitNod || pindMap.Contains(pind))
768 meshDS->SetNodeOnEdge(nodeVec.at(pind), aEdge, param);
773 if (nbp < 3) // second order ?
774 edge = meshDS->AddEdge(nodeVec.at(pinds[0]), nodeVec.at(pinds[1]));
776 edge = meshDS->AddEdge(nodeVec.at(pinds[0]), nodeVec.at(pinds[1]),
777 nodeVec.at(pinds[2]));
780 if ( !comment.size() ) comment << "Cannot create a mesh edge";
781 MESSAGE("Cannot create a mesh edge");
782 nbSeg = nbFac = nbVol = isOK = 0;
786 meshDS->SetMeshElementOnShape(edge, aEdge);
789 // create mesh faces along geometric faces
790 for (i = nbInitFac+1; i <= nbFac/* && isOK*/; ++i )
792 const netgen::Element2d& elem = ngMesh->SurfaceElement(i);
793 int aGeomFaceInd = elem.GetIndex();
795 if (aGeomFaceInd > 0 && aGeomFaceInd <= occgeo.fmap.Extent())
796 aFace = TopoDS::Face(occgeo.fmap(aGeomFaceInd));
797 vector<SMDS_MeshNode*> nodes;
798 for (int j=1; j <= elem.GetNP(); ++j)
800 int pind = elem.PNum(j);
801 SMDS_MeshNode* node = nodeVec.at(pind);
802 nodes.push_back(node);
803 if (pind <= nbInitNod || pindMap.Contains(pind))
807 const netgen::PointGeomInfo& pgi = elem.GeomInfoPi(j);
808 meshDS->SetNodeOnFace(node, aFace, pgi.u, pgi.v);
812 SMDS_MeshFace* face = NULL;
813 switch (elem.GetType())
816 face = meshDS->AddFace(nodes[0],nodes[1],nodes[2]);
819 face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3]);
822 face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[5],nodes[3],nodes[4]);
825 face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3],
826 nodes[4],nodes[7],nodes[5],nodes[6]);
829 MESSAGE("NETGEN created a face of unexpected type, ignoring");
834 if ( !comment.size() ) comment << "Cannot create a mesh face";
835 MESSAGE("Cannot create a mesh face");
836 nbSeg = nbFac = nbVol = isOK = 0;
840 meshDS->SetMeshElementOnShape(face, aFace);
844 for (i = 1; i <= nbVol/* && isOK*/; ++i)
846 const netgen::Element& elem = ngMesh->VolumeElement(i);
847 int aSolidInd = elem.GetIndex();
849 if (aSolidInd > 0 && aSolidInd <= occgeo.somap.Extent())
850 aSolid = TopoDS::Solid(occgeo.somap(aSolidInd));
851 vector<SMDS_MeshNode*> nodes;
852 for (int j=1; j <= elem.GetNP(); ++j)
854 int pind = elem.PNum(j);
855 SMDS_MeshNode* node = nodeVec.at(pind);
856 nodes.push_back(node);
857 if (pind <= nbInitNod || pindMap.Contains(pind))
859 if (!aSolid.IsNull())
862 meshDS->SetNodeInVolume(node, aSolid);
866 SMDS_MeshVolume* vol = NULL;
867 switch (elem.GetType())
870 vol = meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3]);
873 vol = meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3],
874 nodes[4],nodes[7],nodes[5],nodes[6],nodes[8],nodes[9]);
877 MESSAGE("NETGEN created a volume of unexpected type, ignoring");
882 if ( !comment.size() ) comment << "Cannot create a mesh volume";
883 MESSAGE("Cannot create a mesh volume");
884 nbSeg = nbFac = nbVol = isOK = 0;
887 if (!aSolid.IsNull())
888 meshDS->SetMeshElementOnShape(vol, aSolid);
892 if ( error->IsOK() && ( !isOK || comment.size() > 0 ))
893 error->myName = COMPERR_ALGO_FAILED;
894 if ( !comment.empty() )
895 error->myComment = comment;
897 // set bad compute error to subshapes of all failed subshapes shapes
898 if ( !error->IsOK() && err )
900 for (int i = 1; i <= occgeo.fmap.Extent(); i++) {
901 int status = occgeo.facemeshstatus[i-1];
902 if (status == 1 ) continue;
903 if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( occgeo.fmap( i ))) {
904 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
905 if ( !smError || smError->IsOK() ) {
907 smError.reset( new SMESH_ComputeError( error->myName, error->myComment ));
909 smError.reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED, "Ignored" ));
915 nglib::Ng_DeleteMesh((nglib::Ng_Mesh*)ngMesh);
920 return error->IsOK();
923 //================================================================================
925 * \brief Remove "test.out" and "problemfaces" files in current directory
927 //================================================================================
929 void NETGENPlugin_Mesher::RemoveTmpFiles()
931 TCollection_AsciiString str("test.out");
932 OSD_Path path1( str );
933 OSD_File file1( path1 );
935 str = "problemfaces";
936 OSD_Path path2( str );
937 OSD_File file2( path2 );