libMesh
mesh_netgen_interface.C
Go to the documentation of this file.
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  WrappedNgMesh() {
47  _ngmesh = nglib::Ng_NewMesh();
48  }
49 
50  ~WrappedNgMesh() {
51  nglib::Ng_DeleteMesh(_ngmesh);
52  }
53 
54  void clear() {
55  nglib::Ng_DeleteMesh(_ngmesh);
56  _ngmesh = nglib::Ng_NewMesh();
57  }
58 
59  operator nglib::Ng_Mesh* () {
60  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
76  _serializer(mesh)
77 {
78 }
79 
80 
81 
83 {
84  using namespace nglib;
85 
86  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 =
94 
95  std::vector<MeshSerializer> hole_serializers;
96  if (_holes)
97  for (std::unique_ptr<UnstructuredMesh> & hole : *_holes)
98  {
99  const BoundingBox hole_bb =
101 
102  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  (*hole, /* need_serial */ true,
109  /* 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  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  hole_serializers.clear();
120  if (_holes)
121  _holes->clear();
122 
123  // Receive the mesh data rank 0 will send later, then fix it up
124  // together
126  this->_mesh.prepare_for_use();
127  return;
128  }
129 
130  this->check_hull_integrity();
131 
132  Ng_Meshing_Parameters params;
133 
134  // Override any default parameters we might need to, to avoid
135  // inserting nodes we don't want.
136  params.uselocalh = false;
137  params.minh = 0;
138  params.elementsperedge = 1;
139  params.elementspercurve = 1;
140  params.closeedgeenable = false;
141  params.closeedgefact = 0;
142  params.minedgelenenable = false;
143  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  if (_desired_volume == 0) // shorthand for "no refinement"
152  {
153  params.maxh = std::numeric_limits<double>::max();
154  params.fineness = 0; // "coarse" in the docs
155  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  params.optsteps_3d = 0;
160  }
161  else
162  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  std::unordered_map<int, dof_id_type> ng_to_libmesh_id;
167 
168  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  "Netgen surface failure", "Netgen file not found"};
173 
174  if (result+1 >= 0 &&
175  std::size_t(result+1) < result_types.size())
176  libmesh_error_msg_if
177  (result, "Ng_GenerateVolumeMesh failed: " <<
178  result_types[result+1]);
179  else
180  libmesh_error_msg
181  ("Ng_GenerateVolumeMesh failed with an unknown error code");
182  };
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  boundary_id_type, libMesh::hash> side_boundary_id;
190 
191  auto insert_id = []
192  (std::array<dof_id_type,3> & array,
193  dof_id_type n_id)
194  {
195  libmesh_assert_less(n_id, DofObject::invalid_id);
196  unsigned int i=0;
197  while (array[i] < n_id)
198  ++i;
199  while (i < 3)
200  std::swap(array[i++], n_id);
201  };
202 
203  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  int ng_id = 1;
211 
212  auto create_surface_component =
213  [this, &ng_id, &ng_to_libmesh_id, &ngmesh, &side_boundary_id, &insert_id]
214  (UnstructuredMesh & srcmesh, bool hole_mesh,
215  boundary_id_type bcid)
216  {
217  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  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  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  for (const auto * elem : srcmesh.element_ptr_range())
238  {
239  // If someone has triangles we can't triangulate, we have a
240  // problem
241  if (elem->type() == TRI6 ||
242  elem->type() == TRI7)
243  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  if (elem->type() != TRI3)
248  continue;
249 
250  std::array<dof_id_type,3> sorted_ids =
253 
254  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  auto & elem_node = hole_mesh ? elem_nodes[2-ni] : elem_nodes[ni];
260 
261  const Node & n = elem->node_ref(ni);
262  auto n_id = n.id();
263  if (hole_mesh)
264  {
265  if (auto it = hole_to_main_mesh_id.find(n_id);
266  it != hole_to_main_mesh_id.end())
267  {
268  n_id = it->second;
269  }
270  else
271  {
272  Node * n_new = this->_mesh.add_point(n);
273  const dof_id_type n_new_id = n_new->id();
274  hole_to_main_mesh_id.emplace(n_id, n_new_id);
275  n_id = n_new_id;
276  }
277  }
278 
279  if (auto it = libmesh_to_ng_id.find(n_id);
280  it != libmesh_to_ng_id.end())
281  {
282  const int existing_ng_id = it->second;
283  elem_node = existing_ng_id;
284  }
285  else
286  {
287  for (auto i : make_range(3))
288  point_val[i] = double(n(i));
289 
290  Ng_AddPoint(ngmesh, point_val.data());
291 
292  ng_to_libmesh_id[ng_id] = n_id;
293  libmesh_to_ng_id[n_id] = ng_id;
294  elem_node = ng_id;
295  ++ng_id;
296  }
297 
298  insert_id(sorted_ids, n_id);
299  }
300 
301  side_boundary_id[sorted_ids] = bcid;
302 
303  Ng_AddSurfaceElement(ngmesh, NG_TRIG, elem_nodes.data());
304  }
305  };
306 
307  // Number the outer boundary 0, and the holes starting from 1
308  boundary_id_type bcid = 0;
309 
310  create_surface_component(this->_mesh, false, bcid);
311 
312  if (_holes)
313  for (const std::unique_ptr<UnstructuredMesh> & h : *_holes)
314  create_surface_component(*h, true, ++bcid);
315  }
316 
317  {
318  LOG_SCOPE("Ng_GenerateVolumeMesh()", "NetGenMeshInterface");
319 
320  auto result = Ng_GenerateVolumeMesh(ngmesh, &params);
321  handle_ng_result(result);
322  }
323 
324  const int n_elem = Ng_GetNE(ngmesh);
325 
326  libmesh_error_msg_if (n_elem <= 0,
327  "NetGen failed to generate any tetrahedra");
328 
329  const dof_id_type n_points = Ng_GetNP(ngmesh);
330  const dof_id_type old_nodes = this->_mesh.n_nodes();
331 
332  // Netgen may have generated new interior nodes
333  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  if (!_desired_volume)
339  {
340  std::cout <<
341  "NetGen output " << n_points <<
342  " points when we gave it " <<
343  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  "(non-advancing-front?) mesh generator." << std::endl;
348  libmesh_error();
349  }
350  else
351  for (auto i : make_range(old_nodes, n_points))
352  {
353  // i+1 since ng uses ONE-BASED numbering
354  Ng_GetPoint (ngmesh, i+1, point_val.data());
355  const Point p(point_val[0], point_val[1], point_val[2]);
356  Node * n_new = this->_mesh.add_point(p);
357  const dof_id_type n_new_id = n_new->id();
358  ng_to_libmesh_id[i+1] = n_new_id;
359  }
360  }
361 
362  for (auto * elem : this->_mesh.element_ptr_range())
363  this->_mesh.delete_elem(elem);
364 
365  BoundaryInfo * bi = & this->_mesh.get_boundary_info();
366 
367  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  Ng_GetVolumeElement(ngmesh, i+1, ngnodes);
376 
377  // But really nglib shouldn't go nuts
378  libmesh_assert(ngtype == NG_TET);
379  libmesh_ignore(ngtype);
380 
381  auto elem = this->_mesh.add_elem(Elem::build_with_id(TET4, i));
382  for (auto n : make_range(4))
383  {
384  const dof_id_type node_id =
385  libmesh_map_find(ng_to_libmesh_id, ngnodes[n]);
386  elem->set_node(n, this->_mesh.node_ptr(node_id));
387  }
388 
389  // NetGen and we disagree about node numbering orientation
390  elem->orient(bi);
391 
392  for (auto s : make_range(4))
393  {
394  std::array<dof_id_type,3> sorted_ids =
397 
398  std::vector<unsigned int> nos = elem->nodes_on_side(s);
399  for (auto n : nos)
400  insert_id(sorted_ids, elem->node_id(n));
401 
402  if (auto it = side_boundary_id.find(sorted_ids);
403  it != side_boundary_id.end())
404  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  hole_serializers.clear();
411  if (_holes)
412  _holes->clear();
413 
414  // Send our data to other ranks, then fix it up together
416  this->_mesh.prepare_for_use();
417 }
418 
419 
420 
421 } // namespace libMesh
422 
423 
424 #endif // #ifdef LIBMESH_HAVE_NETGEN
NetGenMeshInterface(UnstructuredMesh &mesh)
Constructor.
std::unique_ptr< std::vector< std::unique_ptr< UnstructuredMesh > > > _holes
A pointer to a vector of meshes each defining a hole.
A Node is like a Point, but with more information.
Definition: node.h:52
unsigned check_hull_integrity()
This function checks the integrity of the current set of elements in the Mesh to see if they comprise...
dof_id_type n_elem(const MeshBase::const_element_iterator &begin, const MeshBase::const_element_iterator &end)
Count up the number of elements of a specific type (as defined by an iterator range).
Definition: mesh_tools.C:969
void prepare_for_use(const bool skip_renumber_nodes_and_elements, const bool skip_find_neighbors)
Prepare a newly ecreated (or read) mesh for use.
Definition: mesh_base.C:759
MeshBase & mesh
virtual void triangulate() override
Method invokes NetGen library to compute a tetrahedralization.
The libMesh namespace provides an interface to certain functionality in the library.
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:165
virtual Node * add_point(const Point &p, const dof_id_type id=DofObject::invalid_id, const processor_id_type proc_id=DofObject::invalid_processor_id)=0
Add a new Node at Point p to the end of the vertex array, with processor_id procid.
void libmesh_ignore(const Args &...)
T pow(const T &x)
Definition: utility.h:328
int8_t boundary_id_type
Definition: id_types.h:51
virtual void delete_elem(Elem *e)=0
Removes element e from the mesh.
This is the MeshCommunication class.
dof_id_type id() const
Definition: dof_object.h:828
The UnstructuredMesh class is derived from the MeshBase class.
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
Definition: boundary_info.h:57
libmesh_assert(ctx)
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
Definition: dof_object.h:482
static BoundingBox volume_to_surface_mesh(UnstructuredMesh &mesh)
Remove volume elements from the given mesh, after converting their outer boundary faces to surface el...
Definition: libmesh.C:90
Defines a Cartesian bounding box by the two corner extremum.
Definition: bounding_box.h:40
Real _desired_volume
The desired volume for the elements in the resulting mesh.
static std::unique_ptr< Elem > build_with_id(const ElemType type, dof_id_type id)
Calls the build() method above with a nullptr parent, and additionally sets the newly-created Elem&#39;s ...
Definition: elem.C:558
void broadcast(MeshBase &) const
Finds all the processors that may contain elements that neighbor my elements.
void add_side(const dof_id_type elem, const unsigned short int side, const boundary_id_type id)
Add side side of element number elem with boundary id id to the boundary information data structure...
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
bool contains(const BoundingBox &) const
Definition: bounding_box.h:209
virtual const Node * node_ptr(const dof_id_type i) const =0
processor_id_type processor_id() const
UnstructuredMesh & _mesh
Local reference to the mesh we are working with.
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
virtual dof_id_type n_nodes() const =0
Class MeshTetInterface provides an abstract interface for tetrahedralization of meshes by subclasses...
uint8_t dof_id_type
Definition: id_types.h:67