Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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 :
27 : // C++ headers
28 : #include <numeric> // std::iota
29 :
30 : namespace libMesh
31 : {
32 :
33 : // ------------------------------------------------------------
34 : // C0Polyhedron class member functions
35 :
36 3991 : C0Polyhedron::C0Polyhedron
37 3991 : (const std::vector<std::shared_ptr<Polygon>> & sides, Elem * p) :
38 3991 : Polyhedron(sides, p)
39 : {
40 3991 : this->retriangulate();
41 3991 : }
42 :
43 :
44 :
45 7526 : C0Polyhedron::~C0Polyhedron() = default;
46 :
47 :
48 :
49 15 : std::unique_ptr<Elem> C0Polyhedron::disconnected_clone() const
50 : {
51 19 : auto sides = this->side_clones();
52 :
53 : std::unique_ptr<Elem> returnval =
54 15 : std::make_unique<C0Polyhedron>(sides);
55 :
56 15 : returnval->set_id() = this->id();
57 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
58 15 : if (this->valid_unique_id())
59 2 : returnval->set_unique_id(this->unique_id());
60 : #endif
61 :
62 15 : const auto n_elem_ints = this->n_extra_integers();
63 15 : returnval->add_extra_integers(n_elem_ints);
64 15 : for (unsigned int i = 0; i != n_elem_ints; ++i)
65 0 : returnval->set_extra_integer(i, this->get_extra_integer(i));
66 :
67 13 : returnval->inherit_data_from(*this);
68 :
69 17 : return returnval;
70 11 : }
71 :
72 :
73 :
74 10440 : bool C0Polyhedron::is_vertex(const unsigned int libmesh_dbg_var(i)) const
75 : {
76 1208 : libmesh_assert (i < this->n_nodes());
77 :
78 10440 : return true;
79 : }
80 :
81 :
82 :
83 1704 : bool C0Polyhedron::is_edge(const unsigned int libmesh_dbg_var(i)) const
84 : {
85 48 : libmesh_assert_less (i, this->n_nodes());
86 :
87 1704 : return false;
88 : }
89 :
90 :
91 :
92 1704 : bool C0Polyhedron::is_face(const unsigned int libmesh_dbg_var(i)) const
93 : {
94 48 : libmesh_assert_less (i, this->n_nodes());
95 :
96 1704 : return false;
97 : }
98 :
99 :
100 :
101 2568 : bool C0Polyhedron::is_node_on_side(const unsigned int n,
102 : const unsigned int s) const
103 : {
104 120 : libmesh_assert_less (n, this->n_nodes());
105 120 : libmesh_assert_less (s, this->n_sides());
106 :
107 : const std::vector<unsigned int> & node_map =
108 2568 : std::get<2>(_sidelinks_data[s]);
109 :
110 2568 : const auto it = std::find(node_map.begin(), node_map.end(), n);
111 :
112 2688 : return (it != node_map.end());
113 : }
114 :
115 : std::vector<unsigned int>
116 1218 : C0Polyhedron::nodes_on_side(const unsigned int s) const
117 : {
118 78 : libmesh_assert_less(s, this->n_sides());
119 1296 : return std::get<2>(_sidelinks_data[s]);
120 : }
121 :
122 : std::vector<unsigned int>
123 720 : C0Polyhedron::nodes_on_edge(const unsigned int e) const
124 : {
125 720 : auto [s, se] = _edge_lookup[e];
126 720 : const Polygon & face = *std::get<0>(_sidelinks_data[s]);
127 : const std::vector<unsigned int> & node_map =
128 60 : std::get<2>(_sidelinks_data[s]);
129 : std::vector<unsigned int> nodes_on_edge =
130 720 : face.nodes_on_side(se);
131 2160 : for (auto i : index_range(nodes_on_edge))
132 1560 : nodes_on_edge[i] = node_map[nodes_on_edge[i]];
133 780 : return nodes_on_edge;
134 : }
135 :
136 :
137 :
138 288 : bool C0Polyhedron::is_node_on_edge(const unsigned int n,
139 : const unsigned int e) const
140 : {
141 24 : libmesh_assert_less(e, _edge_lookup.size());
142 288 : auto [s, se] = _edge_lookup[e];
143 :
144 288 : const Polygon & face = *std::get<0>(_sidelinks_data[s]);
145 : const std::vector<unsigned int> & node_map =
146 24 : std::get<2>(_sidelinks_data[s]);
147 : std::vector<unsigned int> nodes_on_edge =
148 312 : face.nodes_on_side(se);
149 432 : for (auto noe : nodes_on_edge)
150 468 : if (node_map[noe] == n)
151 24 : return true;
152 :
153 0 : return false;
154 : }
155 :
156 :
157 :
158 0 : void C0Polyhedron::connectivity(const unsigned int /*sf*/,
159 : const IOPackage /*iop*/,
160 : std::vector<dof_id_type> & /*conn*/) const
161 : {
162 0 : libmesh_not_implemented();
163 : }
164 :
165 :
166 :
167 83 : Real C0Polyhedron::volume () const
168 : {
169 : // This specialization is good for Lagrange mappings only
170 83 : if (this->mapping_type() != LAGRANGE_MAP)
171 0 : return this->Elem::volume();
172 :
173 : // We use a triangulation to calculate here.
174 :
175 3 : Real six_vol = 0;
176 498 : for (const auto & subtet : this->_triangulation)
177 : {
178 415 : const Point p0 = this->point(subtet[0]);
179 415 : const Point p1 = this->point(subtet[1]);
180 415 : const Point p2 = this->point(subtet[2]);
181 430 : const Point p3 = this->point(subtet[3]);
182 :
183 15 : const Point v01 = p1 - p0;
184 15 : const Point v02 = p2 - p0;
185 15 : const Point v03 = p3 - p0;
186 :
187 415 : six_vol += triple_product(v01, v02, v03);
188 : }
189 :
190 83 : return six_vol * (1./6.);
191 : }
192 :
193 :
194 :
195 24 : Point C0Polyhedron::true_centroid () const
196 : {
197 : // This specialization is good for Lagrange mappings only
198 24 : if (this->mapping_type() != LAGRANGE_MAP)
199 0 : return this->Elem::true_centroid();
200 :
201 2 : Real six_vol = 0;
202 2 : Point six_vol_weighted_centroid;
203 144 : for (const auto & subtet : this->_triangulation)
204 : {
205 120 : const Point p0 = this->point(subtet[0]);
206 120 : const Point p1 = this->point(subtet[1]);
207 120 : const Point p2 = this->point(subtet[2]);
208 130 : const Point p3 = this->point(subtet[3]);
209 :
210 10 : const Point v01 = p1 - p0;
211 10 : const Point v02 = p2 - p0;
212 10 : const Point v03 = p3 - p0;
213 :
214 110 : const Real six_tet_vol = triple_product(v01, v02, v03);
215 :
216 10 : const Point tet_centroid = (p0 + p1 + p2 + p3) * 0.25;
217 :
218 120 : six_vol += six_tet_vol;
219 :
220 10 : six_vol_weighted_centroid += six_tet_vol * tet_centroid;
221 : }
222 :
223 2 : return six_vol_weighted_centroid / six_vol;
224 : }
225 :
226 :
227 :
228 : std::pair<unsigned short int, unsigned short int>
229 0 : C0Polyhedron::second_order_child_vertex (const unsigned int /*n*/) const
230 : {
231 0 : libmesh_not_implemented();
232 : return std::pair<unsigned short int, unsigned short int> (0,0);
233 : }
234 :
235 :
236 :
237 216 : ElemType C0Polyhedron::side_type (const unsigned int libmesh_dbg_var(s)) const
238 : {
239 18 : libmesh_assert_less (s, this->n_sides());
240 216 : return C0POLYGON;
241 : }
242 :
243 :
244 :
245 : Point
246 852 : C0Polyhedron::side_vertex_average_normal(const unsigned int s) const
247 : {
248 24 : libmesh_assert_less (s, this->n_sides());
249 24 : libmesh_assert_equal_to(this->mapping_type(), LAGRANGE_MAP);
250 852 : const auto poly_side_ptr = this->side_ptr(s);
251 852 : const auto n_side_edges = poly_side_ptr->n_sides();
252 :
253 : // At the side vertex average, things simplify a bit
254 : // We get the side "plane" normal at all vertices, then average them
255 24 : Point normal;
256 48 : Point current_edge = poly_side_ptr->point(1) - poly_side_ptr->point(0);
257 852 : for (auto i : make_range(n_side_edges))
258 : {
259 852 : const Point next_edge = poly_side_ptr->point((i + 2) % n_side_edges) -
260 876 : poly_side_ptr->point((i + 1) % n_side_edges);
261 828 : const Point normal_at_vertex = current_edge.cross(next_edge);
262 24 : normal += normal_at_vertex;
263 : // Note: the sides are planar, we don't need to test them all
264 852 : if (normal.norm_sq() > TOLERANCE)
265 24 : break;
266 0 : current_edge = next_edge;
267 : }
268 852 : bool outward_normal = std::get<1>(_sidelinks_data[s]);
269 1708 : return (outward_normal ? 1. : -1.) * normal.unit();
270 804 : }
271 :
272 :
273 :
274 3991 : void C0Polyhedron::retriangulate()
275 : {
276 114 : this->_triangulation.clear();
277 :
278 : // We should already have a triangulation for each side. We'll turn
279 : // that into a triangulation of the entire polyhedral surface
280 : // (stored as a mesh, so we can use `find_neighbors()`, then turn
281 : // that into a tetrahedralization by shrinking the volume one tet at
282 : // a time, which should work fine for convex polyhedra.
283 :
284 228 : Parallel::Communicator comm_self; // the default communicator is 1 rank
285 4219 : ReplicatedMesh surface(comm_self);
286 :
287 : // poly_node_to_local_id[poly_node] is the local id of
288 : // poly_node in the Polyhedron, which is also the global id of the
289 : // corresponding node in the surface mesh.
290 228 : std::unordered_map<Node *, unsigned int> poly_node_to_id;
291 :
292 27937 : for (unsigned int s : make_range(this->n_sides()))
293 : {
294 23946 : const auto & [side, inward_normal, node_map] = this->_sidelinks_data[s];
295 :
296 71838 : for (auto t : make_range(side->n_subtriangles()))
297 : {
298 49260 : Elem * tri = surface.add_elem(Elem::build(TRI3));
299 :
300 47892 : const std::array<int, 3> subtri = side->subtriangle(t);
301 :
302 191568 : for (int i : make_range(3))
303 : {
304 143676 : const int side_id = subtri[i];
305 143676 : const Node * poly_node = side->node_ptr(side_id);
306 :
307 4104 : libmesh_assert_less(side_id, node_map.size());
308 143676 : const unsigned int local_id = node_map[side_id];
309 :
310 143676 : Node * surf_node = surface.query_node_ptr(local_id);
311 143676 : if (surf_node)
312 3192 : libmesh_assert_equal_to(*(const Point*)poly_node,
313 : *(const Point*)surf_node);
314 : else
315 31928 : surf_node = surface.add_point(*poly_node, local_id);
316 :
317 : // Get a consistent orientation for surface triangles,
318 : // facing their zeta directions outward
319 143676 : const int tri_node = inward_normal ? i : 2-i;
320 143676 : tri->set_node(tri_node, surf_node);
321 : }
322 : }
323 : }
324 :
325 114 : surface.allow_renumbering(false);
326 3991 : surface.prepare_for_use();
327 :
328 : // We should have a watertight surface with consistent triangle
329 : // orientations. That's expensive to check.
330 : #ifdef DEBUG
331 14250 : auto verify_surface = [& surface] ()
332 : {
333 5130 : for (const Elem * elem : surface.element_ptr_range())
334 : {
335 18240 : for (auto s : make_range(3))
336 : {
337 13680 : const Elem * neigh = elem->neighbor_ptr(s);
338 13680 : libmesh_assert(neigh);
339 13680 : libmesh_assert_equal_to(neigh,
340 : surface.elem_ptr(neigh->id()));
341 13680 : const unsigned int ns = neigh->which_neighbor_am_i(elem);
342 13680 : libmesh_assert_less(ns, 3);
343 13680 : libmesh_assert_equal_to(elem->node_ptr(s),
344 : neigh->node_ptr((ns+1)%3));
345 13680 : libmesh_assert_equal_to(elem->node_ptr((s+1)%3),
346 : neigh->node_ptr(ns));
347 : }
348 : }
349 570 : };
350 :
351 114 : verify_surface();
352 : #endif
353 :
354 : // We'll have to edit this as we change the surface elements, but we
355 : // have a method to initialize it easily.
356 342 : std::vector<std::vector<dof_id_type>> nodes_to_elem_vec_map;
357 3991 : MeshTools::build_nodes_to_elem_map(surface, nodes_to_elem_vec_map);
358 :
359 : // We'll be inserting and deleting entries heavily, so we'll use
360 : // sets rather than vectors. We want to get the same results in
361 : // parallel, so we'll use element ids rather than Elem pointers
362 342 : std::vector<std::set<dof_id_type>> nodes_to_elem_map;
363 35919 : for (auto i : index_range(nodes_to_elem_vec_map))
364 : nodes_to_elem_map.emplace_back
365 62944 : (nodes_to_elem_vec_map[i].begin(),
366 34664 : nodes_to_elem_vec_map[i].end());
367 :
368 : // Now start meshing the volume enclosed by surface, one tet at a
369 : // time, with a greedy heuristic: find the vertex node with the most
370 : // acutely convex (solid) angle and strip it out as tetrahedra.
371 :
372 : // We'll want a vector of surrounding nodes for multiple uses,
373 : // sometimes with a similarly-sorted vector of surrounding elements
374 : // to go with it.
375 : auto surroundings_of =
376 48919 : [&nodes_to_elem_map, & surface]
377 : (const Node & node,
378 62143 : std::vector<Elem *> * surrounding_elems)
379 : {
380 : const std::set<dof_id_type> & elems_by_node =
381 3705 : nodes_to_elem_map[node.id()];
382 :
383 51883 : const unsigned int n_surrounding = elems_by_node.size();
384 1482 : libmesh_assert_greater_equal(n_surrounding, 3);
385 :
386 51883 : if (surrounding_elems)
387 : {
388 570 : libmesh_assert(surrounding_elems->empty());
389 19955 : surrounding_elems->resize(n_surrounding);
390 : }
391 :
392 51883 : std::vector<Node *> surrounding_nodes(n_surrounding);
393 :
394 51883 : Elem * elem = surface.elem_ptr(*elems_by_node.begin());
395 255424 : for (auto i : make_range(n_surrounding))
396 : {
397 197727 : const unsigned int n = elem->get_node_index(&node);
398 5814 : libmesh_assert_not_equal_to(n, invalid_uint);
399 203541 : Node * next_node = elem->node_ptr((n+1)%3);
400 203541 : surrounding_nodes[i] = next_node;
401 203541 : if (surrounding_elems)
402 61575 : (*surrounding_elems)[i] = elem;
403 203541 : elem = elem->neighbor_ptr((n+2)%3);
404 5814 : libmesh_assert(elem);
405 5814 : libmesh_assert_equal_to(elem, surface.elem_ptr(elem->id()));
406 :
407 : // We should have a manifold here, but verifying that is
408 : // expensive
409 : #ifdef DEBUG
410 5814 : libmesh_assert_equal_to
411 : (std::count(surrounding_nodes.begin(),
412 : surrounding_nodes.end(), next_node),
413 : 1);
414 : #endif
415 : }
416 :
417 : // We should have finished a loop
418 1482 : libmesh_assert_equal_to
419 : (elem, surface.elem_ptr(*elems_by_node.begin()));
420 :
421 53365 : return surrounding_nodes;
422 3991 : };
423 :
424 31928 : auto geometry_at = [&surroundings_of](const Node & node)
425 : {
426 : const std::vector<Node *> surrounding_nodes =
427 32840 : surroundings_of(node, nullptr);
428 :
429 : // Now sum up solid angles from tetrahedra created from the
430 : // trivial triangulation of the surrounding nodes loop
431 912 : Real total_solid_angle = 0;
432 : const int n_surrounding =
433 1824 : cast_int<int>(surrounding_nodes.size());
434 :
435 111748 : for (auto n : make_range(n_surrounding-2))
436 : {
437 : const Point
438 82100 : v01 = static_cast<Point>(*surrounding_nodes[n]) - node,
439 82100 : v02 = static_cast<Point>(*surrounding_nodes[n+1]) - node,
440 82100 : v03 = static_cast<Point>(*surrounding_nodes[n+2]) - node;
441 :
442 79820 : total_solid_angle += solid_angle(v01, v02, v03);
443 : }
444 :
445 32840 : return std::make_pair(n_surrounding, total_solid_angle);
446 3877 : };
447 :
448 : // We'll keep track of solid angles and node valences so we don't
449 : // waste time recomputing them when they haven't changed. We need
450 : // to be able to search efficiently for the smallest angles of each
451 : // valence, but also search efficiently for a particular node to
452 : // remove and reinsert it when its connectivity changes.
453 : //
454 : // Since C++11 multimap has guaranteed that pairs with matching keys
455 : // are kept in insertion order, so we can use Node * for values even
456 : // in parallel.
457 : typedef std::multimap<std::pair<int, Real>, Node*> node_map_type;
458 228 : node_map_type nodes_by_geometry;
459 228 : std::map<Node *, node_map_type::iterator> node_index;
460 :
461 39796 : for (auto node : surface.node_ptr_range())
462 31928 : node_index[node] =
463 67619 : nodes_by_geometry.emplace(geometry_at(*node), node);
464 :
465 : // In 3D, this will require nested loops: an outer loop to remove
466 : // each vertex, and an inner loop to remove multiple tetrahedra in
467 : // cases where the vertex has more than 3 neighboring triangles.
468 :
469 : // We'll be done when there are only three "unremoved" nodes left,
470 : // so they don't actually enclose any volume.
471 23946 : for (auto i : make_range(nodes_by_geometry.size()-3))
472 : {
473 570 : auto geometry_it = nodes_by_geometry.begin();
474 19955 : auto geometry_key = geometry_it->first;
475 19955 : auto [valence, angle] = geometry_key;
476 19955 : Node * node = geometry_it->second;
477 570 : libmesh_ignore(i);
478 :
479 : // If our lowest-valence nodes are all points of non-convexity,
480 : // skip to a higher valence.
481 19955 : while (angle > 2*pi-TOLERANCE)
482 : {
483 : geometry_it =
484 : nodes_by_geometry.upper_bound
485 0 : (std::make_pair(valence, Real(100)));
486 0 : libmesh_assert(geometry_it != nodes_by_geometry.end());
487 :
488 0 : std::tie(geometry_key, node) = *geometry_it;
489 0 : std::tie(valence, angle) = geometry_key;
490 : }
491 :
492 1140 : std::vector<Elem *> surrounding_elems;
493 : std::vector<Node *> surrounding_nodes =
494 20525 : surroundings_of(*node, &surrounding_elems);
495 :
496 1140 : const std::size_t n_surrounding = surrounding_nodes.size();
497 :
498 : // As we separate surrounding nodes from our center node, we'll
499 : // be marking them as nullptr; we still need to be able to find
500 : // predecessor and successor nodes in order afterward.
501 : auto find_valid_nodes_around =
502 75260 : [n_surrounding, & surrounding_nodes]
503 166480 : (unsigned int j)
504 : {
505 79820 : unsigned int jnext = (j+1)%n_surrounding;
506 82100 : while (!surrounding_nodes[jnext])
507 0 : jnext = (jnext+1)%n_surrounding;
508 :
509 79820 : unsigned int jprev = (j+n_surrounding-1)%n_surrounding;
510 82100 : while (!surrounding_nodes[jprev])
511 0 : jprev = (jprev+n_surrounding-1)%n_surrounding;
512 :
513 82100 : return std::make_pair(jprev, jnext);
514 19955 : };
515 :
516 : // We may have too many surrounding nodes to handle with
517 : // just one tet. In that case we'll keep a cache of the
518 : // element qualities that we'd get by making a tet with the
519 : // edge from the center node to each surrounding node, so we
520 : // can build the best tets first.
521 : //
522 : // In the case where we just have 3 nodes, we'll just pretend
523 : // they all have the same positive quality, so we can still
524 : // search this vector.
525 20525 : std::vector<Real> local_tet_quality(n_surrounding, 1);
526 :
527 : // From our center node with N surrounding nodes we can make N-2
528 : // tetrahedra. The first N-3 each replace two surface tets with
529 : // two new surface tets.
530 : //
531 : // My first idea was to greedily pick nodes with the smallest
532 : // local (solid) angles to get the best quality. This works in
533 : // 2D, but such a node can give a pancake tet in 3D.
534 : //
535 : // My second idea was to greedily pick nodes with the highest
536 : // prospective tet quality. This works for the first tets, but
537 : // can leave a pancake tet behind.
538 : //
539 : // My third idea is to try to fix the lowest quality tets first,
540 : // by picking cases where they have higher quality neighbors,
541 : // and creating those neighbors so as to change them.
542 :
543 : auto find_new_tet_nodes =
544 18815 : [& local_tet_quality, & find_valid_nodes_around]
545 35345 : ()
546 : {
547 570 : unsigned int jbest = 0;
548 19955 : auto [jminus, jplus] = find_valid_nodes_around(jbest);
549 19955 : Real qneighbest = std::min(local_tet_quality[jminus],
550 20525 : local_tet_quality[jplus]);
551 19955 : for (auto j : make_range(std::size_t(1),
552 61005 : local_tet_quality.size()))
553 : {
554 : // We don't want to build a bad tet
555 41050 : if (local_tet_quality[j] <= 0)
556 0 : continue;
557 :
558 39910 : std::tie(jminus, jplus) = find_valid_nodes_around(j);
559 39910 : Real qneighj = std::min(local_tet_quality[jminus],
560 41050 : local_tet_quality[jplus]);
561 :
562 : // We don't want to build a tet that can't fix a neighbor
563 : // if we can build one that can.
564 39910 : if (qneighbest <= 0 &&
565 : qneighj > 0)
566 0 : continue;
567 :
568 : // We want to try for the best possible fix.
569 41050 : if ((local_tet_quality[j] - qneighj) >
570 41050 : (local_tet_quality[jbest] - qneighj))
571 : {
572 0 : jbest = j;
573 0 : qneighbest = qneighj;
574 : }
575 : }
576 :
577 20525 : libmesh_error_msg_if
578 : (local_tet_quality[jbest] <= 0,
579 : "Cannot build non-singular non-inverted tet");
580 :
581 19955 : std::tie(jminus, jplus) = find_valid_nodes_around(jbest);
582 :
583 20525 : return std::make_tuple(jbest, jminus, jplus);
584 19955 : };
585 :
586 19955 : if (n_surrounding > 3)
587 : {
588 : // We'll be searching local_tet_quality even after tet
589 : // extraction disconnects us from some nodes; when we do we
590 : // don't want to get one.
591 0 : constexpr Real far_node = -1e6;
592 :
593 : // Vectors from the center node to each of its surrounding
594 : // nodes are helpful for calculating prospective tet
595 : // quality.
596 0 : std::vector<Point> v0s(n_surrounding);
597 0 : for (auto j : make_range(n_surrounding))
598 0 : v0s[j] = *(Point *)surrounding_nodes[j] - *node;
599 :
600 : // Find the tet quality we'd potentially get from each
601 : // possible choice of tet
602 : auto local_tet_quality_of =
603 0 : [& surrounding_nodes, & v0s, & find_valid_nodes_around]
604 0 : (unsigned int j)
605 : {
606 0 : auto [jminus, jplus] = find_valid_nodes_around(j);
607 :
608 : // Anything proportional to the ratio of volume to
609 : // total-edge-length-cubed should peak for perfect tets
610 : // but hit 0 for pancakes and slivers.
611 :
612 : const Real total_len =
613 0 : v0s[j].norm() + v0s[jminus].norm() + v0s[jplus].norm() +
614 0 : (*(Point *)surrounding_nodes[jplus] -
615 0 : *(Point *)surrounding_nodes[j]).norm() +
616 0 : (*(Point *)surrounding_nodes[j] -
617 0 : *(Point *)surrounding_nodes[jminus]).norm() +
618 0 : (*(Point *)surrounding_nodes[jminus] -
619 0 : *(Point *)surrounding_nodes[jplus]).norm();
620 :
621 : // Orientation here is tricky. Think of the triple
622 : // product as (v1 cross v2) dot v3, with right hand rule.
623 : const Real six_vol =
624 0 : triple_product(v0s[jminus], v0s[jplus], v0s[j]);
625 :
626 0 : return six_vol / (total_len * total_len * total_len);
627 0 : };
628 :
629 0 : for (auto j : make_range(n_surrounding))
630 0 : local_tet_quality[j] = local_tet_quality_of(j);
631 :
632 : // If we have N surrounding nodes, we can make N tets and
633 : // that'll bring us back to the 3-surrounding-node case to
634 : // finish.
635 0 : for (auto t : make_range(n_surrounding-3))
636 : {
637 0 : libmesh_ignore(t);
638 :
639 0 : auto [jbest, jminus, jplus] = find_new_tet_nodes();
640 :
641 : // Turn these four nodes into a tet
642 0 : Node * nbest = surrounding_nodes[jbest],
643 0 : * nminus = surrounding_nodes[jminus],
644 0 : * nplus = surrounding_nodes[jplus];
645 0 : this->add_tet(nminus->id(), nbest->id(), nplus->id(),
646 0 : node->id());
647 :
648 : // Replace the old two triangles around these nodes with the
649 : // proper two new triangles.
650 0 : Elem * oldtri1 = surrounding_elems[jminus],
651 0 : * oldtri2 = surrounding_elems[jbest],
652 0 : * newtri1 = surface.add_elem(Elem::build(TRI3)),
653 0 : * newtri2 = surface.add_elem(Elem::build(TRI3));
654 :
655 0 : const unsigned int c1 = oldtri1->get_node_index(node),
656 0 : c2 = oldtri2->get_node_index(node);
657 :
658 0 : newtri1->set_node(0, node);
659 0 : newtri1->set_node(1, nminus);
660 0 : newtri1->set_node(2, nplus);
661 :
662 0 : surrounding_elems[jminus] = newtri1;
663 :
664 0 : Elem * neigh10 = oldtri1->neighbor_ptr(c1);
665 0 : Elem * neigh12 = oldtri2->neighbor_ptr((c2+2)%3);
666 0 : newtri1->set_neighbor(0, neigh10);
667 0 : neigh10->set_neighbor(neigh10->which_neighbor_am_i(oldtri1), newtri1);
668 0 : newtri1->set_neighbor(1, newtri2);
669 0 : newtri1->set_neighbor(2, neigh12);
670 0 : neigh12->set_neighbor(neigh12->which_neighbor_am_i(oldtri2), newtri1);
671 :
672 0 : newtri2->set_node(0, nplus);
673 0 : newtri2->set_node(1, nminus);
674 0 : newtri2->set_node(2, nbest);
675 :
676 0 : Elem * neigh21 = oldtri1->neighbor_ptr((c1+1)%3);
677 0 : Elem * neigh22 = oldtri2->neighbor_ptr((c2+1)%3);
678 0 : newtri2->set_neighbor(0, newtri1);
679 0 : newtri2->set_neighbor(1, neigh21);
680 0 : neigh21->set_neighbor(neigh21->which_neighbor_am_i(oldtri1), newtri2);
681 0 : newtri2->set_neighbor(2, neigh22);
682 0 : neigh22->set_neighbor(neigh22->which_neighbor_am_i(oldtri2), newtri2);
683 :
684 0 : for (unsigned int p : make_range(3))
685 : {
686 0 : nodes_to_elem_map[oldtri1->node_id(p)].erase(oldtri1->id());
687 0 : nodes_to_elem_map[oldtri2->node_id(p)].erase(oldtri2->id());
688 0 : nodes_to_elem_map[newtri1->node_id(p)].insert(newtri1->id());
689 0 : nodes_to_elem_map[newtri2->node_id(p)].insert(newtri2->id());
690 : }
691 :
692 : // We've finished replacing the old triangles.
693 0 : surface.delete_elem(oldtri1);
694 0 : surface.delete_elem(oldtri2);
695 :
696 : // The solid angle for the far node should now stay
697 : // unchanged until we're out of this inner loop; let's
698 : // recalculate it here, and then we'll be done with it.
699 0 : Node * & nbestref = surrounding_nodes[jbest];
700 0 : nodes_by_geometry.erase(node_index[nbestref]);
701 0 : node_index[nbestref] =
702 0 : nodes_by_geometry.emplace(geometry_at(*nbestref), nbestref);
703 :
704 : // The far node is no longer sharing an edge with our center
705 : // node. Make sure we don't use it again with the center
706 : // node.
707 0 : local_tet_quality[jbest] = far_node;
708 0 : nbestref = nullptr;
709 :
710 : // The potential tet qualities using the side nodes have
711 : // changed now that they're directly connected to each
712 : // other.
713 0 : local_tet_quality[jminus] =
714 0 : local_tet_quality_of(jminus);
715 :
716 0 : local_tet_quality[jplus] =
717 0 : local_tet_quality_of(jplus);
718 : }
719 : }
720 :
721 : // Now we should have just 3 surrounding nodes, with which to
722 : // make one tetrahedron. Put them in a counterclockwise
723 : // (looking from outside) orientation, not the "best, clockwise,
724 : // counterclockwise" we get from the lambda.
725 19955 : auto [j2, j1, j3] = find_new_tet_nodes();
726 :
727 : // Turn these last four nodes into a tet
728 19955 : Node * n1 = surrounding_nodes[j1],
729 19955 : * n2 = surrounding_nodes[j2],
730 19955 : * n3 = surrounding_nodes[j3];
731 19955 : this->add_tet(n1->id(), n2->id(), n3->id(), node->id());
732 :
733 : // Replace the three surface triangles of that tet with the new
734 : // fourth triangle.
735 19955 : Elem * oldtri1 = surrounding_elems[j1],
736 19955 : * oldtri2 = surrounding_elems[j2],
737 19955 : * oldtri3 = surrounding_elems[j3],
738 20525 : * newtri = surface.add_elem(Elem::build(TRI3));
739 :
740 19385 : const unsigned int c1 = oldtri1->get_node_index(node),
741 19385 : c2 = oldtri2->get_node_index(node),
742 19385 : c3 = oldtri3->get_node_index(node);
743 :
744 19955 : newtri->set_node(0, n1);
745 19955 : newtri->set_node(1, n2);
746 19955 : newtri->set_node(2, n3);
747 19955 : Elem * neigh0 = oldtri1->neighbor_ptr((c1+1)%3);
748 1140 : newtri->set_neighbor(0, neigh0);
749 19955 : neigh0->set_neighbor(neigh0->which_neighbor_am_i(oldtri1), newtri);
750 19955 : Elem * neigh1 = oldtri2->neighbor_ptr((c2+1)%3);
751 1140 : newtri->set_neighbor(1, neigh1);
752 19955 : neigh1->set_neighbor(neigh1->which_neighbor_am_i(oldtri2), newtri);
753 19955 : Elem * neigh2 = oldtri3->neighbor_ptr((c3+1)%3);
754 1140 : newtri->set_neighbor(2, neigh2);
755 19955 : neigh2->set_neighbor(neigh2->which_neighbor_am_i(oldtri3), newtri);
756 :
757 79820 : for (unsigned int p : make_range(3))
758 : {
759 63285 : nodes_to_elem_map[oldtri1->node_id(p)].erase(oldtri1->id());
760 63285 : nodes_to_elem_map[oldtri2->node_id(p)].erase(oldtri2->id());
761 63285 : nodes_to_elem_map[oldtri3->node_id(p)].erase(oldtri3->id());
762 63285 : nodes_to_elem_map[newtri->node_id(p)].insert(newtri->id());
763 : }
764 :
765 : // We've finished replacing the old triangles.
766 19955 : surface.delete_elem(oldtri1);
767 19955 : surface.delete_elem(oldtri2);
768 19955 : surface.delete_elem(oldtri3);
769 :
770 : // We should have used up all our surrounding nodes now, and we
771 : // shouldn't have messed up our surface in the process, and our
772 : // center node should no longer be part of the surface.
773 : #ifdef DEBUG
774 570 : surrounding_nodes[j1] = nullptr;
775 570 : surrounding_nodes[j2] = nullptr;
776 570 : surrounding_nodes[j3] = nullptr;
777 :
778 2280 : for (auto ltq : surrounding_nodes)
779 1710 : libmesh_assert(!ltq);
780 :
781 570 : if (surface.n_elem() > 3)
782 456 : verify_surface();
783 :
784 3990 : for (const Elem * elem : surface.element_ptr_range())
785 13680 : for (auto p : make_range(3))
786 10260 : libmesh_assert_not_equal_to
787 : (elem->node_ptr(p), node);
788 : #endif
789 :
790 : // We've used up our center node, so it's not something we can
791 : // eliminate again.
792 19385 : nodes_by_geometry.erase(geometry_it);
793 : }
794 :
795 : // At this point our surface should just have two triangles left.
796 114 : libmesh_assert_equal_to(surface.n_elem(), 2);
797 11517 : }
798 :
799 :
800 19955 : void C0Polyhedron::add_tet(int n1,
801 : int n2,
802 : int n3,
803 : int n4)
804 : {
805 : #ifndef NDEBUG
806 570 : const auto nn = this->n_nodes();
807 570 : libmesh_assert_less(n1, nn);
808 570 : libmesh_assert_less(n2, nn);
809 570 : libmesh_assert_less(n3, nn);
810 570 : libmesh_assert_less(n4, nn);
811 :
812 570 : const Point v12 = this->point(n2) - this->point(n1);
813 570 : const Point v13 = this->point(n3) - this->point(n1);
814 570 : const Point v14 = this->point(n4) - this->point(n1);
815 570 : const Real six_vol = triple_product(v12, v13, v14);
816 570 : libmesh_assert_greater(six_vol, Real(0));
817 : #endif
818 :
819 19955 : this->_triangulation.push_back({n1, n2, n3, n4});
820 19955 : }
821 :
822 :
823 : } // namespace libMesh
|
