1 // NETGENPlugin : C++ implementation
3 // Copyright (C) 2006 OPEN CASCADE, CEA/DEN, EDF R&D
5 // This library is free software; you can redistribute it and/or
6 // modify it under the terms of the GNU Lesser General Public
7 // License as published by the Free Software Foundation; either
8 // version 2.1 of the License.
10 // This library is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 // Lesser General Public License for more details.
15 // You should have received a copy of the GNU Lesser General Public
16 // License along with this library; if not, write to the Free Software
17 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
22 // File : NETGENPlugin_Mesher.cxx
23 // 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>
45 #include <BRep_Tool.hxx>
47 #include <TopExp_Explorer.hxx>
49 #include <NCollection_Map.hxx>
50 #include <OSD_Path.hxx>
51 #include <OSD_File.hxx>
52 #include <TCollection_AsciiString.hxx>
53 #include <TopTools_ListIteratorOfListOfShape.hxx>
55 // Netgen include files
60 #include <occgeom.hpp>
61 #include <meshing.hpp>
62 //#include <ngexception.hpp>
64 extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, int, int, char*);
65 extern MeshingParameters mparam;
70 //=============================================================================
74 //=============================================================================
76 NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh* mesh,
77 const TopoDS_Shape& aShape,
88 //================================================================================
90 * \brief Initialize global NETGEN parameters with default values
92 //================================================================================
94 void NETGENPlugin_Mesher::defaultParameters()
97 netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters();
99 netgen::MeshingParameters& mparams = netgen::mparam;
101 // maximal mesh edge size
102 mparams.maxh = NETGENPlugin_Hypothesis::GetDefaultMaxSize();
103 // minimal number of segments per edge
104 mparams.segmentsperedge = NETGENPlugin_Hypothesis::GetDefaultNbSegPerEdge();
105 // rate of growth of size between elements
106 mparams.grading = NETGENPlugin_Hypothesis::GetDefaultGrowthRate();
107 // safety factor for curvatures (elements per radius)
108 mparams.curvaturesafety = NETGENPlugin_Hypothesis::GetDefaultNbSegPerRadius();
109 // create elements of second order
110 mparams.secondorder = NETGENPlugin_Hypothesis::GetDefaultSecondOrder() ? 1 : 0;
111 // quad-dominated surface meshing
115 mparams.quad = NETGENPlugin_Hypothesis_2D::GetDefaultQuadAllowed() ? 1 : 0;
118 //=============================================================================
120 * Pass parameters to NETGEN
122 //=============================================================================
123 void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp)
128 netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters();
130 netgen::MeshingParameters& mparams = netgen::mparam;
132 // Initialize global NETGEN parameters:
133 // maximal mesh segment size
134 mparams.maxh = hyp->GetMaxSize();
135 // minimal number of segments per edge
136 mparams.segmentsperedge = hyp->GetNbSegPerEdge();
137 // rate of growth of size between elements
138 mparams.grading = hyp->GetGrowthRate();
139 // safety factor for curvatures (elements per radius)
140 mparams.curvaturesafety = hyp->GetNbSegPerRadius();
141 // create elements of second order
142 mparams.secondorder = hyp->GetSecondOrder() ? 1 : 0;
143 // quad-dominated surface meshing
144 // only triangles are allowed for volumic mesh
146 mparams.quad = static_cast<const NETGENPlugin_Hypothesis_2D*>
147 (hyp)->GetQuadAllowed() ? 1 : 0;
148 _optimize = hyp->GetOptimize();
153 //=============================================================================
155 * Pass simple parameters to NETGEN
157 //=============================================================================
159 void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_SimpleHypothesis_2D* hyp)
166 //=============================================================================
168 * Link - a pair of integer numbers
170 //=============================================================================
174 Link(int _n1, int _n2) : n1(_n1), n2(_n2) {}
175 Link() : n1(0), n2(0) {}
178 int HashCode(const Link& aLink, int aLimit)
180 return HashCode(aLink.n1 + aLink.n2, aLimit);
183 Standard_Boolean IsEqual(const Link& aLink1, const Link& aLink2)
185 return (aLink1.n1 == aLink2.n1 && aLink1.n2 == aLink2.n2 ||
186 aLink1.n1 == aLink2.n2 && aLink1.n2 == aLink2.n1);
189 //================================================================================
191 * \brief Initialize netgen::OCCGeometry with OCCT shape
193 //================================================================================
195 void NETGENPlugin_Mesher::PrepareOCCgeometry(netgen::OCCGeometry& occgeo,
196 const TopoDS_Shape& shape,
198 list< SMESH_subMesh* > * meshedSM)
200 BRepTools::Clean (shape);
201 BRepMesh_IncrementalMesh::BRepMesh_IncrementalMesh (shape, 0.01, true);
203 BRepBndLib::Add (shape, bb);
204 double x1,y1,z1,x2,y2,z2;
205 bb.Get (x1,y1,z1,x2,y2,z2);
206 MESSAGE("shape bounding box:\n" <<
207 "(" << x1 << " " << y1 << " " << z1 << ") " <<
208 "(" << x2 << " " << y2 << " " << z2 << ")");
209 netgen::Point<3> p1 = netgen::Point<3> (x1,y1,z1);
210 netgen::Point<3> p2 = netgen::Point<3> (x2,y2,z2);
211 occgeo.boundingbox = netgen::Box<3> (p1,p2);
213 occgeo.shape = shape;
215 //occgeo.BuildFMap();
217 // fill maps of shapes of occgeo with not yet meshed subshapes
219 // get root submeshes
220 list< SMESH_subMesh* > rootSM;
221 if ( SMESH_subMesh* sm = mesh.GetSubMeshContaining( shape )) {
222 rootSM.push_back( sm );
225 for ( TopoDS_Iterator it( shape ); it.More(); it.Next() )
226 rootSM.push_back( mesh.GetSubMesh( it.Value() ));
229 // add subshapes of empty submeshes
230 list< SMESH_subMesh* >::iterator rootIt = rootSM.begin(), rootEnd = rootSM.end();
231 for ( ; rootIt != rootEnd; ++rootIt ) {
232 SMESH_subMesh * root = *rootIt;
233 SMESH_subMeshIteratorPtr smIt = root->getDependsOnIterator(/*includeSelf=*/true,
234 /*complexShapeFirst=*/true);
235 while ( smIt->more() ) {
236 SMESH_subMesh* sm = smIt->next();
237 if ( sm->IsEmpty() ) {
238 switch ( sm->GetSubShape().ShapeType() ) {
239 case TopAbs_FACE : occgeo.fmap.Add( sm->GetSubShape() ); break;
240 case TopAbs_EDGE : occgeo.emap.Add( sm->GetSubShape() ); break;
241 case TopAbs_VERTEX: occgeo.vmap.Add( sm->GetSubShape() ); break;
242 case TopAbs_SOLID :occgeo.somap.Add( sm->GetSubShape() ); break;
246 // collect submeshes of meshed shapes
248 meshedSM->push_back( sm );
252 occgeo.facemeshstatus.SetSize (occgeo.fmap.Extent());
253 occgeo.facemeshstatus = 0;
257 //================================================================================
259 * \brief return id of netgen point corresponding to SMDS node
261 //================================================================================
263 static int ngNodeId( const SMDS_MeshNode* node,
264 netgen::Mesh& ngMesh,
265 map< const SMDS_MeshNode*, int >& nodeNgIdMap)
267 int newNgId = ngMesh.GetNP() + 1;
269 pair< map< const SMDS_MeshNode*, int >::iterator, bool > it_isNew =
270 nodeNgIdMap.insert( make_pair( node, newNgId ));
272 if ( it_isNew.second ) {
273 netgen::MeshPoint p( netgen::Point<3> (node->X(), node->Y(), node->Z()) );
274 ngMesh.AddPoint( p );
276 return it_isNew.first->second;
279 //================================================================================
281 * \brief fill ngMesh with nodes and elements of computed submeshes
283 //================================================================================
285 bool NETGENPlugin_Mesher::fillNgMesh(netgen::OCCGeometry& occgeom,
286 netgen::Mesh& ngMesh,
287 vector<SMDS_MeshNode*>& nodeVec,
288 const list< SMESH_subMesh* > & meshedSM)
290 map< const SMDS_MeshNode*, int > nodeNgIdMap;
292 TopTools_MapOfShape visitedShapes;
294 SMESH_MesherHelper helper (*_mesh);
296 int faceID = occgeom.fmap.Extent();
297 _faceDescriptors.clear();
299 list< SMESH_subMesh* >::const_iterator smIt, smEnd = meshedSM.end();
300 for ( smIt = meshedSM.begin(); smIt != smEnd; ++smIt )
302 SMESH_subMesh* sm = *smIt;
303 if ( !visitedShapes.Add( sm->GetSubShape() ))
306 SMESHDS_SubMesh * smDS = sm->GetSubMeshDS();
308 switch ( sm->GetSubShape().ShapeType() )
310 case TopAbs_EDGE: { // EDGE
311 // ----------------------
312 const TopoDS_Edge& geomEdge = TopoDS::Edge( sm->GetSubShape() );
314 // Add ng segments for each not meshed face the edge bounds
315 TopTools_MapOfShape visitedAncestors;
316 const TopTools_ListOfShape& ancestors = _mesh->GetAncestors( geomEdge );
317 TopTools_ListIteratorOfListOfShape ancestorIt ( ancestors );
318 for ( ; ancestorIt.More(); ancestorIt.Next() )
320 const TopoDS_Shape & ans = ancestorIt.Value();
321 if ( ans.ShapeType() != TopAbs_FACE || !visitedAncestors.Add( ans ))
323 const TopoDS_Face& face = TopoDS::Face( ans );
325 int faceID = occgeom.fmap.FindIndex( face );
327 continue; // meshed face
329 // find out orientation of geomEdge within face
330 bool isForwad = false;
331 for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next() ) {
332 if ( geomEdge.IsSame( exp.Current() )) {
333 isForwad = ( exp.Current().Orientation() == geomEdge.Orientation() );
337 bool isQuad = smDS->GetElements()->next()->IsQuadratic();
339 // get all nodes from geomEdge
340 StdMeshers_FaceSide fSide( face, geomEdge, _mesh, isForwad, isQuad );
341 const vector<UVPtStruct>& points = fSide.GetUVPtStruct();
342 int i, nbSeg = fSide.NbSegments();
344 double otherSeamParam = 0;
345 helper.SetSubShape( face );
346 bool isSeam = helper.IsRealSeam( geomEdge );
349 helper.GetOtherParam( helper.GetPeriodicIndex() == 1 ? points[0].u : points[0].v );
353 int prevNgId = ngNodeId( points[0].node, ngMesh, nodeNgIdMap );
355 for ( i = 0; i < nbSeg; ++i )
357 const UVPtStruct& p1 = points[ i ];
358 const UVPtStruct& p2 = points[ i+1 ];
363 seg.p2 = prevNgId = ngNodeId( p2.node, ngMesh, nodeNgIdMap );
364 // node param on curve
365 seg.epgeominfo[ 0 ].dist = p1.param;
366 seg.epgeominfo[ 1 ].dist = p2.param;
368 seg.epgeominfo[ 0 ].u = p1.u;
369 seg.epgeominfo[ 0 ].v = p1.v;
370 seg.epgeominfo[ 1 ].u = p2.u;
371 seg.epgeominfo[ 1 ].v = p2.v;
373 //seg.epgeominfo[ iEnd ].edgenr = edgeID; // = geom.emap.FindIndex(edge);
374 seg.si = faceID; // = geom.fmap.FindIndex (face);
375 seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
376 ngMesh.AddSegment (seg);
380 if ( helper.GetPeriodicIndex() == 1 ) {
381 seg.epgeominfo[ 0 ].u = otherSeamParam;
382 seg.epgeominfo[ 1 ].u = otherSeamParam;
383 swap (seg.epgeominfo[0].v, seg.epgeominfo[1].v);
385 seg.epgeominfo[ 0 ].v = otherSeamParam;
386 seg.epgeominfo[ 1 ].v = otherSeamParam;
387 swap (seg.epgeominfo[0].u, seg.epgeominfo[1].u);
389 swap (seg.p1, seg.p2);
390 swap (seg.epgeominfo[0].dist, seg.epgeominfo[1].dist);
391 seg.edgenr = ngMesh.GetNSeg() + 1; // segment id
392 ngMesh.AddSegment (seg);
395 } // loop on geomEdge ancestors
398 } // case TopAbs_EDGE
400 case TopAbs_FACE: { // FACE
401 // ----------------------
402 const TopoDS_Face& geomFace = TopoDS::Face( sm->GetSubShape() );
403 helper.SetSubShape( geomFace );
405 // find solids geomFace bounds
406 int solidID1 = 0, solidID2 = 0;
407 const TopTools_ListOfShape& ancestors = _mesh->GetAncestors( geomFace );
408 TopTools_ListIteratorOfListOfShape ancestorIt ( ancestors );
409 for ( ; ancestorIt.More(); ancestorIt.Next() )
411 const TopoDS_Shape & solid = ancestorIt.Value();
412 if ( solid.ShapeType() == TopAbs_SOLID ) {
413 int id = occgeom.somap.FindIndex ( solid );
414 if ( solidID1 && id != solidID1 ) solidID2 = id;
419 _faceDescriptors[ faceID ].first = solidID1;
420 _faceDescriptors[ faceID ].second = solidID2;
422 // add surface elements
423 SMDS_ElemIteratorPtr faces = smDS->GetElements();
424 while ( faces->more() ) {
426 const SMDS_MeshElement* f = faces->next();
427 if ( f->NbNodes() % 3 != 0 ) { // not triangle
428 for ( ancestorIt.Initialize(ancestors); ancestorIt.More(); ancestorIt.Next() )
429 if ( ancestorIt.Value().ShapeType() == TopAbs_SOLID ) {
430 sm = _mesh->GetSubMesh( ancestorIt.Value() );
433 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
434 smError.reset( new SMESH_ComputeError(COMPERR_BAD_INPUT_MESH,"Not triangle submesh"));
435 smError->myBadElements.push_back( f );
439 netgen::Element2d tri(3);
440 tri.SetIndex ( faceID );
442 for ( int i = 0; i < 3; ++i ) {
443 const SMDS_MeshNode* node = f->GetNode( i ), * inFaceNode=0;
444 if ( helper.IsSeamShape( node->GetPosition()->GetShapeId() ))
445 if ( helper.IsSeamShape( f->GetNode( i+1 )->GetPosition()->GetShapeId() ))
446 inFaceNode = f->GetNode( i-1 );
448 inFaceNode = f->GetNode( i+1 );
450 gp_XY uv = helper.GetNodeUV( geomFace, node, inFaceNode );
451 tri.GeomInfoPi(i+1).u = uv.X();
452 tri.GeomInfoPi(i+1).v = uv.Y();
453 tri.PNum(i+1) = ngNodeId( node, ngMesh, nodeNgIdMap );
456 ngMesh.AddSurfaceElement (tri);
462 case TopAbs_VERTEX: { // VERTEX
463 // --------------------------
464 SMDS_NodeIteratorPtr nodeIt = smDS->GetNodes();
465 if ( nodeIt->more() )
466 ngNodeId( nodeIt->next(), ngMesh, nodeNgIdMap );
471 } // loop on submeshes
474 nodeVec.resize( ngMesh.GetNP() + 1 );
475 map< const SMDS_MeshNode*, int >::iterator node_NgId, nodeNgIdEnd = nodeNgIdMap.end();
476 for ( node_NgId = nodeNgIdMap.begin(); node_NgId != nodeNgIdEnd; ++node_NgId)
477 nodeVec[ node_NgId->second ] = (SMDS_MeshNode*) node_NgId->first;
482 //=============================================================================
484 * Here we are going to use the NETGEN mesher
486 //=============================================================================
487 bool NETGENPlugin_Mesher::Compute()
490 netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters();
492 netgen::MeshingParameters& mparams = netgen::mparam;
494 MESSAGE("Compute with:\n"
495 " max size = " << mparams.maxh << "\n"
496 " segments per edge = " << mparams.segmentsperedge);
498 " growth rate = " << mparams.grading << "\n"
499 " elements per radius = " << mparams.curvaturesafety << "\n"
500 " second order = " << mparams.secondorder << "\n"
501 " quad allowed = " << mparams.quad);
503 SMESH_ComputeErrorPtr error = SMESH_ComputeError::New();
506 // -------------------------
507 // Prepare OCC geometry
508 // -------------------------
510 netgen::OCCGeometry occgeo;
511 list< SMESH_subMesh* > meshedSM;
512 PrepareOCCgeometry( occgeo, _shape, *_mesh, &meshedSM );
514 // -------------------------
516 // -------------------------
518 netgen::Mesh *ngMesh = NULL;
520 SMESH_Comment comment;
525 // vector of nodes in which node index == netgen ID
526 vector< SMDS_MeshNode* > nodeVec;
532 // pass 1D simple parameters to NETGEN
534 if ( int nbSeg = _simpleHyp->GetNumberOfSegments() ) {
536 mparams.segmentsperedge = nbSeg + 0.1;
537 mparams.maxh = occgeo.boundingbox.Diam();
542 mparams.segmentsperedge = 1;
543 mparams.maxh = _simpleHyp->GetLocalLength();
546 // let netgen create ngMesh and calculate element size on not meshed shapes
548 int startWith = netgen::MESHCONST_ANALYSE;
549 int endWith = netgen::MESHCONST_ANALYSE;
550 err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
551 if (err) comment << "Error in netgen::OCCGenerateMesh() at MESHCONST_ANALYSE step";
553 // fill ngMesh with nodes and elements of computed submeshes
554 err = ! fillNgMesh(occgeo, *ngMesh, nodeVec, meshedSM);
555 nbInitNod = ngMesh->GetNP();
556 nbInitSeg = ngMesh->GetNSeg();
557 nbInitFac = ngMesh->GetNSE();
562 startWith = endWith = netgen::MESHCONST_MESHEDGES;
563 err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
564 if (err) comment << "Error in netgen::OCCGenerateMesh() at 1D mesh generation";
566 // ---------------------
567 // compute surface mesh
568 // ---------------------
571 // pass 2D simple parameters to NETGEN
573 if ( double area = _simpleHyp->GetMaxElementArea() ) {
575 mparams.maxh = sqrt(2. * area/sqrt(3.0));
576 mparams.grading = 0.4; // moderate size growth
581 for ( TopExp_Explorer exp( _shape, TopAbs_EDGE ); exp.More(); exp.Next() )
582 length += SMESH_Algo::EdgeLength( TopoDS::Edge( exp.Current() ));
583 if ( ngMesh->GetNSeg() )
584 mparams.maxh = length / ngMesh->GetNSeg();
587 mparams.grading = 0.2; // slow size growth
589 mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
590 ngMesh->SetGlobalH (mparams.maxh);
591 netgen::Box<3> bb = occgeo.GetBoundingBox();
592 bb.Increase (bb.Diam()/20);
593 ngMesh->SetLocalH (bb.PMin(), bb.PMax(), mparams.grading);
595 // let netgen compute 2D mesh
596 startWith = netgen::MESHCONST_MESHSURFACE;
597 endWith = _optimize ? netgen::MESHCONST_OPTSURFACE : netgen::MESHCONST_MESHSURFACE;
598 err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
599 if (err) comment << "Error in netgen::OCCGenerateMesh() at surface mesh generation";
601 // ---------------------
602 // generate volume mesh
603 // ---------------------
604 if (!err && _isVolume)
606 // add ng face descriptors of meshed faces
607 std::map< int, std::pair<int,int> >::iterator fId_soIds = _faceDescriptors.begin();
608 for ( ; fId_soIds != _faceDescriptors.end(); ++fId_soIds ) {
609 int faceID = fId_soIds->first;
610 int solidID1 = fId_soIds->second.first;
611 int solidID2 = fId_soIds->second.second;
612 ngMesh->AddFaceDescriptor (netgen::FaceDescriptor(faceID, solidID1, solidID2, 0));
614 // pass 3D simple parameters to NETGEN
615 const NETGENPlugin_SimpleHypothesis_3D* simple3d =
616 dynamic_cast< const NETGENPlugin_SimpleHypothesis_3D* > ( _simpleHyp );
618 if ( double vol = simple3d->GetMaxElementVolume() ) {
620 mparams.maxh = pow( 72 * vol * vol, 1/3. );
621 mparams.maxh = min( mparams.maxh, occgeo.boundingbox.Diam()/2 );
625 mparams.maxh = ngMesh->AverageH();
627 // netgen::ARRAY<double> maxhdom;
628 // maxhdom.SetSize (occgeo.NrSolids());
629 // maxhdom = mparams.maxh;
630 // ngMesh->SetMaxHDomain (maxhdom);
631 ngMesh->SetGlobalH (mparams.maxh);
632 mparams.grading = 0.4;
633 ngMesh->CalcLocalH();
635 // let netgen compute 3D mesh
636 startWith = netgen::MESHCONST_MESHVOLUME;
637 endWith = _optimize ? netgen::MESHCONST_OPTVOLUME : netgen::MESHCONST_MESHVOLUME;
638 err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
639 if (err) comment << "Error in netgen::OCCGenerateMesh()";
641 if (!err && mparams.secondorder > 0)
643 netgen::OCCRefinementSurfaces ref (occgeo);
644 ref.MakeSecondOrder (*ngMesh);
647 catch (netgen::NgException exc)
649 error->myName = err = COMPERR_ALGO_FAILED;
650 comment << exc.What();
653 int nbNod = ngMesh->GetNP();
654 int nbSeg = ngMesh->GetNSeg();
655 int nbFac = ngMesh->GetNSE();
656 int nbVol = ngMesh->GetNE();
658 MESSAGE((err ? "Mesh Generation failure" : "End of Mesh Generation") <<
659 ", nb nodes: " << nbNod <<
660 ", nb segments: " << nbSeg <<
661 ", nb faces: " << nbFac <<
662 ", nb volumes: " << nbVol);
664 // -----------------------------------------------------------
665 // Feed back the SMESHDS with the generated Nodes and Elements
666 // -----------------------------------------------------------
668 SMESHDS_Mesh* meshDS = _mesh->GetMeshDS();
669 bool isOK = ( !err && (_isVolume ? (nbVol > 0) : (nbFac > 0)) );
670 if ( true /*isOK*/ ) // get whatever built
672 // map of nodes assigned to submeshes
673 NCollection_Map<int> pindMap;
674 // create and insert nodes into nodeVec
675 nodeVec.resize( nbNod + 1 );
677 for (i = nbInitNod+1; i <= nbNod /*&& isOK*/; ++i )
679 const netgen::MeshPoint& ngPoint = ngMesh->Point(i);
680 SMDS_MeshNode* node = NULL;
681 bool newNodeOnVertex = false;
683 if (i-nbInitNod <= occgeo.vmap.Extent())
686 aVert = TopoDS::Vertex(occgeo.vmap(i-nbInitNod));
687 SMESHDS_SubMesh * submesh = meshDS->MeshElements(aVert);
690 SMDS_NodeIteratorPtr it = submesh->GetNodes();
693 node = const_cast<SMDS_MeshNode*> (it->next());
698 newNodeOnVertex = true;
701 node = meshDS->AddNode(ngPoint.X(), ngPoint.Y(), ngPoint.Z());
704 MESSAGE("Cannot create a mesh node");
705 if ( !comment.size() ) comment << "Cannot create a mesh node";
706 nbSeg = nbFac = nbVol = isOK = 0;
709 nodeVec.at(i) = node;
713 meshDS->SetNodeOnVertex(node, aVert);
718 // create mesh segments along geometric edges
719 NCollection_Map<Link> linkMap;
720 for (i = nbInitSeg+1; i <= nbSeg/* && isOK*/; ++i )
722 const netgen::Segment& seg = ngMesh->LineSegment(i);
723 Link link(seg.p1, seg.p2);
724 if (linkMap.Contains(link))
728 int pinds[3] = { seg.p1, seg.p2, seg.pmid };
731 for (int j=0; j < 3; ++j)
734 if (pind <= 0) continue;
741 int aGeomEdgeInd = seg.epgeominfo[j].edgenr;
742 if (aGeomEdgeInd > 0 && aGeomEdgeInd <= occgeo.emap.Extent())
743 aEdge = TopoDS::Edge(occgeo.emap(aGeomEdgeInd));
745 param = seg.epgeominfo[j].dist;
749 param = param2 * 0.5;
750 if (pind <= nbInitNod || pindMap.Contains(pind))
754 meshDS->SetNodeOnEdge(nodeVec.at(pind), aEdge, param);
759 if (nbp < 3) // second order ?
760 edge = meshDS->AddEdge(nodeVec.at(pinds[0]), nodeVec.at(pinds[1]));
762 edge = meshDS->AddEdge(nodeVec.at(pinds[0]), nodeVec.at(pinds[1]),
763 nodeVec.at(pinds[2]));
766 if ( !comment.size() ) comment << "Cannot create a mesh edge";
767 MESSAGE("Cannot create a mesh edge");
768 nbSeg = nbFac = nbVol = isOK = 0;
772 meshDS->SetMeshElementOnShape(edge, aEdge);
775 // create mesh faces along geometric faces
776 for (i = nbInitFac+1; i <= nbFac/* && isOK*/; ++i )
778 const netgen::Element2d& elem = ngMesh->SurfaceElement(i);
779 int aGeomFaceInd = elem.GetIndex();
781 if (aGeomFaceInd > 0 && aGeomFaceInd <= occgeo.fmap.Extent())
782 aFace = TopoDS::Face(occgeo.fmap(aGeomFaceInd));
783 vector<SMDS_MeshNode*> nodes;
784 for (int j=1; j <= elem.GetNP(); ++j)
786 int pind = elem.PNum(j);
787 SMDS_MeshNode* node = nodeVec.at(pind);
788 nodes.push_back(node);
789 if (pind <= nbInitNod || pindMap.Contains(pind))
793 const netgen::PointGeomInfo& pgi = elem.GeomInfoPi(j);
794 meshDS->SetNodeOnFace(node, aFace, pgi.u, pgi.v);
798 SMDS_MeshFace* face = NULL;
799 switch (elem.GetType())
802 face = meshDS->AddFace(nodes[0],nodes[1],nodes[2]);
805 face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3]);
808 face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[5],nodes[3],nodes[4]);
811 face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3],
812 nodes[4],nodes[7],nodes[5],nodes[6]);
815 MESSAGE("NETGEN created a face of unexpected type, ignoring");
820 if ( !comment.size() ) comment << "Cannot create a mesh face";
821 MESSAGE("Cannot create a mesh face");
822 nbSeg = nbFac = nbVol = isOK = 0;
826 meshDS->SetMeshElementOnShape(face, aFace);
830 for (i = 1; i <= nbVol/* && isOK*/; ++i)
832 const netgen::Element& elem = ngMesh->VolumeElement(i);
833 int aSolidInd = elem.GetIndex();
835 if (aSolidInd > 0 && aSolidInd <= occgeo.somap.Extent())
836 aSolid = TopoDS::Solid(occgeo.somap(aSolidInd));
837 vector<SMDS_MeshNode*> nodes;
838 for (int j=1; j <= elem.GetNP(); ++j)
840 int pind = elem.PNum(j);
841 SMDS_MeshNode* node = nodeVec.at(pind);
842 nodes.push_back(node);
843 if (pind <= nbInitNod || pindMap.Contains(pind))
845 if (!aSolid.IsNull())
848 meshDS->SetNodeInVolume(node, aSolid);
852 SMDS_MeshVolume* vol = NULL;
853 switch (elem.GetType())
856 vol = meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3]);
859 vol = meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3],
860 nodes[4],nodes[7],nodes[5],nodes[6],nodes[8],nodes[9]);
863 MESSAGE("NETGEN created a volume of unexpected type, ignoring");
868 if ( !comment.size() ) comment << "Cannot create a mesh volume";
869 MESSAGE("Cannot create a mesh volume");
870 nbSeg = nbFac = nbVol = isOK = 0;
873 if (!aSolid.IsNull())
874 meshDS->SetMeshElementOnShape(vol, aSolid);
878 if ( error->IsOK() && ( !isOK || comment.size() > 0 ))
879 error->myName = COMPERR_ALGO_FAILED;
880 if ( !comment.empty() )
881 error->myComment = comment;
883 // set bad compute error to subshapes of all failed subshapes shapes
884 if ( !error->IsOK() && err )
886 for (int i = 1; i <= occgeo.fmap.Extent(); i++) {
887 int status = occgeo.facemeshstatus[i-1];
888 if (status == 1 ) continue;
889 if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( occgeo.fmap( i ))) {
890 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
891 if ( !smError || smError->IsOK() ) {
893 smError.reset( new SMESH_ComputeError( error->myName, error->myComment ));
895 smError.reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED, "Ignored" ));
901 nglib::Ng_DeleteMesh((nglib::Ng_Mesh*)ngMesh);
906 return error->IsOK();
909 //================================================================================
911 * \brief Remove "test.out" and "problemfaces" files in current directory
913 //================================================================================
915 void NETGENPlugin_Mesher::RemoveTmpFiles()
917 TCollection_AsciiString str("test.out");
918 OSD_Path path1( str );
919 OSD_File file1( path1 );
921 str = "problemfaces";
922 OSD_Path path2( str );
923 OSD_File file2( path2 );