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/cell_tet10.h"
31 : #include "libmesh/elem.h"
32 : #include "libmesh/face_tri3.h"
33 : #include "libmesh/face_tri6.h"
34 : #include "libmesh/libmesh_logging.h"
35 : #include "libmesh/mesh_communication.h"
36 : #include "libmesh/threads.h"
37 : #include "libmesh/unstructured_mesh.h"
38 : #include "libmesh/utility.h" // libmesh_map_find
39 :
40 : namespace nglib {
41 : #include "netgen/nglib/nglib.h"
42 : }
43 :
44 : namespace {
45 :
46 : // RAII for exception safety
47 : class WrappedNgMesh
48 : {
49 : public:
50 120 : WrappedNgMesh() {
51 120 : _ngmesh = nglib::Ng_NewMesh();
52 10 : }
53 :
54 30 : ~WrappedNgMesh() {
55 120 : nglib::Ng_DeleteMesh(_ngmesh);
56 110 : }
57 :
58 : void clear() {
59 : nglib::Ng_DeleteMesh(_ngmesh);
60 : _ngmesh = nglib::Ng_NewMesh();
61 : }
62 :
63 1396 : operator nglib::Ng_Mesh* () {
64 15362 : return _ngmesh;
65 : }
66 :
67 : private:
68 : nglib::Ng_Mesh * _ngmesh;
69 : };
70 :
71 : }
72 :
73 : namespace libMesh
74 : {
75 :
76 : //----------------------------------------------------------------------
77 : // NetGenMeshInterface class members
78 781 : NetGenMeshInterface::NetGenMeshInterface (UnstructuredMesh & mesh) :
79 : MeshTetInterface(mesh),
80 781 : _serializer(mesh)
81 : {
82 781 : }
83 :
84 :
85 :
86 781 : void NetGenMeshInterface::triangulate ()
87 : {
88 : using namespace nglib;
89 :
90 24 : LOG_SCOPE("triangulate()", "NetGenMeshInterface");
91 :
92 803 : if (_elem_type != TET4 &&
93 767 : _elem_type != TET10 &&
94 4 : _elem_type != TET14)
95 0 : libmesh_not_implemented();
96 :
97 : // We're hoping to do volume_to_surface_mesh in parallel at least,
98 : // but then we'll need to serialize any hole meshes to rank 0 so it
99 : // can use them in serial.
100 :
101 : // If the user wants higher-order output, record midpoints from any
102 : // quadratic boundary elements before stripping them to TRI3 for
103 : // NetGen. We key by sorted position pairs (not node IDs) so that
104 : // outer mesh and hole mesh node namespaces cannot conflict.
105 24 : std::map<std::pair<Point,Point>, Point> edge_midpoints;
106 :
107 : // TRI7 boundary faces additionally carry a face-centroid node (local index
108 : // 6) that TET14 output reproduces (its faces are TRI7). Record those too,
109 : // keyed by the sorted triple of the face's three corner positions so the
110 : // key is independent of node id and rotation.
111 24 : std::map<std::array<Point,3>, Point> face_centroids;
112 :
113 : auto record_and_strip =
114 3992 : [&edge_midpoints, &face_centroids](UnstructuredMesh & m)
115 : {
116 8 : bool has_quadratic = false;
117 4976 : for (const auto & elem : m.element_ptr_range())
118 : {
119 2272 : if (elem->type() != TRI6 && elem->type() != TRI7) continue;
120 32 : has_quadratic = true;
121 : // TRI6/TRI7: edge e runs node[e]→node[(e+1)%3], midpoint=node[e+3]
122 4544 : for (auto e : make_range(3u))
123 : {
124 3408 : const Point & pa = elem->point(e);
125 3408 : const Point & pb = elem->point((e+1)%3);
126 3408 : auto key = pa < pb
127 3408 : ? std::make_pair(pa,pb) : std::make_pair(pb,pa);
128 3504 : edge_midpoints[key] = elem->point(e+3);
129 : }
130 : // TRI7: node 6 is the face-centroid node.
131 1136 : if (elem->type() == TRI7)
132 : {
133 : std::array<Point,3> corners =
134 616 : {elem->point(0), elem->point(1), elem->point(2)};
135 16 : std::sort(corners.begin(), corners.end());
136 584 : face_centroids[corners] = elem->point(6);
137 : }
138 268 : }
139 284 : if (has_quadratic)
140 142 : m.all_first_order();
141 306 : };
142 :
143 : const BoundingBox mesh_bb =
144 781 : MeshTetInterface::volume_to_surface_mesh(this->_mesh);
145 :
146 710 : if (_elem_type != TET4)
147 284 : record_and_strip(this->_mesh);
148 :
149 30 : std::vector<MeshSerializer> hole_serializers;
150 690 : if (_holes)
151 284 : for (std::unique_ptr<UnstructuredMesh> & hole : *_holes)
152 : {
153 : const BoundingBox hole_bb =
154 142 : MeshTetInterface::volume_to_surface_mesh(*hole);
155 :
156 142 : if (_elem_type != TET4)
157 0 : record_and_strip(*hole);
158 :
159 142 : libmesh_error_msg_if
160 : (!mesh_bb.contains(hole_bb),
161 : "Found hole with bounding box " << hole_bb <<
162 : "\nextending outside of mesh bounding box " << mesh_bb);
163 :
164 : hole_serializers.emplace_back
165 276 : (*hole, /* need_serial */ true,
166 146 : /* serial_only_needed_on_proc_0 */ true);
167 : }
168 :
169 : // Increasing the element order (all_second_order()/all_complete_order()
170 : // via increase_tet_order()) performs collective MPI communication, so it
171 : // must run on every rank in lockstep -- it cannot run on rank 0 alone
172 : // while the other ranks wait in broadcast(). We therefore broadcast
173 : // NetGen's TET4 result first, then increase the order and restore the
174 : // curved-boundary midpoints identically on all ranks. This mirrors what
175 : // the 2D interfaces do (see TriangulatorInterface::increase_triangle_order,
176 : // which is likewise called on all ranks). edge_midpoints was built while
177 : // the mesh was serialized on every rank, so it is identical everywhere and
178 : // this post-broadcast fixup is deterministic.
179 : auto increase_order_and_restore_midpoints =
180 4448 : [this, &edge_midpoints, &face_centroids]()
181 : {
182 710 : if (_elem_type == TET4)
183 12 : return;
184 :
185 : // find_neighbors() is needed before all_second_order() can place
186 : // shared edge midpoints correctly.
187 284 : this->_mesh.find_neighbors();
188 :
189 : // Refresh the cached element dimensions. We just replaced the 2D
190 : // TRI3 surface elements with 3D TET4 volume elements, but the mesh's
191 : // cached dimension is still that of the original surface (2).
192 : // all_second_order()/all_complete_order() derive the per-element
193 : // unique_id reservation width from mesh_dimension(): a stale value of
194 : // 2 reserves only 9-4=5 slots per element, too few for the 6 new edge
195 : // nodes of a TET10, so adjacent elements' unique_id ranges overlap and
196 : // collide. cache_elem_data() recomputes the dimension to 3 (reserving
197 : // 27-8=19 slots) before the order increase runs.
198 284 : this->_mesh.cache_elem_data();
199 :
200 284 : this->increase_tet_order();
201 :
202 : // Move auto-placed geometric midpoints to the recorded positions,
203 : // preserving any curvature from the original quadratic boundary.
204 284 : if (!edge_midpoints.empty() || !face_centroids.empty())
205 1936 : for (Elem * elem : _mesh.element_ptr_range())
206 2840 : for (auto s : elem->side_index_range())
207 : {
208 2336 : if (elem->neighbor_ptr(s)) continue;
209 :
210 : // build_side_ptr() returns a TRI6 (TET10) or TRI7 (TET14)
211 : // whose node pointers reference the actual mesh nodes; point()
212 : // assignments update mesh node coordinates in place.
213 1168 : auto side = elem->build_side_ptr(s);
214 4544 : for (auto e : make_range(3u))
215 : {
216 192 : const Point & pa = side->point(e);
217 3408 : const Point & pb = side->point((e+1)%3);
218 3408 : auto key = pa < pb
219 3408 : ? std::make_pair(pa,pb) : std::make_pair(pb,pa);
220 3408 : if (auto it = edge_midpoints.find(key);
221 96 : it != edge_midpoints.end())
222 3504 : side->point(e+3) = it->second;
223 : }
224 :
225 : // TRI7 faces (TET14 output) also carry a face-centroid node
226 : // at local index 6; restore its recorded curved position.
227 1136 : if (side->type() == TRI7)
228 : {
229 : std::array<Point,3> corners =
230 616 : {side->point(0), side->point(1), side->point(2)};
231 16 : std::sort(corners.begin(), corners.end());
232 568 : if (auto it = face_centroids.find(corners);
233 16 : it != face_centroids.end())
234 584 : side->point(6) = it->second;
235 : }
236 1206 : }
237 710 : };
238 :
239 : // This should probably only be done on rank 0, but the API is
240 : // designed with the hope that we'll parallelize it eventually
241 710 : auto integrity = this->improve_hull_integrity();
242 710 : this->process_hull_integrity_result(integrity);
243 :
244 : // If we're not rank 0, we're just going to wait for rank 0 to call
245 : // Netgen, then receive its data afterward, we're not going to hope
246 : // that Netgen does the exact same thing on every processor.
247 730 : if (this->_mesh.processor_id() != 0)
248 : {
249 : // We don't need our holes anymore. Delete their serializers
250 : // first to avoid dereferencing dangling pointers.
251 10 : hole_serializers.clear();
252 590 : if (_holes)
253 2 : _holes->clear();
254 :
255 : // Receive the TET4 mesh data rank 0 will send later.
256 590 : MeshCommunication().broadcast(this->_mesh);
257 :
258 : // If we got an empty mesh here then our tetrahedralization
259 : // failed.
260 590 : libmesh_error_msg_if (!this->_mesh.n_elem(),
261 : "NetGen failed to generate any tetrahedra");
262 :
263 : // Increase the element order collectively, in lockstep with rank 0.
264 590 : increase_order_and_restore_midpoints();
265 :
266 590 : this->_mesh.prepare_for_use();
267 10 : return;
268 : }
269 :
270 120 : Ng_Meshing_Parameters params;
271 :
272 120 : Ng_SetNumThreads(cast_int<int>(libMesh::n_threads()));
273 :
274 : // Override any default parameters we might need to, to avoid
275 : // inserting nodes we don't want.
276 120 : params.uselocalh = false;
277 120 : params.minh = 0;
278 120 : params.elementsperedge = 1;
279 120 : params.elementspercurve = 1;
280 120 : params.closeedgeenable = false;
281 120 : params.closeedgefact = 0;
282 120 : params.minedgelenenable = false;
283 120 : params.minedgelen = 0;
284 :
285 : // Try to get a no-extra-nodes mesh if we're asked to, or try to
286 : // translate our desired volume into NetGen terms otherwise.
287 : //
288 : // Spoiler alert: all we can do is try; NetGen uses a marching front
289 : // algorithm that can insert extra nodes despite all my best
290 : // efforts.
291 120 : if (_desired_volume == 0) // shorthand for "no refinement"
292 : {
293 108 : params.maxh = std::numeric_limits<double>::max();
294 108 : params.fineness = 0; // "coarse" in the docs
295 108 : params.grading = 1; // "aggressive local grading" to avoid smoothing??
296 :
297 : // Turning off optimization steps avoids another opportunity for
298 : // Netgen to try to add more nodes.
299 108 : params.optsteps_3d = 0;
300 : }
301 : else
302 12 : params.maxh = double(std::pow(_desired_volume, 1./3.));
303 :
304 : // Keep track of how NetGen copies of nodes map back to our original
305 : // nodes, so we can connect new elements to nodes correctly.
306 20 : std::unordered_map<int, dof_id_type> ng_to_libmesh_id;
307 :
308 120 : auto handle_ng_result = [](Ng_Result result) {
309 : static const std::vector<std::string> result_types =
310 : {"Netgen error", "Netgen success", "Netgen surface input error",
311 : "Netgen volume failure", "Netgen STL input error",
312 214 : "Netgen surface failure", "Netgen file not found"};
313 :
314 130 : if (result+1 >= 0 &&
315 120 : std::size_t(result+1) < result_types.size())
316 120 : libmesh_error_msg_if
317 : (result, "Ng_GenerateVolumeMesh failed: " <<
318 : result_types[result+1]);
319 : else
320 0 : libmesh_error_msg
321 : ("Ng_GenerateVolumeMesh failed with an unknown error code");
322 140 : };
323 :
324 : // Keep track of what boundary ids we want to assign to each new
325 : // triangle. We'll give the outer boundary BC 0, and give holes ids
326 : // starting from 1.
327 : // We key on sorted tuples of node ids to identify a side.
328 : std::unordered_map<std::array<dof_id_type,3>,
329 20 : boundary_id_type, libMesh::hash> side_boundary_id;
330 :
331 90900 : auto insert_id = []
332 : (std::array<dof_id_type,3> & array,
333 9084 : dof_id_type n_id)
334 : {
335 9084 : libmesh_assert_less(n_id, DofObject::invalid_id);
336 9084 : unsigned int i=0;
337 163762 : while (array[i] < n_id)
338 54694 : ++i;
339 381578 : while (i < 3)
340 272510 : std::swap(array[i++], n_id);
341 99984 : };
342 :
343 20 : WrappedNgMesh ngmesh;
344 :
345 : // Create surface mesh in the WrappedNgMesh
346 : {
347 : // NetGen appears to use ONE-BASED numbering for its nodes, and
348 : // since it doesn't return an id when adding nodes we'll have to
349 : // track the numbering ourselves.
350 120 : int ng_id = 1;
351 :
352 : auto create_surface_component =
353 120 : [this, &ng_id, &ng_to_libmesh_id, &ngmesh, &side_boundary_id, &insert_id]
354 : (UnstructuredMesh & srcmesh, bool hole_mesh,
355 9114 : boundary_id_type bcid)
356 : {
357 24 : LOG_SCOPE("create_surface_component()", "NetGenMeshInterface");
358 :
359 : // Keep track of what nodes we've already added to the Netgen
360 : // mesh vs what nodes we need to add. We'll keep track by id,
361 : // not by point location. I don't know if Netgen can handle
362 : // multiple nodes with the same point location, but if they can
363 : // it's not going to be *us* who breaks that feature.
364 24 : std::unordered_map<dof_id_type, int> libmesh_to_ng_id;
365 :
366 : // Keep track of what nodes we've already added to the main
367 : // mesh from a hole mesh.
368 24 : std::unordered_map<dof_id_type, dof_id_type> hole_to_main_mesh_id;
369 :
370 : // Use a separate array for passing points to NetGen, just in case
371 : // we're not using double-precision ourselves.
372 : std::array<double, 3> point_val;
373 :
374 : // And an array for element vertices
375 : std::array<int, 3> elem_nodes;
376 :
377 10572 : for (const auto * elem : srcmesh.element_ptr_range())
378 : {
379 : // If someone has triangles we can't triangulate, we have a
380 : // problem
381 11232 : if (elem->type() == TRI6 ||
382 5616 : elem->type() == TRI7)
383 0 : libmesh_not_implemented_msg
384 : ("Netgen tetrahedralization currently only supports TRI3 boundaries");
385 :
386 : // If someone has non-triangles, let's just ignore them.
387 5616 : if (elem->type() != TRI3)
388 0 : continue;
389 :
390 5616 : std::array<dof_id_type,3> sorted_ids =
391 : {DofObject::invalid_id, DofObject::invalid_id,
392 : DofObject::invalid_id};
393 :
394 22464 : for (int ni : make_range(3))
395 : {
396 : // Just using the "invert_trigs" option in NetGen params
397 : // doesn't work for me, so we'll have to have properly
398 : // oriented the tris earlier.
399 16848 : auto & elem_node = hole_mesh ? elem_nodes[2-ni] : elem_nodes[ni];
400 :
401 2808 : const Node & n = elem->node_ref(ni);
402 16848 : auto n_id = n.id();
403 16848 : if (hole_mesh)
404 : {
405 7200 : if (auto it = hole_to_main_mesh_id.find(n_id);
406 600 : it != hole_to_main_mesh_id.end())
407 : {
408 5952 : n_id = it->second;
409 : }
410 : else
411 : {
412 1248 : Node * n_new = this->_mesh.add_point(n);
413 1248 : const dof_id_type n_new_id = n_new->id();
414 104 : hole_to_main_mesh_id.emplace(n_id, n_new_id);
415 1248 : n_id = n_new_id;
416 : }
417 : }
418 :
419 16848 : if (auto it = libmesh_to_ng_id.find(n_id);
420 1404 : it != libmesh_to_ng_id.end())
421 : {
422 13752 : const int existing_ng_id = it->second;
423 13752 : elem_node = existing_ng_id;
424 : }
425 : else
426 : {
427 12384 : for (auto i : make_range(3))
428 9288 : point_val[i] = double(n(i));
429 :
430 3096 : Ng_AddPoint(ngmesh, point_val.data());
431 :
432 3096 : ng_to_libmesh_id[ng_id] = n_id;
433 3096 : libmesh_to_ng_id[n_id] = ng_id;
434 3096 : elem_node = ng_id;
435 3096 : ++ng_id;
436 : }
437 :
438 16848 : insert_id(sorted_ids, n_id);
439 : }
440 :
441 5616 : side_boundary_id[sorted_ids] = bcid;
442 :
443 5616 : Ng_AddSurfaceElement(ngmesh, NG_TRIG, elem_nodes.data());
444 120 : }
445 144 : };
446 :
447 : // Number the outer boundary 0, and the holes starting from 1
448 10 : boundary_id_type bcid = 0;
449 :
450 120 : create_surface_component(this->_mesh, false, bcid);
451 :
452 120 : if (_holes)
453 48 : for (const std::unique_ptr<UnstructuredMesh> & h : *_holes)
454 26 : create_surface_component(*h, true, ++bcid);
455 : }
456 :
457 : {
458 20 : LOG_SCOPE("Ng_GenerateVolumeMesh()", "NetGenMeshInterface");
459 :
460 120 : auto result = Ng_GenerateVolumeMesh(ngmesh, ¶ms);
461 120 : handle_ng_result(result);
462 : }
463 :
464 120 : const int n_elem = Ng_GetNE(ngmesh);
465 :
466 : // If Netgen fails us, we're likely to get n_elem <= 0. This is a
467 : // common enough failure from bad setups that I want to make sure
468 : // it's thrown in parallel so as to not desynchronize any unit tests
469 : // that trigger it. So we'll broadcast the empty mesh to indicate
470 : // the problem and enable throwing exceptions in parallel.
471 120 : if (n_elem <= 0)
472 : {
473 0 : this->_mesh.clear();
474 0 : MeshCommunication().broadcast(this->_mesh);
475 0 : libmesh_error_msg ("NetGen failed to generate any tetrahedra");
476 : }
477 :
478 120 : const dof_id_type n_points = Ng_GetNP(ngmesh);
479 120 : const dof_id_type old_nodes = this->_mesh.n_nodes();
480 :
481 : // Netgen may have generated new interior nodes
482 120 : if (n_points != old_nodes)
483 : {
484 : std::array<double, 3> point_val;
485 :
486 : // We should only be getting new nodes if we asked for them
487 1 : if (!_desired_volume)
488 : {
489 : std::cout <<
490 0 : "NetGen output " << n_points <<
491 0 : " points when we gave it " <<
492 0 : old_nodes << " and disabled refinement\n" <<
493 : "If new interior points are acceptable in your mesh, please set\n" <<
494 : "a non-zero desired_volume to indicate that. If new interior\n" <<
495 : "points are not acceptable in your mesh, you may need a different\n" <<
496 0 : "(non-advancing-front?) mesh generator." << std::endl;
497 0 : libmesh_error();
498 : }
499 : else
500 2 : for (auto i : make_range(old_nodes, n_points))
501 : {
502 : // i+1 since ng uses ONE-BASED numbering
503 1 : Ng_GetPoint (ngmesh, i+1, point_val.data());
504 1 : const Point p(point_val[0], point_val[1], point_val[2]);
505 1 : Node * n_new = this->_mesh.add_point(p);
506 0 : const dof_id_type n_new_id = n_new->id();
507 1 : ng_to_libmesh_id[i+1] = n_new_id;
508 : }
509 : }
510 :
511 6126 : for (auto * elem : this->_mesh.element_ptr_range())
512 3316 : this->_mesh.delete_elem(elem);
513 :
514 120 : BoundaryInfo * bi = & this->_mesh.get_boundary_info();
515 :
516 7805 : for (auto i : make_range(n_elem))
517 : {
518 : // Enough data to return even a Tet10 without a segfault if nglib
519 : // went nuts
520 : int ngnodes[11];
521 :
522 : // i+1 since we must be 1-based with these ids too...
523 : Ng_Volume_Element_Type ngtype =
524 7685 : Ng_GetVolumeElement(ngmesh, i+1, ngnodes);
525 :
526 : // But really nglib shouldn't go nuts
527 640 : libmesh_assert(ngtype == NG_TET);
528 640 : libmesh_ignore(ngtype);
529 :
530 8325 : auto elem = this->_mesh.add_elem(Elem::build_with_id(TET4, i));
531 38425 : for (auto n : make_range(4))
532 : {
533 : const dof_id_type node_id =
534 30740 : libmesh_map_find(ng_to_libmesh_id, ngnodes[n]);
535 30740 : elem->set_node(n, this->_mesh.node_ptr(node_id));
536 : }
537 :
538 : // NetGen and we disagree about node numbering orientation
539 7045 : elem->orient(bi);
540 :
541 38425 : for (auto s : make_range(4))
542 : {
543 30740 : std::array<dof_id_type,3> sorted_ids =
544 : {DofObject::invalid_id, DofObject::invalid_id,
545 : DofObject::invalid_id};
546 :
547 33300 : std::vector<unsigned int> nos = elem->nodes_on_side(s);
548 122960 : for (auto n : nos)
549 92220 : insert_id(sorted_ids, elem->node_id(n));
550 :
551 30740 : if (auto it = side_boundary_id.find(sorted_ids);
552 2560 : it != side_boundary_id.end())
553 5616 : bi->add_side(elem, s, it->second);
554 : }
555 : }
556 :
557 : // We don't need our holes anymore. Delete their serializers
558 : // first to avoid dereferencing dangling pointers.
559 10 : hole_serializers.clear();
560 120 : if (_holes)
561 2 : _holes->clear();
562 :
563 : // Send NetGen's TET4 result to the other ranks, then increase the element
564 : // order collectively on all ranks together. increase_tet_order() performs
565 : // collective communication, so it must not run on rank 0 alone -- doing so
566 : // would deadlock against the other ranks waiting in broadcast() above.
567 120 : MeshCommunication().broadcast(this->_mesh);
568 120 : increase_order_and_restore_midpoints();
569 120 : this->_mesh.prepare_for_use();
570 670 : }
571 :
572 :
573 :
574 : } // namespace libMesh
575 :
576 :
577 : #endif // #ifdef LIBMESH_HAVE_NETGEN
|
