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 4472 : [this, &edge_midpoints, &face_centroids]()
181 : {
182 710 : if (_elem_type == TET4)
183 12 : return;
184 :
185 : // NetGen (and volume_to_surface_mesh) can leave gaps in rank 0's
186 : // element and node id spaces, while broadcast() compacts them on the
187 : // other ranks. Renumber here -- on every rank, in lockstep -- so all
188 : // ranks agree on n_elem()/max_elem_id() before the collective order
189 : // increase runs its cross-rank consistency checks. Note that
190 : // renumber_nodes_and_elements() itself performs a collective reduction
191 : // for the unique-id counter, so it must be called on all ranks together,
192 : // never on rank 0 alone (that would deadlock the others in broadcast()).
193 284 : this->_mesh.renumber_nodes_and_elements();
194 :
195 : // find_neighbors() is needed before all_second_order() can place
196 : // shared edge midpoints correctly.
197 284 : this->_mesh.find_neighbors();
198 :
199 : // Refresh the cached element dimensions. We just replaced the 2D
200 : // TRI3 surface elements with 3D TET4 volume elements, but the mesh's
201 : // cached dimension is still that of the original surface (2).
202 : // all_second_order()/all_complete_order() derive the per-element
203 : // unique_id reservation width from mesh_dimension(): a stale value of
204 : // 2 reserves only 9-4=5 slots per element, too few for the 6 new edge
205 : // nodes of a TET10, so adjacent elements' unique_id ranges overlap and
206 : // collide. cache_elem_data() recomputes the dimension to 3 (reserving
207 : // 27-8=19 slots) before the order increase runs.
208 284 : this->_mesh.cache_elem_data();
209 :
210 284 : this->increase_tet_order();
211 :
212 : // Move auto-placed geometric midpoints to the recorded positions,
213 : // preserving any curvature from the original quadratic boundary.
214 284 : if (!edge_midpoints.empty() || !face_centroids.empty())
215 1936 : for (Elem * elem : _mesh.element_ptr_range())
216 2840 : for (auto s : elem->side_index_range())
217 : {
218 2336 : if (elem->neighbor_ptr(s)) continue;
219 :
220 : // build_side_ptr() returns a TRI6 (TET10) or TRI7 (TET14)
221 : // whose node pointers reference the actual mesh nodes; point()
222 : // assignments update mesh node coordinates in place.
223 1168 : auto side = elem->build_side_ptr(s);
224 4544 : for (auto e : make_range(3u))
225 : {
226 192 : const Point & pa = side->point(e);
227 3408 : const Point & pb = side->point((e+1)%3);
228 3408 : auto key = pa < pb
229 3408 : ? std::make_pair(pa,pb) : std::make_pair(pb,pa);
230 3408 : if (auto it = edge_midpoints.find(key);
231 96 : it != edge_midpoints.end())
232 3504 : side->point(e+3) = it->second;
233 : }
234 :
235 : // TRI7 faces (TET14 output) also carry a face-centroid node
236 : // at local index 6; restore its recorded curved position.
237 1136 : if (side->type() == TRI7)
238 : {
239 : std::array<Point,3> corners =
240 616 : {side->point(0), side->point(1), side->point(2)};
241 16 : std::sort(corners.begin(), corners.end());
242 568 : if (auto it = face_centroids.find(corners);
243 16 : it != face_centroids.end())
244 584 : side->point(6) = it->second;
245 : }
246 1206 : }
247 710 : };
248 :
249 : // This should probably only be done on rank 0, but the API is
250 : // designed with the hope that we'll parallelize it eventually
251 710 : auto integrity = this->improve_hull_integrity();
252 710 : this->process_hull_integrity_result(integrity);
253 :
254 : // If we're not rank 0, we're just going to wait for rank 0 to call
255 : // Netgen, then receive its data afterward, we're not going to hope
256 : // that Netgen does the exact same thing on every processor.
257 730 : if (this->_mesh.processor_id() != 0)
258 : {
259 : // We don't need our holes anymore. Delete their serializers
260 : // first to avoid dereferencing dangling pointers.
261 10 : hole_serializers.clear();
262 590 : if (_holes)
263 2 : _holes->clear();
264 :
265 : // Receive the TET4 mesh data rank 0 will send later.
266 590 : MeshCommunication().broadcast(this->_mesh);
267 :
268 : // If we got an empty mesh here then our tetrahedralization
269 : // failed.
270 590 : libmesh_error_msg_if (!this->_mesh.n_elem(),
271 : "NetGen failed to generate any tetrahedra");
272 :
273 : // Increase the element order collectively, in lockstep with rank 0.
274 590 : increase_order_and_restore_midpoints();
275 :
276 590 : this->_mesh.prepare_for_use();
277 10 : return;
278 : }
279 :
280 120 : Ng_Meshing_Parameters params;
281 :
282 120 : Ng_SetNumThreads(cast_int<int>(libMesh::n_threads()));
283 :
284 : // Override any default parameters we might need to, to avoid
285 : // inserting nodes we don't want.
286 120 : params.uselocalh = false;
287 120 : params.minh = 0;
288 120 : params.elementsperedge = 1;
289 120 : params.elementspercurve = 1;
290 120 : params.closeedgeenable = false;
291 120 : params.closeedgefact = 0;
292 120 : params.minedgelenenable = false;
293 120 : params.minedgelen = 0;
294 :
295 : // Try to get a no-extra-nodes mesh if we're asked to, or try to
296 : // translate our desired volume into NetGen terms otherwise.
297 : //
298 : // Spoiler alert: all we can do is try; NetGen uses a marching front
299 : // algorithm that can insert extra nodes despite all my best
300 : // efforts.
301 120 : if (_desired_volume == 0) // shorthand for "no refinement"
302 : {
303 108 : params.maxh = std::numeric_limits<double>::max();
304 108 : params.fineness = 0; // "coarse" in the docs
305 108 : params.grading = 1; // "aggressive local grading" to avoid smoothing??
306 :
307 : // Turning off optimization steps avoids another opportunity for
308 : // Netgen to try to add more nodes.
309 108 : params.optsteps_3d = 0;
310 : }
311 : else
312 12 : params.maxh = double(std::pow(_desired_volume, 1./3.));
313 :
314 : // Keep track of how NetGen copies of nodes map back to our original
315 : // nodes, so we can connect new elements to nodes correctly.
316 20 : std::unordered_map<int, dof_id_type> ng_to_libmesh_id;
317 :
318 120 : auto handle_ng_result = [](Ng_Result result) {
319 : static const std::vector<std::string> result_types =
320 : {"Netgen error", "Netgen success", "Netgen surface input error",
321 : "Netgen volume failure", "Netgen STL input error",
322 214 : "Netgen surface failure", "Netgen file not found"};
323 :
324 130 : if (result+1 >= 0 &&
325 120 : std::size_t(result+1) < result_types.size())
326 120 : libmesh_error_msg_if
327 : (result, "Ng_GenerateVolumeMesh failed: " <<
328 : result_types[result+1]);
329 : else
330 0 : libmesh_error_msg
331 : ("Ng_GenerateVolumeMesh failed with an unknown error code");
332 140 : };
333 :
334 : // Keep track of what boundary ids we want to assign to each new
335 : // triangle. We'll give the outer boundary BC 0, and give holes ids
336 : // starting from 1.
337 : // We key on sorted tuples of node ids to identify a side.
338 : std::unordered_map<std::array<dof_id_type,3>,
339 20 : boundary_id_type, libMesh::hash> side_boundary_id;
340 :
341 90900 : auto insert_id = []
342 : (std::array<dof_id_type,3> & array,
343 9084 : dof_id_type n_id)
344 : {
345 9084 : libmesh_assert_less(n_id, DofObject::invalid_id);
346 9084 : unsigned int i=0;
347 163762 : while (array[i] < n_id)
348 54694 : ++i;
349 381578 : while (i < 3)
350 272510 : std::swap(array[i++], n_id);
351 99984 : };
352 :
353 20 : WrappedNgMesh ngmesh;
354 :
355 : // Create surface mesh in the WrappedNgMesh
356 : {
357 : // NetGen appears to use ONE-BASED numbering for its nodes, and
358 : // since it doesn't return an id when adding nodes we'll have to
359 : // track the numbering ourselves.
360 120 : int ng_id = 1;
361 :
362 : auto create_surface_component =
363 120 : [this, &ng_id, &ng_to_libmesh_id, &ngmesh, &side_boundary_id, &insert_id]
364 : (UnstructuredMesh & srcmesh, bool hole_mesh,
365 9114 : boundary_id_type bcid)
366 : {
367 24 : LOG_SCOPE("create_surface_component()", "NetGenMeshInterface");
368 :
369 : // Keep track of what nodes we've already added to the Netgen
370 : // mesh vs what nodes we need to add. We'll keep track by id,
371 : // not by point location. I don't know if Netgen can handle
372 : // multiple nodes with the same point location, but if they can
373 : // it's not going to be *us* who breaks that feature.
374 24 : std::unordered_map<dof_id_type, int> libmesh_to_ng_id;
375 :
376 : // Keep track of what nodes we've already added to the main
377 : // mesh from a hole mesh.
378 24 : std::unordered_map<dof_id_type, dof_id_type> hole_to_main_mesh_id;
379 :
380 : // Use a separate array for passing points to NetGen, just in case
381 : // we're not using double-precision ourselves.
382 : std::array<double, 3> point_val;
383 :
384 : // And an array for element vertices
385 : std::array<int, 3> elem_nodes;
386 :
387 10572 : for (const auto * elem : srcmesh.element_ptr_range())
388 : {
389 : // If someone has triangles we can't triangulate, we have a
390 : // problem
391 11232 : if (elem->type() == TRI6 ||
392 5616 : elem->type() == TRI7)
393 0 : libmesh_not_implemented_msg
394 : ("Netgen tetrahedralization currently only supports TRI3 boundaries");
395 :
396 : // If someone has non-triangles, let's just ignore them.
397 5616 : if (elem->type() != TRI3)
398 0 : continue;
399 :
400 5616 : std::array<dof_id_type,3> sorted_ids =
401 : {DofObject::invalid_id, DofObject::invalid_id,
402 : DofObject::invalid_id};
403 :
404 22464 : for (int ni : make_range(3))
405 : {
406 : // Just using the "invert_trigs" option in NetGen params
407 : // doesn't work for me, so we'll have to have properly
408 : // oriented the tris earlier.
409 16848 : auto & elem_node = hole_mesh ? elem_nodes[2-ni] : elem_nodes[ni];
410 :
411 2808 : const Node & n = elem->node_ref(ni);
412 16848 : auto n_id = n.id();
413 16848 : if (hole_mesh)
414 : {
415 7200 : if (auto it = hole_to_main_mesh_id.find(n_id);
416 600 : it != hole_to_main_mesh_id.end())
417 : {
418 5952 : n_id = it->second;
419 : }
420 : else
421 : {
422 1248 : Node * n_new = this->_mesh.add_point(n);
423 1248 : const dof_id_type n_new_id = n_new->id();
424 104 : hole_to_main_mesh_id.emplace(n_id, n_new_id);
425 1248 : n_id = n_new_id;
426 : }
427 : }
428 :
429 16848 : if (auto it = libmesh_to_ng_id.find(n_id);
430 1404 : it != libmesh_to_ng_id.end())
431 : {
432 13752 : const int existing_ng_id = it->second;
433 13752 : elem_node = existing_ng_id;
434 : }
435 : else
436 : {
437 12384 : for (auto i : make_range(3))
438 9288 : point_val[i] = double(n(i));
439 :
440 3096 : Ng_AddPoint(ngmesh, point_val.data());
441 :
442 3096 : ng_to_libmesh_id[ng_id] = n_id;
443 3096 : libmesh_to_ng_id[n_id] = ng_id;
444 3096 : elem_node = ng_id;
445 3096 : ++ng_id;
446 : }
447 :
448 16848 : insert_id(sorted_ids, n_id);
449 : }
450 :
451 5616 : side_boundary_id[sorted_ids] = bcid;
452 :
453 5616 : Ng_AddSurfaceElement(ngmesh, NG_TRIG, elem_nodes.data());
454 120 : }
455 144 : };
456 :
457 : // Number the outer boundary 0, and the holes starting from 1
458 10 : boundary_id_type bcid = 0;
459 :
460 120 : create_surface_component(this->_mesh, false, bcid);
461 :
462 120 : if (_holes)
463 48 : for (const std::unique_ptr<UnstructuredMesh> & h : *_holes)
464 26 : create_surface_component(*h, true, ++bcid);
465 : }
466 :
467 : {
468 20 : LOG_SCOPE("Ng_GenerateVolumeMesh()", "NetGenMeshInterface");
469 :
470 120 : auto result = Ng_GenerateVolumeMesh(ngmesh, ¶ms);
471 120 : handle_ng_result(result);
472 : }
473 :
474 120 : const int n_elem = Ng_GetNE(ngmesh);
475 :
476 : // If Netgen fails us, we're likely to get n_elem <= 0. This is a
477 : // common enough failure from bad setups that I want to make sure
478 : // it's thrown in parallel so as to not desynchronize any unit tests
479 : // that trigger it. So we'll broadcast the empty mesh to indicate
480 : // the problem and enable throwing exceptions in parallel.
481 120 : if (n_elem <= 0)
482 : {
483 0 : this->_mesh.clear();
484 0 : MeshCommunication().broadcast(this->_mesh);
485 0 : libmesh_error_msg ("NetGen failed to generate any tetrahedra");
486 : }
487 :
488 120 : const dof_id_type n_points = Ng_GetNP(ngmesh);
489 120 : const dof_id_type old_nodes = this->_mesh.n_nodes();
490 :
491 : // Netgen may have generated new interior nodes
492 120 : if (n_points != old_nodes)
493 : {
494 : std::array<double, 3> point_val;
495 :
496 : // We should only be getting new nodes if we asked for them
497 1 : if (!_desired_volume)
498 : {
499 : std::cout <<
500 0 : "NetGen output " << n_points <<
501 0 : " points when we gave it " <<
502 0 : old_nodes << " and disabled refinement\n" <<
503 : "If new interior points are acceptable in your mesh, please set\n" <<
504 : "a non-zero desired_volume to indicate that. If new interior\n" <<
505 : "points are not acceptable in your mesh, you may need a different\n" <<
506 0 : "(non-advancing-front?) mesh generator." << std::endl;
507 0 : libmesh_error();
508 : }
509 : else
510 2 : for (auto i : make_range(old_nodes, n_points))
511 : {
512 : // i+1 since ng uses ONE-BASED numbering
513 1 : Ng_GetPoint (ngmesh, i+1, point_val.data());
514 1 : const Point p(point_val[0], point_val[1], point_val[2]);
515 1 : Node * n_new = this->_mesh.add_point(p);
516 0 : const dof_id_type n_new_id = n_new->id();
517 1 : ng_to_libmesh_id[i+1] = n_new_id;
518 : }
519 : }
520 :
521 6126 : for (auto * elem : this->_mesh.element_ptr_range())
522 3316 : this->_mesh.delete_elem(elem);
523 :
524 120 : BoundaryInfo * bi = & this->_mesh.get_boundary_info();
525 :
526 7805 : for (auto i : make_range(n_elem))
527 : {
528 : // Enough data to return even a Tet10 without a segfault if nglib
529 : // went nuts
530 : int ngnodes[11];
531 :
532 : // i+1 since we must be 1-based with these ids too...
533 : Ng_Volume_Element_Type ngtype =
534 7685 : Ng_GetVolumeElement(ngmesh, i+1, ngnodes);
535 :
536 : // But really nglib shouldn't go nuts
537 640 : libmesh_assert(ngtype == NG_TET);
538 640 : libmesh_ignore(ngtype);
539 :
540 8325 : auto elem = this->_mesh.add_elem(Elem::build_with_id(TET4, i));
541 38425 : for (auto n : make_range(4))
542 : {
543 : const dof_id_type node_id =
544 30740 : libmesh_map_find(ng_to_libmesh_id, ngnodes[n]);
545 30740 : elem->set_node(n, this->_mesh.node_ptr(node_id));
546 : }
547 :
548 : // NetGen and we disagree about node numbering orientation
549 7045 : elem->orient(bi);
550 :
551 38425 : for (auto s : make_range(4))
552 : {
553 30740 : std::array<dof_id_type,3> sorted_ids =
554 : {DofObject::invalid_id, DofObject::invalid_id,
555 : DofObject::invalid_id};
556 :
557 33300 : std::vector<unsigned int> nos = elem->nodes_on_side(s);
558 122960 : for (auto n : nos)
559 92220 : insert_id(sorted_ids, elem->node_id(n));
560 :
561 30740 : if (auto it = side_boundary_id.find(sorted_ids);
562 2560 : it != side_boundary_id.end())
563 5616 : bi->add_side(elem, s, it->second);
564 : }
565 : }
566 :
567 : // We don't need our holes anymore. Delete their serializers
568 : // first to avoid dereferencing dangling pointers.
569 10 : hole_serializers.clear();
570 120 : if (_holes)
571 2 : _holes->clear();
572 :
573 : // Send NetGen's TET4 result to the other ranks, then increase the element
574 : // order collectively on all ranks together. increase_tet_order() performs
575 : // collective communication, so it must not run on rank 0 alone -- doing so
576 : // would deadlock against the other ranks waiting in broadcast() above.
577 120 : MeshCommunication().broadcast(this->_mesh);
578 120 : increase_order_and_restore_midpoints();
579 120 : this->_mesh.prepare_for_use();
580 670 : }
581 :
582 :
583 :
584 : } // namespace libMesh
585 :
586 :
587 : #endif // #ifdef LIBMESH_HAVE_NETGEN
|
