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 3849 : C0Polyhedron::C0Polyhedron
37 3849 : (const std::vector<std::shared_ptr<Polygon>> & sides, Elem * p) :
38 3849 : Polyhedron(sides, p)
39 : {
40 3849 : this->retriangulate();
41 3849 : }
42 :
43 :
44 :
45 7258 : 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 3849 : void C0Polyhedron::retriangulate()
246 : {
247 110 : this->_triangulation.clear();
248 :
249 : // We should already have a triangulation for each side. We'll turn
250 : // that into a triangulation of the entire polyhedral surface
251 : // (stored as a mesh, so we can use `find_neighbors()`, then turn
252 : // that into a tetrahedralization by shrinking the volume one tet at
253 : // a time, which should work fine for convex polyhedra.
254 :
255 220 : Parallel::Communicator comm_self; // the default communicator is 1 rank
256 4069 : ReplicatedMesh surface(comm_self);
257 :
258 : // poly_node_to_local_id[poly_node] is the local id of
259 : // poly_node in the Polyhedron, which is also the global id of the
260 : // corresponding node in the surface mesh.
261 220 : std::unordered_map<Node *, unsigned int> poly_node_to_id;
262 :
263 26943 : for (unsigned int s : make_range(this->n_sides()))
264 : {
265 23094 : const auto & [side, inward_normal, node_map] = this->_sidelinks_data[s];
266 :
267 69282 : for (auto t : make_range(side->n_subtriangles()))
268 : {
269 47508 : Elem * tri = surface.add_elem(Elem::build(TRI3));
270 :
271 46188 : const std::array<int, 3> subtri = side->subtriangle(t);
272 :
273 184752 : for (int i : make_range(3))
274 : {
275 138564 : const int side_id = subtri[i];
276 138564 : const Node * poly_node = side->node_ptr(side_id);
277 :
278 3960 : libmesh_assert_less(side_id, node_map.size());
279 138564 : const unsigned int local_id = node_map[side_id];
280 :
281 138564 : Node * surf_node = surface.query_node_ptr(local_id);
282 138564 : if (surf_node)
283 3080 : libmesh_assert_equal_to(*(const Point*)poly_node,
284 : *(const Point*)surf_node);
285 : else
286 30792 : surf_node = surface.add_point(*poly_node, local_id);
287 :
288 : // Get a consistent orientation for surface triangles,
289 : // facing their zeta directions outward
290 138564 : const int tri_node = inward_normal ? i : 2-i;
291 138564 : tri->set_node(tri_node, surf_node);
292 : }
293 : }
294 : }
295 :
296 110 : surface.allow_renumbering(false);
297 3849 : surface.prepare_for_use();
298 :
299 : // We should have a watertight surface with consistent triangle
300 : // orientations. That's expensive to check.
301 : #ifdef DEBUG
302 13750 : auto verify_surface = [& surface] ()
303 : {
304 4950 : for (const Elem * elem : surface.element_ptr_range())
305 : {
306 17600 : for (auto s : make_range(3))
307 : {
308 13200 : const Elem * neigh = elem->neighbor_ptr(s);
309 13200 : libmesh_assert(neigh);
310 13200 : libmesh_assert_equal_to(neigh,
311 : surface.elem_ptr(neigh->id()));
312 13200 : const unsigned int ns = neigh->which_neighbor_am_i(elem);
313 13200 : libmesh_assert_less(ns, 3);
314 13200 : libmesh_assert_equal_to(elem->node_ptr(s),
315 : neigh->node_ptr((ns+1)%3));
316 13200 : libmesh_assert_equal_to(elem->node_ptr((s+1)%3),
317 : neigh->node_ptr(ns));
318 : }
319 : }
320 550 : };
321 :
322 110 : verify_surface();
323 : #endif
324 :
325 : // We'll have to edit this as we change the surface elements, but we
326 : // have a method to initialize it easily.
327 330 : std::vector<std::vector<dof_id_type>> nodes_to_elem_vec_map;
328 3849 : MeshTools::build_nodes_to_elem_map(surface, nodes_to_elem_vec_map);
329 :
330 : // We'll be inserting and deleting entries heavily, so we'll use
331 : // sets rather than vectors. We want to get the same results in
332 : // parallel, so we'll use element ids rather than Elem pointers
333 330 : std::vector<std::set<dof_id_type>> nodes_to_elem_map;
334 34641 : for (auto i : index_range(nodes_to_elem_vec_map))
335 : nodes_to_elem_map.emplace_back
336 60704 : (nodes_to_elem_vec_map[i].begin(),
337 33432 : nodes_to_elem_vec_map[i].end());
338 :
339 : // Now start meshing the volume enclosed by surface, one tet at a
340 : // time, with a greedy heuristic: find the vertex node with the most
341 : // acutely convex (solid) angle and strip it out as tetrahedra.
342 :
343 : // We'll want a vector of surrounding nodes for multiple uses,
344 : // sometimes with a similarly-sorted vector of surrounding elements
345 : // to go with it.
346 : auto surroundings_of =
347 47177 : [&nodes_to_elem_map, & surface]
348 : (const Node & node,
349 59937 : std::vector<Elem *> * surrounding_elems)
350 : {
351 : const std::set<dof_id_type> & elems_by_node =
352 3575 : nodes_to_elem_map[node.id()];
353 :
354 50037 : const unsigned int n_surrounding = elems_by_node.size();
355 1430 : libmesh_assert_greater_equal(n_surrounding, 3);
356 :
357 50037 : if (surrounding_elems)
358 : {
359 550 : libmesh_assert(surrounding_elems->empty());
360 19245 : surrounding_elems->resize(n_surrounding);
361 : }
362 :
363 50037 : std::vector<Node *> surrounding_nodes(n_surrounding);
364 :
365 50037 : Elem * elem = surface.elem_ptr(*elems_by_node.begin());
366 246336 : for (auto i : make_range(n_surrounding))
367 : {
368 190689 : const unsigned int n = elem->get_node_index(&node);
369 5610 : libmesh_assert_not_equal_to(n, invalid_uint);
370 196299 : Node * next_node = elem->node_ptr((n+1)%3);
371 196299 : surrounding_nodes[i] = next_node;
372 196299 : if (surrounding_elems)
373 59385 : (*surrounding_elems)[i] = elem;
374 196299 : elem = elem->neighbor_ptr((n+2)%3);
375 5610 : libmesh_assert(elem);
376 5610 : libmesh_assert_equal_to(elem, surface.elem_ptr(elem->id()));
377 :
378 : // We should have a manifold here, but verifying that is
379 : // expensive
380 : #ifdef DEBUG
381 5610 : libmesh_assert_equal_to
382 : (std::count(surrounding_nodes.begin(),
383 : surrounding_nodes.end(), next_node),
384 : 1);
385 : #endif
386 : }
387 :
388 : // We should have finished a loop
389 1430 : libmesh_assert_equal_to
390 : (elem, surface.elem_ptr(*elems_by_node.begin()));
391 :
392 51467 : return surrounding_nodes;
393 3849 : };
394 :
395 30792 : auto geometry_at = [&surroundings_of](const Node & node)
396 : {
397 : const std::vector<Node *> surrounding_nodes =
398 31672 : surroundings_of(node, nullptr);
399 :
400 : // Now sum up solid angles from tetrahedra created from the
401 : // trivial triangulation of the surrounding nodes loop
402 880 : Real total_solid_angle = 0;
403 : const int n_surrounding =
404 1760 : cast_int<int>(surrounding_nodes.size());
405 :
406 107772 : for (auto n : make_range(n_surrounding-2))
407 : {
408 : const Point
409 79180 : v01 = static_cast<Point>(*surrounding_nodes[n]) - node,
410 79180 : v02 = static_cast<Point>(*surrounding_nodes[n+1]) - node,
411 79180 : v03 = static_cast<Point>(*surrounding_nodes[n+2]) - node;
412 :
413 76980 : total_solid_angle += solid_angle(v01, v02, v03);
414 : }
415 :
416 31672 : return std::make_pair(n_surrounding, total_solid_angle);
417 3739 : };
418 :
419 : // We'll keep track of solid angles and node valences so we don't
420 : // waste time recomputing them when they haven't changed. We need
421 : // to be able to search efficiently for the smallest angles of each
422 : // valence, but also search efficiently for a particular node to
423 : // remove and reinsert it when its connectivity changes.
424 : //
425 : // Since C++11 multimap has guaranteed that pairs with matching keys
426 : // are kept in insertion order, so we can use Node * for values even
427 : // in parallel.
428 : typedef std::multimap<std::pair<int, Real>, Node*> node_map_type;
429 220 : node_map_type nodes_by_geometry;
430 220 : std::map<Node *, node_map_type::iterator> node_index;
431 :
432 38380 : for (auto node : surface.node_ptr_range())
433 30792 : node_index[node] =
434 65213 : nodes_by_geometry.emplace(geometry_at(*node), node);
435 :
436 : // In 3D, this will require nested loops: an outer loop to remove
437 : // each vertex, and an inner loop to remove multiple tetrahedra in
438 : // cases where the vertex has more than 3 neighboring triangles.
439 :
440 : // We'll be done when there are only three "unremoved" nodes left,
441 : // so they don't actually enclose any volume.
442 23094 : for (auto i : make_range(nodes_by_geometry.size()-3))
443 : {
444 550 : auto geometry_it = nodes_by_geometry.begin();
445 19245 : auto geometry_key = geometry_it->first;
446 19245 : auto [valence, angle] = geometry_key;
447 19245 : Node * node = geometry_it->second;
448 550 : libmesh_ignore(i);
449 :
450 : // If our lowest-valence nodes are all points of non-convexity,
451 : // skip to a higher valence.
452 19245 : while (angle > 2*pi-TOLERANCE)
453 : {
454 : geometry_it =
455 : nodes_by_geometry.upper_bound
456 0 : (std::make_pair(valence, Real(100)));
457 0 : libmesh_assert(geometry_it != nodes_by_geometry.end());
458 :
459 0 : std::tie(geometry_key, node) = *geometry_it;
460 0 : std::tie(valence, angle) = geometry_key;
461 : }
462 :
463 1100 : std::vector<Elem *> surrounding_elems;
464 : std::vector<Node *> surrounding_nodes =
465 19795 : surroundings_of(*node, &surrounding_elems);
466 :
467 1100 : const std::size_t n_surrounding = surrounding_nodes.size();
468 :
469 : // As we separate surrounding nodes from our center node, we'll
470 : // be marking them as nullptr; we still need to be able to find
471 : // predecessor and successor nodes in order afterward.
472 : auto find_valid_nodes_around =
473 72580 : [n_surrounding, & surrounding_nodes]
474 160560 : (unsigned int j)
475 : {
476 76980 : unsigned int jnext = (j+1)%n_surrounding;
477 79180 : while (!surrounding_nodes[jnext])
478 0 : jnext = (jnext+1)%n_surrounding;
479 :
480 76980 : unsigned int jprev = (j+n_surrounding-1)%n_surrounding;
481 79180 : while (!surrounding_nodes[jprev])
482 0 : jprev = (jprev+n_surrounding-1)%n_surrounding;
483 :
484 79180 : return std::make_pair(jprev, jnext);
485 19245 : };
486 :
487 : // We may have too many surrounding nodes to handle with
488 : // just one tet. In that case we'll keep a cache of the
489 : // element qualities that we'd get by making a tet with the
490 : // edge from the center node to each surrounding node, so we
491 : // can build the best tets first.
492 : //
493 : // In the case where we just have 3 nodes, we'll just pretend
494 : // they all have the same positive quality, so we can still
495 : // search this vector.
496 19795 : std::vector<Real> local_tet_quality(n_surrounding, 1);
497 :
498 : // From our center node with N surrounding nodes we can make N-2
499 : // tetrahedra. The first N-3 each replace two surface tets with
500 : // two new surface tets.
501 : //
502 : // My first idea was to greedily pick nodes with the smallest
503 : // local (solid) angles to get the best quality. This works in
504 : // 2D, but such a node can give a pancake tet in 3D.
505 : //
506 : // My second idea was to greedily pick nodes with the highest
507 : // prospective tet quality. This works for the first tets, but
508 : // can leave a pancake tet behind.
509 : //
510 : // My third idea is to try to fix the lowest quality tets first,
511 : // by picking cases where they have higher quality neighbors,
512 : // and creating those neighbors so as to change them.
513 :
514 : auto find_new_tet_nodes =
515 18145 : [& local_tet_quality, & find_valid_nodes_around]
516 34095 : ()
517 : {
518 550 : unsigned int jbest = 0;
519 19245 : auto [jminus, jplus] = find_valid_nodes_around(jbest);
520 19245 : Real qneighbest = std::min(local_tet_quality[jminus],
521 19795 : local_tet_quality[jplus]);
522 19245 : for (auto j : make_range(std::size_t(1),
523 58835 : local_tet_quality.size()))
524 : {
525 : // We don't want to build a bad tet
526 39590 : if (local_tet_quality[j] <= 0)
527 0 : continue;
528 :
529 38490 : std::tie(jminus, jplus) = find_valid_nodes_around(j);
530 38490 : Real qneighj = std::min(local_tet_quality[jminus],
531 39590 : local_tet_quality[jplus]);
532 :
533 : // We don't want to build a tet that can't fix a neighbor
534 : // if we can build one that can.
535 38490 : if (qneighbest <= 0 &&
536 : qneighj > 0)
537 0 : continue;
538 :
539 : // We want to try for the best possible fix.
540 39590 : if ((local_tet_quality[j] - qneighj) >
541 39590 : (local_tet_quality[jbest] - qneighj))
542 : {
543 0 : jbest = j;
544 0 : qneighbest = qneighj;
545 : }
546 : }
547 :
548 19795 : libmesh_error_msg_if
549 : (local_tet_quality[jbest] <= 0,
550 : "Cannot build non-singular non-inverted tet");
551 :
552 19245 : std::tie(jminus, jplus) = find_valid_nodes_around(jbest);
553 :
554 19795 : return std::make_tuple(jbest, jminus, jplus);
555 19245 : };
556 :
557 19245 : if (n_surrounding > 3)
558 : {
559 : // We'll be searching local_tet_quality even after tet
560 : // extraction disconnects us from some nodes; when we do we
561 : // don't want to get one.
562 0 : constexpr Real far_node = -1e6;
563 :
564 : // Vectors from the center node to each of its surrounding
565 : // nodes are helpful for calculating prospective tet
566 : // quality.
567 0 : std::vector<Point> v0s(n_surrounding);
568 0 : for (auto j : make_range(n_surrounding))
569 0 : v0s[j] = *(Point *)surrounding_nodes[j] - *node;
570 :
571 : // Find the tet quality we'd potentially get from each
572 : // possible choice of tet
573 : auto local_tet_quality_of =
574 0 : [& surrounding_nodes, & v0s, & find_valid_nodes_around]
575 0 : (unsigned int j)
576 : {
577 0 : auto [jminus, jplus] = find_valid_nodes_around(j);
578 :
579 : // Anything proportional to the ratio of volume to
580 : // total-edge-length-cubed should peak for perfect tets
581 : // but hit 0 for pancakes and slivers.
582 :
583 : const Real total_len =
584 0 : v0s[j].norm() + v0s[jminus].norm() + v0s[jplus].norm() +
585 0 : (*(Point *)surrounding_nodes[jplus] -
586 0 : *(Point *)surrounding_nodes[j]).norm() +
587 0 : (*(Point *)surrounding_nodes[j] -
588 0 : *(Point *)surrounding_nodes[jminus]).norm() +
589 0 : (*(Point *)surrounding_nodes[jminus] -
590 0 : *(Point *)surrounding_nodes[jplus]).norm();
591 :
592 : // Orientation here is tricky. Think of the triple
593 : // product as (v1 cross v2) dot v3, with right hand rule.
594 : const Real six_vol =
595 0 : triple_product(v0s[jminus], v0s[jplus], v0s[j]);
596 :
597 0 : return six_vol / (total_len * total_len * total_len);
598 0 : };
599 :
600 0 : for (auto j : make_range(n_surrounding))
601 0 : local_tet_quality[j] = local_tet_quality_of(j);
602 :
603 : // If we have N surrounding nodes, we can make N tets and
604 : // that'll bring us back to the 3-surrounding-node case to
605 : // finish.
606 0 : for (auto t : make_range(n_surrounding-3))
607 : {
608 0 : libmesh_ignore(t);
609 :
610 0 : auto [jbest, jminus, jplus] = find_new_tet_nodes();
611 :
612 : // Turn these four nodes into a tet
613 0 : Node * nbest = surrounding_nodes[jbest],
614 0 : * nminus = surrounding_nodes[jminus],
615 0 : * nplus = surrounding_nodes[jplus];
616 0 : this->add_tet(nminus->id(), nbest->id(), nplus->id(),
617 0 : node->id());
618 :
619 : // Replace the old two triangles around these nodes with the
620 : // proper two new triangles.
621 0 : Elem * oldtri1 = surrounding_elems[jminus],
622 0 : * oldtri2 = surrounding_elems[jbest],
623 0 : * newtri1 = surface.add_elem(Elem::build(TRI3)),
624 0 : * newtri2 = surface.add_elem(Elem::build(TRI3));
625 :
626 0 : const unsigned int c1 = oldtri1->get_node_index(node),
627 0 : c2 = oldtri2->get_node_index(node);
628 :
629 0 : newtri1->set_node(0, node);
630 0 : newtri1->set_node(1, nminus);
631 0 : newtri1->set_node(2, nplus);
632 :
633 0 : surrounding_elems[jminus] = newtri1;
634 :
635 0 : Elem * neigh10 = oldtri1->neighbor_ptr(c1);
636 0 : Elem * neigh12 = oldtri2->neighbor_ptr((c2+2)%3);
637 0 : newtri1->set_neighbor(0, neigh10);
638 0 : neigh10->set_neighbor(neigh10->which_neighbor_am_i(oldtri1), newtri1);
639 0 : newtri1->set_neighbor(1, newtri2);
640 0 : newtri1->set_neighbor(2, neigh12);
641 0 : neigh12->set_neighbor(neigh12->which_neighbor_am_i(oldtri2), newtri1);
642 :
643 0 : newtri2->set_node(0, nplus);
644 0 : newtri2->set_node(1, nminus);
645 0 : newtri2->set_node(2, nbest);
646 :
647 0 : Elem * neigh21 = oldtri1->neighbor_ptr((c1+1)%3);
648 0 : Elem * neigh22 = oldtri2->neighbor_ptr((c2+1)%3);
649 0 : newtri2->set_neighbor(0, newtri1);
650 0 : newtri2->set_neighbor(1, neigh21);
651 0 : neigh21->set_neighbor(neigh21->which_neighbor_am_i(oldtri1), newtri2);
652 0 : newtri2->set_neighbor(2, neigh22);
653 0 : neigh22->set_neighbor(neigh22->which_neighbor_am_i(oldtri2), newtri2);
654 :
655 0 : for (unsigned int p : make_range(3))
656 : {
657 0 : nodes_to_elem_map[oldtri1->node_id(p)].erase(oldtri1->id());
658 0 : nodes_to_elem_map[oldtri2->node_id(p)].erase(oldtri2->id());
659 0 : nodes_to_elem_map[newtri1->node_id(p)].insert(newtri1->id());
660 0 : nodes_to_elem_map[newtri2->node_id(p)].insert(newtri2->id());
661 : }
662 :
663 : // We've finished replacing the old triangles.
664 0 : surface.delete_elem(oldtri1);
665 0 : surface.delete_elem(oldtri2);
666 :
667 : // The solid angle for the far node should now stay
668 : // unchanged until we're out of this inner loop; let's
669 : // recalculate it here, and then we'll be done with it.
670 0 : Node * & nbestref = surrounding_nodes[jbest];
671 0 : nodes_by_geometry.erase(node_index[nbestref]);
672 0 : node_index[nbestref] =
673 0 : nodes_by_geometry.emplace(geometry_at(*nbestref), nbestref);
674 :
675 : // The far node is no longer sharing an edge with our center
676 : // node. Make sure we don't use it again with the center
677 : // node.
678 0 : local_tet_quality[jbest] = far_node;
679 0 : nbestref = nullptr;
680 :
681 : // The potential tet qualities using the side nodes have
682 : // changed now that they're directly connected to each
683 : // other.
684 0 : local_tet_quality[jminus] =
685 0 : local_tet_quality_of(jminus);
686 :
687 0 : local_tet_quality[jplus] =
688 0 : local_tet_quality_of(jplus);
689 : }
690 : }
691 :
692 : // Now we should have just 3 surrounding nodes, with which to
693 : // make one tetrahedron. Put them in a counterclockwise
694 : // (looking from outside) orientation, not the "best, clockwise,
695 : // counterclockwise" we get from the lambda.
696 19245 : auto [j2, j1, j3] = find_new_tet_nodes();
697 :
698 : // Turn these last four nodes into a tet
699 19245 : Node * n1 = surrounding_nodes[j1],
700 19245 : * n2 = surrounding_nodes[j2],
701 19245 : * n3 = surrounding_nodes[j3];
702 19245 : this->add_tet(n1->id(), n2->id(), n3->id(), node->id());
703 :
704 : // Replace the three surface triangles of that tet with the new
705 : // fourth triangle.
706 19245 : Elem * oldtri1 = surrounding_elems[j1],
707 19245 : * oldtri2 = surrounding_elems[j2],
708 19245 : * oldtri3 = surrounding_elems[j3],
709 19795 : * newtri = surface.add_elem(Elem::build(TRI3));
710 :
711 18695 : const unsigned int c1 = oldtri1->get_node_index(node),
712 18695 : c2 = oldtri2->get_node_index(node),
713 18695 : c3 = oldtri3->get_node_index(node);
714 :
715 19245 : newtri->set_node(0, n1);
716 19245 : newtri->set_node(1, n2);
717 19245 : newtri->set_node(2, n3);
718 19245 : Elem * neigh0 = oldtri1->neighbor_ptr((c1+1)%3);
719 1100 : newtri->set_neighbor(0, neigh0);
720 19245 : neigh0->set_neighbor(neigh0->which_neighbor_am_i(oldtri1), newtri);
721 19245 : Elem * neigh1 = oldtri2->neighbor_ptr((c2+1)%3);
722 1100 : newtri->set_neighbor(1, neigh1);
723 19245 : neigh1->set_neighbor(neigh1->which_neighbor_am_i(oldtri2), newtri);
724 19245 : Elem * neigh2 = oldtri3->neighbor_ptr((c3+1)%3);
725 1100 : newtri->set_neighbor(2, neigh2);
726 19245 : neigh2->set_neighbor(neigh2->which_neighbor_am_i(oldtri3), newtri);
727 :
728 76980 : for (unsigned int p : make_range(3))
729 : {
730 61035 : nodes_to_elem_map[oldtri1->node_id(p)].erase(oldtri1->id());
731 61035 : nodes_to_elem_map[oldtri2->node_id(p)].erase(oldtri2->id());
732 61035 : nodes_to_elem_map[oldtri3->node_id(p)].erase(oldtri3->id());
733 61035 : nodes_to_elem_map[newtri->node_id(p)].insert(newtri->id());
734 : }
735 :
736 : // We've finished replacing the old triangles.
737 19245 : surface.delete_elem(oldtri1);
738 19245 : surface.delete_elem(oldtri2);
739 19245 : surface.delete_elem(oldtri3);
740 :
741 : // We should have used up all our surrounding nodes now, and we
742 : // shouldn't have messed up our surface in the process, and our
743 : // center node should no longer be part of the surface.
744 : #ifdef DEBUG
745 550 : surrounding_nodes[j1] = nullptr;
746 550 : surrounding_nodes[j2] = nullptr;
747 550 : surrounding_nodes[j3] = nullptr;
748 :
749 2200 : for (auto ltq : surrounding_nodes)
750 1650 : libmesh_assert(!ltq);
751 :
752 550 : if (surface.n_elem() > 3)
753 440 : verify_surface();
754 :
755 3850 : for (const Elem * elem : surface.element_ptr_range())
756 13200 : for (auto p : make_range(3))
757 9900 : libmesh_assert_not_equal_to
758 : (elem->node_ptr(p), node);
759 : #endif
760 :
761 : // We've used up our center node, so it's not something we can
762 : // eliminate again.
763 18695 : nodes_by_geometry.erase(geometry_it);
764 : }
765 :
766 : // At this point our surface should just have two triangles left.
767 110 : libmesh_assert_equal_to(surface.n_elem(), 2);
768 11107 : }
769 :
770 :
771 19245 : void C0Polyhedron::add_tet(int n1,
772 : int n2,
773 : int n3,
774 : int n4)
775 : {
776 : #ifndef NDEBUG
777 550 : const auto nn = this->n_nodes();
778 550 : libmesh_assert_less(n1, nn);
779 550 : libmesh_assert_less(n2, nn);
780 550 : libmesh_assert_less(n3, nn);
781 550 : libmesh_assert_less(n4, nn);
782 :
783 550 : const Point v12 = this->point(n2) - this->point(n1);
784 550 : const Point v13 = this->point(n3) - this->point(n1);
785 550 : const Point v14 = this->point(n4) - this->point(n1);
786 550 : const Real six_vol = triple_product(v12, v13, v14);
787 550 : libmesh_assert_greater(six_vol, Real(0));
788 : #endif
789 :
790 19245 : this->_triangulation.push_back({n1, n2, n3, n4});
791 19245 : }
792 :
793 :
794 : } // namespace libMesh
|
