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/elem_range.h"
30 : #include "libmesh/mesh_modification.h"
31 : #include "libmesh/mesh_serializer.h"
32 : #include "libmesh/mesh_tools.h"
33 : #include "libmesh/remote_elem.h"
34 : #include "libmesh/unstructured_mesh.h"
35 :
36 : #include "timpi/parallel_implementation.h"
37 :
38 : namespace {
39 : using namespace libMesh;
40 :
41 : std::unordered_set<Elem *>
42 639 : flood_component (std::unordered_set<Elem *> & all_components,
43 : Elem * elem)
44 : {
45 18 : libmesh_assert(!all_components.count(elem));
46 :
47 18 : std::unordered_set<Elem *> current_component;
48 :
49 744 : std::unordered_set<Elem *> elements_to_consider = {elem};
50 :
51 4898 : while (!elements_to_consider.empty())
52 : {
53 236 : std::unordered_set<Elem *> next_elements_to_consider;
54 :
55 43094 : for (Elem * considering : elements_to_consider)
56 : {
57 1003 : all_components.insert(considering);
58 1003 : current_component.insert(considering);
59 :
60 155198 : for (auto s : make_range(considering->n_sides()))
61 : {
62 116434 : Elem * neigh = considering->neighbor_ptr(s);
63 :
64 116505 : libmesh_error_msg_if
65 : (!neigh,
66 : "Tet generation encountered a 2D element with a null neighbor, but a\n"
67 : "boundary must be a 2D closed manifold (surface).\n");
68 :
69 116363 : if (all_components.find(neigh) == all_components.end())
70 1358 : next_elements_to_consider.insert(neigh);
71 : }
72 : }
73 116 : elements_to_consider = next_elements_to_consider;
74 : }
75 :
76 584 : return current_component;
77 : }
78 :
79 : // Returns six times the signed volume of a tet formed by the given
80 : // 3 points and the origin
81 30956 : Real six_times_signed_tet_volume (const Point & p1,
82 : const Point & p2,
83 : const Point & p3)
84 : {
85 30956 : return p1(0)*p2(1)*p3(2)
86 30956 : - p1(0)*p3(1)*p2(2)
87 30956 : - p2(0)*p1(1)*p3(2)
88 30956 : + p2(0)*p3(1)*p1(2)
89 30956 : + p3(0)*p1(1)*p2(2)
90 30956 : - p3(0)*p2(1)*p1(2);
91 : }
92 :
93 : // Returns six times the signed volume of the space defined by the
94 : // manifold of surface elements in component
95 568 : Real six_times_signed_volume (const std::unordered_set<Elem *> component)
96 : {
97 16 : Real six_vol = 0;
98 :
99 31524 : for (const Elem * elem: component)
100 : {
101 872 : libmesh_assert_equal_to(elem->dim(), 2);
102 61912 : for (auto n : make_range(elem->n_vertices()-2))
103 32700 : six_vol += six_times_signed_tet_volume(elem->point(0),
104 : elem->point(n+1),
105 : elem->point(n+2));
106 : }
107 :
108 568 : return six_vol;
109 : }
110 : }
111 :
112 : namespace libMesh
113 : {
114 :
115 : //----------------------------------------------------------------------
116 : // MeshTetInterface class members
117 525 : MeshTetInterface::MeshTetInterface (UnstructuredMesh & mesh) :
118 469 : _verbosity(0), _desired_volume(0), _smooth_after_generating(false),
119 525 : _elem_type(TET4), _mesh(mesh)
120 : {
121 525 : }
122 :
123 :
124 469 : MeshTetInterface::~MeshTetInterface() = default;
125 :
126 :
127 142 : void MeshTetInterface::attach_hole_list
128 : (std::unique_ptr<std::vector<std::unique_ptr<UnstructuredMesh>>> holes)
129 : {
130 4 : _holes = std::move(holes);
131 142 : }
132 :
133 :
134 639 : BoundingBox MeshTetInterface::volume_to_surface_mesh(UnstructuredMesh & mesh)
135 : {
136 : // If we've been handed an unprepared mesh then we need to be made
137 : // aware of that and fix that; we're relying on neighbor pointers.
138 18 : libmesh_assert(MeshTools::valid_is_prepared(mesh));
139 :
140 639 : if (!mesh.is_prepared())
141 71 : mesh.prepare_for_use();
142 :
143 : // We'll return a bounding box for use by subclasses in basic sanity checks.
144 621 : BoundingBox surface_bb;
145 :
146 : // First convert all volume boundaries to surface elements; this
147 : // gives us a manifold bounding the mesh, though it may not be a
148 : // connected manifold even if the volume mesh was connected.
149 : {
150 : // Make sure ids are in sync and valid on a DistributedMesh
151 639 : const dof_id_type max_orig_id = mesh.max_elem_id();
152 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
153 639 : const unique_id_type max_unique_id = mesh.parallel_max_unique_id();
154 : #endif
155 :
156 : // Change this if we add arbitrary polyhedra...
157 18 : const dof_id_type max_sides = 6;
158 :
159 36 : std::unordered_set<Elem *> elems_to_delete;
160 :
161 54 : std::vector<std::unique_ptr<Elem>> elems_to_add;
162 :
163 : // Convert all faces to surface elements
164 510132 : for (auto * elem : mesh.active_element_ptr_range())
165 : {
166 265276 : libmesh_error_msg_if (elem->dim() < 2,
167 : "Cannot use meshes with 0D or 1D elements to define a volume");
168 :
169 : // If we've already got 2D elements then those are (part of)
170 : // our surface.
171 265276 : if (elem->dim() == 2)
172 2523 : continue;
173 :
174 : // 3D elements will be removed after we've extracted their
175 : // surface faces.
176 10762 : elems_to_delete.insert(elem);
177 :
178 1313517 : for (auto s : make_range(elem->n_sides()))
179 : {
180 : // If there's a neighbor on this side then there's not a
181 : // boundary
182 1093894 : if (elem->neighbor_ptr(s))
183 : {
184 : // We're not supporting AMR meshes here yet
185 1031260 : if (elem->level() != elem->neighbor_ptr(s)->level())
186 0 : libmesh_not_implemented_msg
187 : ("Tetrahedralizaton of adapted meshes is not currently supported");
188 989004 : continue;
189 946748 : }
190 :
191 38368 : elems_to_add.push_back(elem->build_side_ptr(s));
192 796 : Elem * side_elem = elems_to_add.back().get();
193 :
194 : // Wipe the interior_parent before it can become a
195 : // dangling pointer later
196 19582 : side_elem->set_interior_parent(nullptr);
197 :
198 : // If the mesh is replicated then its automatic id
199 : // setting is fine. If not, then we need unambiguous ids
200 : // independent of element traversal.
201 19582 : if (!mesh.is_replicated())
202 : {
203 16796 : side_elem->set_id(max_orig_id + max_sides*elem->id() + s);
204 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
205 16796 : side_elem->set_unique_id(max_unique_id + max_sides*elem->id() + s);
206 : #endif
207 : }
208 : }
209 603 : }
210 :
211 : // If the mesh is replicated then its automatic neighbor finding
212 : // is fine. If not, then we need to insert them ourselves, but
213 : // it's easy because we can use the fact (from our implementation
214 : // above) that our new elements have no parents or children, plus
215 : // the fact (from the tiny fraction of homology I understand) that
216 : // a manifold boundary is a manifold with no boundary.
217 : //
218 : // See UnstructuredMesh::find_neighbors() for more explanation of
219 : // (a more complicated version of) the algorithm here.
220 639 : if (!mesh.is_replicated())
221 : {
222 : typedef dof_id_type key_type;
223 : typedef std::pair<Elem *, unsigned char> val_type;
224 : typedef std::unordered_multimap<key_type, val_type> map_type;
225 0 : map_type side_to_elem_map;
226 :
227 576 : std::unique_ptr<Elem> my_side, their_side;
228 :
229 17372 : for (auto & elem : elems_to_add)
230 : {
231 67568 : for (auto s : elem->side_index_range())
232 : {
233 50772 : if (elem->neighbor_ptr(s))
234 0 : continue;
235 50772 : const dof_id_type key = elem->low_order_key(s);
236 0 : auto bounds = side_to_elem_map.equal_range(key);
237 50772 : if (bounds.first != bounds.second)
238 : {
239 0 : elem->side_ptr(my_side, s);
240 0 : while (bounds.first != bounds.second)
241 : {
242 0 : Elem * potential_neighbor = bounds.first->second.first;
243 0 : const unsigned int ns = bounds.first->second.second;
244 0 : potential_neighbor->side_ptr(their_side, ns);
245 0 : if (*my_side == *their_side)
246 : {
247 0 : elem->set_neighbor(s, potential_neighbor);
248 0 : potential_neighbor->set_neighbor(ns, elem.get());
249 0 : side_to_elem_map.erase (bounds.first);
250 0 : break;
251 : }
252 0 : ++bounds.first;
253 : }
254 :
255 0 : if (!elem->neighbor_ptr(s))
256 : side_to_elem_map.emplace
257 0 : (key, std::make_pair(elem.get(), cast_int<unsigned char>(s)));
258 : }
259 : }
260 : }
261 :
262 : // At this point we *should* have a match for everything, so
263 : // anything we don't have a match for is remote.
264 17372 : for (auto & elem : elems_to_add)
265 67568 : for (auto s : elem->side_index_range())
266 50772 : if (!elem->neighbor_ptr(s))
267 50772 : elem->set_neighbor(s, const_cast<RemoteElem*>(remote_elem));
268 576 : }
269 :
270 : // Remove volume and edge elements
271 263314 : for (Elem * elem : elems_to_delete)
272 262675 : mesh.delete_elem(elem);
273 :
274 : // Add the new elements outside the loop so we don't risk
275 : // invalidating iterators.
276 20221 : for (auto & elem : elems_to_add)
277 21174 : mesh.add_elem(std::move(elem));
278 603 : }
279 :
280 : // Fix up neighbor pointers, element counts, etc.
281 639 : mesh.prepare_for_use();
282 :
283 : // We're making tets; we need to start with tris
284 639 : MeshTools::Modification::all_tri(mesh);
285 :
286 : // Partition surface into connected components. At this point I'm
287 : // finally going to give up and serialize, because at least we got
288 : // from 3D down to 2D first, and because I don't want to have to
289 : // turn flood_component into a while loop with a parallel sync in
290 : // the middle, and because we do have to serialize *eventually*
291 : // anyways unless we get a parallel tetrahedralizer backend someday.
292 675 : MeshSerializer mesh_serializer(mesh);
293 :
294 54 : std::vector<std::unordered_set<Elem *>> components;
295 36 : std::unordered_set<Elem *> in_component;
296 :
297 91587 : for (auto * elem : mesh.element_ptr_range())
298 1748 : if (!in_component.count(elem))
299 1794 : components.emplace_back(flood_component(in_component, elem));
300 :
301 16 : const std::unordered_set<Elem *> * biggest_component = nullptr;
302 16 : Real biggest_six_vol = 0;
303 1136 : for (const auto & component : components)
304 : {
305 1189 : Real six_vol = six_times_signed_volume(component);
306 568 : if (std::abs(six_vol) > std::abs(biggest_six_vol))
307 : {
308 16 : biggest_six_vol = six_vol;
309 16 : biggest_component = &component;
310 : }
311 : }
312 :
313 568 : if (!biggest_component)
314 0 : libmesh_error_msg("No non-zero-volume component found among " <<
315 : components.size() << " boundary components");
316 :
317 1136 : for (const auto & component : components)
318 568 : if (&component != biggest_component)
319 : {
320 0 : for (Elem * elem: component)
321 0 : mesh.delete_elem(elem);
322 : }
323 : else
324 : {
325 31524 : for (Elem * elem: component)
326 : {
327 30956 : if (biggest_six_vol < 0)
328 1168 : elem->flip(&mesh.get_boundary_info());
329 :
330 123824 : for (auto & node : elem->node_ref_range())
331 90252 : surface_bb.union_with(node);
332 : }
333 : }
334 :
335 568 : mesh.prepare_for_use();
336 :
337 584 : return surface_bb;
338 670 : }
339 :
340 :
341 432 : std::set<MeshTetInterface::SurfaceIntegrity> MeshTetInterface::check_hull_integrity() const
342 : {
343 : // Check for easy return: if the Mesh is empty (i.e. if
344 : // somebody called triangulate_conformingDelaunayMesh on
345 : // a Mesh with no elements, then hull integrity check must
346 : // fail...
347 432 : if (_mesh.n_elem() == 0)
348 0 : return {EMPTY_MESH};
349 :
350 32 : std::set<MeshTetInterface::SurfaceIntegrity> returnval;
351 :
352 432 : const BoundingBox bb = MeshTools::create_bounding_box(this->_mesh);
353 16 : const Point extents = bb.max() - bb.min();
354 432 : if (extents(0) == 0 ||
355 864 : extents(1) == 0 ||
356 16 : extents(2) == 0)
357 0 : returnval.insert(DEGENERATE_MESH);
358 :
359 : // Figure a area to use for relative tolerances when detecting
360 : // degenerate elements
361 432 : const Real ref_area = std::abs(extents(0) * extents(1)) +
362 432 : std::abs(extents(0) * extents(2)) +
363 432 : std::abs(extents(1) * extents(2));
364 :
365 416 : struct TriChecker {
366 : std::set<MeshTetInterface::SurfaceIntegrity> my_returnval;
367 : const Real my_ref_area;
368 : const unsigned int my_verbosity;
369 :
370 432 : TriChecker (Real ref_area, unsigned int verbosity) :
371 402 : my_returnval(), my_ref_area(ref_area),
372 432 : my_verbosity(verbosity) {}
373 0 : TriChecker (TriChecker & other, Threads::split) :
374 0 : my_returnval(), my_ref_area(other.my_ref_area),
375 0 : my_verbosity(other.my_verbosity) {}
376 :
377 432 : void operator()(const ConstElemRange & range) {
378 :
379 3272 : for (const Elem * elem : range)
380 : {
381 : // Check for proper element type
382 2840 : if (elem->type() != TRI3)
383 : {
384 0 : if (my_verbosity >= 50)
385 0 : std::cerr << "Non-Tri3: " << elem->get_info() << std::endl;
386 0 : my_returnval.insert(NON_TRI3);
387 : }
388 :
389 : // Make sure it's a decent element.
390 2840 : if (elem->volume() < my_ref_area * TOLERANCE * TOLERANCE)
391 : {
392 0 : if (my_verbosity >= 50)
393 0 : std::cerr << "Degenerate element: " << elem->get_info() << std::endl;
394 0 : my_returnval.insert(DEGENERATE_ELEMENT);
395 : }
396 :
397 11360 : for (auto s : elem->side_index_range())
398 : {
399 8520 : const Elem * const neigh = elem->neighbor_ptr(s);
400 :
401 8520 : if (neigh == nullptr)
402 : {
403 0 : if (my_verbosity >= 50)
404 0 : std::cerr << "Element missing neighbor " << s << ": " << elem->get_info() << std::endl;
405 0 : my_returnval.insert(MISSING_NEIGHBOR);
406 0 : continue;
407 : }
408 :
409 : // Make sure our neighbor points back to us
410 8520 : const unsigned int nn = neigh->which_neighbor_am_i(elem);
411 :
412 8520 : if (nn >= 3)
413 : {
414 0 : if (my_verbosity >= 50)
415 0 : std::cerr << "Element missing backlink " << s << ": " << elem->get_info() << std::endl;
416 0 : my_returnval.insert(MISSING_BACKLINK);
417 0 : continue;
418 : }
419 :
420 : // Our neighbor should have the same the edge nodes we do on
421 : // the neighboring edgei
422 1440 : const Node * const n1 = elem->node_ptr(s);
423 8520 : const Node * const n2 = elem->node_ptr((s+1)%3);
424 :
425 8520 : const unsigned int i1 = neigh->local_node(n1->id());
426 8520 : const unsigned int i2 = neigh->local_node(n2->id());
427 8520 : if (i1 >= 3 || i2 >= 3)
428 : {
429 0 : if (my_verbosity >= 50)
430 0 : std::cerr << "Element with bad neighbor " << s << " nodes: " << elem->get_info() << std::endl;
431 0 : my_returnval.insert(BAD_NEIGHBOR_NODES);
432 0 : continue;
433 : }
434 :
435 : // It should have those edge nodes in the opposite order
436 : // (because they have the same orientation we do)
437 8520 : if ((i2 + 1)%3 != i1)
438 : {
439 144 : if (my_verbosity >= 50)
440 0 : std::cerr << "Element orientation mismatch with neighbor " << s << ": " << elem->get_info() << std::endl;
441 144 : my_returnval.insert(NON_ORIENTED);
442 144 : continue;
443 : }
444 :
445 : // And it should have those edge nodes in the expected
446 : // places relative to its neighbor link
447 8376 : if (i2 != nn)
448 : {
449 0 : if (my_verbosity >= 50)
450 0 : std::cerr << "Element with bad links on neighbor " << s << ": " << elem->get_info() << std::endl;
451 0 : my_returnval.insert(BAD_NEIGHBOR_LINKS);
452 0 : continue;
453 : }
454 : }
455 : }
456 432 : }
457 :
458 0 : void join(TriChecker & other) {
459 0 : my_returnval.merge(other.my_returnval);
460 0 : }
461 : };
462 :
463 448 : TriChecker checker (ref_area, this->_verbosity);
464 :
465 : Threads::parallel_reduce
466 432 : (this->_mesh.active_local_element_stored_range(), checker);
467 :
468 : // Join problems found in threaded loop
469 16 : returnval.merge(checker.my_returnval);
470 :
471 : // Join problems found on other ranks
472 32 : std::set<char> int_set;
473 : std::transform
474 402 : (returnval.begin(), returnval.end(),
475 : std::inserter<std::set<char>>(int_set, int_set.end()),
476 32 : [](SurfaceIntegrity i){return int(i);});
477 432 : _mesh.comm().set_union(int_set);
478 : std::transform
479 402 : (int_set.begin(), int_set.end(),
480 : std::inserter<std::set<SurfaceIntegrity>>(returnval, returnval.end()),
481 99 : [](int i){return SurfaceIntegrity(i);});
482 :
483 16 : return returnval;
484 : }
485 :
486 :
487 :
488 430 : std::set<MeshTetInterface::SurfaceIntegrity> MeshTetInterface::improve_hull_integrity()
489 : {
490 : // We don't really do anything parallel here, but we aspire to.
491 14 : libmesh_parallel_only(this->_mesh.comm());
492 :
493 : std::set<MeshTetInterface::SurfaceIntegrity> integrityproblems =
494 444 : this->check_hull_integrity();
495 :
496 : // If we have no problem, or a problem we can't fix, we're done.
497 14 : if (integrityproblems.empty() ||
498 432 : integrityproblems.count(NON_TRI3) ||
499 418 : integrityproblems.count(EMPTY_MESH))
500 12 : return integrityproblems;
501 :
502 : // Possibly the user gave us an unprepared mesh with missing or bad
503 : // neighbor links?
504 69 : if (integrityproblems.count(MISSING_NEIGHBOR) ||
505 71 : integrityproblems.count(MISSING_BACKLINK) ||
506 71 : integrityproblems.count(BAD_NEIGHBOR_LINKS))
507 : {
508 0 : this->_mesh.find_neighbors();
509 0 : integrityproblems = this->check_hull_integrity();
510 : }
511 :
512 : // If find_neighbors() doesn't fix these, I give up.
513 69 : if (integrityproblems.count(MISSING_NEIGHBOR) ||
514 71 : integrityproblems.count(MISSING_BACKLINK) ||
515 71 : integrityproblems.count(BAD_NEIGHBOR_LINKS))
516 0 : return integrityproblems;
517 :
518 : // find_neighbors() might have fixed everything
519 71 : if (integrityproblems.empty())
520 0 : return integrityproblems;
521 :
522 : // A non-oriented (but orientable!) surface is the only thing we
523 : // shouldn't have fixed or given up on by now.
524 2 : libmesh_assert_equal_to(integrityproblems.size(), 1);
525 2 : libmesh_assert_equal_to(integrityproblems.count(NON_ORIENTED), 1);
526 :
527 : // We need one known-good triangle to start from. We'll pick the
528 : // most-negative-x normal among the triangles on the most-negative-x
529 : // point.
530 :
531 : // We'll just implement this in serial for now.
532 75 : MeshSerializer mesh_serializer(this->_mesh);
533 :
534 : // I don't see why we'd need boundary info here, but maybe we'll
535 : // want to preserve edge/node conditions eventually?
536 71 : BoundaryInfo & bi = this->_mesh.get_boundary_info();
537 :
538 140 : const Node * lowest_point = (*this->_mesh.elements_begin())->node_ptr(0);
539 :
540 : // Index by ids, not pointers, for consistency in parallel
541 4 : std::unordered_set<dof_id_type> attached_elements;
542 :
543 1260 : for (Elem * elem : this->_mesh.element_ptr_range())
544 : {
545 2272 : for (const Node & node : elem->node_ref_range())
546 : {
547 1704 : if (node(0) < (*lowest_point)(0))
548 : {
549 2 : lowest_point = &node;
550 2 : attached_elements.clear();
551 : }
552 1704 : if (&node == lowest_point)
553 390 : attached_elements.insert(elem->id());
554 : }
555 67 : }
556 :
557 2 : Elem * best_elem = nullptr;
558 2 : Real best_abs_normal_0 = 0;
559 :
560 355 : for (dof_id_type id : attached_elements)
561 : {
562 284 : Elem * elem = this->_mesh.elem_ptr(id);
563 16 : const Point e01 = elem->point(1) - elem->point(0);
564 8 : const Point e02 = elem->point(2) - elem->point(0);
565 284 : const Point normal = e01.cross(e02).unit();
566 276 : const Real abs_normal_0 = std::abs(normal(0));
567 :
568 284 : if (!best_elem || abs_normal_0 > best_abs_normal_0)
569 : {
570 2 : best_elem = elem;
571 2 : best_abs_normal_0 = abs_normal_0;
572 :
573 : // Make sure that element is actually a good one, by
574 : // flipping it if it's not.
575 71 : if (abs_normal_0 == normal(0))
576 0 : elem->flip(&bi);
577 : }
578 : }
579 :
580 : // Now flood-fill from that element to get a consistent orientation
581 : // for the others.
582 73 : std::unordered_set<dof_id_type> frontier_elements{best_elem->id()},
583 73 : finished_elements{};
584 :
585 639 : while (!frontier_elements.empty())
586 : {
587 568 : const dof_id_type elem_id = *frontier_elements.begin();
588 568 : Elem & elem = this->_mesh.elem_ref(elem_id);
589 2272 : for (auto s : elem.side_index_range())
590 : {
591 1704 : Elem * neigh = elem.neighbor_ptr(s);
592 48 : libmesh_assert(neigh);
593 48 : libmesh_assert_less(neigh->which_neighbor_am_i(&elem), 3);
594 :
595 96 : const Node * const n1 = elem.node_ptr(s);
596 1704 : const Node * const n2 = elem.node_ptr((s+1)%3);
597 1704 : const unsigned int i1 = neigh->local_node(n1->id());
598 1704 : const unsigned int i2 = neigh->local_node(n2->id());
599 48 : libmesh_assert_less(i1, 3);
600 48 : libmesh_assert_less(i2, 3);
601 :
602 1704 : const dof_id_type neigh_id = neigh->id();
603 :
604 48 : const bool frontier_neigh = frontier_elements.count(neigh_id);
605 48 : const bool finished_neigh = finished_elements.count(neigh_id);
606 :
607 : // Are we flipped?
608 1704 : if ((i2 + 1)%3 != i1)
609 : {
610 : // Are we a Moebius strip??? We give up.
611 142 : if (frontier_neigh || finished_neigh)
612 0 : return integrityproblems;
613 :
614 142 : neigh->flip(&bi);
615 : }
616 :
617 1704 : if (!frontier_neigh && !finished_neigh)
618 14 : frontier_elements.insert(neigh_id);
619 : }
620 :
621 16 : finished_elements.insert(elem_id);
622 16 : frontier_elements.erase(elem_id);
623 : }
624 :
625 71 : this->_mesh.find_neighbors();
626 :
627 2 : libmesh_assert(this->check_hull_integrity().empty());
628 :
629 71 : return {};
630 67 : }
631 :
632 :
633 430 : void MeshTetInterface::process_hull_integrity_result
634 : (const std::set<MeshTetInterface::SurfaceIntegrity> & result) const
635 : {
636 444 : std::ostringstream err_msg;
637 :
638 430 : if (result.empty()) // success
639 444 : return;
640 :
641 0 : err_msg << "Error! Conforming Delaunay mesh tetrahedralization requires a convex hull." << std::endl;
642 :
643 0 : if (result.count(NON_TRI3))
644 : {
645 0 : err_msg << "At least one non-Tri3 element was found in the input boundary mesh. ";
646 0 : err_msg << "Our constrained Delaunay tetrahedralization boundary must be a triangulation of Tri3 elements." << std::endl;
647 : }
648 0 : if (result.count(MISSING_NEIGHBOR))
649 : {
650 0 : err_msg << "At least one triangle without three neighbors was found in the input boundary mesh. ";
651 0 : err_msg << "A constrained Delaunay tetrahedralization boundary must be a triangular manifold without boundary." << std::endl;
652 : }
653 0 : if (result.count(EMPTY_MESH))
654 : {
655 0 : err_msg << "The input boundary mesh was empty!" << std::endl;
656 0 : err_msg << "Our constrained Delaunay tetrahedralization boundary must be a triangulation of Tri3 elements." << std::endl;
657 : }
658 0 : if (result.count(MISSING_BACKLINK))
659 : {
660 0 : err_msg << "At least one triangle neighbor without a return neighbor link was found in the input boundary mesh. ";
661 0 : err_msg << "A constrained Delaunay tetrahedralization boundary must be a conforming and non-adaptively-refined mesh." << std::endl;
662 : }
663 0 : if (result.count(BAD_NEIGHBOR_NODES))
664 : {
665 0 : err_msg << "At least one triangle neighbor without expected node links was found in the input boundary mesh. ";
666 0 : err_msg << "A constrained Delaunay tetrahedralization boundary must be a conforming and non-adaptively-refined mesh." << std::endl;
667 : }
668 0 : if (result.count(NON_ORIENTED))
669 : {
670 0 : err_msg << "At least one triangle neighbor with an inconsistent orientation was found in the input boundary mesh. ";
671 0 : err_msg << "A constrained Delaunay tetrahedralization boundary must be an oriented Tri3 mesh." << std::endl;
672 : }
673 0 : if (result.count(BAD_NEIGHBOR_LINKS))
674 0 : err_msg << "At least one triangle neighbor with inconsistent node and neighbor links was found in the input boundary mesh." << std::endl;
675 0 : if (result.count(DEGENERATE_ELEMENT))
676 0 : err_msg << "At least one input triangle is degenerate, with near-zero area relative to the manifold." << std::endl;
677 0 : if (result.count(DEGENERATE_MESH))
678 0 : err_msg << "The input mesh is degenerate, with zero thickness in at least one direction." << std::endl;
679 :
680 0 : libmesh_error_msg(err_msg.str());
681 402 : }
682 :
683 :
684 :
685 4 : void MeshTetInterface::delete_2D_hull_elements()
686 : {
687 10078 : for (auto & elem : this->_mesh.element_ptr_range())
688 : {
689 : // Check for proper element type. Yes, we legally delete elements while
690 : // iterating over them because no entries from the underlying container
691 : // are actually erased.
692 10072 : if (elem->type() == TRI3)
693 96 : _mesh.delete_elem(elem);
694 0 : }
695 :
696 : // We just removed any boundary info associated with hull element
697 : // edges, so let's update the boundary id caches.
698 6 : this->_mesh.get_boundary_info().regenerate_id_sets();
699 4 : }
700 :
701 :
702 :
703 : } // namespace libMesh
|
