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