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 :
20 :
21 : // C++ includes
22 : #include <sstream>
23 :
24 : // Local includes
25 : #include "libmesh/mesh_tet_interface.h"
26 :
27 : #include "libmesh/boundary_info.h"
28 : #include "libmesh/elem.h"
29 : #include "libmesh/mesh_modification.h"
30 : #include "libmesh/mesh_serializer.h"
31 : #include "libmesh/mesh_tools.h"
32 : #include "libmesh/remote_elem.h"
33 : #include "libmesh/unstructured_mesh.h"
34 :
35 : namespace {
36 : using namespace libMesh;
37 :
38 : std::unordered_set<Elem *>
39 568 : flood_component (std::unordered_set<Elem *> & all_components,
40 : Elem * elem)
41 : {
42 16 : libmesh_assert(!all_components.count(elem));
43 :
44 16 : std::unordered_set<Elem *> current_component;
45 :
46 669 : std::unordered_set<Elem *> elements_to_consider = {elem};
47 :
48 4533 : while (!elements_to_consider.empty())
49 : {
50 218 : std::unordered_set<Elem *> next_elements_to_consider;
51 :
52 42196 : for (Elem * considering : elements_to_consider)
53 : {
54 994 : all_components.insert(considering);
55 994 : current_component.insert(considering);
56 :
57 152782 : for (auto s : make_range(considering->n_sides()))
58 : {
59 114622 : Elem * neigh = considering->neighbor_ptr(s);
60 :
61 114693 : libmesh_error_msg_if
62 : (!neigh,
63 : "Tet generation encountered a 2D element with a null neighbor, but a\n"
64 : "boundary must be a 2D closed manifold (surface).\n");
65 :
66 114551 : if (all_components.find(neigh) == all_components.end())
67 1340 : next_elements_to_consider.insert(neigh);
68 : }
69 : }
70 107 : elements_to_consider = next_elements_to_consider;
71 : }
72 :
73 511 : return current_component;
74 : }
75 :
76 : // Returns six times the signed volume of a tet formed by the given
77 : // 3 points and the origin
78 30388 : Real six_times_signed_tet_volume (const Point & p1,
79 : const Point & p2,
80 : const Point & p3)
81 : {
82 30388 : return p1(0)*p2(1)*p3(2)
83 30388 : - p1(0)*p3(1)*p2(2)
84 30388 : - p2(0)*p1(1)*p3(2)
85 30388 : + p2(0)*p3(1)*p1(2)
86 30388 : + p3(0)*p1(1)*p2(2)
87 30388 : - p3(0)*p2(1)*p1(2);
88 : }
89 :
90 : // Returns six times the signed volume of the space defined by the
91 : // manifold of surface elements in component
92 497 : Real six_times_signed_volume (const std::unordered_set<Elem *> component)
93 : {
94 14 : Real six_vol = 0;
95 :
96 30885 : for (const Elem * elem: component)
97 : {
98 856 : libmesh_assert_equal_to(elem->dim(), 2);
99 60776 : for (auto n : make_range(elem->n_vertices()-2))
100 32100 : six_vol += six_times_signed_tet_volume(elem->point(0),
101 : elem->point(n+1),
102 : elem->point(n+2));
103 : }
104 :
105 497 : return six_vol;
106 : }
107 : }
108 :
109 : namespace libMesh
110 : {
111 :
112 : //----------------------------------------------------------------------
113 : // MeshTetInterface class members
114 454 : MeshTetInterface::MeshTetInterface (UnstructuredMesh & mesh) :
115 402 : _desired_volume(0), _smooth_after_generating(false),
116 454 : _elem_type(TET4), _mesh(mesh)
117 : {
118 454 : }
119 :
120 :
121 402 : MeshTetInterface::~MeshTetInterface() = default;
122 :
123 :
124 142 : void MeshTetInterface::attach_hole_list
125 : (std::unique_ptr<std::vector<std::unique_ptr<UnstructuredMesh>>> holes)
126 : {
127 4 : _holes = std::move(holes);
128 142 : }
129 :
130 :
131 568 : BoundingBox MeshTetInterface::volume_to_surface_mesh(UnstructuredMesh & mesh)
132 : {
133 : // If we've been handed an unprepared mesh then we need to be made
134 : // aware of that and fix that; we're relying on neighbor pointers.
135 16 : libmesh_assert(MeshTools::valid_is_prepared(mesh));
136 :
137 568 : if (!mesh.is_prepared())
138 0 : mesh.prepare_for_use();
139 :
140 : // We'll return a bounding box for use by subclasses in basic sanity checks.
141 552 : BoundingBox surface_bb;
142 :
143 : // First convert all volume boundaries to surface elements; this
144 : // gives us a manifold bounding the mesh, though it may not be a
145 : // connected manifold even if the volume mesh was connected.
146 : {
147 : // Make sure ids are in sync and valid on a DistributedMesh
148 568 : const dof_id_type max_orig_id = mesh.max_elem_id();
149 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
150 568 : const unique_id_type max_unique_id = mesh.parallel_max_unique_id();
151 : #endif
152 :
153 : // Change this if we add arbitrary polyhedra...
154 16 : const dof_id_type max_sides = 6;
155 :
156 32 : std::unordered_set<Elem *> elems_to_delete;
157 :
158 48 : std::vector<std::unique_ptr<Elem>> elems_to_add;
159 :
160 : // Convert all faces to surface elements
161 508546 : for (auto * elem : mesh.active_element_ptr_range())
162 : {
163 264537 : libmesh_error_msg_if (elem->dim() < 2,
164 : "Cannot use meshes with 0D or 1D elements to define a volume");
165 :
166 : // If we've already got 2D elements then those are (part of)
167 : // our surface.
168 264537 : if (elem->dim() == 2)
169 1971 : continue;
170 :
171 : // 3D elements will be removed after we've extracted their
172 : // surface faces.
173 10762 : elems_to_delete.insert(elem);
174 :
175 1312662 : for (auto s : make_range(elem->n_sides()))
176 : {
177 : // If there's a neighbor on this side then there's not a
178 : // boundary
179 1093210 : if (elem->neighbor_ptr(s))
180 : {
181 : // We're not supporting AMR meshes here yet
182 1030584 : if (elem->level() != elem->neighbor_ptr(s)->level())
183 0 : libmesh_not_implemented_msg
184 : ("Tetrahedralizaton of adapted meshes is not currently supported");
185 988328 : continue;
186 946072 : }
187 :
188 38352 : elems_to_add.push_back(elem->build_side_ptr(s));
189 796 : Elem * side_elem = elems_to_add.back().get();
190 :
191 : // Wipe the interior_parent before it can become a
192 : // dangling pointer later
193 19574 : side_elem->set_interior_parent(nullptr);
194 :
195 : // If the mesh is replicated then its automatic id
196 : // setting is fine. If not, then we need unambiguous ids
197 : // independent of element traversal.
198 19574 : if (!mesh.is_replicated())
199 : {
200 16788 : side_elem->set_id(max_orig_id + max_sides*elem->id() + s);
201 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
202 0 : side_elem->set_unique_id(max_unique_id + max_sides*elem->id() + s);
203 : #endif
204 : }
205 : }
206 536 : }
207 :
208 : // If the mesh is replicated then its automatic neighbor finding
209 : // is fine. If not, then we need to insert them ourselves, but
210 : // it's easy because we can use the fact (from our implementation
211 : // above) that our new elements have no parents or children, plus
212 : // the fact (from the tiny fraction of homology I understand) that
213 : // a manifold boundary is a manifold with no boundary.
214 : //
215 : // See UnstructuredMesh::find_neighbors() for more explanation of
216 : // (a more complicated version of) the algorithm here.
217 568 : if (!mesh.is_replicated())
218 : {
219 : typedef dof_id_type key_type;
220 : typedef std::pair<Elem *, unsigned char> val_type;
221 : typedef std::unordered_multimap<key_type, val_type> map_type;
222 0 : map_type side_to_elem_map;
223 :
224 512 : std::unique_ptr<Elem> my_side, their_side;
225 :
226 17300 : for (auto & elem : elems_to_add)
227 : {
228 67536 : for (auto s : elem->side_index_range())
229 : {
230 50748 : if (elem->neighbor_ptr(s))
231 0 : continue;
232 50748 : const dof_id_type key = elem->low_order_key(s);
233 0 : auto bounds = side_to_elem_map.equal_range(key);
234 50748 : if (bounds.first != bounds.second)
235 : {
236 0 : elem->side_ptr(my_side, s);
237 0 : while (bounds.first != bounds.second)
238 : {
239 0 : Elem * potential_neighbor = bounds.first->second.first;
240 0 : const unsigned int ns = bounds.first->second.second;
241 0 : potential_neighbor->side_ptr(their_side, ns);
242 0 : if (*my_side == *their_side)
243 : {
244 0 : elem->set_neighbor(s, potential_neighbor);
245 0 : potential_neighbor->set_neighbor(ns, elem.get());
246 0 : side_to_elem_map.erase (bounds.first);
247 0 : break;
248 : }
249 0 : ++bounds.first;
250 : }
251 :
252 0 : if (!elem->neighbor_ptr(s))
253 : side_to_elem_map.emplace
254 0 : (key, std::make_pair(elem.get(), cast_int<unsigned char>(s)));
255 : }
256 : }
257 : }
258 :
259 : // At this point we *should* have a match for everything, so
260 : // anything we don't have a match for is remote.
261 17300 : for (auto & elem : elems_to_add)
262 67536 : for (auto s : elem->side_index_range())
263 50748 : if (!elem->neighbor_ptr(s))
264 50748 : elem->set_neighbor(s, const_cast<RemoteElem*>(remote_elem));
265 512 : }
266 :
267 : // Remove volume and edge elements
268 263072 : for (Elem * elem : elems_to_delete)
269 262504 : mesh.delete_elem(elem);
270 :
271 : // Add the new elements outside the loop so we don't risk
272 : // invalidating iterators.
273 20142 : for (auto & elem : elems_to_add)
274 21166 : mesh.add_elem(std::move(elem));
275 536 : }
276 :
277 : // Fix up neighbor pointers, element counts, etc.
278 568 : mesh.prepare_for_use();
279 :
280 : // We're making tets; we need to start with tris
281 568 : MeshTools::Modification::all_tri(mesh);
282 :
283 : // Partition surface into connected components. At this point I'm
284 : // finally going to give up and serialize, because at least we got
285 : // from 3D down to 2D first, and because I don't want to have to
286 : // turn flood_component into a while loop with a parallel sync in
287 : // the middle, and because we do have to serialize *eventually*
288 : // anyways unless we get a parallel tetrahedralizer backend someday.
289 600 : MeshSerializer mesh_serializer(mesh);
290 :
291 48 : std::vector<std::unordered_set<Elem *>> components;
292 32 : std::unordered_set<Elem *> in_component;
293 :
294 89791 : for (auto * elem : mesh.element_ptr_range())
295 1716 : if (!in_component.count(elem))
296 1587 : components.emplace_back(flood_component(in_component, elem));
297 :
298 14 : const std::unordered_set<Elem *> * biggest_component = nullptr;
299 14 : Real biggest_six_vol = 0;
300 994 : for (const auto & component : components)
301 : {
302 1049 : Real six_vol = six_times_signed_volume(component);
303 497 : if (std::abs(six_vol) > std::abs(biggest_six_vol))
304 : {
305 14 : biggest_six_vol = six_vol;
306 14 : biggest_component = &component;
307 : }
308 : }
309 :
310 497 : if (!biggest_component)
311 0 : libmesh_error_msg("No non-zero-volume component found among " <<
312 : components.size() << " boundary components");
313 :
314 994 : for (const auto & component : components)
315 497 : if (&component != biggest_component)
316 : {
317 0 : for (Elem * elem: component)
318 0 : mesh.delete_elem(elem);
319 : }
320 : else
321 : {
322 30885 : for (Elem * elem: component)
323 : {
324 30388 : if (biggest_six_vol < 0)
325 0 : elem->flip(&mesh.get_boundary_info());
326 :
327 121552 : for (auto & node : elem->node_ref_range())
328 88596 : surface_bb.union_with(node);
329 : }
330 : }
331 :
332 497 : mesh.prepare_for_use();
333 :
334 511 : return surface_bb;
335 603 : }
336 :
337 :
338 64 : unsigned MeshTetInterface::check_hull_integrity()
339 : {
340 : // Check for easy return: if the Mesh is empty (i.e. if
341 : // somebody called triangulate_conformingDelaunayMesh on
342 : // a Mesh with no elements, then hull integrity check must
343 : // fail...
344 64 : if (_mesh.n_elem() == 0)
345 0 : return 3;
346 :
347 5233 : for (auto & elem : this->_mesh.element_ptr_range())
348 : {
349 : // Check for proper element type
350 2832 : if (elem->type() != TRI3)
351 : {
352 : //libmesh_error_msg("ERROR: Some of the elements in the original mesh were not TRI3!");
353 0 : return 1;
354 : }
355 :
356 11604 : for (auto neigh : elem->neighbor_ptr_range())
357 : {
358 8496 : if (neigh == nullptr)
359 : {
360 : // libmesh_error_msg("ERROR: Non-convex hull, cannot be tetrahedralized.");
361 0 : return 2;
362 : }
363 : }
364 50 : }
365 :
366 : // If we made it here, return success!
367 64 : return 0;
368 : }
369 :
370 :
371 :
372 4 : void MeshTetInterface::process_hull_integrity_result(unsigned result)
373 : {
374 4 : if (result != 0)
375 : {
376 0 : libMesh::err << "Error! Conforming Delaunay mesh tetrahedralization requires a convex hull." << std::endl;
377 :
378 0 : if (result==1)
379 : {
380 0 : libMesh::err << "Non-TRI3 elements were found in the input Mesh. ";
381 0 : libMesh::err << "A constrained Delaunay triangulation requires a convex hull of TRI3 elements." << std::endl;
382 : }
383 :
384 0 : libmesh_error_msg("Consider calling TetGenMeshInterface::pointset_convexhull() followed by Mesh::find_neighbors() first.");
385 : }
386 4 : }
387 :
388 :
389 :
390 4 : void MeshTetInterface::delete_2D_hull_elements()
391 : {
392 10078 : for (auto & elem : this->_mesh.element_ptr_range())
393 : {
394 : // Check for proper element type. Yes, we legally delete elements while
395 : // iterating over them because no entries from the underlying container
396 : // are actually erased.
397 10072 : if (elem->type() == TRI3)
398 96 : _mesh.delete_elem(elem);
399 0 : }
400 :
401 : // We just removed any boundary info associated with hull element
402 : // edges, so let's update the boundary id caches.
403 6 : this->_mesh.get_boundary_info().regenerate_id_sets();
404 4 : }
405 :
406 :
407 :
408 : } // namespace libMesh
|
