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 : #include "libmesh/cell_polyhedron.h"
19 :
20 : // Local includes
21 : #include "libmesh/face_polygon.h"
22 : #include "libmesh/enum_elem_quality.h"
23 : #include "libmesh/hashword.h"
24 :
25 : // C++ includes
26 : #include <algorithm>
27 : #include <array>
28 : #include <unordered_map>
29 : #include <unordered_set>
30 :
31 :
32 : namespace libMesh
33 : {
34 :
35 : // ------------------------------------------------------------
36 : // Polyhedron class static member initializations
37 : const int Polyhedron::num_children;
38 :
39 : // ------------------------------------------------------------
40 : // Polyhedron class member functions
41 :
42 :
43 9932 : Polyhedron::Polyhedron (const std::vector<std::shared_ptr<Polygon>> & sides,
44 9932 : Elem * p) :
45 568 : Cell(/* unused here */ 0, sides.size(), p, nullptr, nullptr),
46 568 : _elemlinks_data(sides.size()+2), // neighbors + parent + interior_parent
47 9364 : _nodelinks_data(0), // We'll have to resize *this* later too!
48 19580 : _sidelinks_data(sides.size())
49 : {
50 : // Set our sides, and while we're at it figure out our node maps
51 : // and side normal directions and edge lookup table, and count our
52 : // sides' nodes. If we have internal nodes too then the subclass
53 : // will append those afterward.
54 9932 : unsigned int nn = 0;
55 568 : std::unordered_map<Node *, unsigned int> local_node_number;
56 : std::unordered_set<std::pair<const Node *, const Node *>,
57 568 : libMesh::hash> edges_seen;
58 9932 : std::unique_ptr<const Elem> edge;
59 70046 : for (unsigned int s : index_range(sides))
60 : {
61 1724 : libmesh_assert(sides[s].get());
62 60114 : auto & side_tuple = _sidelinks_data[s];
63 3448 : std::get<0>(side_tuple) = sides[s];
64 :
65 3448 : Polygon & side = *sides[s]; // not const, for writeable nodes
66 301614 : for (auto n : make_range(side.n_nodes()))
67 : {
68 248436 : Node * node = side.node_ptr(n);
69 241500 : if (auto it = local_node_number.find(node);
70 6936 : it != local_node_number.end())
71 : {
72 161000 : std::get<2>(side_tuple).push_back(it->second);
73 : }
74 : else
75 : {
76 80500 : std::get<2>(side_tuple).push_back(nn);
77 80500 : local_node_number[node] = nn++;
78 80500 : _nodelinks_data.push_back(node);
79 : }
80 : }
81 :
82 301614 : for (unsigned int e : make_range(side.n_edges()))
83 : {
84 241500 : side.build_edge_ptr(edge, e);
85 241500 : auto edge_vertices = std::make_pair(edge->node_ptr(0), edge->node_ptr(1));
86 241500 : if (edge_vertices.first > edge_vertices.second)
87 3464 : std::swap(edge_vertices.first, edge_vertices.second);
88 :
89 13872 : if (!edges_seen.count(edge_vertices))
90 : {
91 3468 : edges_seen.insert(edge_vertices);
92 120750 : _edge_lookup.emplace_back(s, e);
93 : }
94 : }
95 : }
96 :
97 : // Plan room for an extra node so we don't invalidate this->_nodes when we add to
98 : // nodelinks_data in the derived class
99 10216 : _nodelinks_data.reserve(_nodelinks_data.size() + 1);
100 :
101 : // Do the manual initialization that Elem::Elem and Cell::Cell
102 : // couldn't, now that we've resized both our vectors. No need to
103 : // manually set nullptr, though, since std::vector does that.
104 9932 : this->_elemlinks = _elemlinks_data.data();
105 9932 : this->_nodes = _nodelinks_data.data();
106 9932 : this->_elemlinks[0] = p;
107 :
108 284 : libmesh_assert_equal_to(nn, this->n_nodes());
109 :
110 : // Figure out the orientation of our sides, now that we've got our
111 : // nodes organized enough to find our center. The algorithm below
112 : // only works for convex polyhedra, but that's all we're
113 : // supporting for now.
114 284 : Point center;
115 90432 : for (auto n : make_range(nn))
116 2312 : center.add (this->point(n));
117 9932 : center /= static_cast<Real>(nn);
118 :
119 70046 : for (unsigned int s : index_range(sides))
120 : {
121 3448 : const Polygon & side = *sides[s];
122 61838 : const Point x_i = side.point(0);
123 : const Point n_i =
124 58390 : (side.point(1) - side.point(0)).cross
125 61838 : (side.point(0) - side.point(side.n_sides()-1)).unit();
126 :
127 3448 : bool & inward_normal = std::get<1>(_sidelinks_data[s]);
128 60114 : inward_normal = (n_i * (center - x_i) > TOLERANCE);
129 : }
130 :
131 : // We're betting a lot on "our polyhedra are all convex", so let's
132 : // check that if we have time.
133 : #ifdef DEBUG
134 2008 : for (unsigned int s : index_range(sides))
135 : {
136 1724 : const Polygon & side = *sides[s];
137 1724 : const Point x_i = side.point(0);
138 1724 : const bool inward_normal = std::get<1>(this->_sidelinks_data[s]);
139 :
140 : const Point n_i =
141 1724 : (side.point(1) - side.point(0)).cross
142 1724 : (side.point(0) - side.point(side.n_sides()-1)).unit() *
143 3448 : (inward_normal ? -1 : 1);
144 :
145 15900 : for (const Point & node : this->node_ref_range())
146 : {
147 14176 : const Point d_n = node - x_i;
148 14176 : if (d_n * n_i > TOLERANCE * d_n.norm())
149 0 : libmesh_not_implemented_msg
150 : ("Cannot create a non-convex polyhedron");
151 : }
152 : }
153 : #endif
154 :
155 : // Is this likely to ever be used? We may do refinement with
156 : // polyhedra but it's probably not going to have a hierarchy...
157 9932 : if (p)
158 : {
159 0 : this->subdomain_id() = p->subdomain_id();
160 0 : this->processor_id() = p->processor_id();
161 0 : _map_type = p->mapping_type();
162 0 : _map_data = p->mapping_data();
163 :
164 : #ifdef LIBMESH_ENABLE_AMR
165 0 : this->set_p_level(p->p_level());
166 : #endif
167 : }
168 :
169 : // Make sure the interior parent isn't undefined
170 9932 : this->set_interior_parent(nullptr);
171 19296 : }
172 :
173 :
174 :
175 888573036 : Point Polyhedron::master_point (const unsigned int i) const
176 : {
177 961647264 : return this->point(i);
178 : }
179 :
180 :
181 :
182 0 : bool Polyhedron::convex()
183 : {
184 0 : for (unsigned int s : make_range(this->n_sides()))
185 : {
186 0 : const Polygon & side = *std::get<0>(this->_sidelinks_data[s]);
187 0 : const Point x_i = side.point(0);
188 0 : const bool inward_normal = std::get<1>(this->_sidelinks_data[s]);
189 :
190 : const Point n_i =
191 0 : (side.point(1) - side.point(0)).cross
192 0 : (side.point(0) - side.point(side.n_sides()-1)).unit() *
193 0 : (inward_normal ? -1 : 1);
194 :
195 0 : for (const Point & node : this->node_ref_range())
196 : {
197 0 : const Point d_n = node - x_i;
198 0 : if (d_n * n_i > TOLERANCE * d_n.norm())
199 0 : return false;
200 : }
201 : }
202 0 : return true;
203 : }
204 :
205 :
206 :
207 442986 : bool Polyhedron::on_reference_element(const Point & p,
208 : const Real eps) const
209 : {
210 57898 : const unsigned int ns = this->n_sides();
211 :
212 : // Check that the point is on the same side of all the faces by
213 : // testing whether:
214 : //
215 : // n_i.(p - x_i) <= 0
216 : //
217 : // for each i, where:
218 : // n_i is the outward normal of face i,
219 : // x_i is a point on face i.
220 :
221 2349102 : for (auto i : make_range(ns))
222 : {
223 2051436 : const Polygon & face = *std::get<0>(this->_sidelinks_data[i]);
224 2051436 : const bool inward_normal = std::get<1>(this->_sidelinks_data[i]);
225 :
226 2051436 : const Point x_i = face.point(0);
227 :
228 : const Point n_i =
229 1891928 : (face.point(1) - face.point(0)).cross
230 2210944 : (face.point(0) - face.point(face.n_sides()-1)).unit() *
231 2348284 : (inward_normal ? -1 : 1);
232 :
233 : // This only works for polyhedra with flat sides.
234 : #ifdef DEBUG
235 1187392 : for (auto j : make_range(face.n_sides()-1))
236 : {
237 890544 : const Point x_j = face.point(j+1);
238 890544 : const Point d_j = x_j - x_i;
239 890544 : if (std::abs(d_j * n_i) > eps * d_j.norm())
240 0 : libmesh_not_implemented_msg
241 : ("Polyhedra with non-flat sides are not fully supported.");
242 : }
243 : #endif
244 :
245 2051436 : if (n_i * (p - x_i) > eps)
246 133210 : return false;
247 : }
248 :
249 68686 : return true;
250 : }
251 :
252 :
253 :
254 0 : dof_id_type Polyhedron::key (const unsigned int s) const
255 : {
256 0 : libmesh_assert_less (s, this->n_sides());
257 :
258 0 : const Polygon & face = *std::get<0>(this->_sidelinks_data[s]);
259 :
260 0 : return face.key();
261 : }
262 :
263 :
264 :
265 61348 : dof_id_type Polyhedron::low_order_key (const unsigned int s) const
266 : {
267 1672 : libmesh_assert_less (s, this->n_sides());
268 :
269 61348 : const Polygon & face = *std::get<0>(this->_sidelinks_data[s]);
270 :
271 1672 : const unsigned int nv = face.n_vertices();
272 63020 : std::vector<dof_id_type> vertex_ids(nv);
273 307744 : for (unsigned int v : make_range(nv))
274 259836 : vertex_ids[v] = face.node_id(v);
275 :
276 : // Sort the vertex ids before hashing so the key is independent of the
277 : // order in which the face's vertices are stored, just like
278 : // Elem::compute_key() does for every other element type. This is
279 : // essential for find_neighbors(): two polyhedra sharing a face
280 : // generally wind that face in opposite orders (each orienting it
281 : // outward from its own cell), so without sorting the two sides would
282 : // hash to different keys, never get compared, and never be linked as
283 : // neighbors.
284 61348 : std::sort(vertex_ids.begin(), vertex_ids.end());
285 :
286 63020 : return Utility::hashword(vertex_ids);
287 : }
288 :
289 :
290 :
291 0 : unsigned int Polyhedron::local_side_node(unsigned int side,
292 : unsigned int side_node) const
293 : {
294 0 : libmesh_assert_less (side, this->n_sides());
295 :
296 : const std::vector<unsigned int> & node_map =
297 0 : std::get<2>(this->_sidelinks_data[side]);
298 0 : libmesh_assert_less (side_node, node_map.size());
299 :
300 0 : return node_map[side_node];
301 : }
302 :
303 :
304 :
305 4608 : unsigned int Polyhedron::local_edge_node(unsigned int edge,
306 : unsigned int edge_node) const
307 : {
308 384 : libmesh_assert_less (edge, this->n_edges());
309 384 : libmesh_assert_less (edge, _edge_lookup.size());
310 :
311 4608 : auto [side, edge_of_side] = _edge_lookup[edge];
312 :
313 4608 : const Polygon & face = *std::get<0>(this->_sidelinks_data[side]);
314 :
315 : const std::vector<unsigned int> & node_map =
316 384 : std::get<2>(this->_sidelinks_data[side]);
317 :
318 4608 : return node_map[face.local_edge_node(edge_of_side, edge_node)];
319 : }
320 :
321 :
322 :
323 0 : dof_id_type Polyhedron::key () const
324 : {
325 0 : std::vector<dof_id_type> node_ids;
326 0 : for (const auto & n : this->node_ref_range())
327 0 : node_ids.push_back(n.id());
328 :
329 0 : return Utility::hashword(node_ids);
330 : }
331 :
332 :
333 :
334 37898 : std::unique_ptr<Elem> Polyhedron::side_ptr (const unsigned int i)
335 : {
336 37898 : return const_cast<const Polyhedron *>(this)->side_ptr(i);
337 : }
338 :
339 :
340 :
341 38750 : std::unique_ptr<Elem> Polyhedron::side_ptr (const unsigned int i) const
342 : {
343 1506 : libmesh_assert_less (i, this->n_sides());
344 :
345 38750 : Polygon & face = *std::get<0>(this->_sidelinks_data[i]);
346 38750 : std::unique_ptr<Elem> face_copy = face.disconnected_clone();
347 193424 : for (auto n : face.node_index_range())
348 160674 : face_copy->set_node(n, face.node_ptr(n));
349 :
350 38750 : return face_copy;
351 0 : }
352 :
353 :
354 :
355 32344 : void Polyhedron::side_ptr (std::unique_ptr<Elem> & side,
356 : const unsigned int i)
357 : {
358 1242 : libmesh_assert_less (i, this->n_sides());
359 :
360 : // Polyhedra are irregular enough that we're not even going to try
361 : // and bother optimizing heap access here.
362 62204 : side = this->side_ptr(i);
363 32344 : }
364 :
365 :
366 :
367 5554 : std::unique_ptr<Elem> Polyhedron::build_side_ptr (const unsigned int i)
368 : {
369 5554 : auto returnval = this->side_ptr(i);
370 5554 : returnval->set_interior_parent(this);
371 5314 : returnval->inherit_data_from(*this);
372 5554 : return returnval;
373 0 : }
374 :
375 :
376 :
377 22554 : void Polyhedron::build_side_ptr (std::unique_ptr<Elem> & side,
378 : const unsigned int i)
379 : {
380 22554 : this->side_ptr(side, i);
381 22554 : side->set_interior_parent(this);
382 21556 : side->inherit_data_from(*this);
383 22554 : }
384 :
385 :
386 :
387 0 : std::unique_ptr<Elem> Polyhedron::build_edge_ptr (const unsigned int i)
388 : {
389 0 : auto [s, se] = _edge_lookup[i];
390 0 : Polygon & face = *std::get<0>(_sidelinks_data[s]);
391 0 : return face.build_edge_ptr(se);
392 : }
393 :
394 :
395 :
396 0 : void Polyhedron::build_edge_ptr (std::unique_ptr<Elem> & elem,
397 : const unsigned int i)
398 : {
399 0 : auto [s, se] = _edge_lookup[i];
400 0 : Polygon & face = *std::get<0>(_sidelinks_data[s]);
401 0 : face.build_edge_ptr(elem, se);
402 0 : }
403 :
404 :
405 :
406 0 : bool Polyhedron::is_child_on_side(const unsigned int /*c*/,
407 : const unsigned int /*s*/) const
408 : {
409 0 : libmesh_not_implemented();
410 : return false;
411 : }
412 :
413 :
414 :
415 0 : unsigned int Polyhedron::opposite_side(const unsigned int /* side_in */) const
416 : {
417 : // This is too ambiguous in general.
418 0 : libmesh_not_implemented();
419 : return libMesh::invalid_uint;
420 : }
421 :
422 :
423 :
424 0 : unsigned int Polyhedron::opposite_node(const unsigned int /* n */,
425 : const unsigned int /* s */) const
426 : {
427 : // This is too ambiguous in general.
428 0 : libmesh_not_implemented();
429 : return libMesh::invalid_uint;
430 : }
431 :
432 :
433 :
434 205 : bool Polyhedron::is_flipped() const
435 : {
436 205 : if (this->_triangulation.empty())
437 0 : return false;
438 :
439 10 : auto & tet = this->_triangulation[0];
440 :
441 215 : const Point v01 = this->point(tet[1]) - this->point(tet[0]);
442 205 : const Point v02 = this->point(tet[2]) - this->point(tet[0]);
443 205 : const Point v03 = this->point(tet[3]) - this->point(tet[0]);
444 :
445 205 : return (triple_product(v01, v02, v03) < 0);
446 : }
447 :
448 :
449 : std::vector<unsigned int>
450 420 : Polyhedron::edges_adjacent_to_node(const unsigned int n) const
451 : {
452 35 : libmesh_assert_less(n, this->n_nodes());
453 :
454 : // For mid-edge or mid-face nodes, the subclass had better have
455 : // overridden this.
456 35 : libmesh_assert_less(n, this->n_vertices());
457 :
458 35 : libmesh_assert_equal_to(this->n_edges(), _edge_lookup.size());
459 :
460 35 : std::vector<unsigned int> adjacent_edges;
461 :
462 5460 : for (const auto edge : index_range(_edge_lookup))
463 : {
464 5040 : const auto [side, side_edge] = _edge_lookup[edge];
465 5040 : const Polygon & face = *std::get<0>(_sidelinks_data[side]);
466 : const std::vector<unsigned int> & node_map =
467 420 : std::get<2>(_sidelinks_data[side]);
468 :
469 420 : const unsigned int fnv = face.n_vertices();
470 420 : libmesh_assert_equal_to(fnv, face.n_edges());
471 420 : libmesh_assert_equal_to(fnv, face.n_sides());
472 420 : libmesh_assert_less_equal(fnv, node_map.size());
473 420 : libmesh_assert_less(side_edge, fnv);
474 :
475 5040 : const unsigned int vn = node_map[side_edge];
476 5040 : const unsigned int vnp = node_map[(side_edge + 1) % fnv];
477 420 : libmesh_assert_less(vn, this->n_vertices());
478 420 : libmesh_assert_less(vnp, this->n_vertices());
479 5040 : if (vn == n || vnp == n)
480 1260 : adjacent_edges.push_back(edge);
481 : }
482 :
483 420 : return adjacent_edges;
484 : }
485 :
486 :
487 60 : std::pair<Real, Real> Polyhedron::qual_bounds (const ElemQuality q) const
488 : {
489 5 : std::pair<Real, Real> bounds;
490 :
491 60 : switch (q)
492 : {
493 1 : case EDGE_LENGTH_RATIO:
494 1 : bounds.first = 1.;
495 1 : bounds.second = 4.;
496 1 : break;
497 :
498 12 : case MIN_ANGLE:
499 1 : bounds.first = 30.;
500 1 : bounds.second = 180.;
501 12 : break;
502 :
503 12 : case MAX_ANGLE:
504 1 : bounds.first = 60.;
505 1 : bounds.second = 180.;
506 12 : break;
507 :
508 : // The generic dihedral metrics measure the unoriented angle
509 : // between adjacent face planes, so their values are in [0, 90].
510 : // Use the same recommended lower bounds as MIN,MAX_ANGLE.
511 12 : case MIN_DIHEDRAL_ANGLE:
512 1 : bounds.first = 30.;
513 1 : bounds.second = 90.;
514 12 : break;
515 :
516 12 : case MAX_DIHEDRAL_ANGLE:
517 1 : bounds.first = 60.;
518 1 : bounds.second = 90.;
519 12 : break;
520 :
521 0 : case JACOBIAN:
522 : case SCALED_JACOBIAN:
523 0 : bounds.first = 0.5;
524 0 : bounds.second = 1.;
525 0 : break;
526 :
527 0 : default:
528 0 : libMesh::out << "Warning: Invalid quality measure chosen." << std::endl;
529 0 : bounds.first = -1;
530 0 : bounds.second = -1;
531 : }
532 :
533 60 : return bounds;
534 : }
535 :
536 :
537 :
538 : std::vector<std::shared_ptr<Polygon>>
539 15 : Polyhedron::side_clones() const
540 : {
541 2 : const auto ns = this->n_sides();
542 :
543 2 : libmesh_assert_equal_to(ns, _sidelinks_data.size());
544 :
545 15 : std::vector<std::shared_ptr<Polygon>> cloned_sides(ns);
546 :
547 105 : for (auto i : make_range(ns))
548 : {
549 90 : const Polygon & face = *std::get<0>(this->_sidelinks_data[i]);
550 :
551 90 : Elem * clone = face.disconnected_clone().release();
552 12 : Polygon * polygon_clone = cast_ptr<Polygon *>(clone);
553 90 : cloned_sides[i] = std::shared_ptr<Polygon>(polygon_clone);
554 :
555 : // We can't actually use a *disconnected* clone to reconstruct
556 : // links between sides, so we'll temporarily give the clone our
557 : // own nodes; user code that typically replaces the usual
558 : // nullptr with permanent nodes will then instead place our
559 : // nodes with permanent nodes.
560 528 : for (auto n : make_range(face.n_nodes()))
561 96 : cloned_sides[i]->set_node
562 408 : (n, const_cast<Node *>(face.node_ptr(n)));
563 : }
564 :
565 17 : return cloned_sides;
566 0 : }
567 :
568 :
569 :
570 3384 : bool Polyhedron::side_has_edge_nodes(unsigned int s,
571 : unsigned int min_node,
572 : unsigned int max_node) const
573 : {
574 3384 : const Polygon & face = *std::get<0>(_sidelinks_data[s]);
575 : const std::vector<unsigned int> & node_map =
576 282 : std::get<2>(this->_sidelinks_data[s]);
577 :
578 14484 : for (unsigned int e : make_range(face.n_sides()))
579 : {
580 : std::vector<unsigned int> nodes_on_edge =
581 12108 : face.nodes_on_side(e);
582 1009 : libmesh_assert_equal_to(nodes_on_edge.size(), 2);
583 12108 : nodes_on_edge[0] = node_map[nodes_on_edge[0]];
584 12108 : nodes_on_edge[1] = node_map[nodes_on_edge[1]];
585 12217 : if ((nodes_on_edge[0] == min_node) &&
586 109 : (nodes_on_edge[1] == max_node))
587 70 : return true;
588 11811 : if ((nodes_on_edge[1] == min_node) &&
589 123 : (nodes_on_edge[0] == max_node))
590 98 : return true;
591 : }
592 :
593 396 : return false;
594 : }
595 :
596 :
597 :
598 : std::vector<unsigned int>
599 864 : Polyhedron::sides_on_edge(const unsigned int e) const
600 : {
601 864 : std::vector<unsigned int> returnval(2);
602 864 : auto [s1, s1e] = _edge_lookup[e];
603 864 : returnval[0] = s1;
604 :
605 864 : const Polygon & face1 = *std::get<0>(_sidelinks_data[s1]);
606 : const std::vector<unsigned int> & node_map =
607 72 : std::get<2>(this->_sidelinks_data[s1]);
608 :
609 : std::vector<unsigned int> nodes_on_edge =
610 864 : face1.nodes_on_side(s1e);
611 72 : libmesh_assert_equal_to(nodes_on_edge.size(), 2);
612 864 : nodes_on_edge[0] = node_map[nodes_on_edge[0]];
613 864 : nodes_on_edge[1] = node_map[nodes_on_edge[1]];
614 :
615 864 : if (nodes_on_edge[0] > nodes_on_edge[1])
616 18 : std::swap(nodes_on_edge[0], nodes_on_edge[1]);
617 :
618 4104 : for (unsigned int s2 : make_range(this->n_sides()))
619 : {
620 4104 : if (s2 == s1)
621 792 : continue;
622 :
623 3240 : if (this->side_has_edge_nodes(s2, nodes_on_edge[0],
624 540 : nodes_on_edge[1]))
625 : {
626 864 : returnval[1] = s2;
627 936 : return returnval;
628 : }
629 : }
630 :
631 0 : libmesh_error();
632 :
633 : return returnval;
634 : }
635 :
636 :
637 :
638 288 : bool Polyhedron::is_edge_on_side(const unsigned int e,
639 : const unsigned int s) const
640 : {
641 288 : auto [s1, s1e] = _edge_lookup[e];
642 :
643 : // Did we get lucky with our cache?
644 288 : if (s1 == s)
645 12 : return true;
646 :
647 144 : const Polygon & face1 = *std::get<0>(_sidelinks_data[s1]);
648 : const std::vector<unsigned int> & node_map =
649 12 : std::get<2>(this->_sidelinks_data[s1]);
650 : std::vector<unsigned int> nodes_on_edge1 =
651 156 : face1.nodes_on_side(s1e);
652 12 : libmesh_assert_equal_to(nodes_on_edge1.size(), 2);
653 :
654 144 : nodes_on_edge1[0] = node_map[nodes_on_edge1[0]];
655 144 : nodes_on_edge1[1] = node_map[nodes_on_edge1[1]];
656 144 : if (nodes_on_edge1[0] > nodes_on_edge1[1])
657 3 : std::swap(nodes_on_edge1[0], nodes_on_edge1[1]);
658 :
659 144 : return this->side_has_edge_nodes(s,
660 12 : nodes_on_edge1[0],
661 24 : nodes_on_edge1[1]);
662 : }
663 :
664 :
665 :
666 222143235 : std::array<Point, 4> Polyhedron::master_subelement (unsigned int i) const
667 : {
668 19108205 : libmesh_assert_less(i, this->_triangulation.size());
669 :
670 222143235 : const auto & tet = this->_triangulation[i];
671 :
672 222143235 : return { this->master_point(tet[0]),
673 222143235 : this->master_point(tet[1]),
674 222143235 : this->master_point(tet[2]),
675 222143235 : this->master_point(tet[3]) };
676 : }
677 :
678 :
679 :
680 : std::tuple<unsigned int, Real, Real, Real>
681 59113506 : Polyhedron::subelement_coordinates (const Point & p, Real tol) const
682 : {
683 5085353 : std::tuple<unsigned int, Real, Real, Real> returnval =
684 : {libMesh::invalid_uint, -1, -1, -1};
685 :
686 5085353 : Real best_bad_coord = -1;
687 :
688 205708960 : for (auto s : make_range(this->n_subelements()))
689 : {
690 : const std::array<Point, 4> subtet =
691 200808184 : this->master_subelement(s);
692 :
693 : // Find barycentric coordinates in subelem
694 17336257 : const Point v0 = p - subtet[0];
695 : // const Point v1 = p - subtet[1];
696 :
697 17336257 : const Point v01 = subtet[1] - subtet[0];
698 17336257 : const Point v02 = subtet[2] - subtet[0];
699 17336257 : const Point v03 = subtet[3] - subtet[0];
700 :
701 : // const Point v12 = subtet[2] - subtet[1];
702 : // const Point v13 = subtet[3] - subtet[1];
703 :
704 : // const Real tp0 = triple_product(v1, v13, v12);
705 184311577 : const Real tp1 = triple_product(v0, v02, v03);
706 184311577 : const Real tp2 = triple_product(v0, v03, v01);
707 184311577 : const Real tp3 = triple_product(v0, v01, v02);
708 :
709 184311577 : const Real six_vol = triple_product(v01, v02, v03);
710 :
711 200808184 : const Real xi = tp1 / six_vol;
712 200808184 : const Real eta = tp2 / six_vol;
713 200808184 : const Real zeta = tp3 / six_vol;
714 :
715 200808184 : if (xi>=0 && eta>=0 && zeta>=0 && xi+eta+zeta<=1)
716 54212730 : return { s, xi, eta, zeta };
717 :
718 : const Real my_best_bad_coord =
719 146595454 : std::min(std::min(std::min(xi, eta), zeta), 1-xi-eta-zeta);
720 :
721 146595454 : if (my_best_bad_coord > best_bad_coord)
722 : {
723 5230582 : best_bad_coord = my_best_bad_coord;
724 5230582 : returnval = { s, xi, eta, zeta };
725 : }
726 : }
727 :
728 4900776 : if (best_bad_coord > -tol)
729 408398 : return returnval;
730 :
731 0 : return {libMesh::invalid_uint, -1, -1, -1};
732 : }
733 :
734 :
735 : } // namespace libMesh
|
