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 910 : Polyhedron::Polyhedron (const std::vector<std::shared_ptr<Polygon>> & sides,
43 910 : Elem * p) :
44 520 : Cell(/* unused here */ 0, sides.size(), p, nullptr, nullptr),
45 520 : _elemlinks_data(sides.size()+2), // neighbors + parent + interior_parent
46 390 : _nodelinks_data(0), // We'll have to resize *this* later too!
47 1560 : _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 910 : unsigned int nn = 0;
54 520 : std::unordered_map<Node *, unsigned int> local_node_number;
55 : std::unordered_set<std::pair<const Node *, const Node *>,
56 520 : libMesh::hash> edges_seen;
57 910 : std::unique_ptr<const Elem> edge;
58 6440 : for (unsigned int s : index_range(sides))
59 : {
60 1580 : libmesh_assert(sides[s].get());
61 5530 : auto & side_tuple = _sidelinks_data[s];
62 3160 : std::get<0>(side_tuple) = sides[s];
63 :
64 3160 : Polygon & side = *sides[s]; // not const, for writeable nodes
65 27790 : for (auto n : make_range(side.n_nodes()))
66 : {
67 28620 : Node * node = side.node_ptr(n);
68 22260 : if (auto it = local_node_number.find(node);
69 6360 : it != local_node_number.end())
70 : {
71 14840 : std::get<2>(side_tuple).push_back(it->second);
72 : }
73 : else
74 : {
75 7420 : std::get<2>(side_tuple).push_back(nn);
76 7420 : local_node_number[node] = nn++;
77 7420 : _nodelinks_data.push_back(node);
78 : }
79 : }
80 :
81 27790 : for (unsigned int e : make_range(side.n_edges()))
82 : {
83 22260 : side.build_edge_ptr(edge, e);
84 22260 : auto edge_vertices = std::make_pair(edge->node_ptr(0), edge->node_ptr(1));
85 22260 : if (edge_vertices.first > edge_vertices.second)
86 3184 : std::swap(edge_vertices.first, edge_vertices.second);
87 :
88 12720 : if (!edges_seen.count(edge_vertices))
89 : {
90 3180 : edges_seen.insert(edge_vertices);
91 11130 : _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 1170 : _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 910 : this->_elemlinks = _elemlinks_data.data();
104 910 : this->_nodes = _nodelinks_data.data();
105 910 : this->_elemlinks[0] = p;
106 :
107 260 : 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 260 : Point center;
114 8330 : for (auto n : make_range(nn))
115 2120 : center.add (this->point(n));
116 910 : center /= static_cast<Real>(nn);
117 :
118 6440 : for (unsigned int s : index_range(sides))
119 : {
120 3160 : const Polygon & side = *sides[s];
121 7110 : const Point x_i = side.point(0);
122 : const Point n_i =
123 3950 : (side.point(1) - side.point(0)).cross
124 7110 : (side.point(0) - side.point(side.n_sides()-1)).unit();
125 :
126 3160 : bool & inward_normal = std::get<1>(_sidelinks_data[s]);
127 5530 : 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 1840 : for (unsigned int s : index_range(sides))
134 : {
135 1580 : const Polygon & side = *sides[s];
136 1580 : const Point x_i = side.point(0);
137 1580 : const bool inward_normal = std::get<1>(this->_sidelinks_data[s]);
138 :
139 : const Point n_i =
140 1580 : (side.point(1) - side.point(0)).cross
141 1580 : (side.point(0) - side.point(side.n_sides()-1)).unit() *
142 3160 : (inward_normal ? -1 : 1);
143 :
144 14540 : for (const Point & node : this->node_ref_range())
145 : {
146 12960 : const Point d_n = node - x_i;
147 12960 : 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 910 : 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 910 : this->set_interior_parent(nullptr);
170 1300 : }
171 :
172 :
173 :
174 294136668 : Point Polyhedron::master_point (const unsigned int i) const
175 : {
176 366875460 : 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 162922 : 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 889854 : for (auto i : make_range(ns))
221 : {
222 775372 : const Polygon & face = *std::get<0>(this->_sidelinks_data[i]);
223 775372 : const bool inward_normal = std::get<1>(this->_sidelinks_data[i]);
224 :
225 775372 : const Point x_i = face.point(0);
226 :
227 : const Point n_i =
228 615864 : (face.point(1) - face.point(0)).cross
229 934880 : (face.point(0) - face.point(face.n_sides()-1)).unit() *
230 1072220 : (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 775372 : if (n_i * (p - x_i) > eps)
245 36330 : 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 5404 : dof_id_type Polyhedron::low_order_key (const unsigned int s) const
265 : {
266 1544 : libmesh_assert_less (s, this->n_sides());
267 :
268 5404 : const Polygon & face = *std::get<0>(this->_sidelinks_data[s]);
269 :
270 1544 : const unsigned int nv = face.n_vertices();
271 6948 : std::vector<dof_id_type> vertex_ids(nv);
272 27160 : for (unsigned int v : make_range(nv))
273 34188 : vertex_ids[v] = face.node_id(v);
274 :
275 6948 : 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 2208 : unsigned int Polyhedron::local_edge_node(unsigned int edge,
295 : unsigned int edge_node) const
296 : {
297 552 : libmesh_assert_less (edge, this->n_edges());
298 552 : libmesh_assert_less (edge, _edge_lookup.size());
299 :
300 2208 : auto [side, edge_of_side] = _edge_lookup[edge];
301 :
302 2208 : const Polygon & face = *std::get<0>(this->_sidelinks_data[side]);
303 :
304 : const std::vector<unsigned int> & node_map =
305 552 : std::get<2>(this->_sidelinks_data[side]);
306 :
307 2208 : 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 4156 : std::unique_ptr<Elem> Polyhedron::side_ptr (const unsigned int i)
324 : {
325 4156 : return const_cast<const Polyhedron *>(this)->side_ptr(i);
326 : }
327 :
328 :
329 :
330 4240 : std::unique_ptr<Elem> Polyhedron::side_ptr (const unsigned int i) const
331 : {
332 1322 : libmesh_assert_less (i, this->n_sides());
333 :
334 4240 : Polygon & face = *std::get<0>(this->_sidelinks_data[i]);
335 4240 : std::unique_ptr<Elem> face_copy = face.disconnected_clone();
336 21256 : for (auto n : face.node_index_range())
337 22320 : face_copy->set_node(n, face.node_ptr(n));
338 :
339 4240 : return face_copy;
340 0 : }
341 :
342 :
343 :
344 3350 : void Polyhedron::side_ptr (std::unique_ptr<Elem> & side,
345 : const unsigned int i)
346 : {
347 1051 : 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 4598 : side = this->side_ptr(i);
352 3350 : }
353 :
354 :
355 :
356 806 : std::unique_ptr<Elem> Polyhedron::build_side_ptr (const unsigned int i)
357 : {
358 806 : auto returnval = this->side_ptr(i);
359 806 : returnval->set_interior_parent(this);
360 559 : returnval->inherit_data_from(*this);
361 806 : return returnval;
362 0 : }
363 :
364 :
365 :
366 3350 : void Polyhedron::build_side_ptr (std::unique_ptr<Elem> & side,
367 : const unsigned int i)
368 : {
369 3350 : this->side_ptr(side, i);
370 3350 : side->set_interior_parent(this);
371 2299 : side->inherit_data_from(*this);
372 3350 : }
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 37 : bool Polyhedron::is_flipped() const
424 : {
425 37 : if (this->_triangulation.empty())
426 0 : return false;
427 :
428 10 : auto & tet = this->_triangulation[0];
429 :
430 47 : const Point v01 = this->point(tet[1]) - this->point(tet[0]);
431 37 : const Point v02 = this->point(tet[2]) - this->point(tet[0]);
432 37 : const Point v03 = this->point(tet[3]) - this->point(tet[0]);
433 :
434 37 : return (triple_product(v01, v02, v03) < 0);
435 : }
436 :
437 :
438 : std::vector<unsigned int>
439 252 : Polyhedron::edges_adjacent_to_node(const unsigned int n) const
440 : {
441 63 : 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 63 : libmesh_assert_less(n, this->n_vertices());
446 :
447 63 : libmesh_assert_equal_to(this->n_edges(), _edge_lookup.size());
448 :
449 63 : std::vector<unsigned int> adjacent_edges;
450 :
451 3276 : for (const auto edge : index_range(_edge_lookup))
452 : {
453 3024 : const auto [side, side_edge] = _edge_lookup[edge];
454 3024 : const Polygon & face = *std::get<0>(_sidelinks_data[side]);
455 : const std::vector<unsigned int> & node_map =
456 756 : std::get<2>(_sidelinks_data[side]);
457 :
458 756 : const unsigned int fnv = face.n_vertices();
459 756 : libmesh_assert_equal_to(fnv, face.n_edges());
460 756 : libmesh_assert_equal_to(fnv, face.n_sides());
461 756 : libmesh_assert_less_equal(fnv, node_map.size());
462 756 : libmesh_assert_less(side_edge, fnv);
463 :
464 3024 : const unsigned int vn = node_map[side_edge];
465 3024 : const unsigned int vnp = node_map[(side_edge + 1) % fnv];
466 756 : libmesh_assert_less(vn, this->n_vertices());
467 756 : libmesh_assert_less(vnp, this->n_vertices());
468 3024 : if (vn == n || vnp == n)
469 756 : adjacent_edges.push_back(edge);
470 : }
471 :
472 252 : return adjacent_edges;
473 : }
474 :
475 :
476 20 : std::pair<Real, Real> Polyhedron::qual_bounds (const ElemQuality q) const
477 : {
478 5 : std::pair<Real, Real> bounds;
479 :
480 20 : switch (q)
481 : {
482 1 : case EDGE_LENGTH_RATIO:
483 1 : bounds.first = 1.;
484 1 : bounds.second = 4.;
485 1 : break;
486 :
487 4 : case MIN_ANGLE:
488 1 : bounds.first = 30.;
489 1 : bounds.second = 180.;
490 4 : break;
491 :
492 4 : case MAX_ANGLE:
493 1 : bounds.first = 60.;
494 1 : bounds.second = 180.;
495 4 : break;
496 :
497 8 : case JACOBIAN:
498 : case SCALED_JACOBIAN:
499 2 : bounds.first = 0.5;
500 2 : bounds.second = 1.;
501 8 : break;
502 :
503 0 : default:
504 0 : libMesh::out << "Warning: Invalid quality measure chosen." << std::endl;
505 0 : bounds.first = -1;
506 0 : bounds.second = -1;
507 : }
508 :
509 20 : return bounds;
510 : }
511 :
512 :
513 :
514 : std::vector<std::shared_ptr<Polygon>>
515 7 : Polyhedron::side_clones() const
516 : {
517 2 : const auto ns = this->n_sides();
518 :
519 2 : libmesh_assert_equal_to(ns, _sidelinks_data.size());
520 :
521 7 : std::vector<std::shared_ptr<Polygon>> cloned_sides(ns);
522 :
523 49 : for (auto i : make_range(ns))
524 : {
525 42 : const Polygon & face = *std::get<0>(this->_sidelinks_data[i]);
526 :
527 42 : Elem * clone = face.disconnected_clone().release();
528 12 : Polygon * polygon_clone = cast_ptr<Polygon *>(clone);
529 42 : cloned_sides[i] = std::shared_ptr<Polygon>(polygon_clone);
530 :
531 : // We can't actually use a *disconnected* clone to reconstruct
532 : // links between sides, so we'll temporarily give the clone our
533 : // own nodes; user code that typically replaces the usual
534 : // nullptr with permanent nodes will then instead place our
535 : // nodes with permanent nodes.
536 240 : for (auto n : make_range(face.n_nodes()))
537 96 : cloned_sides[i]->set_node
538 216 : (n, const_cast<Node *>(face.node_ptr(n)));
539 : }
540 :
541 9 : return cloned_sides;
542 0 : }
543 :
544 :
545 :
546 768 : bool Polyhedron::side_has_edge_nodes(unsigned int s,
547 : unsigned int min_node,
548 : unsigned int max_node) const
549 : {
550 768 : const Polygon & face = *std::get<0>(_sidelinks_data[s]);
551 : const std::vector<unsigned int> & node_map =
552 192 : std::get<2>(this->_sidelinks_data[s]);
553 :
554 3260 : for (unsigned int e : make_range(face.n_sides()))
555 : {
556 : std::vector<unsigned int> nodes_on_edge =
557 2732 : face.nodes_on_side(e);
558 683 : libmesh_assert_equal_to(nodes_on_edge.size(), 2);
559 2732 : nodes_on_edge[0] = node_map[nodes_on_edge[0]];
560 2732 : nodes_on_edge[1] = node_map[nodes_on_edge[1]];
561 2807 : if ((nodes_on_edge[0] == min_node) &&
562 75 : (nodes_on_edge[1] == max_node))
563 50 : return true;
564 2717 : if ((nodes_on_edge[1] == min_node) &&
565 85 : (nodes_on_edge[0] == max_node))
566 70 : return true;
567 : }
568 :
569 264 : return false;
570 : }
571 :
572 :
573 :
574 : std::vector<unsigned int>
575 192 : Polyhedron::sides_on_edge(const unsigned int e) const
576 : {
577 192 : std::vector<unsigned int> returnval(2);
578 192 : auto [s1, s1e] = _edge_lookup[e];
579 192 : returnval[0] = s1;
580 :
581 192 : const Polygon & face1 = *std::get<0>(_sidelinks_data[s1]);
582 : const std::vector<unsigned int> & node_map =
583 48 : std::get<2>(this->_sidelinks_data[s1]);
584 :
585 : std::vector<unsigned int> nodes_on_edge =
586 192 : face1.nodes_on_side(s1e);
587 48 : libmesh_assert_equal_to(nodes_on_edge.size(), 2);
588 192 : nodes_on_edge[0] = node_map[nodes_on_edge[0]];
589 192 : nodes_on_edge[1] = node_map[nodes_on_edge[1]];
590 :
591 192 : if (nodes_on_edge[0] > nodes_on_edge[1])
592 12 : std::swap(nodes_on_edge[0], nodes_on_edge[1]);
593 :
594 912 : for (unsigned int s2 : make_range(this->n_sides()))
595 : {
596 912 : if (s2 == s1)
597 144 : continue;
598 :
599 720 : if (this->side_has_edge_nodes(s2, nodes_on_edge[0],
600 360 : nodes_on_edge[1]))
601 : {
602 192 : returnval[1] = s2;
603 240 : return returnval;
604 : }
605 : }
606 :
607 0 : libmesh_error();
608 :
609 : return returnval;
610 : }
611 :
612 :
613 :
614 96 : bool Polyhedron::is_edge_on_side(const unsigned int e,
615 : const unsigned int s) const
616 : {
617 96 : auto [s1, s1e] = _edge_lookup[e];
618 :
619 : // Did we get lucky with our cache?
620 96 : if (s1 == s)
621 12 : return true;
622 :
623 48 : const Polygon & face1 = *std::get<0>(_sidelinks_data[s1]);
624 : const std::vector<unsigned int> & node_map =
625 12 : std::get<2>(this->_sidelinks_data[s1]);
626 : std::vector<unsigned int> nodes_on_edge1 =
627 60 : face1.nodes_on_side(s1e);
628 12 : libmesh_assert_equal_to(nodes_on_edge1.size(), 2);
629 :
630 48 : nodes_on_edge1[0] = node_map[nodes_on_edge1[0]];
631 48 : nodes_on_edge1[1] = node_map[nodes_on_edge1[1]];
632 48 : if (nodes_on_edge1[0] > nodes_on_edge1[1])
633 3 : std::swap(nodes_on_edge1[0], nodes_on_edge1[1]);
634 :
635 48 : return this->side_has_edge_nodes(s,
636 12 : nodes_on_edge1[0],
637 24 : nodes_on_edge1[1]);
638 : }
639 :
640 :
641 :
642 73534159 : std::array<Point, 4> Polyhedron::master_subelement (unsigned int i) const
643 : {
644 19024346 : libmesh_assert_less(i, this->_triangulation.size());
645 :
646 73534159 : const auto & tet = this->_triangulation[i];
647 :
648 73534159 : return { this->master_point(tet[0]),
649 73534159 : this->master_point(tet[1]),
650 73534159 : this->master_point(tet[2]),
651 73534159 : this->master_point(tet[3]) };
652 : }
653 :
654 :
655 :
656 : std::tuple<unsigned int, Real, Real, Real>
657 19701290 : Polyhedron::subelement_coordinates (const Point & p, Real tol) const
658 : {
659 5071145 : std::tuple<unsigned int, Real, Real, Real> returnval =
660 : {libMesh::invalid_uint, -1, -1, -1};
661 :
662 5071145 : Real best_bad_coord = -1;
663 :
664 68101272 : for (auto s : make_range(this->n_subelements()))
665 : {
666 : const std::array<Point, 4> subtet =
667 66467680 : this->master_subelement(s);
668 :
669 : // Find barycentric coordinates in subelem
670 17257345 : const Point v0 = p - subtet[0];
671 : // const Point v1 = p - subtet[1];
672 :
673 17257345 : const Point v01 = subtet[1] - subtet[0];
674 17257345 : const Point v02 = subtet[2] - subtet[0];
675 17257345 : const Point v03 = subtet[3] - subtet[0];
676 :
677 : // const Point v12 = subtet[2] - subtet[1];
678 : // const Point v13 = subtet[3] - subtet[1];
679 :
680 : // const Real tp0 = triple_product(v1, v13, v12);
681 50049985 : const Real tp1 = triple_product(v0, v02, v03);
682 50049985 : const Real tp2 = triple_product(v0, v03, v01);
683 50049985 : const Real tp3 = triple_product(v0, v01, v02);
684 :
685 50049985 : const Real six_vol = triple_product(v01, v02, v03);
686 :
687 66467680 : const Real xi = tp1 / six_vol;
688 66467680 : const Real eta = tp2 / six_vol;
689 66467680 : const Real zeta = tp3 / six_vol;
690 :
691 66467680 : if (xi>=0 && eta>=0 && zeta>=0 && xi+eta+zeta<=1)
692 18067698 : return { s, xi, eta, zeta };
693 :
694 : const Real my_best_bad_coord =
695 48399982 : std::min(std::min(std::min(xi, eta), zeta), 1-xi-eta-zeta);
696 :
697 48399982 : if (my_best_bad_coord > best_bad_coord)
698 : {
699 5202878 : best_bad_coord = my_best_bad_coord;
700 5202878 : returnval = { s, xi, eta, zeta };
701 : }
702 : }
703 :
704 1633592 : if (best_bad_coord > -tol)
705 408398 : return returnval;
706 :
707 0 : return {libMesh::invalid_uint, -1, -1, -1};
708 : }
709 :
710 :
711 : } // namespace libMesh
|
