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 3991 : Polyhedron::Polyhedron (const std::vector<std::shared_ptr<Polygon>> & sides,
43 3991 : Elem * p) :
44 228 : Cell(/* unused here */ 0, sides.size(), p, nullptr, nullptr),
45 228 : _elemlinks_data(sides.size()+2), // neighbors + parent + interior_parent
46 3763 : _nodelinks_data(0), // We'll have to resize *this* later too!
47 7868 : _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 3991 : unsigned int nn = 0;
54 228 : std::unordered_map<Node *, unsigned int> local_node_number;
55 : std::unordered_set<std::pair<const Node *, const Node *>,
56 228 : libMesh::hash> edges_seen;
57 3991 : std::unique_ptr<const Elem> edge;
58 27937 : for (unsigned int s : index_range(sides))
59 : {
60 684 : libmesh_assert(sides[s].get());
61 23946 : auto & side_tuple = _sidelinks_data[s];
62 1368 : std::get<0>(side_tuple) = sides[s];
63 :
64 1368 : Polygon & side = *sides[s]; // not const, for writeable nodes
65 119730 : for (auto n : make_range(side.n_nodes()))
66 : {
67 98520 : Node * node = side.node_ptr(n);
68 95784 : if (auto it = local_node_number.find(node);
69 2736 : it != local_node_number.end())
70 : {
71 63856 : std::get<2>(side_tuple).push_back(it->second);
72 : }
73 : else
74 : {
75 31928 : std::get<2>(side_tuple).push_back(nn);
76 31928 : local_node_number[node] = nn++;
77 31928 : _nodelinks_data.push_back(node);
78 : }
79 : }
80 :
81 119730 : for (unsigned int e : make_range(side.n_edges()))
82 : {
83 95784 : side.build_edge_ptr(edge, e);
84 95784 : auto edge_vertices = std::make_pair(edge->node_ptr(0), edge->node_ptr(1));
85 95784 : if (edge_vertices.first > edge_vertices.second)
86 1444 : std::swap(edge_vertices.first, edge_vertices.second);
87 :
88 5472 : if (!edges_seen.count(edge_vertices))
89 : {
90 1368 : edges_seen.insert(edge_vertices);
91 47892 : _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 3991 : this->_elemlinks = _elemlinks_data.data();
100 3991 : this->_nodes = _nodelinks_data.data();
101 3991 : this->_elemlinks[0] = p;
102 :
103 114 : 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 114 : Point center;
110 35919 : for (auto n : make_range(nn))
111 912 : center.add (this->point(n));
112 3991 : center /= static_cast<Real>(nn);
113 :
114 27937 : for (unsigned int s : index_range(sides))
115 : {
116 1368 : const Polygon & side = *sides[s];
117 24630 : const Point x_i = side.point(0);
118 : const Point n_i =
119 23262 : (side.point(1) - side.point(0)).cross
120 24630 : (side.point(0) - side.point(side.n_sides()-1)).unit();
121 :
122 1368 : bool & inward_normal = std::get<1>(_sidelinks_data[s]);
123 23946 : 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 798 : for (unsigned int s : index_range(sides))
130 : {
131 684 : const Polygon & side = *sides[s];
132 684 : const Point x_i = side.point(0);
133 684 : const bool inward_normal = std::get<1>(this->_sidelinks_data[s]);
134 :
135 : const Point n_i =
136 684 : (side.point(1) - side.point(0)).cross
137 684 : (side.point(0) - side.point(side.n_sides()-1)).unit() *
138 1368 : (inward_normal ? -1 : 1);
139 :
140 6156 : for (const Point & node : this->node_ref_range())
141 : {
142 5472 : const Point d_n = node - x_i;
143 5472 : 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 3991 : 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 3991 : this->set_interior_parent(nullptr);
166 7754 : }
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 930 : return const_cast<const Polyhedron *>(this)->side_ptr(i);
322 : }
323 :
324 :
325 :
326 1782 : std::unique_ptr<Elem> Polyhedron::side_ptr (const unsigned int i) const
327 : {
328 78 : libmesh_assert_less (i, this->n_sides());
329 :
330 1782 : Polygon & face = *std::get<0>(this->_sidelinks_data[i]);
331 1782 : std::unique_ptr<Elem> face_copy = face.disconnected_clone();
332 8910 : for (auto n : face.node_index_range())
333 7440 : face_copy->set_node(n, face.node_ptr(n));
334 :
335 1782 : return face_copy;
336 0 : }
337 :
338 :
339 :
340 276 : void Polyhedron::side_ptr (std::unique_ptr<Elem> & side,
341 : const unsigned int i)
342 : {
343 23 : libmesh_assert_less (i, this->n_sides());
344 :
345 : // Polyhedra are irregular enough that we're not even going to try
346 : // and bother optimizing heap access here.
347 506 : side = this->side_ptr(i);
348 276 : }
349 :
350 :
351 :
352 654 : std::unique_ptr<Elem> Polyhedron::build_side_ptr (const unsigned int i)
353 : {
354 654 : auto returnval = this->side_ptr(i);
355 654 : returnval->set_interior_parent(this);
356 623 : returnval->inherit_data_from(*this);
357 654 : return returnval;
358 0 : }
359 :
360 :
361 :
362 276 : void Polyhedron::build_side_ptr (std::unique_ptr<Elem> & side,
363 : const unsigned int i)
364 : {
365 276 : this->side_ptr(side, i);
366 276 : side->set_interior_parent(this);
367 253 : side->inherit_data_from(*this);
368 276 : }
369 :
370 :
371 :
372 0 : std::unique_ptr<Elem> Polyhedron::build_edge_ptr (const unsigned int i)
373 : {
374 0 : auto [s, se] = _edge_lookup[i];
375 0 : Polygon & face = *std::get<0>(_sidelinks_data[s]);
376 0 : return face.build_edge_ptr(se);
377 : }
378 :
379 :
380 :
381 0 : void Polyhedron::build_edge_ptr (std::unique_ptr<Elem> & elem,
382 : const unsigned int i)
383 : {
384 0 : auto [s, se] = _edge_lookup[i];
385 0 : Polygon & face = *std::get<0>(_sidelinks_data[s]);
386 0 : face.build_edge_ptr(elem, se);
387 0 : }
388 :
389 :
390 :
391 0 : bool Polyhedron::is_child_on_side(const unsigned int /*c*/,
392 : const unsigned int /*s*/) const
393 : {
394 0 : libmesh_not_implemented();
395 : return false;
396 : }
397 :
398 :
399 :
400 0 : unsigned int Polyhedron::opposite_side(const unsigned int /* side_in */) const
401 : {
402 : // This is too ambiguous in general.
403 0 : libmesh_not_implemented();
404 : return libMesh::invalid_uint;
405 : }
406 :
407 :
408 :
409 0 : unsigned int Polyhedron::opposite_node(const unsigned int /* n */,
410 : const unsigned int /* s */) const
411 : {
412 : // This is too ambiguous in general.
413 0 : libmesh_not_implemented();
414 : return libMesh::invalid_uint;
415 : }
416 :
417 :
418 :
419 134 : bool Polyhedron::is_flipped() const
420 : {
421 134 : if (this->_triangulation.empty())
422 0 : return false;
423 :
424 8 : auto & tet = this->_triangulation[0];
425 :
426 142 : const Point v01 = this->point(tet[1]) - this->point(tet[0]);
427 134 : const Point v02 = this->point(tet[2]) - this->point(tet[0]);
428 134 : const Point v03 = this->point(tet[3]) - this->point(tet[0]);
429 :
430 134 : return (triple_product(v01, v02, v03) < 0);
431 : }
432 :
433 :
434 : std::vector<unsigned int>
435 480 : Polyhedron::edges_adjacent_to_node(const unsigned int n) const
436 : {
437 40 : libmesh_assert_less(n, this->n_nodes());
438 :
439 : // For mid-edge or mid-face nodes, the subclass had better have
440 : // overridden this.
441 40 : libmesh_assert_less(n, this->n_vertices());
442 :
443 40 : const unsigned int ne = this->n_edges();
444 40 : libmesh_assert_equal_to(ne, _edge_lookup.size());
445 :
446 40 : std::vector<unsigned int> adjacent_edges;
447 :
448 480 : unsigned int next_edge = 0;
449 :
450 : // Look for any adjacent edges on each side. Make use of the fact
451 : // that we number our edges in order as we encounter them from each
452 : // side.
453 2400 : for (auto t : index_range(_sidelinks_data))
454 : {
455 400 : const Polygon & face = *std::get<0>(_sidelinks_data[t]);
456 : const std::vector<unsigned int> & node_map =
457 200 : std::get<2>(_sidelinks_data[t]);
458 :
459 3640 : while (_edge_lookup[next_edge].first < t)
460 : {
461 960 : ++next_edge;
462 960 : if (next_edge == ne)
463 40 : return adjacent_edges;
464 : }
465 :
466 : // If we haven't seen the next edge on this or an earlier side
467 : // then we might as well go on to the next.
468 2400 : if (_edge_lookup[next_edge].first > t)
469 0 : continue;
470 :
471 200 : const unsigned int fnv = face.n_vertices();
472 200 : libmesh_assert_equal_to(fnv, face.n_edges());
473 200 : libmesh_assert_equal_to(fnv, face.n_sides());
474 200 : libmesh_assert_less_equal(fnv, node_map.size());
475 :
476 : // Polygon faces have one edge per vertex
477 9600 : for (auto v : make_range(fnv))
478 : {
479 720 : libmesh_assert_equal_to (_edge_lookup[next_edge].first, t);
480 :
481 9360 : if (_edge_lookup[next_edge].second > v)
482 1320 : continue;
483 :
484 13360 : while (_edge_lookup[next_edge].first == t &&
485 10560 : _edge_lookup[next_edge].second < v)
486 : {
487 4800 : ++next_edge;
488 4800 : if (next_edge == ne)
489 40 : return adjacent_edges;
490 : }
491 :
492 6720 : if (_edge_lookup[next_edge].first > t)
493 80 : break;
494 :
495 5760 : const unsigned int vn = node_map[v];
496 5760 : const unsigned int vnp = node_map[(v+1)%fnv];
497 480 : libmesh_assert_less(vn, this->n_vertices());
498 480 : libmesh_assert_less(vnp, this->n_vertices());
499 5760 : if (vn == n || vnp == n)
500 1440 : adjacent_edges.push_back(next_edge);
501 : }
502 : }
503 :
504 0 : return adjacent_edges;
505 : }
506 :
507 :
508 0 : std::pair<Real, Real> Polyhedron::qual_bounds (const ElemQuality q) const
509 : {
510 0 : std::pair<Real, Real> bounds;
511 :
512 0 : switch (q)
513 : {
514 0 : case EDGE_LENGTH_RATIO:
515 0 : bounds.first = 1.;
516 0 : bounds.second = 4.;
517 0 : break;
518 :
519 0 : case MIN_ANGLE:
520 0 : bounds.first = 30.;
521 0 : bounds.second = 180.;
522 0 : break;
523 :
524 0 : case MAX_ANGLE:
525 0 : bounds.first = 60.;
526 0 : bounds.second = 180.;
527 0 : break;
528 :
529 0 : case JACOBIAN:
530 : case SCALED_JACOBIAN:
531 0 : bounds.first = 0.5;
532 0 : bounds.second = 1.;
533 0 : break;
534 :
535 0 : default:
536 0 : libMesh::out << "Warning: Invalid quality measure chosen." << std::endl;
537 0 : bounds.first = -1;
538 0 : bounds.second = -1;
539 : }
540 :
541 0 : return bounds;
542 : }
543 :
544 :
545 :
546 : std::vector<std::shared_ptr<Polygon>>
547 15 : Polyhedron::side_clones() const
548 : {
549 2 : const auto ns = this->n_sides();
550 :
551 2 : libmesh_assert_equal_to(ns, _sidelinks_data.size());
552 :
553 15 : std::vector<std::shared_ptr<Polygon>> cloned_sides(ns);
554 :
555 105 : for (auto i : make_range(ns))
556 : {
557 90 : const Polygon & face = *std::get<0>(this->_sidelinks_data[i]);
558 :
559 90 : Elem * clone = face.disconnected_clone().release();
560 12 : Polygon * polygon_clone = cast_ptr<Polygon *>(clone);
561 90 : cloned_sides[i] = std::shared_ptr<Polygon>(polygon_clone);
562 :
563 : // We can't actually use a *disconnected* clone to reconstruct
564 : // links between sides, so we'll temporarily give the clone our
565 : // own nodes; user code that typically replaces the usual
566 : // nullptr with permanent nodes will then instead place our
567 : // nodes with permanent nodes.
568 528 : for (auto n : make_range(face.n_nodes()))
569 96 : cloned_sides[i]->set_node
570 408 : (n, const_cast<Node *>(face.node_ptr(n)));
571 : }
572 :
573 17 : return cloned_sides;
574 0 : }
575 :
576 :
577 :
578 2208 : bool Polyhedron::side_has_edge_nodes(unsigned int s,
579 : unsigned int min_node,
580 : unsigned int max_node) const
581 : {
582 2208 : const Polygon & face = *std::get<0>(_sidelinks_data[s]);
583 : const std::vector<unsigned int> & node_map =
584 184 : std::get<2>(this->_sidelinks_data[s]);
585 :
586 9300 : for (unsigned int e : make_range(face.n_sides()))
587 : {
588 : std::vector<unsigned int> nodes_on_edge =
589 7812 : face.nodes_on_side(e);
590 651 : libmesh_assert_equal_to(nodes_on_edge.size(), 2);
591 7812 : nodes_on_edge[0] = node_map[nodes_on_edge[0]];
592 7812 : nodes_on_edge[1] = node_map[nodes_on_edge[1]];
593 7901 : if ((nodes_on_edge[0] == min_node) &&
594 89 : (nodes_on_edge[1] == max_node))
595 60 : return true;
596 7536 : if ((nodes_on_edge[1] == min_node) &&
597 84 : (nodes_on_edge[0] == max_node))
598 60 : return true;
599 : }
600 :
601 248 : return false;
602 : }
603 :
604 :
605 :
606 : std::vector<unsigned int>
607 576 : Polyhedron::sides_on_edge(const unsigned int e) const
608 : {
609 576 : std::vector<unsigned int> returnval(2);
610 576 : auto [s1, s1e] = _edge_lookup[e];
611 576 : returnval[0] = s1;
612 :
613 576 : const Polygon & face1 = *std::get<0>(_sidelinks_data[s1]);
614 : const std::vector<unsigned int> & node_map =
615 48 : std::get<2>(this->_sidelinks_data[s1]);
616 :
617 : std::vector<unsigned int> nodes_on_edge =
618 576 : face1.nodes_on_side(s1e);
619 48 : libmesh_assert_equal_to(nodes_on_edge.size(), 2);
620 576 : nodes_on_edge[0] = node_map[nodes_on_edge[0]];
621 576 : nodes_on_edge[1] = node_map[nodes_on_edge[1]];
622 :
623 576 : if (nodes_on_edge[0] > nodes_on_edge[1])
624 16 : std::swap(nodes_on_edge[0], nodes_on_edge[1]);
625 :
626 2640 : for (unsigned int s2 : make_range(this->n_sides()))
627 : {
628 2640 : if (s2 == s1)
629 528 : continue;
630 :
631 2064 : if (this->side_has_edge_nodes(s2, nodes_on_edge[0],
632 344 : nodes_on_edge[1]))
633 : {
634 576 : returnval[1] = s2;
635 624 : return returnval;
636 : }
637 : }
638 :
639 0 : libmesh_error();
640 :
641 : return returnval;
642 : }
643 :
644 :
645 :
646 288 : bool Polyhedron::is_edge_on_side(const unsigned int e,
647 : const unsigned int s) const
648 : {
649 288 : auto [s1, s1e] = _edge_lookup[e];
650 :
651 : // Did we get lucky with our cache?
652 288 : if (s1 == s)
653 12 : return true;
654 :
655 144 : const Polygon & face1 = *std::get<0>(_sidelinks_data[s1]);
656 : const std::vector<unsigned int> & node_map =
657 12 : std::get<2>(this->_sidelinks_data[s1]);
658 : std::vector<unsigned int> nodes_on_edge1 =
659 156 : face1.nodes_on_side(s1e);
660 12 : libmesh_assert_equal_to(nodes_on_edge1.size(), 2);
661 :
662 144 : nodes_on_edge1[0] = node_map[nodes_on_edge1[0]];
663 144 : nodes_on_edge1[1] = node_map[nodes_on_edge1[1]];
664 144 : if (nodes_on_edge1[0] > nodes_on_edge1[1])
665 4 : std::swap(nodes_on_edge1[0], nodes_on_edge1[1]);
666 :
667 144 : return this->side_has_edge_nodes(s,
668 12 : nodes_on_edge1[0],
669 24 : nodes_on_edge1[1]);
670 : }
671 :
672 :
673 :
674 92183699 : std::array<Point, 4> Polyhedron::master_subelement (unsigned int i) const
675 : {
676 7990237 : libmesh_assert_less(i, this->_triangulation.size());
677 :
678 92183699 : const auto & tet = this->_triangulation[i];
679 :
680 92183699 : return { this->master_point(tet[0]),
681 92183699 : this->master_point(tet[1]),
682 92183699 : this->master_point(tet[2]),
683 92183699 : this->master_point(tet[3]) };
684 : }
685 :
686 :
687 : std::tuple<unsigned int, Real, Real, Real>
688 31795656 : Polyhedron::subelement_coordinates (const Point & p, Real tol) const
689 : {
690 2743208 : std::tuple<unsigned int, Real, Real, Real> returnval =
691 : {libMesh::invalid_uint, -1, -1, -1};
692 :
693 2743208 : Real best_bad_coord = -1;
694 :
695 82709560 : for (auto s : make_range(this->n_subelements()))
696 : {
697 : const std::array<Point, 4> subtet =
698 80051512 : this->master_subelement(s);
699 :
700 : // Find barycentric coordinates in subelem
701 6980416 : const Point v0 = p - subtet[0];
702 : // const Point v1 = p - subtet[1];
703 :
704 6980416 : const Point v01 = subtet[1] - subtet[0];
705 6980416 : const Point v02 = subtet[2] - subtet[0];
706 6980416 : const Point v03 = subtet[3] - subtet[0];
707 :
708 : // const Point v12 = subtet[2] - subtet[1];
709 : // const Point v13 = subtet[3] - subtet[1];
710 :
711 : // const Real tp0 = triple_product(v1, v13, v12);
712 73416376 : const Real tp1 = triple_product(v0, v02, v03);
713 73416376 : const Real tp2 = triple_product(v0, v03, v01);
714 73416376 : const Real tp3 = triple_product(v0, v01, v02);
715 :
716 73416376 : const Real six_vol = triple_product(v01, v02, v03);
717 :
718 80051512 : const Real xi = tp1 / six_vol;
719 80051512 : const Real eta = tp2 / six_vol;
720 80051512 : const Real zeta = tp3 / six_vol;
721 :
722 80051512 : if (xi>=0 && eta>=0 && zeta>=0 && xi+eta+zeta<=1)
723 29137608 : return { s, xi, eta, zeta };
724 :
725 : const Real my_best_bad_coord =
726 50913904 : std::min(std::min(std::min(xi, eta), zeta), 1-xi-eta-zeta);
727 :
728 50913904 : if (my_best_bad_coord > best_bad_coord)
729 : {
730 2173880 : best_bad_coord = my_best_bad_coord;
731 2173880 : returnval = { s, xi, eta, zeta };
732 : }
733 : }
734 :
735 2658048 : if (best_bad_coord > -tol)
736 221504 : return returnval;
737 :
738 0 : return {libMesh::invalid_uint, -1, -1, -1};
739 : }
740 :
741 :
742 : } // namespace libMesh
|
