Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2023 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 : #include "libmesh/libmesh_config.h"
19 : #ifdef LIBMESH_HAVE_NETGEN
20 :
21 :
22 : // C++ includes
23 : #include <sstream>
24 :
25 : // Local includes
26 : #include "libmesh/mesh_netgen_interface.h"
27 :
28 : #include "libmesh/boundary_info.h"
29 : #include "libmesh/cell_tet4.h"
30 : #include "libmesh/face_tri3.h"
31 : #include "libmesh/libmesh_logging.h"
32 : #include "libmesh/mesh_communication.h"
33 : #include "libmesh/unstructured_mesh.h"
34 : #include "libmesh/utility.h" // libmesh_map_find
35 :
36 : namespace nglib {
37 : #include "netgen/nglib/nglib.h"
38 : }
39 :
40 : namespace {
41 :
42 : // RAII for exception safety
43 : class WrappedNgMesh
44 : {
45 : public:
46 60 : WrappedNgMesh() {
47 60 : _ngmesh = nglib::Ng_NewMesh();
48 5 : }
49 :
50 15 : ~WrappedNgMesh() {
51 60 : nglib::Ng_DeleteMesh(_ngmesh);
52 55 : }
53 :
54 : void clear() {
55 : nglib::Ng_DeleteMesh(_ngmesh);
56 : _ngmesh = nglib::Ng_NewMesh();
57 : }
58 :
59 1291 : operator nglib::Ng_Mesh* () {
60 14283 : return _ngmesh;
61 : }
62 :
63 : private:
64 : nglib::Ng_Mesh * _ngmesh;
65 : };
66 :
67 : }
68 :
69 : namespace libMesh
70 : {
71 :
72 : //----------------------------------------------------------------------
73 : // NetGenMeshInterface class members
74 426 : NetGenMeshInterface::NetGenMeshInterface (UnstructuredMesh & mesh) :
75 : MeshTetInterface(mesh),
76 426 : _serializer(mesh)
77 : {
78 426 : }
79 :
80 :
81 :
82 426 : void NetGenMeshInterface::triangulate ()
83 : {
84 : using namespace nglib;
85 :
86 14 : LOG_SCOPE("triangulate()", "NetGenMeshInterface");
87 :
88 : // We're hoping to do volume_to_surface_mesh in parallel at least,
89 : // but then we'll need to serialize any hole meshes to rank 0 so it
90 : // can use them in serial.
91 :
92 : const BoundingBox mesh_bb =
93 426 : MeshTetInterface::volume_to_surface_mesh(this->_mesh);
94 :
95 15 : std::vector<MeshSerializer> hole_serializers;
96 355 : if (_holes)
97 284 : for (std::unique_ptr<UnstructuredMesh> & hole : *_holes)
98 : {
99 : const BoundingBox hole_bb =
100 142 : MeshTetInterface::volume_to_surface_mesh(*hole);
101 :
102 142 : libmesh_error_msg_if
103 : (!mesh_bb.contains(hole_bb),
104 : "Found hole with bounding box " << hole_bb <<
105 : "\nextending outside of mesh bounding box " << mesh_bb);
106 :
107 : hole_serializers.emplace_back
108 276 : (*hole, /* need_serial */ true,
109 146 : /* serial_only_needed_on_proc_0 */ true);
110 : }
111 :
112 : // If we're not rank 0, we're just going to wait for rank 0 to call
113 : // Netgen, then receive its data afterward, we're not going to hope
114 : // that Netgen does the exact same thing on every processor.
115 710 : if (this->_mesh.processor_id() != 0)
116 : {
117 : // We don't need our holes anymore. Delete their serializers
118 : // first to avoid dereferencing dangling pointers.
119 5 : hole_serializers.clear();
120 295 : if (_holes)
121 2 : _holes->clear();
122 :
123 : // Receive the mesh data rank 0 will send later, then fix it up
124 : // together
125 295 : MeshCommunication().broadcast(this->_mesh);
126 295 : this->_mesh.prepare_for_use();
127 10 : return;
128 : }
129 :
130 60 : this->check_hull_integrity();
131 :
132 60 : Ng_Meshing_Parameters params;
133 :
134 : // Override any default parameters we might need to, to avoid
135 : // inserting nodes we don't want.
136 60 : params.uselocalh = false;
137 60 : params.minh = 0;
138 60 : params.elementsperedge = 1;
139 60 : params.elementspercurve = 1;
140 60 : params.closeedgeenable = false;
141 60 : params.closeedgefact = 0;
142 60 : params.minedgelenenable = false;
143 60 : params.minedgelen = 0;
144 :
145 : // Try to get a no-extra-nodes mesh if we're asked to, or try to
146 : // translate our desired volume into NetGen terms otherwise.
147 : //
148 : // Spoiler alert: all we can do is try; NetGen uses a marching front
149 : // algorithm that can insert extra nodes despite all my best
150 : // efforts.
151 60 : if (_desired_volume == 0) // shorthand for "no refinement"
152 : {
153 48 : params.maxh = std::numeric_limits<double>::max();
154 48 : params.fineness = 0; // "coarse" in the docs
155 48 : params.grading = 1; // "aggressive local grading" to avoid smoothing??
156 :
157 : // Turning off optimization steps avoids another opportunity for
158 : // Netgen to try to add more nodes.
159 48 : params.optsteps_3d = 0;
160 : }
161 : else
162 12 : params.maxh = double(std::pow(_desired_volume, 1./3.));
163 :
164 : // Keep track of how NetGen copies of nodes map back to our original
165 : // nodes, so we can connect new elements to nodes correctly.
166 10 : std::unordered_map<int, dof_id_type> ng_to_libmesh_id;
167 :
168 60 : auto handle_ng_result = [](Ng_Result result) {
169 : static const std::vector<std::string> result_types =
170 : {"Netgen error", "Netgen success", "Netgen surface input error",
171 : "Netgen volume failure", "Netgen STL input error",
172 154 : "Netgen surface failure", "Netgen file not found"};
173 :
174 65 : if (result+1 >= 0 &&
175 60 : std::size_t(result+1) < result_types.size())
176 60 : libmesh_error_msg_if
177 : (result, "Ng_GenerateVolumeMesh failed: " <<
178 : result_types[result+1]);
179 : else
180 0 : libmesh_error_msg
181 : ("Ng_GenerateVolumeMesh failed with an unknown error code");
182 80 : };
183 :
184 : // Keep track of what boundary ids we want to assign to each new
185 : // triangle. We'll give the outer boundary BC 0, and give holes ids
186 : // starting from 1.
187 : // We key on sorted tuples of node ids to identify a side.
188 : std::unordered_map<std::array<dof_id_type,3>,
189 10 : boundary_id_type, libMesh::hash> side_boundary_id;
190 :
191 88056 : auto insert_id = []
192 : (std::array<dof_id_type,3> & array,
193 8724 : dof_id_type n_id)
194 : {
195 8724 : libmesh_assert_less(n_id, DofObject::invalid_id);
196 8724 : unsigned int i=0;
197 157506 : while (array[i] < n_id)
198 52002 : ++i;
199 370014 : while (i < 3)
200 264510 : std::swap(array[i++], n_id);
201 96780 : };
202 :
203 10 : WrappedNgMesh ngmesh;
204 :
205 : // Create surface mesh in the WrappedNgMesh
206 : {
207 : // NetGen appears to use ONE-BASED numbering for its nodes, and
208 : // since it doesn't return an id when adding nodes we'll have to
209 : // track the numbering ourselves.
210 60 : int ng_id = 1;
211 :
212 : auto create_surface_component =
213 70 : [this, &ng_id, &ng_to_libmesh_id, &ngmesh, &side_boundary_id, &insert_id]
214 : (UnstructuredMesh & srcmesh, bool hole_mesh,
215 8314 : boundary_id_type bcid)
216 : {
217 14 : LOG_SCOPE("create_surface_component()", "NetGenMeshInterface");
218 :
219 : // Keep track of what nodes we've already added to the Netgen
220 : // mesh vs what nodes we need to add. We'll keep track by id,
221 : // not by point location. I don't know if Netgen can handle
222 : // multiple nodes with the same point location, but if they can
223 : // it's not going to be *us* who breaks that feature.
224 14 : std::unordered_map<dof_id_type, int> libmesh_to_ng_id;
225 :
226 : // Keep track of what nodes we've already added to the main
227 : // mesh from a hole mesh.
228 14 : std::unordered_map<dof_id_type, dof_id_type> hole_to_main_mesh_id;
229 :
230 : // Use a separate array for passing points to NetGen, just in case
231 : // we're not using double-precision ourselves.
232 : std::array<double, 3> point_val;
233 :
234 : // And an array for element vertices
235 : std::array<int, 3> elem_nodes;
236 :
237 9577 : for (const auto * elem : srcmesh.element_ptr_range())
238 : {
239 : // If someone has triangles we can't triangulate, we have a
240 : // problem
241 10272 : if (elem->type() == TRI6 ||
242 5136 : elem->type() == TRI7)
243 0 : libmesh_not_implemented_msg
244 : ("Netgen tetrahedralization currently only supports TRI3 boundaries");
245 :
246 : // If someone has non-triangles, let's just ignore them.
247 5136 : if (elem->type() != TRI3)
248 0 : continue;
249 :
250 5136 : std::array<dof_id_type,3> sorted_ids =
251 : {DofObject::invalid_id, DofObject::invalid_id,
252 : DofObject::invalid_id};
253 :
254 20544 : for (int ni : make_range(3))
255 : {
256 : // Just using the "invert_trigs" option in NetGen params
257 : // doesn't work for me, so we'll have to have properly
258 : // oriented the tris earlier.
259 15408 : auto & elem_node = hole_mesh ? elem_nodes[2-ni] : elem_nodes[ni];
260 :
261 2568 : const Node & n = elem->node_ref(ni);
262 15408 : auto n_id = n.id();
263 15408 : if (hole_mesh)
264 : {
265 7200 : if (auto it = hole_to_main_mesh_id.find(n_id);
266 600 : it != hole_to_main_mesh_id.end())
267 : {
268 5952 : n_id = it->second;
269 : }
270 : else
271 : {
272 1248 : Node * n_new = this->_mesh.add_point(n);
273 1248 : const dof_id_type n_new_id = n_new->id();
274 104 : hole_to_main_mesh_id.emplace(n_id, n_new_id);
275 1248 : n_id = n_new_id;
276 : }
277 : }
278 :
279 15408 : if (auto it = libmesh_to_ng_id.find(n_id);
280 1284 : it != libmesh_to_ng_id.end())
281 : {
282 12672 : const int existing_ng_id = it->second;
283 12672 : elem_node = existing_ng_id;
284 : }
285 : else
286 : {
287 10944 : for (auto i : make_range(3))
288 8208 : point_val[i] = double(n(i));
289 :
290 2736 : Ng_AddPoint(ngmesh, point_val.data());
291 :
292 2736 : ng_to_libmesh_id[ng_id] = n_id;
293 2736 : libmesh_to_ng_id[n_id] = ng_id;
294 2736 : elem_node = ng_id;
295 2736 : ++ng_id;
296 : }
297 :
298 15408 : insert_id(sorted_ids, n_id);
299 : }
300 :
301 5136 : side_boundary_id[sorted_ids] = bcid;
302 :
303 5136 : Ng_AddSurfaceElement(ngmesh, NG_TRIG, elem_nodes.data());
304 70 : }
305 84 : };
306 :
307 : // Number the outer boundary 0, and the holes starting from 1
308 5 : boundary_id_type bcid = 0;
309 :
310 60 : create_surface_component(this->_mesh, false, bcid);
311 :
312 60 : if (_holes)
313 48 : for (const std::unique_ptr<UnstructuredMesh> & h : *_holes)
314 26 : create_surface_component(*h, true, ++bcid);
315 : }
316 :
317 : {
318 10 : LOG_SCOPE("Ng_GenerateVolumeMesh()", "NetGenMeshInterface");
319 :
320 60 : auto result = Ng_GenerateVolumeMesh(ngmesh, ¶ms);
321 60 : handle_ng_result(result);
322 : }
323 :
324 60 : const int n_elem = Ng_GetNE(ngmesh);
325 :
326 60 : libmesh_error_msg_if (n_elem <= 0,
327 : "NetGen failed to generate any tetrahedra");
328 :
329 60 : const dof_id_type n_points = Ng_GetNP(ngmesh);
330 60 : const dof_id_type old_nodes = this->_mesh.n_nodes();
331 :
332 : // Netgen may have generated new interior nodes
333 60 : if (n_points != old_nodes)
334 : {
335 : std::array<double, 3> point_val;
336 :
337 : // We should only be getting new nodes if we asked for them
338 3 : if (!_desired_volume)
339 : {
340 : std::cout <<
341 0 : "NetGen output " << n_points <<
342 0 : " points when we gave it " <<
343 0 : old_nodes << " and disabled refinement\n" <<
344 : "If new interior points are acceptable in your mesh, please set\n" <<
345 : "a non-zero desired_volume to indicate that. If new interior\n" <<
346 : "points are not acceptable in your mesh, you may need a different\n" <<
347 0 : "(non-advancing-front?) mesh generator." << std::endl;
348 0 : libmesh_error();
349 : }
350 : else
351 17 : for (auto i : make_range(old_nodes, n_points))
352 : {
353 : // i+1 since ng uses ONE-BASED numbering
354 14 : Ng_GetPoint (ngmesh, i+1, point_val.data());
355 14 : const Point p(point_val[0], point_val[1], point_val[2]);
356 14 : Node * n_new = this->_mesh.add_point(p);
357 0 : const dof_id_type n_new_id = n_new->id();
358 14 : ng_to_libmesh_id[i+1] = n_new_id;
359 : }
360 : }
361 :
362 5131 : for (auto * elem : this->_mesh.element_ptr_range())
363 2786 : this->_mesh.delete_elem(elem);
364 :
365 60 : BoundaryInfo * bi = & this->_mesh.get_boundary_info();
366 :
367 7568 : for (auto i : make_range(n_elem))
368 : {
369 : // Enough data to return even a Tet10 without a segfault if nglib
370 : // went nuts
371 : int ngnodes[11];
372 :
373 : // i+1 since we must be 1-based with these ids too...
374 : Ng_Volume_Element_Type ngtype =
375 7508 : Ng_GetVolumeElement(ngmesh, i+1, ngnodes);
376 :
377 : // But really nglib shouldn't go nuts
378 620 : libmesh_assert(ngtype == NG_TET);
379 620 : libmesh_ignore(ngtype);
380 :
381 8128 : auto elem = this->_mesh.add_elem(Elem::build_with_id(TET4, i));
382 37540 : for (auto n : make_range(4))
383 : {
384 : const dof_id_type node_id =
385 30032 : libmesh_map_find(ng_to_libmesh_id, ngnodes[n]);
386 30032 : elem->set_node(n, this->_mesh.node_ptr(node_id));
387 : }
388 :
389 : // NetGen and we disagree about node numbering orientation
390 6888 : elem->orient(bi);
391 :
392 37540 : for (auto s : make_range(4))
393 : {
394 30032 : std::array<dof_id_type,3> sorted_ids =
395 : {DofObject::invalid_id, DofObject::invalid_id,
396 : DofObject::invalid_id};
397 :
398 32512 : std::vector<unsigned int> nos = elem->nodes_on_side(s);
399 120128 : for (auto n : nos)
400 90096 : insert_id(sorted_ids, elem->node_id(n));
401 :
402 30032 : if (auto it = side_boundary_id.find(sorted_ids);
403 2480 : it != side_boundary_id.end())
404 5136 : bi->add_side(elem, s, it->second);
405 : }
406 : }
407 :
408 : // We don't need our holes anymore. Delete their serializers
409 : // first to avoid dereferencing dangling pointers.
410 5 : hole_serializers.clear();
411 60 : if (_holes)
412 2 : _holes->clear();
413 :
414 : // Send our data to other ranks, then fix it up together
415 60 : MeshCommunication().broadcast(this->_mesh);
416 60 : this->_mesh.prepare_for_use();
417 335 : }
418 :
419 :
420 :
421 : } // namespace libMesh
422 :
423 :
424 : #endif // #ifdef LIBMESH_HAVE_NETGEN
|
