18 #include "libmesh/libmesh_config.h" 19 #ifdef LIBMESH_HAVE_NETGEN 26 #include "libmesh/mesh_netgen_interface.h" 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" 37 #include "netgen/nglib/nglib.h" 47 _ngmesh = nglib::Ng_NewMesh();
51 nglib::Ng_DeleteMesh(_ngmesh);
55 nglib::Ng_DeleteMesh(_ngmesh);
56 _ngmesh = nglib::Ng_NewMesh();
59 operator nglib::Ng_Mesh* () {
64 nglib::Ng_Mesh * _ngmesh;
84 using namespace nglib;
86 LOG_SCOPE(
"triangulate()",
"NetGenMeshInterface");
95 std::vector<MeshSerializer> hole_serializers;
97 for (std::unique_ptr<UnstructuredMesh> & hole : *
_holes)
104 "Found hole with bounding box " << hole_bb <<
105 "\nextending outside of mesh bounding box " << mesh_bb);
107 hole_serializers.emplace_back
119 hole_serializers.clear();
132 Ng_Meshing_Parameters params;
136 params.uselocalh =
false;
138 params.elementsperedge = 1;
139 params.elementspercurve = 1;
140 params.closeedgeenable =
false;
141 params.closeedgefact = 0;
142 params.minedgelenenable =
false;
143 params.minedgelen = 0;
153 params.maxh = std::numeric_limits<double>::max();
159 params.optsteps_3d = 0;
166 std::unordered_map<int, dof_id_type> ng_to_libmesh_id;
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"};
175 std::size_t(result+1) < result_types.size())
177 (result,
"Ng_GenerateVolumeMesh failed: " <<
178 result_types[result+1]);
181 (
"Ng_GenerateVolumeMesh failed with an unknown error code");
188 std::unordered_map<std::array<dof_id_type,3>,
192 (std::array<dof_id_type,3> & array,
197 while (array[i] < n_id)
200 std::swap(array[i++], n_id);
203 WrappedNgMesh ngmesh;
212 auto create_surface_component =
213 [
this, &ng_id, &ng_to_libmesh_id, &ngmesh, &side_boundary_id, &insert_id]
217 LOG_SCOPE(
"create_surface_component()",
"NetGenMeshInterface");
224 std::unordered_map<dof_id_type, int> libmesh_to_ng_id;
228 std::unordered_map<dof_id_type, dof_id_type> hole_to_main_mesh_id;
232 std::array<double, 3> point_val;
235 std::array<int, 3> elem_nodes;
237 for (
const auto * elem : srcmesh.element_ptr_range())
241 if (elem->type() ==
TRI6 ||
242 elem->type() ==
TRI7)
243 libmesh_not_implemented_msg
244 (
"Netgen tetrahedralization currently only supports TRI3 boundaries");
247 if (elem->type() !=
TRI3)
250 std::array<dof_id_type,3> sorted_ids =
259 auto & elem_node = hole_mesh ? elem_nodes[2-ni] : elem_nodes[ni];
261 const Node & n = elem->node_ref(ni);
265 if (
auto it = hole_to_main_mesh_id.find(n_id);
266 it != hole_to_main_mesh_id.end())
274 hole_to_main_mesh_id.emplace(n_id, n_new_id);
279 if (
auto it = libmesh_to_ng_id.find(n_id);
280 it != libmesh_to_ng_id.end())
282 const int existing_ng_id = it->second;
283 elem_node = existing_ng_id;
288 point_val[i] =
double(n(i));
290 Ng_AddPoint(ngmesh, point_val.data());
292 ng_to_libmesh_id[ng_id] = n_id;
293 libmesh_to_ng_id[n_id] = ng_id;
298 insert_id(sorted_ids, n_id);
301 side_boundary_id[sorted_ids] = bcid;
303 Ng_AddSurfaceElement(ngmesh, NG_TRIG, elem_nodes.data());
310 create_surface_component(this->
_mesh,
false, bcid);
313 for (
const std::unique_ptr<UnstructuredMesh> & h : *
_holes)
314 create_surface_component(*h,
true, ++bcid);
318 LOG_SCOPE(
"Ng_GenerateVolumeMesh()",
"NetGenMeshInterface");
320 auto result = Ng_GenerateVolumeMesh(ngmesh, ¶ms);
321 handle_ng_result(result);
324 const int n_elem = Ng_GetNE(ngmesh);
326 libmesh_error_msg_if (
n_elem <= 0,
327 "NetGen failed to generate any tetrahedra");
333 if (n_points != old_nodes)
335 std::array<double, 3> point_val;
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;
351 for (
auto i :
make_range(old_nodes, n_points))
354 Ng_GetPoint (ngmesh, i+1, point_val.data());
355 const Point p(point_val[0], point_val[1], point_val[2]);
358 ng_to_libmesh_id[i+1] = n_new_id;
362 for (
auto * elem : this->
_mesh.element_ptr_range())
374 Ng_Volume_Element_Type ngtype =
375 Ng_GetVolumeElement(ngmesh, i+1, ngnodes);
385 libmesh_map_find(ng_to_libmesh_id, ngnodes[n]);
394 std::array<dof_id_type,3> sorted_ids =
398 std::vector<unsigned int> nos = elem->nodes_on_side(s);
400 insert_id(sorted_ids, elem->node_id(n));
402 if (
auto it = side_boundary_id.find(sorted_ids);
403 it != side_boundary_id.end())
410 hole_serializers.clear();
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.
unsigned check_hull_integrity()
This function checks the integrity of the current set of elements in the Mesh to see if they comprise...
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.
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.
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 &...)
virtual void delete_elem(Elem *e)=0
Removes element e from the mesh.
This is the MeshCommunication class.
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...
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
static BoundingBox volume_to_surface_mesh(UnstructuredMesh &mesh)
Remove volume elements from the given mesh, after converting their outer boundary faces to surface el...
Defines a Cartesian bounding box by the two corner extremum.
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's ...
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...
bool contains(const BoundingBox &) const
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.
virtual dof_id_type n_nodes() const =0
Class MeshTetInterface provides an abstract interface for tetrahedralization of meshes by subclasses...