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.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
22 // File : NETGENPlugin_Mesher.cxx
23 // Author : Michael Sazonov (OCN)
27 //=============================================================================
30 #include "NETGENPlugin_Mesher.hxx"
31 #include "NETGENPlugin_Hypothesis_2D.hxx"
33 #include <SMESHDS_Mesh.hxx>
34 #include <SMDS_MeshElement.hxx>
35 #include <SMDS_MeshNode.hxx>
36 #include <utilities.h>
40 #include <BRep_Tool.hxx>
42 #include <TopExp_Explorer.hxx>
44 #include <NCollection_Map.hxx>
46 // Netgen include files
51 #include <occgeom.hpp>
52 #include <meshing.hpp>
53 //#include <ngexception.hpp>
55 extern int OCCGenerateMesh (OCCGeometry&, Mesh*&, int, int, char*);
56 extern MeshingParameters mparam;
59 //=============================================================================
63 //=============================================================================
65 NETGENPlugin_Mesher::NETGENPlugin_Mesher (SMESHDS_Mesh* meshDS,
66 const TopoDS_Shape& aShape,
73 // Initialize global NETGEN parameters by default values:
74 // maximal mesh edge size
75 netgen::mparam.maxh = NETGENPlugin_Hypothesis::GetDefaultMaxSize();
76 // minimal number of segments per edge
77 netgen::mparam.segmentsperedge = NETGENPlugin_Hypothesis::GetDefaultNbSegPerEdge();
78 // rate of growth of size between elements
79 netgen::mparam.grading = NETGENPlugin_Hypothesis::GetDefaultGrowthRate();
80 // safety factor for curvatures (elements per radius)
81 netgen::mparam.curvaturesafety = NETGENPlugin_Hypothesis::GetDefaultNbSegPerRadius();
82 // create elements of second order
83 netgen::mparam.secondorder = NETGENPlugin_Hypothesis::GetDefaultSecondOrder() ? 1 : 0;
84 // quad-dominated surface meshing
86 netgen::mparam.quad = 0;
88 netgen::mparam.quad = NETGENPlugin_Hypothesis_2D::GetDefaultQuadAllowed() ? 1 : 0;
91 //=============================================================================
93 * Pass parameters to NETGEN
95 //=============================================================================
96 void NETGENPlugin_Mesher::SetParameters(const NETGENPlugin_Hypothesis* hyp)
100 // Initialize global NETGEN parameters:
101 // maximal mesh segment size
102 netgen::mparam.maxh = hyp->GetMaxSize();
103 // minimal number of segments per edge
104 netgen::mparam.segmentsperedge = hyp->GetNbSegPerEdge();
105 // rate of growth of size between elements
106 netgen::mparam.grading = hyp->GetGrowthRate();
107 // safety factor for curvatures (elements per radius)
108 netgen::mparam.curvaturesafety = hyp->GetNbSegPerRadius();
109 // create elements of second order
110 netgen::mparam.secondorder = hyp->GetSecondOrder() ? 1 : 0;
111 // quad-dominated surface meshing
112 // only triangles are allowed for volumic mesh
114 netgen::mparam.quad = static_cast<const NETGENPlugin_Hypothesis_2D*>
115 (hyp)->GetQuadAllowed() ? 1 : 0;
116 _optimize = hyp->GetOptimize();
120 //=============================================================================
122 * Link - a pair of integer numbers
124 //=============================================================================
128 Link(int _n1, int _n2) : n1(_n1), n2(_n2) {}
129 Link() : n1(0), n2(0) {}
132 int HashCode(const Link& aLink, int aLimit)
134 return HashCode(aLink.n1 + aLink.n2, aLimit);
137 Standard_Boolean IsEqual(const Link& aLink1, const Link& aLink2)
139 return (aLink1.n1 == aLink2.n1 && aLink1.n2 == aLink2.n2 ||
140 aLink1.n1 == aLink2.n2 && aLink1.n2 == aLink2.n1);
143 //=============================================================================
145 * Here we are going to use the NETGEN mesher
147 //=============================================================================
148 bool NETGENPlugin_Mesher::Compute()
150 MESSAGE("Compute with:\n"
151 " max size = " << netgen::mparam.maxh << "\n"
152 " segments per edge = " << netgen::mparam.segmentsperedge);
154 " growth rate = " << netgen::mparam.grading << "\n"
155 " elements per radius = " << netgen::mparam.curvaturesafety << "\n"
156 " second order = " << netgen::mparam.secondorder << "\n"
157 " quad allowed = " << netgen::mparam.quad);
161 // -------------------------
162 // Prepare OCC geometry
163 // -------------------------
165 netgen::OCCGeometry occgeo;
166 occgeo.shape = _shape;
169 BRepTools::Clean (_shape);
170 BRepMesh_IncrementalMesh::BRepMesh_IncrementalMesh (_shape, 0.01, true);
172 BRepBndLib::Add (_shape, bb);
173 double x1,y1,z1,x2,y2,z2;
174 bb.Get (x1,y1,z1,x2,y2,z2);
175 MESSAGE("shape bounding box:\n" <<
176 "(" << x1 << " " << y1 << " " << z1 << ") " <<
177 "(" << x2 << " " << y2 << " " << z2 << ")");
178 netgen::Point<3> p1 = netgen::Point<3> (x1,y1,z1);
179 netgen::Point<3> p2 = netgen::Point<3> (x2,y2,z2);
180 occgeo.boundingbox = netgen::Box<3> (p1,p2);
182 // -------------------------
184 // -------------------------
186 netgen::Mesh *ngMesh = NULL;
187 // we start always with ANALYSE,
188 // but end depending on _optimize and _isVolume
189 int startWith = netgen::MESHCONST_ANALYSE;
190 int endWith = (_optimize
191 ? (_isVolume ? netgen::MESHCONST_OPTVOLUME : netgen::MESHCONST_OPTSURFACE)
192 : netgen::MESHCONST_MESHSURFACE);
198 err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
199 if (!err && !_optimize)
201 // we have got surface mesh only, so generate volume mesh
202 startWith = endWith = netgen::MESHCONST_MESHVOLUME;
203 err = netgen::OCCGenerateMesh(occgeo, ngMesh, startWith, endWith, optstr);
205 if (!err && netgen::mparam.secondorder > 0)
207 netgen::OCCRefinementSurfaces ref (occgeo);
208 ref.MakeSecondOrder (*ngMesh);
211 catch (netgen::NgException exc)
213 MESSAGE ("Exception in NETGEN: " << exc.What());
217 int nbNod = ngMesh->GetNP();
218 int nbSeg = ngMesh->GetNSeg();
219 int nbFac = ngMesh->GetNSE();
220 int nbVol = ngMesh->GetNE();
222 MESSAGE((err ? "Mesh Generation failure" : "End of Mesh Generation") <<
223 ", nb nodes: " << nbNod <<
224 ", nb segments: " << nbSeg <<
225 ", nb faces: " << nbFac <<
226 ", nb volumes: " << nbVol);
228 // -----------------------------------------------------------
229 // Feed back the SMESHDS with the generated Nodes and Elements
230 // -----------------------------------------------------------
232 bool isOK = ( !err && (_isVolume ? (nbVol > 0) : (nbFac > 0)) );
235 // vector of nodes in which node index == netgen ID
236 vector< SMDS_MeshNode* > nodeVec ( nbNod + 1 );
237 // map of nodes assigned to submeshes
238 NCollection_Map<int> pindMap;
239 // create and insert nodes into nodeVec
241 for (i = 1; i <= nbNod && isOK; ++i )
243 const netgen::MeshPoint& ngPoint = ngMesh->Point(i);
244 SMDS_MeshNode* node = NULL;
245 bool newNodeOnVertex = false;
247 if (i <= occgeo.vmap.Extent())
250 aVert = TopoDS::Vertex(occgeo.vmap(i));
251 SMESHDS_SubMesh * submesh = _meshDS->MeshElements(aVert);
254 SMDS_NodeIteratorPtr it = submesh->GetNodes();
257 node = const_cast<SMDS_MeshNode*> (it->next());
262 newNodeOnVertex = true;
265 node = _meshDS->AddNode(ngPoint.X(), ngPoint.Y(), ngPoint.Z());
268 MESSAGE("Cannot create a mesh node");
272 nodeVec.at(i) = node;
276 _meshDS->SetNodeOnVertex(node, aVert);
281 // create mesh segments along geometric edges
282 NCollection_Map<Link> linkMap;
283 for (i = 1; i <= nbSeg && isOK; ++i )
285 const netgen::Segment& seg = ngMesh->LineSegment(i);
286 Link link(seg.p1, seg.p2);
287 if (linkMap.Contains(link))
291 int pinds[3] = { seg.p1, seg.p2, seg.pmid };
294 for (int j=0; j < 3; ++j)
297 if (pind <= 0) continue;
304 int aGeomEdgeInd = seg.epgeominfo[j].edgenr;
305 if (aGeomEdgeInd > 0 && aGeomEdgeInd <= occgeo.emap.Extent())
306 aEdge = TopoDS::Edge(occgeo.emap(aGeomEdgeInd));
308 param = seg.epgeominfo[j].dist;
312 param = param2 * 0.5;
313 if (pindMap.Contains(pind))
317 _meshDS->SetNodeOnEdge(nodeVec.at(pind), aEdge, param);
322 if (nbp < 3) // second order ?
323 edge = _meshDS->AddEdge(nodeVec.at(pinds[0]), nodeVec.at(pinds[1]));
325 edge = _meshDS->AddEdge(nodeVec.at(pinds[0]), nodeVec.at(pinds[1]),
326 nodeVec.at(pinds[2]));
329 MESSAGE("Cannot create a mesh edge");
334 _meshDS->SetMeshElementOnShape(edge, aEdge);
337 // create mesh faces along geometric faces
338 for (i = 1; i <= nbFac && isOK; ++i )
340 const netgen::Element2d& elem = ngMesh->SurfaceElement(i);
341 int aGeomFaceInd = elem.GetIndex();
343 if (aGeomFaceInd > 0 && aGeomFaceInd <= occgeo.fmap.Extent())
344 aFace = TopoDS::Face(occgeo.fmap(aGeomFaceInd));
345 vector<SMDS_MeshNode*> nodes;
346 for (int j=1; j <= elem.GetNP(); ++j)
348 int pind = elem.PNum(j);
349 SMDS_MeshNode* node = nodeVec.at(pind);
350 nodes.push_back(node);
351 if (pindMap.Contains(pind))
355 const netgen::PointGeomInfo& pgi = elem.GeomInfoPi(j);
356 _meshDS->SetNodeOnFace(node, aFace, pgi.u, pgi.v);
360 SMDS_MeshFace* face = NULL;
361 switch (elem.GetType())
364 face = _meshDS->AddFace(nodes[0],nodes[1],nodes[2]);
367 face = _meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3]);
370 face = _meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[5],nodes[3],nodes[4]);
373 face = _meshDS->AddFace(nodes[0],nodes[1],nodes[2],nodes[3],
374 nodes[4],nodes[7],nodes[5],nodes[6]);
377 MESSAGE("NETGEN created a face of unexpected type, ignoring");
382 MESSAGE("Cannot create a mesh face");
387 _meshDS->SetMeshElementOnShape(face, aFace);
391 for (i = 1; i <= nbVol && isOK; ++i)
393 const netgen::Element& elem = ngMesh->VolumeElement(i);
394 int aSolidInd = elem.GetIndex();
396 if (aSolidInd > 0 && aSolidInd <= occgeo.somap.Extent())
397 aSolid = TopoDS::Solid(occgeo.somap(aSolidInd));
398 vector<SMDS_MeshNode*> nodes;
399 for (int j=1; j <= elem.GetNP(); ++j)
401 int pind = elem.PNum(j);
402 SMDS_MeshNode* node = nodeVec.at(pind);
403 nodes.push_back(node);
404 if (pindMap.Contains(pind))
406 if (!aSolid.IsNull())
409 _meshDS->SetNodeInVolume(node, aSolid);
413 SMDS_MeshVolume* vol = NULL;
414 switch (elem.GetType())
417 vol = _meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3]);
420 vol = _meshDS->AddVolume(nodes[0],nodes[1],nodes[2],nodes[3],
421 nodes[4],nodes[7],nodes[5],nodes[6],nodes[8],nodes[9]);
424 MESSAGE("NETGEN created a volume of unexpected type, ignoring");
429 MESSAGE("Cannot create a mesh volume");
433 if (!aSolid.IsNull())
434 _meshDS->SetMeshElementOnShape(vol, aSolid);
438 nglib::Ng_DeleteMesh((nglib::Ng_Mesh*)ngMesh);