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"
32 #include <SMESH_Mesh.hxx>
33 #include <SMESH_Comment.hxx>
34 #include <SMESH_ComputeError.hxx>
35 #include <SMESH_subMesh.hxx>
36 #include <SMESHDS_Mesh.hxx>
37 #include <SMDS_MeshElement.hxx>
38 #include <SMDS_MeshNode.hxx>
39 #include <utilities.h>
43 #include <BRep_Tool.hxx>
45 #include <TopExp_Explorer.hxx>
47 #include <NCollection_Map.hxx>
48 #include <OSD_Path.hxx>
49 #include <OSD_File.hxx>
50 #include <TCollection_AsciiString.hxx>
52 // Netgen include files
57 #include <occgeom.hpp>
58 #include <meshing.hpp>
59 //#include <ngexception.hpp>
61 extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, int, int, char*);
62 extern MeshingParameters mparam;
67 //=============================================================================
71 //=============================================================================
73 NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESH_Mesh* mesh,
74 const TopoDS_Shape& aShape,
82 netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters();
84 netgen::MeshingParameters& mparams = netgen::mparam;
86 // Initialize global NETGEN parameters by default values:
87 // maximal mesh edge size
88 mparams.maxh = NETGENPlugin_Hypothesis::GetDefaultMaxSize();
89 // minimal number of segments per edge
90 mparams.segmentsperedge = NETGENPlugin_Hypothesis::GetDefaultNbSegPerEdge();
91 // rate of growth of size between elements
92 mparams.grading = NETGENPlugin_Hypothesis::GetDefaultGrowthRate();
93 // safety factor for curvatures (elements per radius)
94 mparams.curvaturesafety = NETGENPlugin_Hypothesis::GetDefaultNbSegPerRadius();
95 // create elements of second order
96 mparams.secondorder = NETGENPlugin_Hypothesis::GetDefaultSecondOrder() ? 1 : 0;
97 // quad-dominated surface meshing
101 mparams.quad = NETGENPlugin_Hypothesis_2D::GetDefaultQuadAllowed() ? 1 : 0;
104 //=============================================================================
106 * Pass parameters to NETGEN
108 //=============================================================================
109 void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp)
114 netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters();
116 netgen::MeshingParameters& mparams = netgen::mparam;
118 // Initialize global NETGEN parameters:
119 // maximal mesh segment size
120 mparams.maxh = hyp->GetMaxSize();
121 // minimal number of segments per edge
122 mparams.segmentsperedge = hyp->GetNbSegPerEdge();
123 // rate of growth of size between elements
124 mparams.grading = hyp->GetGrowthRate();
125 // safety factor for curvatures (elements per radius)
126 mparams.curvaturesafety = hyp->GetNbSegPerRadius();
127 // create elements of second order
128 mparams.secondorder = hyp->GetSecondOrder() ? 1 : 0;
129 // quad-dominated surface meshing
130 // only triangles are allowed for volumic mesh
132 mparams.quad = static_cast<const NETGENPlugin_Hypothesis_2D*>
133 (hyp)->GetQuadAllowed() ? 1 : 0;
134 _optimize = hyp->GetOptimize();
138 //=============================================================================
140 * Link - a pair of integer numbers
142 //=============================================================================
146 Link(int _n1, int _n2) : n1(_n1), n2(_n2) {}
147 Link() : n1(0), n2(0) {}
150 int HashCode(const Link& aLink, int aLimit)
152 return HashCode(aLink.n1 + aLink.n2, aLimit);
155 Standard_Boolean IsEqual(const Link& aLink1, const Link& aLink2)
157 return (aLink1.n1 == aLink2.n1 && aLink1.n2 == aLink2.n2 ||
158 aLink1.n1 == aLink2.n2 && aLink1.n2 == aLink2.n1);
161 //================================================================================
163 * \brief Initialize netgen::OCCGeometry with OCCT shape
165 //================================================================================
167 void NETGENPlugin_Mesher::PrepareOCCgeometry(netgen::OCCGeometry& occgeo,
168 const TopoDS_Shape& shape)
170 occgeo.shape = shape;
174 BRepTools::Clean (shape);
175 BRepMesh_IncrementalMesh::BRepMesh_IncrementalMesh (shape, 0.01, true);
177 BRepBndLib::Add (shape, bb);
178 double x1,y1,z1,x2,y2,z2;
179 bb.Get (x1,y1,z1,x2,y2,z2);
180 MESSAGE("shape bounding box:\n" <<
181 "(" << x1 << " " << y1 << " " << z1 << ") " <<
182 "(" << x2 << " " << y2 << " " << z2 << ")");
183 netgen::Point<3> p1 = netgen::Point<3> (x1,y1,z1);
184 netgen::Point<3> p2 = netgen::Point<3> (x2,y2,z2);
185 occgeo.boundingbox = netgen::Box<3> (p1,p2);
188 //=============================================================================
190 * Here we are going to use the NETGEN mesher
192 //=============================================================================
193 bool NETGENPlugin_Mesher::Compute()
196 netgen::MeshingParameters& mparams = netgen::GlobalMeshingParameters();
198 netgen::MeshingParameters& mparams = netgen::mparam;
200 MESSAGE("Compute with:\n"
201 " max size = " << mparams.maxh << "\n"
202 " segments per edge = " << mparams.segmentsperedge);
204 " growth rate = " << mparams.grading << "\n"
205 " elements per radius = " << mparams.curvaturesafety << "\n"
206 " second order = " << mparams.secondorder << "\n"
207 " quad allowed = " << mparams.quad);
209 SMESH_ComputeErrorPtr error = SMESH_ComputeError::New();
212 // -------------------------
213 // Prepare OCC geometry
214 // -------------------------
216 netgen::OCCGeometry occgeo;
217 PrepareOCCgeometry( occgeo, _shape );
219 // -------------------------
221 // -------------------------
223 netgen::Mesh *ngMesh = NULL;
224 // we start always with ANALYSE,
225 // but end depending on _optimize and _isVolume
226 int startWith = netgen::MESHCONST_ANALYSE;
227 int endWith = (_optimize
228 ? (_isVolume ? netgen::MESHCONST_OPTVOLUME : netgen::MESHCONST_OPTSURFACE)
229 : netgen::MESHCONST_MESHSURFACE);
233 SMESH_Comment comment;
236 err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
237 if (err) comment << "Error in netgen::OCCGenerateMesh()";
238 if (!err && !_optimize)
240 // we have got surface mesh only, so generate volume mesh
241 startWith = endWith = netgen::MESHCONST_MESHVOLUME;
242 err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
243 if (err) comment << "Error in netgen::OCCGenerateMesh()";
245 if (!err && mparams.secondorder > 0)
247 netgen::OCCRefinementSurfaces ref (occgeo);
248 ref.MakeSecondOrder (*ngMesh);
251 catch (netgen::NgException exc)
253 error->myName = err = COMPERR_ALGO_FAILED;
254 comment << exc.What();
257 int nbNod = ngMesh->GetNP();
258 int nbSeg = ngMesh->GetNSeg();
259 int nbFac = ngMesh->GetNSE();
260 int nbVol = ngMesh->GetNE();
262 MESSAGE((err ? "Mesh Generation failure" : "End of Mesh Generation") <<
263 ", nb nodes: " << nbNod <<
264 ", nb segments: " << nbSeg <<
265 ", nb faces: " << nbFac <<
266 ", nb volumes: " << nbVol);
268 // -----------------------------------------------------------
269 // Feed back the SMESHDS with the generated Nodes and Elements
270 // -----------------------------------------------------------
272 SMESHDS_Mesh* meshDS = _mesh->GetMeshDS();
273 bool isOK = ( !err && (_isVolume ? (nbVol > 0) : (nbFac > 0)) );
274 if ( true /*isOK*/ ) // get whatever built
276 // vector of nodes in which node index == netgen ID
277 vector< SMDS_MeshNode* > nodeVec ( nbNod + 1 );
278 // map of nodes assigned to submeshes
279 NCollection_Map<int> pindMap;
280 // create and insert nodes into nodeVec
282 for (i = 1; i <= nbNod /*&& isOK*/; ++i )
284 const netgen::MeshPoint& ngPoint = ngMesh->Point(i);
285 SMDS_MeshNode* node = NULL;
286 bool newNodeOnVertex = false;
288 if (i <= occgeo.vmap.Extent())
291 aVert = TopoDS::Vertex(occgeo.vmap(i));
292 SMESHDS_SubMesh * submesh = meshDS->MeshElements(aVert);
295 SMDS_NodeIteratorPtr it = submesh->GetNodes();
298 node = const_cast<SMDS_MeshNode*> (it->next());
303 newNodeOnVertex = true;
306 node = meshDS->AddNode(ngPoint.X(), ngPoint.Y(), ngPoint.Z());
309 MESSAGE("Cannot create a mesh node");
310 if ( !comment.size() ) comment << "Cannot create a mesh node";
311 nbSeg = nbFac = nbVol = isOK = 0;
314 nodeVec.at(i) = node;
318 meshDS->SetNodeOnVertex(node, aVert);
323 // create mesh segments along geometric edges
324 NCollection_Map<Link> linkMap;
325 for (i = 1; i <= nbSeg/* && isOK*/; ++i )
327 const netgen::Segment& seg = ngMesh->LineSegment(i);
328 Link link(seg.p1, seg.p2);
329 if (linkMap.Contains(link))
333 int pinds[3] = { seg.p1, seg.p2, seg.pmid };
336 for (int j=0; j < 3; ++j)
339 if (pind <= 0) continue;
346 int aGeomEdgeInd = seg.epgeominfo[j].edgenr;
347 if (aGeomEdgeInd > 0 && aGeomEdgeInd <= occgeo.emap.Extent())
348 aEdge = TopoDS::Edge(occgeo.emap(aGeomEdgeInd));
350 param = seg.epgeominfo[j].dist;
354 param = param2 * 0.5;
355 if (pindMap.Contains(pind))
359 meshDS->SetNodeOnEdge(nodeVec.at(pind), aEdge, param);
364 if (nbp < 3) // second order ?
365 edge = meshDS->AddEdge(nodeVec.at(pinds[0]), nodeVec.at(pinds[1]));
367 edge = meshDS->AddEdge(nodeVec.at(pinds[0]), nodeVec.at(pinds[1]),
368 nodeVec.at(pinds[2]));
371 if ( !comment.size() ) comment << "Cannot create a mesh edge";
372 MESSAGE("Cannot create a mesh edge");
373 nbSeg = nbFac = nbVol = isOK = 0;
377 meshDS->SetMeshElementOnShape(edge, aEdge);
380 // create mesh faces along geometric faces
381 for (i = 1; i <= nbFac/* && isOK*/; ++i )
383 const netgen::Element2d& elem = ngMesh->SurfaceElement(i);
384 int aGeomFaceInd = elem.GetIndex();
386 if (aGeomFaceInd > 0 && aGeomFaceInd <= occgeo.fmap.Extent())
387 aFace = TopoDS::Face(occgeo.fmap(aGeomFaceInd));
388 vector<SMDS_MeshNode*> nodes;
389 for (int j=1; j <= elem.GetNP(); ++j)
391 int pind = elem.PNum(j);
392 SMDS_MeshNode* node = nodeVec.at(pind);
393 nodes.push_back(node);
394 if (pindMap.Contains(pind))
398 const netgen::PointGeomInfo& pgi = elem.GeomInfoPi(j);
399 meshDS->SetNodeOnFace(node, aFace, pgi.u, pgi.v);
403 SMDS_MeshFace* face = NULL;
404 switch (elem.GetType())
407 face = meshDS->AddFace(nodes[0],nodes[1],nodes[2]);
410 face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3]);
413 face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[5],nodes[3],nodes[4]);
416 face = meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3],
417 nodes[4],nodes[7],nodes[5],nodes[6]);
420 MESSAGE("NETGEN created a face of unexpected type, ignoring");
425 if ( !comment.size() ) comment << "Cannot create a mesh face";
426 MESSAGE("Cannot create a mesh face");
427 nbSeg = nbFac = nbVol = isOK = 0;
431 meshDS->SetMeshElementOnShape(face, aFace);
435 for (i = 1; i <= nbVol/* && isOK*/; ++i)
437 const netgen::Element& elem = ngMesh->VolumeElement(i);
438 int aSolidInd = elem.GetIndex();
440 if (aSolidInd > 0 && aSolidInd <= occgeo.somap.Extent())
441 aSolid = TopoDS::Solid(occgeo.somap(aSolidInd));
442 vector<SMDS_MeshNode*> nodes;
443 for (int j=1; j <= elem.GetNP(); ++j)
445 int pind = elem.PNum(j);
446 SMDS_MeshNode* node = nodeVec.at(pind);
447 nodes.push_back(node);
448 if (pindMap.Contains(pind))
450 if (!aSolid.IsNull())
453 meshDS->SetNodeInVolume(node, aSolid);
457 SMDS_MeshVolume* vol = NULL;
458 switch (elem.GetType())
461 vol = meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3]);
464 vol = meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3],
465 nodes[4],nodes[7],nodes[5],nodes[6],nodes[8],nodes[9]);
468 MESSAGE("NETGEN created a volume of unexpected type, ignoring");
473 if ( !comment.size() ) comment << "Cannot create a mesh volume";
474 MESSAGE("Cannot create a mesh volume");
475 nbSeg = nbFac = nbVol = isOK = 0;
478 if (!aSolid.IsNull())
479 meshDS->SetMeshElementOnShape(vol, aSolid);
483 if ( error->IsOK() && ( !isOK || comment.size() > 0 ))
484 error->myName = COMPERR_ALGO_FAILED;
485 if ( !comment.empty() )
486 error->myComment = comment;
488 // set bad compute error to subshapes of all failed subshapes shapes
489 if ( !error->IsOK() && err )
491 for (int i = 1; i <= occgeo.fmap.Extent(); i++) {
492 int status = occgeo.facemeshstatus[i-1];
493 if (status == 1 ) continue;
494 if ( SMESH_subMesh* sm = _mesh->GetSubMeshContaining( occgeo.fmap( i ))) {
495 SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
496 if ( !smError || smError->IsOK() ) {
498 smError.reset( new SMESH_ComputeError( error->myName, error->myComment ));
500 smError.reset( new SMESH_ComputeError( COMPERR_ALGO_FAILED, "Ignored" ));
506 nglib::Ng_DeleteMesh((nglib::Ng_Mesh*)ngMesh);
511 return error->IsOK();
514 //================================================================================
516 * \brief Remove "test.out" and "problemfaces" files in current directory
518 //================================================================================
520 void NETGENPlugin_Mesher::RemoveTmpFiles()
522 TCollection_AsciiString str("test.out");
523 OSD_Path path1( str );
524 OSD_File file1( path1 );
526 str = "problemfaces";
527 OSD_Path path2( str );
528 OSD_File file2( path2 );