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