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