Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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 : // Local includes
19 : #include "libmesh/cell_c0polyhedron.h"
20 :
21 : #include "libmesh/enum_order.h"
22 : #include "libmesh/face_polygon.h"
23 : #include "libmesh/mesh_tools.h"
24 : #include "libmesh/replicated_mesh.h"
25 : #include "libmesh/tensor_value.h"
26 : #include "libmesh/node.h"
27 :
28 : // C++ headers
29 : #include <numeric> // std::iota
30 :
31 : namespace libMesh
32 : {
33 :
34 : // ------------------------------------------------------------
35 : // C0Polyhedron class member functions
36 :
37 6121 : C0Polyhedron::C0Polyhedron
38 6121 : (const std::vector<std::shared_ptr<Polygon>> & sides, std::unique_ptr<Node> & mid_elem_node, Elem * p) :
39 : Polyhedron(sides, p),
40 6121 : _has_mid_elem_node(false)
41 : {
42 6121 : this->retriangulate();
43 :
44 : // Grab the last node
45 6121 : if (_has_mid_elem_node)
46 2059 : mid_elem_node.reset(this->_nodelinks_data.back());
47 6121 : }
48 :
49 :
50 :
51 11546 : C0Polyhedron::~C0Polyhedron() = default;
52 :
53 :
54 :
55 15 : std::unique_ptr<Elem> C0Polyhedron::disconnected_clone() const
56 : {
57 19 : auto sides = this->side_clones();
58 :
59 15 : std::unique_ptr<Node> mid_elem_node;
60 15 : std::unique_ptr<Elem> returnval = std::make_unique<C0Polyhedron>(sides, mid_elem_node);
61 : // Swap out the new mid elem node with the previous one
62 15 : if (mid_elem_node)
63 : {
64 0 : libmesh_assert(returnval->n_nodes() == this->n_nodes());
65 0 : returnval->set_node(returnval->n_nodes() - 1, this->_nodelinks_data.back());
66 : }
67 :
68 15 : returnval->set_id() = this->id();
69 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
70 15 : if (this->valid_unique_id())
71 2 : returnval->set_unique_id(this->unique_id());
72 : #endif
73 :
74 15 : const auto n_elem_ints = this->n_extra_integers();
75 15 : returnval->add_extra_integers(n_elem_ints);
76 15 : for (unsigned int i = 0; i != n_elem_ints; ++i)
77 0 : returnval->set_extra_integer(i, this->get_extra_integer(i));
78 :
79 13 : returnval->inherit_data_from(*this);
80 :
81 17 : return returnval;
82 11 : }
83 :
84 :
85 :
86 : unsigned int
87 12737170 : C0Polyhedron::n_vertices() const
88 : {
89 13798968 : return this->_nodelinks_data.size() - _has_mid_elem_node;
90 : }
91 :
92 :
93 :
94 20286 : bool C0Polyhedron::is_vertex(const unsigned int i) const
95 : {
96 2252 : libmesh_assert_less (i, this->n_nodes());
97 :
98 20286 : if (i < this->n_vertices())
99 2144 : return true;
100 : else
101 810 : return false;
102 : }
103 :
104 :
105 :
106 4260 : bool C0Polyhedron::is_edge(const unsigned int libmesh_dbg_var(i)) const
107 : {
108 120 : libmesh_assert_less (i, this->n_nodes());
109 :
110 4260 : return false;
111 : }
112 :
113 :
114 :
115 4260 : bool C0Polyhedron::is_face(const unsigned int libmesh_dbg_var(i)) const
116 : {
117 120 : libmesh_assert_less (i, this->n_nodes());
118 :
119 4260 : return false;
120 : }
121 :
122 :
123 :
124 5124 : bool C0Polyhedron::is_node_on_side(const unsigned int n,
125 : const unsigned int s) const
126 : {
127 192 : libmesh_assert_less (n, this->n_nodes());
128 192 : libmesh_assert_less (s, this->n_sides());
129 :
130 : const std::vector<unsigned int> & node_map =
131 5124 : std::get<2>(_sidelinks_data[s]);
132 :
133 5124 : const auto it = std::find(node_map.begin(), node_map.end(), n);
134 :
135 5316 : return (it != node_map.end());
136 : }
137 :
138 : std::vector<unsigned int>
139 16270 : C0Polyhedron::nodes_on_side(const unsigned int s) const
140 : {
141 502 : libmesh_assert_less(s, this->n_sides());
142 16772 : return std::get<2>(_sidelinks_data[s]);
143 : }
144 :
145 : std::vector<unsigned int>
146 684 : C0Polyhedron::nodes_on_edge(const unsigned int e) const
147 : {
148 684 : auto [s, se] = _edge_lookup[e];
149 684 : const Polygon & face = *std::get<0>(_sidelinks_data[s]);
150 : const std::vector<unsigned int> & node_map =
151 57 : std::get<2>(_sidelinks_data[s]);
152 : std::vector<unsigned int> nodes_on_edge =
153 684 : face.nodes_on_side(se);
154 2052 : for (auto i : index_range(nodes_on_edge))
155 1482 : nodes_on_edge[i] = node_map[nodes_on_edge[i]];
156 741 : return nodes_on_edge;
157 : }
158 :
159 :
160 :
161 288 : bool C0Polyhedron::is_node_on_edge(const unsigned int n,
162 : const unsigned int e) const
163 : {
164 24 : libmesh_assert_less(e, _edge_lookup.size());
165 288 : auto [s, se] = _edge_lookup[e];
166 :
167 288 : const Polygon & face = *std::get<0>(_sidelinks_data[s]);
168 : const std::vector<unsigned int> & node_map =
169 24 : std::get<2>(_sidelinks_data[s]);
170 : std::vector<unsigned int> nodes_on_edge =
171 312 : face.nodes_on_side(se);
172 432 : for (auto noe : nodes_on_edge)
173 468 : if (node_map[noe] == n)
174 24 : return true;
175 :
176 0 : return false;
177 : }
178 :
179 :
180 :
181 : std::vector<unsigned int>
182 480 : C0Polyhedron::edges_adjacent_to_node(const unsigned int n) const
183 : {
184 480 : if (n == n_nodes() - 1)
185 55 : return {};
186 420 : return Polyhedron::edges_adjacent_to_node(n);
187 : }
188 :
189 :
190 :
191 0 : void C0Polyhedron::connectivity(const unsigned int /*sf*/,
192 : const IOPackage /*iop*/,
193 : std::vector<dof_id_type> & /*conn*/) const
194 : {
195 0 : libmesh_not_implemented();
196 : }
197 :
198 :
199 :
200 296 : Real C0Polyhedron::volume () const
201 : {
202 : // This specialization is ... probably good for general non-Lagrange
203 : // mappings, since we require our polyhedral faces to be planar, but
204 : // I'd rather be slow than wrong.
205 296 : if (this->mapping_type() != LAGRANGE_MAP)
206 0 : return this->Elem::volume();
207 :
208 : // We use a triangulation to calculate here.
209 :
210 9 : Real six_vol = 0;
211 3906 : for (const auto & subtet : this->_triangulation)
212 : {
213 3610 : const Point p0 = this->point(subtet[0]);
214 3610 : const Point p1 = this->point(subtet[1]);
215 3610 : const Point p2 = this->point(subtet[2]);
216 3715 : const Point p3 = this->point(subtet[3]);
217 :
218 105 : const Point v01 = p1 - p0;
219 105 : const Point v02 = p2 - p0;
220 105 : const Point v03 = p3 - p0;
221 :
222 3610 : six_vol += triple_product(v01, v02, v03);
223 : }
224 :
225 296 : return six_vol * (1./6.);
226 : }
227 :
228 :
229 :
230 24 : Point C0Polyhedron::true_centroid () const
231 : {
232 : // This specialization is good for Lagrange mappings only
233 24 : if (this->mapping_type() != LAGRANGE_MAP)
234 0 : return this->Elem::true_centroid();
235 :
236 2 : Real six_vol = 0;
237 2 : Point six_vol_weighted_centroid;
238 144 : for (const auto & subtet : this->_triangulation)
239 : {
240 120 : const Point p0 = this->point(subtet[0]);
241 120 : const Point p1 = this->point(subtet[1]);
242 120 : const Point p2 = this->point(subtet[2]);
243 130 : const Point p3 = this->point(subtet[3]);
244 :
245 10 : const Point v01 = p1 - p0;
246 10 : const Point v02 = p2 - p0;
247 10 : const Point v03 = p3 - p0;
248 :
249 110 : const Real six_tet_vol = triple_product(v01, v02, v03);
250 :
251 10 : const Point tet_centroid = (p0 + p1 + p2 + p3) * 0.25;
252 :
253 120 : six_vol += six_tet_vol;
254 :
255 10 : six_vol_weighted_centroid += six_tet_vol * tet_centroid;
256 : }
257 :
258 2 : return six_vol_weighted_centroid / six_vol;
259 : }
260 :
261 :
262 :
263 : std::pair<unsigned short int, unsigned short int>
264 0 : C0Polyhedron::second_order_child_vertex (const unsigned int /*n*/) const
265 : {
266 0 : libmesh_not_implemented();
267 : return std::pair<unsigned short int, unsigned short int> (0,0);
268 : }
269 :
270 :
271 :
272 216 : ElemType C0Polyhedron::side_type (const unsigned int libmesh_dbg_var(s)) const
273 : {
274 18 : libmesh_assert_less (s, this->n_sides());
275 216 : return C0POLYGON;
276 : }
277 :
278 :
279 :
280 : Point
281 852 : C0Polyhedron::side_vertex_average_normal(const unsigned int s) const
282 : {
283 24 : libmesh_assert_less (s, this->n_sides());
284 24 : libmesh_assert_equal_to(this->mapping_type(), LAGRANGE_MAP);
285 852 : const auto poly_side_ptr = this->side_ptr(s);
286 852 : const auto n_side_edges = poly_side_ptr->n_sides();
287 :
288 : // At the side vertex average, things simplify a bit
289 : // We get the side "plane" normal at all vertices, then average them
290 24 : Point normal;
291 48 : Point current_edge = poly_side_ptr->point(1) - poly_side_ptr->point(0);
292 852 : for (auto i : make_range(n_side_edges))
293 : {
294 852 : const Point next_edge = poly_side_ptr->point((i + 2) % n_side_edges) -
295 876 : poly_side_ptr->point((i + 1) % n_side_edges);
296 828 : const Point normal_at_vertex = current_edge.cross(next_edge);
297 24 : normal += normal_at_vertex;
298 : // Note: the sides are planar, we don't need to test them all
299 852 : if (normal.norm_sq() > TOLERANCE)
300 24 : break;
301 0 : current_edge = next_edge;
302 : }
303 852 : bool outward_normal = std::get<1>(_sidelinks_data[s]);
304 1708 : return (outward_normal ? 1. : -1.) * normal.unit();
305 804 : }
306 :
307 :
308 : std::array<int, 4>
309 1917 : C0Polyhedron::subelement_sides_to_poly_sides(unsigned int tri_i) const
310 : {
311 54 : std::array<int, 4> sides = {invalid_int, invalid_int, invalid_int, invalid_int};
312 54 : libmesh_assert(tri_i < this->n_subelements());
313 1917 : const auto & tet_nodes = _triangulation[tri_i];
314 16401 : for (const auto poly_side : make_range(n_sides()))
315 : {
316 14892 : const auto sd_nodes = nodes_on_side(poly_side);
317 14892 : bool zero_in = std::find(sd_nodes.begin(), sd_nodes.end(), tet_nodes[0]) != sd_nodes.end();
318 14892 : bool one_in = std::find(sd_nodes.begin(), sd_nodes.end(), tet_nodes[1]) != sd_nodes.end();
319 14892 : bool two_in = std::find(sd_nodes.begin(), sd_nodes.end(), tet_nodes[2]) != sd_nodes.end();
320 14892 : bool three_in = std::find(sd_nodes.begin(), sd_nodes.end(), tet_nodes[3]) != sd_nodes.end();
321 :
322 : // Take advantage of the fact that the poly side can only contain one tet side
323 : // Note: we could optimize by replacing the booleans with the searches above
324 14484 : if (zero_in)
325 : {
326 5751 : if (one_in)
327 : {
328 2485 : if (two_in)
329 1491 : sides[0] = poly_side;
330 994 : else if (three_in)
331 355 : sides[1] = poly_side;
332 : }
333 3266 : else if (two_in && three_in)
334 355 : sides[3] = poly_side;
335 : }
336 8733 : else if (three_in && two_in && one_in)
337 355 : sides[2] = poly_side;
338 : }
339 1917 : return sides;
340 : }
341 :
342 :
343 6121 : void C0Polyhedron::retriangulate()
344 : {
345 174 : this->_triangulation.clear();
346 :
347 : // We should already have a triangulation for each side. We'll turn
348 : // that into a triangulation of the entire polyhedral surface
349 : // (stored as a mesh, so we can use `find_neighbors()`, then turn
350 : // that into a tetrahedralization by shrinking the volume one tet at
351 : // a time, which should work fine for convex polyhedra.
352 :
353 348 : Parallel::Communicator comm_self; // the default communicator is 1 rank
354 6469 : ReplicatedMesh surface(comm_self);
355 :
356 : // poly_node_to_local_id[poly_node] is the local id of
357 : // poly_node in the Polyhedron, which is also the global id of the
358 : // corresponding node in the surface mesh.
359 348 : std::unordered_map<Node *, unsigned int> poly_node_to_id;
360 :
361 43131 : for (unsigned int s : make_range(this->n_sides()))
362 : {
363 37010 : const auto & [side, inward_normal, node_map] = this->_sidelinks_data[s];
364 :
365 111598 : for (auto t : make_range(side->n_subtriangles()))
366 : {
367 76708 : Elem * tri = surface.add_elem(Elem::build(TRI3));
368 :
369 74588 : const std::array<int, 3> subtri = side->subtriangle(t);
370 :
371 298352 : for (int i : make_range(3))
372 : {
373 223764 : const int side_id = subtri[i];
374 223764 : const Node * poly_node = side->node_ptr(side_id);
375 :
376 6360 : libmesh_assert_less(side_id, node_map.size());
377 223764 : const unsigned int local_id = node_map[side_id];
378 :
379 223764 : Node * surf_node = surface.query_node_ptr(local_id);
380 223764 : if (surf_node)
381 4952 : libmesh_assert_equal_to(*(const Point*)poly_node,
382 : *(const Point*)surf_node);
383 : else
384 49536 : surf_node = surface.add_point(*poly_node, local_id);
385 :
386 : // Get a consistent orientation for surface triangles,
387 : // facing their zeta directions outward
388 223764 : const int tri_node = inward_normal ? i : 2-i;
389 223764 : tri->set_node(tri_node, surf_node);
390 : }
391 : }
392 : }
393 :
394 174 : surface.allow_renumbering(false);
395 6121 : surface.prepare_for_use();
396 :
397 : // We should have a watertight surface with consistent triangle
398 : // orientations. That's expensive to check.
399 : #ifdef DEBUG
400 22582 : auto verify_surface = [& surface] ()
401 : {
402 8118 : for (const Elem * elem : surface.element_ptr_range())
403 : {
404 28928 : for (auto s : make_range(3))
405 : {
406 21696 : const Elem * neigh = elem->neighbor_ptr(s);
407 21696 : libmesh_assert(neigh);
408 21696 : libmesh_assert_equal_to(neigh,
409 : surface.elem_ptr(neigh->id()));
410 21696 : const unsigned int ns = neigh->which_neighbor_am_i(elem);
411 21696 : libmesh_assert_less(ns, 3);
412 21696 : libmesh_assert_equal_to(elem->node_ptr(s),
413 : neigh->node_ptr((ns+1)%3));
414 21696 : libmesh_assert_equal_to(elem->node_ptr((s+1)%3),
415 : neigh->node_ptr(ns));
416 : }
417 : }
418 886 : };
419 :
420 174 : verify_surface();
421 : #endif
422 :
423 : // First heuristic: try with no interior point
424 : // This might not succeed, not every surface triangulation gives a tetrahedralization
425 : // with no additional interior point
426 : // But if it succeeds, it uses less tetrahedra to cut the polyhedron
427 : #ifdef LIBMESH_ENABLE_EXCEPTIONS
428 : try
429 : {
430 : // We'll have to edit this as we change the surface elements, but we
431 : // have a method to initialize it easily.
432 522 : std::vector<std::vector<dof_id_type>> nodes_to_elem_vec_map;
433 6121 : MeshTools::build_nodes_to_elem_map(surface, nodes_to_elem_vec_map);
434 :
435 : // We'll be inserting and deleting entries heavily, so we'll use
436 : // sets rather than vectors. We want to get the same results in
437 : // parallel, so we'll use element ids rather than Elem pointers
438 464 : std::vector<std::set<dof_id_type>> nodes_to_elem_map;
439 55657 : for (auto i : index_range(nodes_to_elem_vec_map))
440 : nodes_to_elem_map.emplace_back
441 99665 : (nodes_to_elem_vec_map[i].begin(),
442 53818 : nodes_to_elem_vec_map[i].end());
443 :
444 : // Now start meshing the volume enclosed by surface, one tet at a
445 : // time, with a greedy heuristic: find the vertex node with the most
446 : // acutely convex (solid) angle and strip it out as tetrahedra.
447 :
448 : // We'll want a vector of surrounding nodes for multiple uses,
449 : // sometimes with a similarly-sorted vector of surrounding elements
450 : // to go with it.
451 : auto surroundings_of =
452 80275 : [&nodes_to_elem_map, & surface]
453 : (const Node & node,
454 101879 : std::vector<Elem *> * surrounding_elems)
455 : {
456 : const std::set<dof_id_type> & elems_by_node =
457 6045 : nodes_to_elem_map[node.id()];
458 :
459 85111 : const unsigned int n_surrounding = elems_by_node.size();
460 2418 : libmesh_assert_greater_equal(n_surrounding, 3);
461 :
462 85111 : if (surrounding_elems)
463 : {
464 886 : libmesh_assert(surrounding_elems->empty());
465 31173 : surrounding_elems->resize(n_surrounding);
466 : }
467 :
468 85111 : std::vector<Node *> surrounding_nodes(n_surrounding);
469 :
470 85111 : Elem * elem = surface.elem_ptr(*elems_by_node.begin());
471 420002 : for (auto i : make_range(n_surrounding))
472 : {
473 325377 : const unsigned int n = elem->get_node_index(&node);
474 9514 : libmesh_assert_not_equal_to(n, invalid_uint);
475 334891 : Node * next_node = elem->node_ptr((n+1)%3);
476 334891 : surrounding_nodes[i] = next_node;
477 334891 : if (surrounding_elems)
478 100703 : (*surrounding_elems)[i] = elem;
479 334891 : elem = elem->neighbor_ptr((n+2)%3);
480 9514 : libmesh_assert(elem);
481 9514 : libmesh_assert_equal_to(elem, surface.elem_ptr(elem->id()));
482 :
483 : // We should have a manifold here, but verifying that is
484 : // expensive
485 : #ifdef DEBUG
486 9514 : libmesh_assert_equal_to
487 : (std::count(surrounding_nodes.begin(),
488 : surrounding_nodes.end(), next_node),
489 : 1);
490 : #endif
491 : }
492 :
493 : // We should have finished a loop
494 2418 : libmesh_assert_equal_to
495 : (elem, surface.elem_ptr(*elems_by_node.begin()));
496 :
497 87529 : return surrounding_nodes;
498 6121 : };
499 :
500 53938 : auto geometry_at = [&surroundings_of](const Node & node)
501 : {
502 : const std::vector<Node *> surrounding_nodes =
503 55470 : surroundings_of(node, nullptr);
504 :
505 : // Now sum up solid angles from tetrahedra created from the
506 : // trivial triangulation of the surrounding nodes loop
507 1532 : Real total_solid_angle = 0;
508 : const int n_surrounding =
509 3064 : cast_int<int>(surrounding_nodes.size());
510 :
511 183032 : for (auto n : make_range(n_surrounding-2))
512 : {
513 : const Point
514 132762 : v01 = static_cast<Point>(*surrounding_nodes[n]) - node,
515 132762 : v02 = static_cast<Point>(*surrounding_nodes[n+1]) - node,
516 132762 : v03 = static_cast<Point>(*surrounding_nodes[n+2]) - node;
517 :
518 129094 : total_solid_angle += solid_angle(v01, v02, v03);
519 : }
520 :
521 55470 : return std::make_pair(n_surrounding, total_solid_angle);
522 5947 : };
523 :
524 : // We'll keep track of solid angles and node valences so we don't
525 : // waste time recomputing them when they haven't changed. We need
526 : // to be able to search efficiently for the smallest angles of each
527 : // valence, but also search efficiently for a particular node to
528 : // remove and reinsert it when its connectivity changes.
529 : //
530 : // Since C++11 multimap has guaranteed that pairs with matching keys
531 : // are kept in insertion order, so we can use Node * for values even
532 : // in parallel.
533 : typedef std::multimap<std::pair<int, Real>, Node*> node_map_type;
534 348 : node_map_type nodes_by_geometry;
535 232 : std::map<Node *, node_map_type::iterator> node_index;
536 :
537 61604 : for (auto node : surface.node_ptr_range())
538 49536 : node_index[node] =
539 104845 : nodes_by_geometry.emplace(geometry_at(*node), node);
540 :
541 : // In 3D, this will require nested loops: an outer loop to remove
542 : // each vertex, and an inner loop to remove multiple tetrahedra in
543 : // cases where the vertex has more than 3 neighboring triangles.
544 :
545 : // We'll be done when there are only three "unremoved" nodes left,
546 : // so they don't actually enclose any volume.
547 35235 : for (auto i : make_range(nodes_by_geometry.size()-3))
548 : {
549 886 : auto geometry_it = nodes_by_geometry.begin();
550 31173 : auto geometry_key = geometry_it->first;
551 31173 : auto [valence, angle] = geometry_key;
552 31173 : Node * node = geometry_it->second;
553 886 : libmesh_ignore(i);
554 :
555 : // If our lowest-valence nodes are all points of non-convexity,
556 : // skip to a higher valence.
557 31173 : while (angle > 2*pi-TOLERANCE)
558 : {
559 : geometry_it =
560 : nodes_by_geometry.upper_bound
561 0 : (std::make_pair(valence, Real(100)));
562 0 : libmesh_assert(geometry_it != nodes_by_geometry.end());
563 :
564 0 : std::tie(geometry_key, node) = *geometry_it;
565 0 : std::tie(valence, angle) = geometry_key;
566 : }
567 :
568 1772 : std::vector<Elem *> surrounding_elems;
569 : std::vector<Node *> surrounding_nodes =
570 32059 : surroundings_of(*node, &surrounding_elems);
571 :
572 1772 : const std::size_t n_surrounding = surrounding_nodes.size();
573 :
574 : // As we separate surrounding nodes from our center node, we'll
575 : // be marking them as nullptr; we still need to be able to find
576 : // predecessor and successor nodes in order afterward.
577 : auto find_valid_nodes_around =
578 162762 : [n_surrounding, & surrounding_nodes]
579 369794 : (unsigned int j)
580 : {
581 172546 : unsigned int jnext = (j+1)%n_surrounding;
582 186636 : while (!surrounding_nodes[jnext])
583 8946 : jnext = (jnext+1)%n_surrounding;
584 :
585 172546 : unsigned int jprev = (j+n_surrounding-1)%n_surrounding;
586 190578 : while (!surrounding_nodes[jprev])
587 12780 : jprev = (jprev+n_surrounding-1)%n_surrounding;
588 :
589 177438 : return std::make_pair(jprev, jnext);
590 31173 : };
591 :
592 : // We may have too many surrounding nodes to handle with
593 : // just one tet. In that case we'll keep a cache of the
594 : // element qualities that we'd get by making a tet with the
595 : // edge from the center node to each surrounding node, so we
596 : // can build the best tets first.
597 : //
598 : // In the case where we just have 3 nodes, we'll just pretend
599 : // they all have the same positive quality, so we can still
600 : // search this vector.
601 34060 : std::vector<Real> local_tet_quality(n_surrounding, 1);
602 :
603 : // From our center node with N surrounding nodes we can make N-2
604 : // tetrahedra. The first N-3 each replace two surface tets with
605 : // two new surface tets.
606 : //
607 : // My first idea was to greedily pick nodes with the smallest
608 : // local (solid) angles to get the best quality. This works in
609 : // 2D, but such a node can give a pancake tet in 3D.
610 : //
611 : // My second idea was to greedily pick nodes with the highest
612 : // prospective tet quality. This works for the first tets, but
613 : // can leave a pancake tet behind.
614 : //
615 : // My third idea is to try to fix the lowest quality tets first,
616 : // by picking cases where they have higher quality neighbors,
617 : // and creating those neighbors so as to change them.
618 :
619 : auto find_new_tet_nodes =
620 33555 : [& local_tet_quality, & find_valid_nodes_around]
621 63849 : ()
622 : {
623 1010 : unsigned int jbest = 0;
624 35575 : auto [jminus, jplus] = find_valid_nodes_around(jbest);
625 35575 : Real qneighbest = std::min(local_tet_quality[jminus],
626 36719 : local_tet_quality[jplus]);
627 35823 : for (auto j : make_range(std::size_t(1),
628 117797 : local_tet_quality.size()))
629 : {
630 : // We don't want to build a bad tet
631 82222 : if (local_tet_quality[j] <= 0)
632 4830 : continue;
633 :
634 74984 : std::tie(jminus, jplus) = find_valid_nodes_around(j);
635 74984 : Real qneighj = std::min(local_tet_quality[jminus],
636 77112 : local_tet_quality[jplus]);
637 :
638 : // We don't want to build a tet that can't fix a neighbor
639 : // if we can build one that can.
640 74984 : if (qneighbest <= 0 &&
641 : qneighj > 0)
642 4278 : continue;
643 :
644 : // We want to try for the best possible fix.
645 72586 : if ((local_tet_quality[j] - qneighj) >
646 72586 : (local_tet_quality[jbest] - qneighj))
647 : {
648 240 : jbest = j;
649 240 : qneighbest = qneighj;
650 : }
651 : }
652 :
653 36585 : libmesh_error_msg_if
654 : (local_tet_quality[jbest] <= 0,
655 : "Cannot build non-singular non-inverted tet");
656 :
657 35575 : std::tie(jminus, jplus) = find_valid_nodes_around(jbest);
658 :
659 36585 : return std::make_tuple(jbest, jminus, jplus);
660 31173 : };
661 :
662 31173 : if (n_surrounding > 3)
663 : {
664 : // We'll be searching local_tet_quality even after tet
665 : // extraction disconnects us from some nodes; when we do we
666 : // don't want to get one.
667 124 : constexpr Real far_node = -1e6;
668 :
669 : // Vectors from the center node to each of its surrounding
670 : // nodes are helpful for calculating prospective tet
671 : // quality.
672 4526 : std::vector<Point> v0s(n_surrounding);
673 22010 : for (auto j : make_range(n_surrounding))
674 19096 : v0s[j] = *(Point *)surrounding_nodes[j] - *node;
675 :
676 : // Find the tet quality we'd potentially get from each
677 : // possible choice of tet
678 : auto local_tet_quality_of =
679 24924 : [& surrounding_nodes, & v0s, & find_valid_nodes_around]
680 65472 : (unsigned int j)
681 : {
682 26412 : auto [jminus, jplus] = find_valid_nodes_around(j);
683 :
684 : // Anything proportional to the ratio of volume to
685 : // total-edge-length-cubed should peak for perfect tets
686 : // but hit 0 for pancakes and slivers.
687 :
688 : const Real total_len =
689 30132 : v0s[j].norm() + v0s[jminus].norm() + v0s[jplus].norm() +
690 24924 : (*(Point *)surrounding_nodes[jplus] -
691 53568 : *(Point *)surrounding_nodes[j]).norm() +
692 24924 : (*(Point *)surrounding_nodes[j] -
693 27900 : *(Point *)surrounding_nodes[jminus]).norm() +
694 24924 : (*(Point *)surrounding_nodes[jminus] -
695 28644 : *(Point *)surrounding_nodes[jplus]).norm();
696 :
697 : // Orientation here is tricky. Think of the triple
698 : // product as (v1 cross v2) dot v3, with right hand rule.
699 : const Real six_vol =
700 26412 : triple_product(v0s[jminus], v0s[jplus], v0s[j]);
701 :
702 26412 : return six_vol / (total_len * total_len * total_len);
703 4402 : };
704 :
705 22010 : for (auto j : make_range(n_surrounding))
706 17608 : local_tet_quality[j] = local_tet_quality_of(j);
707 :
708 : // If we have N surrounding nodes, we can make N tets and
709 : // that'll bring us back to the 3-surrounding-node case to
710 : // finish.
711 8804 : for (auto t : make_range(n_surrounding-3))
712 : {
713 124 : libmesh_ignore(t);
714 :
715 4402 : auto [jbest, jminus, jplus] = find_new_tet_nodes();
716 :
717 : // Turn these four nodes into a tet
718 4402 : Node * nbest = surrounding_nodes[jbest],
719 4402 : * nminus = surrounding_nodes[jminus],
720 4402 : * nplus = surrounding_nodes[jplus];
721 4402 : this->add_tet(nminus->id(), nbest->id(), nplus->id(),
722 248 : node->id());
723 :
724 : // Replace the old two triangles around these nodes with the
725 : // proper two new triangles.
726 4402 : Elem * oldtri1 = surrounding_elems[jminus],
727 4402 : * oldtri2 = surrounding_elems[jbest],
728 4402 : * newtri1 = surface.add_elem(Elem::build(TRI3)),
729 4526 : * newtri2 = surface.add_elem(Elem::build(TRI3));
730 :
731 4278 : const unsigned int c1 = oldtri1->get_node_index(node),
732 4278 : c2 = oldtri2->get_node_index(node);
733 :
734 4402 : newtri1->set_node(0, node);
735 4402 : newtri1->set_node(1, nminus);
736 4402 : newtri1->set_node(2, nplus);
737 :
738 4402 : surrounding_elems[jminus] = newtri1;
739 :
740 248 : Elem * neigh10 = oldtri1->neighbor_ptr(c1);
741 4402 : Elem * neigh12 = oldtri2->neighbor_ptr((c2+2)%3);
742 248 : newtri1->set_neighbor(0, neigh10);
743 4402 : neigh10->set_neighbor(neigh10->which_neighbor_am_i(oldtri1), newtri1);
744 248 : newtri1->set_neighbor(1, newtri2);
745 124 : newtri1->set_neighbor(2, neigh12);
746 4402 : neigh12->set_neighbor(neigh12->which_neighbor_am_i(oldtri2), newtri1);
747 :
748 4402 : newtri2->set_node(0, nplus);
749 4402 : newtri2->set_node(1, nminus);
750 4402 : newtri2->set_node(2, nbest);
751 :
752 4402 : Elem * neigh21 = oldtri1->neighbor_ptr((c1+1)%3);
753 4402 : Elem * neigh22 = oldtri2->neighbor_ptr((c2+1)%3);
754 248 : newtri2->set_neighbor(0, newtri1);
755 124 : newtri2->set_neighbor(1, neigh21);
756 4402 : neigh21->set_neighbor(neigh21->which_neighbor_am_i(oldtri1), newtri2);
757 248 : newtri2->set_neighbor(2, neigh22);
758 4402 : neigh22->set_neighbor(neigh22->which_neighbor_am_i(oldtri2), newtri2);
759 :
760 17608 : for (unsigned int p : make_range(3))
761 : {
762 13950 : nodes_to_elem_map[oldtri1->node_id(p)].erase(oldtri1->id());
763 13950 : nodes_to_elem_map[oldtri2->node_id(p)].erase(oldtri2->id());
764 13950 : nodes_to_elem_map[newtri1->node_id(p)].insert(newtri1->id());
765 13950 : nodes_to_elem_map[newtri2->node_id(p)].insert(newtri2->id());
766 : }
767 :
768 : // We've finished replacing the old triangles.
769 4402 : surface.delete_elem(oldtri1);
770 4402 : surface.delete_elem(oldtri2);
771 :
772 : // The solid angle for the far node should now stay
773 : // unchanged until we're out of this inner loop; let's
774 : // recalculate it here, and then we'll be done with it.
775 248 : Node * & nbestref = surrounding_nodes[jbest];
776 4402 : nodes_by_geometry.erase(node_index[nbestref]);
777 4402 : node_index[nbestref] =
778 8556 : nodes_by_geometry.emplace(geometry_at(*nbestref), nbestref);
779 :
780 : // The far node is no longer sharing an edge with our center
781 : // node. Make sure we don't use it again with the center
782 : // node.
783 4402 : local_tet_quality[jbest] = far_node;
784 4402 : nbestref = nullptr;
785 :
786 : // The potential tet qualities using the side nodes have
787 : // changed now that they're directly connected to each
788 : // other.
789 4526 : local_tet_quality[jminus] =
790 4402 : local_tet_quality_of(jminus);
791 :
792 4526 : local_tet_quality[jplus] =
793 4402 : local_tet_quality_of(jplus);
794 : }
795 : }
796 :
797 : // Now we should have just 3 surrounding nodes, with which to
798 : // make one tetrahedron. Put them in a counterclockwise
799 : // (looking from outside) orientation, not the "best, clockwise,
800 : // counterclockwise" we get from the lambda.
801 31173 : auto [j2, j1, j3] = find_new_tet_nodes();
802 :
803 : // Turn these last four nodes into a tet
804 31173 : Node * n1 = surrounding_nodes[j1],
805 31173 : * n2 = surrounding_nodes[j2],
806 31173 : * n3 = surrounding_nodes[j3];
807 31173 : this->add_tet(n1->id(), n2->id(), n3->id(), node->id());
808 :
809 : // Replace the three surface triangles of that tet with the new
810 : // fourth triangle.
811 29114 : Elem * oldtri1 = surrounding_elems[j1],
812 29114 : * oldtri2 = surrounding_elems[j2],
813 29114 : * oldtri3 = surrounding_elems[j3],
814 29942 : * newtri = surface.add_elem(Elem::build(TRI3));
815 :
816 28286 : const unsigned int c1 = oldtri1->get_node_index(node),
817 28286 : c2 = oldtri2->get_node_index(node),
818 28286 : c3 = oldtri3->get_node_index(node);
819 :
820 29114 : newtri->set_node(0, n1);
821 29114 : newtri->set_node(1, n2);
822 29114 : newtri->set_node(2, n3);
823 29114 : Elem * neigh0 = oldtri1->neighbor_ptr((c1+1)%3);
824 1656 : newtri->set_neighbor(0, neigh0);
825 29114 : neigh0->set_neighbor(neigh0->which_neighbor_am_i(oldtri1), newtri);
826 29114 : Elem * neigh1 = oldtri2->neighbor_ptr((c2+1)%3);
827 1656 : newtri->set_neighbor(1, neigh1);
828 29114 : neigh1->set_neighbor(neigh1->which_neighbor_am_i(oldtri2), newtri);
829 29114 : Elem * neigh2 = oldtri3->neighbor_ptr((c3+1)%3);
830 1656 : newtri->set_neighbor(2, neigh2);
831 29114 : neigh2->set_neighbor(neigh2->which_neighbor_am_i(oldtri3), newtri);
832 :
833 116456 : for (unsigned int p : make_range(3))
834 : {
835 92310 : nodes_to_elem_map[oldtri1->node_id(p)].erase(oldtri1->id());
836 92310 : nodes_to_elem_map[oldtri2->node_id(p)].erase(oldtri2->id());
837 92310 : nodes_to_elem_map[oldtri3->node_id(p)].erase(oldtri3->id());
838 94311 : nodes_to_elem_map[newtri->node_id(p)].insert(newtri->id());
839 : }
840 :
841 : // We've finished replacing the old triangles.
842 29114 : surface.delete_elem(oldtri1);
843 29114 : surface.delete_elem(oldtri2);
844 29114 : surface.delete_elem(oldtri3);
845 :
846 : // We should have used up all our surrounding nodes now, and we
847 : // shouldn't have messed up our surface in the process, and our
848 : // center node should no longer be part of the surface.
849 : #ifdef DEBUG
850 828 : surrounding_nodes[j1] = nullptr;
851 828 : surrounding_nodes[j2] = nullptr;
852 828 : surrounding_nodes[j3] = nullptr;
853 :
854 3436 : for (auto ltq : surrounding_nodes)
855 2608 : libmesh_assert(!ltq);
856 :
857 828 : if (surface.n_elem() > 3)
858 712 : verify_surface();
859 :
860 6172 : for (const Elem * elem : surface.element_ptr_range())
861 21376 : for (auto p : make_range(3))
862 16032 : libmesh_assert_not_equal_to
863 : (elem->node_ptr(p), node);
864 : #endif
865 :
866 : // We've used up our center node, so it's not something we can
867 : // eliminate again.
868 28286 : nodes_by_geometry.erase(geometry_it);
869 : }
870 : // At this point our surface should just have two triangles left.
871 116 : libmesh_assert_equal_to(surface.n_elem(), 2);
872 4062 : _has_mid_elem_node = false;
873 7716 : }
874 : // Failed without an interior point.
875 : // Use a single vertex-average interior point and tetrahedralize around it
876 2175 : catch ( libMesh::LogicError & )
877 : #endif
878 : {
879 : // Clear the triangulation we started building
880 58 : this->_triangulation.clear();
881 :
882 : // Get the vertex-average, no need to triangulate for this
883 2059 : const auto v_avg = this->vertex_average();
884 2059 : if (!_has_mid_elem_node)
885 : {
886 : // Create the mid element node. Add it to nodelinks
887 : // We'll attach a unique ptr to it later
888 2059 : Node * mid_elem_node = new Node(v_avg);
889 2059 : this->_nodelinks_data.push_back(mid_elem_node);
890 2059 : _has_mid_elem_node = true;
891 : }
892 : else
893 : {
894 : // We are triangulating for a second time, we have already added this mid-element node
895 0 : libmesh_assert(n_vertices() == n_nodes() - 1);
896 : }
897 :
898 : // Build the tetrahedralization with each of the triangles on each side
899 14697 : for (unsigned int s : make_range(this->n_sides()))
900 : {
901 12638 : const auto & [side, inward_normal, node_map] = this->_sidelinks_data[s];
902 :
903 38482 : for (auto t : make_range(side->n_subtriangles()))
904 : {
905 : // Get all the nodes
906 25844 : const auto & n1 = node_map[side->subtriangle(t)[0]];
907 25844 : const auto & n2 = node_map[side->subtriangle(t)[1]];
908 25844 : const auto & n3 = node_map[side->subtriangle(t)[2]];
909 :
910 : // The mid node is the last node in the _nodes array
911 25844 : if (!inward_normal)
912 6106 : this->add_tet((int)n1, (int)n2, (int)n3, this->n_nodes() - 1);
913 : else
914 19738 : this->add_tet((int)n1, (int)n3, (int)n2, this->n_nodes() - 1);
915 : }
916 : }
917 1943 : }
918 11894 : }
919 :
920 :
921 61419 : void C0Polyhedron::add_tet(int n1,
922 : int n2,
923 : int n3,
924 : int n4)
925 : {
926 : #ifndef NDEBUG
927 1738 : const auto nn = this->n_nodes();
928 1738 : libmesh_assert_less(n1, nn);
929 1738 : libmesh_assert_less(n2, nn);
930 1738 : libmesh_assert_less(n3, nn);
931 1738 : libmesh_assert_less(n4, nn);
932 : #endif
933 :
934 63157 : const Point v12 = this->point(n2) - this->point(n1);
935 61419 : const Point v13 = this->point(n3) - this->point(n1);
936 61419 : const Point v14 = this->point(n4) - this->point(n1);
937 59681 : const Real six_vol = triple_product(v12, v13, v14);
938 : // We need to error on this in optimized modes to fall back onto the
939 : // tetrahedralization with a mid node
940 63478 : libmesh_error_msg_if(six_vol <= 0, "Creating flat tet");
941 :
942 59360 : this->_triangulation.push_back({n1, n2, n3, n4});
943 59360 : }
944 :
945 :
946 : } // namespace libMesh
|
