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