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/threads.h"
34 : #include "libmesh/unstructured_mesh.h"
35 : #include "libmesh/utility.h" // libmesh_map_find
36 :
37 : namespace nglib {
38 : #include "netgen/nglib/nglib.h"
39 : }
40 :
41 : namespace {
42 :
43 : // RAII for exception safety
44 : class WrappedNgMesh
45 : {
46 : public:
47 72 : WrappedNgMesh() {
48 72 : _ngmesh = nglib::Ng_NewMesh();
49 6 : }
50 :
51 18 : ~WrappedNgMesh() {
52 72 : nglib::Ng_DeleteMesh(_ngmesh);
53 66 : }
54 :
55 : void clear() {
56 : nglib::Ng_DeleteMesh(_ngmesh);
57 : _ngmesh = nglib::Ng_NewMesh();
58 : }
59 :
60 1312 : operator nglib::Ng_Mesh* () {
61 14438 : return _ngmesh;
62 : }
63 :
64 : private:
65 : nglib::Ng_Mesh * _ngmesh;
66 : };
67 :
68 : }
69 :
70 : namespace libMesh
71 : {
72 :
73 : //----------------------------------------------------------------------
74 : // NetGenMeshInterface class members
75 497 : NetGenMeshInterface::NetGenMeshInterface (UnstructuredMesh & mesh) :
76 : MeshTetInterface(mesh),
77 497 : _serializer(mesh)
78 : {
79 497 : }
80 :
81 :
82 :
83 497 : void NetGenMeshInterface::triangulate ()
84 : {
85 : using namespace nglib;
86 :
87 16 : LOG_SCOPE("triangulate()", "NetGenMeshInterface");
88 :
89 : // We're hoping to do volume_to_surface_mesh in parallel at least,
90 : // but then we'll need to serialize any hole meshes to rank 0 so it
91 : // can use them in serial.
92 :
93 : const BoundingBox mesh_bb =
94 497 : MeshTetInterface::volume_to_surface_mesh(this->_mesh);
95 :
96 18 : std::vector<MeshSerializer> hole_serializers;
97 426 : if (_holes)
98 284 : for (std::unique_ptr<UnstructuredMesh> & hole : *_holes)
99 : {
100 : const BoundingBox hole_bb =
101 142 : MeshTetInterface::volume_to_surface_mesh(*hole);
102 :
103 142 : libmesh_error_msg_if
104 : (!mesh_bb.contains(hole_bb),
105 : "Found hole with bounding box " << hole_bb <<
106 : "\nextending outside of mesh bounding box " << mesh_bb);
107 :
108 : hole_serializers.emplace_back
109 276 : (*hole, /* need_serial */ true,
110 146 : /* serial_only_needed_on_proc_0 */ true);
111 : }
112 :
113 : // This should probably only be done on rank 0, but the API is
114 : // designed with the hope that we'll parallelize it eventually
115 426 : auto integrity = this->improve_hull_integrity();
116 426 : this->process_hull_integrity_result(integrity);
117 :
118 : // If we're not rank 0, we're just going to wait for rank 0 to call
119 : // Netgen, then receive its data afterward, we're not going to hope
120 : // that Netgen does the exact same thing on every processor.
121 438 : if (this->_mesh.processor_id() != 0)
122 : {
123 : // We don't need our holes anymore. Delete their serializers
124 : // first to avoid dereferencing dangling pointers.
125 6 : hole_serializers.clear();
126 354 : if (_holes)
127 2 : _holes->clear();
128 :
129 : // Receive the mesh data rank 0 will send later, then fix it up
130 : // together
131 354 : MeshCommunication().broadcast(this->_mesh);
132 :
133 : // If we got an empty mesh here then our tetrahedralization
134 : // failed.
135 354 : libmesh_error_msg_if (!this->_mesh.n_elem(),
136 : "NetGen failed to generate any tetrahedra");
137 :
138 354 : this->_mesh.prepare_for_use();
139 6 : return;
140 : }
141 :
142 72 : Ng_Meshing_Parameters params;
143 :
144 72 : Ng_SetNumThreads(cast_int<int>(libMesh::n_threads()));
145 :
146 : // Override any default parameters we might need to, to avoid
147 : // inserting nodes we don't want.
148 72 : params.uselocalh = false;
149 72 : params.minh = 0;
150 72 : params.elementsperedge = 1;
151 72 : params.elementspercurve = 1;
152 72 : params.closeedgeenable = false;
153 72 : params.closeedgefact = 0;
154 72 : params.minedgelenenable = false;
155 72 : params.minedgelen = 0;
156 :
157 : // Try to get a no-extra-nodes mesh if we're asked to, or try to
158 : // translate our desired volume into NetGen terms otherwise.
159 : //
160 : // Spoiler alert: all we can do is try; NetGen uses a marching front
161 : // algorithm that can insert extra nodes despite all my best
162 : // efforts.
163 72 : if (_desired_volume == 0) // shorthand for "no refinement"
164 : {
165 60 : params.maxh = std::numeric_limits<double>::max();
166 60 : params.fineness = 0; // "coarse" in the docs
167 60 : params.grading = 1; // "aggressive local grading" to avoid smoothing??
168 :
169 : // Turning off optimization steps avoids another opportunity for
170 : // Netgen to try to add more nodes.
171 60 : params.optsteps_3d = 0;
172 : }
173 : else
174 12 : params.maxh = double(std::pow(_desired_volume, 1./3.));
175 :
176 : // Keep track of how NetGen copies of nodes map back to our original
177 : // nodes, so we can connect new elements to nodes correctly.
178 12 : std::unordered_map<int, dof_id_type> ng_to_libmesh_id;
179 :
180 72 : auto handle_ng_result = [](Ng_Result result) {
181 : static const std::vector<std::string> result_types =
182 : {"Netgen error", "Netgen success", "Netgen surface input error",
183 : "Netgen volume failure", "Netgen STL input error",
184 166 : "Netgen surface failure", "Netgen file not found"};
185 :
186 78 : if (result+1 >= 0 &&
187 72 : std::size_t(result+1) < result_types.size())
188 72 : libmesh_error_msg_if
189 : (result, "Ng_GenerateVolumeMesh failed: " <<
190 : result_types[result+1]);
191 : else
192 0 : libmesh_error_msg
193 : ("Ng_GenerateVolumeMesh failed with an unknown error code");
194 92 : };
195 :
196 : // Keep track of what boundary ids we want to assign to each new
197 : // triangle. We'll give the outer boundary BC 0, and give holes ids
198 : // starting from 1.
199 : // We key on sorted tuples of node ids to identify a side.
200 : std::unordered_map<std::array<dof_id_type,3>,
201 12 : boundary_id_type, libMesh::hash> side_boundary_id;
202 :
203 88020 : auto insert_id = []
204 : (std::array<dof_id_type,3> & array,
205 8796 : dof_id_type n_id)
206 : {
207 8796 : libmesh_assert_less(n_id, DofObject::invalid_id);
208 8796 : unsigned int i=0;
209 157886 : while (array[i] < n_id)
210 52274 : ++i;
211 370174 : while (i < 3)
212 264562 : std::swap(array[i++], n_id);
213 96816 : };
214 :
215 12 : WrappedNgMesh ngmesh;
216 :
217 : // Create surface mesh in the WrappedNgMesh
218 : {
219 : // NetGen appears to use ONE-BASED numbering for its nodes, and
220 : // since it doesn't return an id when adding nodes we'll have to
221 : // track the numbering ourselves.
222 72 : int ng_id = 1;
223 :
224 : auto create_surface_component =
225 80 : [this, &ng_id, &ng_to_libmesh_id, &ngmesh, &side_boundary_id, &insert_id]
226 : (UnstructuredMesh & srcmesh, bool hole_mesh,
227 8474 : boundary_id_type bcid)
228 : {
229 16 : LOG_SCOPE("create_surface_component()", "NetGenMeshInterface");
230 :
231 : // Keep track of what nodes we've already added to the Netgen
232 : // mesh vs what nodes we need to add. We'll keep track by id,
233 : // not by point location. I don't know if Netgen can handle
234 : // multiple nodes with the same point location, but if they can
235 : // it's not going to be *us* who breaks that feature.
236 16 : std::unordered_map<dof_id_type, int> libmesh_to_ng_id;
237 :
238 : // Keep track of what nodes we've already added to the main
239 : // mesh from a hole mesh.
240 16 : std::unordered_map<dof_id_type, dof_id_type> hole_to_main_mesh_id;
241 :
242 : // Use a separate array for passing points to NetGen, just in case
243 : // we're not using double-precision ourselves.
244 : std::array<double, 3> point_val;
245 :
246 : // And an array for element vertices
247 : std::array<int, 3> elem_nodes;
248 :
249 9776 : for (const auto * elem : srcmesh.element_ptr_range())
250 : {
251 : // If someone has triangles we can't triangulate, we have a
252 : // problem
253 10464 : if (elem->type() == TRI6 ||
254 5232 : elem->type() == TRI7)
255 0 : libmesh_not_implemented_msg
256 : ("Netgen tetrahedralization currently only supports TRI3 boundaries");
257 :
258 : // If someone has non-triangles, let's just ignore them.
259 5232 : if (elem->type() != TRI3)
260 0 : continue;
261 :
262 5232 : std::array<dof_id_type,3> sorted_ids =
263 : {DofObject::invalid_id, DofObject::invalid_id,
264 : DofObject::invalid_id};
265 :
266 20928 : for (int ni : make_range(3))
267 : {
268 : // Just using the "invert_trigs" option in NetGen params
269 : // doesn't work for me, so we'll have to have properly
270 : // oriented the tris earlier.
271 15696 : auto & elem_node = hole_mesh ? elem_nodes[2-ni] : elem_nodes[ni];
272 :
273 2616 : const Node & n = elem->node_ref(ni);
274 15696 : auto n_id = n.id();
275 15696 : if (hole_mesh)
276 : {
277 7200 : if (auto it = hole_to_main_mesh_id.find(n_id);
278 600 : it != hole_to_main_mesh_id.end())
279 : {
280 5952 : n_id = it->second;
281 : }
282 : else
283 : {
284 1248 : Node * n_new = this->_mesh.add_point(n);
285 1248 : const dof_id_type n_new_id = n_new->id();
286 104 : hole_to_main_mesh_id.emplace(n_id, n_new_id);
287 1248 : n_id = n_new_id;
288 : }
289 : }
290 :
291 15696 : if (auto it = libmesh_to_ng_id.find(n_id);
292 1308 : it != libmesh_to_ng_id.end())
293 : {
294 12888 : const int existing_ng_id = it->second;
295 12888 : elem_node = existing_ng_id;
296 : }
297 : else
298 : {
299 11232 : for (auto i : make_range(3))
300 8424 : point_val[i] = double(n(i));
301 :
302 2808 : Ng_AddPoint(ngmesh, point_val.data());
303 :
304 2808 : ng_to_libmesh_id[ng_id] = n_id;
305 2808 : libmesh_to_ng_id[n_id] = ng_id;
306 2808 : elem_node = ng_id;
307 2808 : ++ng_id;
308 : }
309 :
310 15696 : insert_id(sorted_ids, n_id);
311 : }
312 :
313 5232 : side_boundary_id[sorted_ids] = bcid;
314 :
315 5232 : Ng_AddSurfaceElement(ngmesh, NG_TRIG, elem_nodes.data());
316 80 : }
317 96 : };
318 :
319 : // Number the outer boundary 0, and the holes starting from 1
320 6 : boundary_id_type bcid = 0;
321 :
322 72 : create_surface_component(this->_mesh, false, bcid);
323 :
324 72 : if (_holes)
325 48 : for (const std::unique_ptr<UnstructuredMesh> & h : *_holes)
326 26 : create_surface_component(*h, true, ++bcid);
327 : }
328 :
329 : {
330 12 : LOG_SCOPE("Ng_GenerateVolumeMesh()", "NetGenMeshInterface");
331 :
332 72 : auto result = Ng_GenerateVolumeMesh(ngmesh, ¶ms);
333 72 : handle_ng_result(result);
334 : }
335 :
336 72 : const int n_elem = Ng_GetNE(ngmesh);
337 :
338 : // If Netgen fails us, we're likely to get n_elem <= 0. This is a
339 : // common enough failure from bad setups that I want to make sure
340 : // it's thrown in parallel so as to not desynchronize any unit tests
341 : // that trigger it. So we'll broadcast the empty mesh to indicate
342 : // the problem and enable throwing exceptions in parallel.
343 72 : if (n_elem <= 0)
344 : {
345 0 : this->_mesh.clear();
346 0 : MeshCommunication().broadcast(this->_mesh);
347 0 : libmesh_error_msg ("NetGen failed to generate any tetrahedra");
348 : }
349 :
350 72 : const dof_id_type n_points = Ng_GetNP(ngmesh);
351 72 : const dof_id_type old_nodes = this->_mesh.n_nodes();
352 :
353 : // Netgen may have generated new interior nodes
354 72 : if (n_points != old_nodes)
355 : {
356 : std::array<double, 3> point_val;
357 :
358 : // We should only be getting new nodes if we asked for them
359 1 : if (!_desired_volume)
360 : {
361 : std::cout <<
362 0 : "NetGen output " << n_points <<
363 0 : " points when we gave it " <<
364 0 : old_nodes << " and disabled refinement\n" <<
365 : "If new interior points are acceptable in your mesh, please set\n" <<
366 : "a non-zero desired_volume to indicate that. If new interior\n" <<
367 : "points are not acceptable in your mesh, you may need a different\n" <<
368 0 : "(non-advancing-front?) mesh generator." << std::endl;
369 0 : libmesh_error();
370 : }
371 : else
372 2 : for (auto i : make_range(old_nodes, n_points))
373 : {
374 : // i+1 since ng uses ONE-BASED numbering
375 1 : Ng_GetPoint (ngmesh, i+1, point_val.data());
376 1 : const Point p(point_val[0], point_val[1], point_val[2]);
377 1 : Node * n_new = this->_mesh.add_point(p);
378 0 : const dof_id_type n_new_id = n_new->id();
379 1 : ng_to_libmesh_id[i+1] = n_new_id;
380 : }
381 : }
382 :
383 5330 : for (auto * elem : this->_mesh.element_ptr_range())
384 2892 : this->_mesh.delete_elem(elem);
385 :
386 72 : BoundaryInfo * bi = & this->_mesh.get_boundary_info();
387 :
388 7565 : for (auto i : make_range(n_elem))
389 : {
390 : // Enough data to return even a Tet10 without a segfault if nglib
391 : // went nuts
392 : int ngnodes[11];
393 :
394 : // i+1 since we must be 1-based with these ids too...
395 : Ng_Volume_Element_Type ngtype =
396 7493 : Ng_GetVolumeElement(ngmesh, i+1, ngnodes);
397 :
398 : // But really nglib shouldn't go nuts
399 624 : libmesh_assert(ngtype == NG_TET);
400 624 : libmesh_ignore(ngtype);
401 :
402 8117 : auto elem = this->_mesh.add_elem(Elem::build_with_id(TET4, i));
403 37465 : for (auto n : make_range(4))
404 : {
405 : const dof_id_type node_id =
406 29972 : libmesh_map_find(ng_to_libmesh_id, ngnodes[n]);
407 29972 : elem->set_node(n, this->_mesh.node_ptr(node_id));
408 : }
409 :
410 : // NetGen and we disagree about node numbering orientation
411 6869 : elem->orient(bi);
412 :
413 37465 : for (auto s : make_range(4))
414 : {
415 29972 : std::array<dof_id_type,3> sorted_ids =
416 : {DofObject::invalid_id, DofObject::invalid_id,
417 : DofObject::invalid_id};
418 :
419 32468 : std::vector<unsigned int> nos = elem->nodes_on_side(s);
420 119888 : for (auto n : nos)
421 89916 : insert_id(sorted_ids, elem->node_id(n));
422 :
423 29972 : if (auto it = side_boundary_id.find(sorted_ids);
424 2496 : it != side_boundary_id.end())
425 5232 : bi->add_side(elem, s, it->second);
426 : }
427 : }
428 :
429 : // We don't need our holes anymore. Delete their serializers
430 : // first to avoid dereferencing dangling pointers.
431 6 : hole_serializers.clear();
432 72 : if (_holes)
433 2 : _holes->clear();
434 :
435 : // Send our data to other ranks, then fix it up together
436 72 : MeshCommunication().broadcast(this->_mesh);
437 72 : this->_mesh.prepare_for_use();
438 402 : }
439 :
440 :
441 :
442 : } // namespace libMesh
443 :
444 :
445 : #endif // #ifdef LIBMESH_HAVE_NETGEN
|
