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 :
19 : // Local includes
20 : #include "libmesh/elem.h"
21 :
22 : #include "libmesh/boundary_info.h"
23 : #include "libmesh/fe_type.h"
24 : #include "libmesh/fe_interface.h"
25 : #include "libmesh/tensor_value.h"
26 : #include "libmesh/node_elem.h"
27 : #include "libmesh/edge_edge2.h"
28 : #include "libmesh/edge_edge3.h"
29 : #include "libmesh/edge_edge4.h"
30 : #include "libmesh/edge_inf_edge2.h"
31 : #include "libmesh/face_c0polygon.h"
32 : #include "libmesh/face_tri3.h"
33 : #include "libmesh/face_tri3_subdivision.h"
34 : #include "libmesh/face_tri3_shell.h"
35 : #include "libmesh/face_tri6.h"
36 : #include "libmesh/face_tri7.h"
37 : #include "libmesh/face_quad4.h"
38 : #include "libmesh/face_quad4_shell.h"
39 : #include "libmesh/face_quad8.h"
40 : #include "libmesh/face_quad8_shell.h"
41 : #include "libmesh/face_quad9.h"
42 : #include "libmesh/face_quad9_shell.h"
43 : #include "libmesh/face_inf_quad4.h"
44 : #include "libmesh/face_inf_quad6.h"
45 : #include "libmesh/cell_tet4.h"
46 : #include "libmesh/cell_tet10.h"
47 : #include "libmesh/cell_tet14.h"
48 : #include "libmesh/cell_hex8.h"
49 : #include "libmesh/cell_hex20.h"
50 : #include "libmesh/cell_hex27.h"
51 : #include "libmesh/cell_inf_hex8.h"
52 : #include "libmesh/cell_inf_hex16.h"
53 : #include "libmesh/cell_inf_hex18.h"
54 : #include "libmesh/cell_prism6.h"
55 : #include "libmesh/cell_prism15.h"
56 : #include "libmesh/cell_prism18.h"
57 : #include "libmesh/cell_prism20.h"
58 : #include "libmesh/cell_prism21.h"
59 : #include "libmesh/cell_inf_prism6.h"
60 : #include "libmesh/cell_inf_prism12.h"
61 : #include "libmesh/cell_pyramid5.h"
62 : #include "libmesh/cell_pyramid13.h"
63 : #include "libmesh/cell_pyramid14.h"
64 : #include "libmesh/cell_pyramid18.h"
65 : #include "libmesh/fe_base.h"
66 : #include "libmesh/mesh_base.h"
67 : #include "libmesh/quadrature_nodal.h"
68 : #include "libmesh/quadrature_gauss.h"
69 : #include "libmesh/remote_elem.h"
70 : #include "libmesh/reference_elem.h"
71 : #include "libmesh/enum_to_string.h"
72 : #include "libmesh/threads.h"
73 : #include "libmesh/enum_elem_quality.h"
74 : #include "libmesh/enum_io_package.h"
75 : #include "libmesh/enum_order.h"
76 : #include "libmesh/elem_internal.h"
77 :
78 : #ifdef LIBMESH_ENABLE_PERIODIC
79 : #include "libmesh/mesh.h"
80 : #include "libmesh/periodic_boundaries.h"
81 : #endif
82 :
83 :
84 : // C++ includes
85 : #include <algorithm> // for std::sort
86 : #include <array>
87 : #include <iterator> // for std::ostream_iterator
88 : #include <sstream>
89 : #include <limits> // for std::numeric_limits<>
90 : #include <cmath> // for std::sqrt()
91 : #include <memory>
92 : #include <regex> // for exceptions in volume()
93 :
94 :
95 : namespace libMesh
96 : {
97 :
98 : Threads::spin_mutex parent_indices_mutex;
99 : Threads::spin_mutex parent_bracketing_nodes_mutex;
100 :
101 : // Initialize static member variables
102 : const unsigned int Elem::type_to_dim_map [] =
103 : {
104 : 1, // EDGE2
105 : 1, // EDGE3
106 : 1, // EDGE4
107 :
108 : 2, // TRI3
109 : 2, // TRI6
110 :
111 : 2, // QUAD4
112 : 2, // QUAD8
113 : 2, // QUAD9
114 :
115 : 3, // TET4
116 : 3, // TET10
117 :
118 : 3, // HEX8
119 : 3, // HEX20
120 : 3, // HEX27
121 :
122 : 3, // PRISM6
123 : 3, // PRISM15
124 : 3, // PRISM18
125 :
126 : 3, // PYRAMID5
127 : 3, // PYRAMID13
128 : 3, // PYRAMID14
129 :
130 : 1, // INFEDGE2
131 :
132 : 2, // INFQUAD4
133 : 2, // INFQUAD6
134 :
135 : 3, // INFHEX8
136 : 3, // INFHEX16
137 : 3, // INFHEX18
138 :
139 : 3, // INFPRISM6
140 : 3, // INFPRISM12
141 :
142 : 0, // NODEELEM
143 :
144 : 0, // REMOTEELEM
145 :
146 : 2, // TRI3SUBDIVISION
147 : 2, // TRISHELL3
148 : 2, // QUADSHELL4
149 : 2, // QUADSHELL8
150 :
151 : 2, // TRI7
152 : 3, // TET14
153 : 3, // PRISM20
154 : 3, // PRISM21
155 : 3, // PYRAMID18
156 :
157 : 2, // QUADSHELL9
158 :
159 : 2, // C0POLYGON
160 : 3, // C0POLYHEDRON
161 : };
162 :
163 : const unsigned int Elem::max_n_nodes;
164 :
165 : const unsigned int Elem::type_to_n_sides_map [] =
166 : {
167 : 2, // EDGE2
168 : 2, // EDGE3
169 : 2, // EDGE4
170 :
171 : 3, // TRI3
172 : 3, // TRI6
173 :
174 : 4, // QUAD4
175 : 4, // QUAD8
176 : 4, // QUAD9
177 :
178 : 4, // TET4
179 : 4, // TET10
180 :
181 : 6, // HEX8
182 : 6, // HEX20
183 : 6, // HEX27
184 :
185 : 5, // PRISM6
186 : 5, // PRISM15
187 : 5, // PRISM18
188 :
189 : 5, // PYRAMID5
190 : 5, // PYRAMID13
191 : 5, // PYRAMID14
192 :
193 : 2, // INFEDGE2
194 :
195 : 3, // INFQUAD4
196 : 3, // INFQUAD6
197 :
198 : 5, // INFHEX8
199 : 5, // INFHEX16
200 : 5, // INFHEX18
201 :
202 : 4, // INFPRISM6
203 : 4, // INFPRISM12
204 :
205 : 0, // NODEELEM
206 :
207 : 0, // REMOTEELEM
208 :
209 : 3, // TRI3SUBDIVISION
210 : 3, // TRISHELL3
211 : 4, // QUADSHELL4
212 : 4, // QUADSHELL8
213 :
214 : 3, // TRI7
215 : 4, // TET14
216 : 5, // PRISM20
217 : 5, // PRISM21
218 : 5, // PYRAMID18
219 :
220 : 4, // QUADSHELL9
221 :
222 : invalid_uint, // C0POLYGON
223 : invalid_uint, // C0POLYHEDRON
224 : };
225 :
226 : const unsigned int Elem::type_to_n_edges_map [] =
227 : {
228 : 0, // EDGE2
229 : 0, // EDGE3
230 : 0, // EDGE4
231 :
232 : 3, // TRI3
233 : 3, // TRI6
234 :
235 : 4, // QUAD4
236 : 4, // QUAD8
237 : 4, // QUAD9
238 :
239 : 6, // TET4
240 : 6, // TET10
241 :
242 : 12, // HEX8
243 : 12, // HEX20
244 : 12, // HEX27
245 :
246 : 9, // PRISM6
247 : 9, // PRISM15
248 : 9, // PRISM18
249 :
250 : 8, // PYRAMID5
251 : 8, // PYRAMID13
252 : 8, // PYRAMID14
253 :
254 : 0, // INFEDGE2
255 :
256 : 3, // INFQUAD4
257 : 3, // INFQUAD6
258 :
259 : 8, // INFHEX8
260 : 8, // INFHEX16
261 : 8, // INFHEX18
262 :
263 : 6, // INFPRISM6
264 : 6, // INFPRISM12
265 :
266 : 0, // NODEELEM
267 :
268 : 0, // REMOTEELEM
269 :
270 : 3, // TRI3SUBDIVISION
271 : 3, // TRISHELL3
272 : 4, // QUADSHELL4
273 : 4, // QUADSHELL8
274 :
275 : 3, // TRI7
276 : 6, // TET14
277 : 9, // PRISM20
278 : 9, // PRISM21
279 : 8, // PYRAMID18
280 :
281 : 4, // QUADSHELL9
282 :
283 : invalid_uint, // C0POLYGON
284 : invalid_uint, // C0POLYHEDRON
285 : };
286 :
287 : const Order Elem::type_to_default_order_map [] =
288 : {
289 : FIRST, // EDGE2
290 : SECOND, // EDGE3
291 : THIRD, // EDGE4
292 :
293 : FIRST, // TRI3
294 : SECOND, // TRI6
295 :
296 : FIRST, // QUAD4
297 : SECOND, // QUAD8
298 : SECOND, // QUAD9
299 :
300 : FIRST, // TET4
301 : SECOND, // TET10
302 :
303 : FIRST, // HEX8
304 : SECOND, // HEX20
305 : SECOND, // HEX27
306 :
307 : FIRST, // PRISM6
308 : SECOND, // PRISM15
309 : SECOND, // PRISM18
310 :
311 : FIRST, // PYRAMID5
312 : SECOND, // PYRAMID13
313 : SECOND, // PYRAMID14
314 :
315 : FIRST, // INFEDGE2
316 :
317 : FIRST, // INFQUAD4
318 : SECOND, // INFQUAD6
319 :
320 : FIRST, // INFHEX8
321 : SECOND, // INFHEX16
322 : SECOND, // INFHEX18
323 :
324 : FIRST, // INFPRISM6
325 : SECOND, // INFPRISM12
326 :
327 : CONSTANT, // NODEELEM
328 :
329 : INVALID_ORDER, // REMOTEELEM
330 :
331 : FIRST, // TRI3SUBDIVISION
332 : FIRST, // TRISHELL3
333 : FIRST, // QUADSHELL4
334 : SECOND, // QUADSHELL8
335 :
336 : THIRD, // TRI7
337 : THIRD, // TET14
338 : THIRD, // PRISM20
339 : THIRD, // PRISM21
340 : THIRD, // PYRAMID18
341 :
342 : SECOND, // QUADSHELL9
343 :
344 : FIRST, // C0POLYGON
345 : FIRST, // C0POLYHEDRON
346 : };
347 :
348 : // ------------------------------------------------------------
349 : // Elem class member functions
350 10113341 : std::unique_ptr<Elem> Elem::disconnected_clone() const
351 : {
352 10113341 : std::unique_ptr<Elem> returnval;
353 :
354 10113341 : switch (this->type())
355 : {
356 31667 : case C0POLYGON:
357 31667 : returnval = std::make_unique<C0Polygon>(this->n_sides());
358 31667 : break;
359 :
360 10081674 : default:
361 19839952 : returnval = Elem::build(this->type());
362 : }
363 :
364 10113341 : returnval->set_id() = this->id();
365 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
366 10113341 : if (this->valid_unique_id())
367 323398 : returnval->set_unique_id(this->unique_id());
368 : #endif
369 :
370 10113341 : const auto n_elem_ints = this->n_extra_integers();
371 10113341 : returnval->add_extra_integers(n_elem_ints);
372 10130594 : for (unsigned int i = 0; i != n_elem_ints; ++i)
373 17253 : returnval->set_extra_integer(i, this->get_extra_integer(i));
374 :
375 9799609 : returnval->inherit_data_from(*this);
376 :
377 10113341 : return returnval;
378 0 : }
379 :
380 :
381 :
382 67721845 : std::unique_ptr<Elem> Elem::build(const ElemType type,
383 : Elem * p)
384 : {
385 67721845 : switch (type)
386 : {
387 : // 0D elements
388 223064 : case NODEELEM:
389 223064 : return std::make_unique<NodeElem>(p);
390 :
391 : // 1D elements
392 510191 : case EDGE2:
393 510191 : return std::make_unique<Edge2>(p);
394 240882 : case EDGE3:
395 240882 : return std::make_unique<Edge3>(p);
396 51640 : case EDGE4:
397 51640 : return std::make_unique<Edge4>(p);
398 :
399 : // 2D elements
400 4517541 : case TRI3:
401 4517541 : return std::make_unique<Tri3>(p);
402 81506 : case TRISHELL3:
403 81506 : return std::make_unique<TriShell3>(p);
404 62956 : case TRI3SUBDIVISION:
405 62956 : return std::make_unique<Tri3Subdivision>(p);
406 1064976 : case TRI6:
407 1064976 : return std::make_unique<Tri6>(p);
408 590308 : case TRI7:
409 590308 : return std::make_unique<Tri7>(p);
410 25443130 : case QUAD4:
411 25443130 : return std::make_unique<Quad4>(p);
412 190772 : case QUADSHELL4:
413 190772 : return std::make_unique<QuadShell4>(p);
414 266209 : case QUAD8:
415 266209 : return std::make_unique<Quad8>(p);
416 185472 : case QUADSHELL8:
417 185472 : return std::make_unique<QuadShell8>(p);
418 3271128 : case QUAD9:
419 3271128 : return std::make_unique<Quad9>(p);
420 70479 : case QUADSHELL9:
421 70479 : return std::make_unique<QuadShell9>(p);
422 :
423 : // Well, a hexagon is *a* polygon...
424 41656 : case C0POLYGON:
425 41656 : return std::make_unique<C0Polygon>(6, p);
426 :
427 : // Building a polyhedron can't currently be done without creating
428 : // its nodes first
429 0 : case C0POLYHEDRON:
430 0 : libmesh_not_implemented_msg
431 : ("Polyhedra cannot be built via Elem::build()");
432 :
433 : // 3D elements
434 16641397 : case TET4:
435 16641397 : return std::make_unique<Tet4>(p);
436 1978527 : case TET10:
437 1978527 : return std::make_unique<Tet10>(p);
438 3879982 : case TET14:
439 3879982 : return std::make_unique<Tet14>(p);
440 3352691 : case HEX8:
441 3352691 : return std::make_unique<Hex8>(p);
442 148124 : case HEX20:
443 148124 : return std::make_unique<Hex20>(p);
444 1971521 : case HEX27:
445 1971521 : return std::make_unique<Hex27>(p);
446 384826 : case PRISM6:
447 384826 : return std::make_unique<Prism6>(p);
448 194262 : case PRISM15:
449 194262 : return std::make_unique<Prism15>(p);
450 640812 : case PRISM18:
451 640812 : return std::make_unique<Prism18>(p);
452 139023 : case PRISM20:
453 139023 : return std::make_unique<Prism20>(p);
454 279817 : case PRISM21:
455 279817 : return std::make_unique<Prism21>(p);
456 637345 : case PYRAMID5:
457 637345 : return std::make_unique<Pyramid5>(p);
458 224803 : case PYRAMID13:
459 224803 : return std::make_unique<Pyramid13>(p);
460 227059 : case PYRAMID14:
461 227059 : return std::make_unique<Pyramid14>(p);
462 202088 : case PYRAMID18:
463 202088 : return std::make_unique<Pyramid18>(p);
464 :
465 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
466 : // 1D infinite elements
467 0 : case INFEDGE2:
468 0 : return std::make_unique<InfEdge2>(p);
469 :
470 : // 2D infinite elements
471 175 : case INFQUAD4:
472 175 : return std::make_unique<InfQuad4>(p);
473 175 : case INFQUAD6:
474 175 : return std::make_unique<InfQuad6>(p);
475 :
476 : // 3D infinite elements
477 3021 : case INFHEX8:
478 3021 : return std::make_unique<InfHex8>(p);
479 255 : case INFHEX16:
480 255 : return std::make_unique<InfHex16>(p);
481 1702 : case INFHEX18:
482 1702 : return std::make_unique<InfHex18>(p);
483 255 : case INFPRISM6:
484 255 : return std::make_unique<InfPrism6>(p);
485 2075 : case INFPRISM12:
486 2075 : return std::make_unique<InfPrism12>(p);
487 : #endif
488 :
489 0 : default:
490 0 : libmesh_error_msg("ERROR: Undefined element type == " << Utility::enum_to_string(type));
491 : }
492 : }
493 :
494 :
495 :
496 7111669 : std::unique_ptr<Elem> Elem::build_with_id (const ElemType type,
497 : dof_id_type id)
498 : {
499 : // Call the other build() method with nullptr parent, then set the
500 : // required id.
501 7111669 : auto temp = Elem::build(type, nullptr);
502 186847 : temp->set_id(id);
503 7111669 : return temp;
504 : }
505 :
506 :
507 :
508 418343 : const Elem * Elem::reference_elem () const
509 : {
510 418343 : return &(ReferenceElem::get(this->type()));
511 : }
512 :
513 :
514 :
515 86472 : Point Elem::true_centroid() const
516 : {
517 : // The base class implementation builds a finite element of the correct
518 : // order and computes the centroid, c=(cx, cy, cz), where:
519 : //
520 : // [cx] [\int x dV]
521 : // [cy] := (1/V) * [\int y dV]
522 : // [cz] [\int z dV]
523 : //
524 : // using quadrature. Note that we can expand "x" in the FE space as:
525 : //
526 : // x = \sum_i x_i \phi_i
527 : //
528 : // where x_i are the nodal positions of the element and \phi_i are the
529 : // associated Lagrange shape functions. This allows us to write the
530 : // integrals above as e.g.:
531 : //
532 : // \int x dV = \sum_i x_i \int \phi_i dV
533 : //
534 : // Defining:
535 : //
536 : // V_i := \int \phi_i dV
537 : //
538 : // we then have:
539 : //
540 : // [cx] [\sum_i x_i V_i]
541 : // [cy] = (1/V) * [\sum_i y_i V_i]
542 : // [cz] [\sum_i z_i V_i]
543 : //
544 : // where:
545 : // V = \sum_i V_i
546 : //
547 : // Derived element types can overload this method to compute
548 : // the centroid more efficiently when possible.
549 :
550 : // If this Elem has an elevated p_level, then we need to generate a
551 : // barebones copy of it with zero p_level and call true_centroid()
552 : // on that instead. This workaround allows us to avoid issues with
553 : // calling FE::reinit() with a default_order() FEType, and then
554 : // having that order incorrectly boosted by p_level.
555 86472 : if (this->p_level())
556 : {
557 73 : auto elem_copy = this->disconnected_clone();
558 : #ifdef LIBMESH_ENABLE_AMR
559 71 : elem_copy->set_p_level(0);
560 : #endif
561 :
562 : // Set node pointers
563 1491 : for (auto n : this->node_index_range())
564 1420 : elem_copy->set_node(n, _nodes[n]);
565 :
566 71 : return elem_copy->true_centroid();
567 67 : }
568 :
569 86401 : const FEFamily mapping_family = FEMap::map_fe_type(*this);
570 86401 : const FEType fe_type(this->default_order(), mapping_family);
571 :
572 : // Build FE and attach quadrature rule. The default quadrature rule
573 : // integrates the mass matrix exactly, thus it is overkill to
574 : // integrate the basis functions, but this is convenient.
575 93276 : std::unique_ptr<FEBase> fe = FEBase::build(this->dim(), fe_type);
576 93276 : QGauss qrule (this->dim(), fe_type.default_quadrature_order());
577 86401 : fe->attach_quadrature_rule(&qrule);
578 :
579 : // Pre-request required data
580 86401 : const auto & JxW = fe->get_JxW();
581 6875 : const auto & phi = fe->get_phi();
582 :
583 : // Re-compute element-specific values
584 86401 : fe->reinit(this);
585 :
586 : // Number of basis functions
587 13750 : auto N = phi.size();
588 6875 : libmesh_assert_equal_to(N, this->n_nodes());
589 :
590 : // Compute V_i
591 86401 : std::vector<Real> V(N);
592 2570250 : for (auto qp : index_range(JxW))
593 37985489 : for (auto i : make_range(N))
594 44331225 : V[i] += JxW[qp] * phi[i][qp];
595 :
596 : // Compute centroid
597 6875 : Point cp;
598 6875 : Real vol = 0.;
599 :
600 1193097 : for (auto i : make_range(N))
601 : {
602 361692 : cp += this->point(i) * V[i];
603 1106696 : vol += V[i];
604 : }
605 :
606 6875 : return cp / vol;
607 72651 : }
608 :
609 536533945 : Point Elem::vertex_average() const
610 : {
611 41917575 : Point cp;
612 :
613 536533945 : const auto n_vertices = this->n_vertices();
614 :
615 3182753884 : for (unsigned int n=0; n<n_vertices; n++)
616 418069090 : cp.add (this->point(n));
617 :
618 536533945 : return (cp /= static_cast<Real>(n_vertices));
619 : }
620 :
621 :
622 :
623 1855014 : Real Elem::hmin() const
624 : {
625 : // Avoid calling a virtual a lot of times
626 1855014 : const auto n_vertices = this->n_vertices();
627 :
628 : // Special case for NodeElem
629 1855014 : if (!n_vertices)
630 0 : return 0;
631 :
632 1855014 : Real h_min=std::numeric_limits<Real>::max();
633 :
634 8842856 : for (unsigned int n_outer=0; n_outer<n_vertices; n_outer++)
635 17431575 : for (unsigned int n_inner=n_outer+1; n_inner<n_vertices; n_inner++)
636 : {
637 2022368 : const auto diff = (this->point(n_outer) - this->point(n_inner));
638 :
639 13609788 : h_min = std::min(h_min, diff.norm_sq());
640 : }
641 :
642 1855014 : return std::sqrt(h_min);
643 : }
644 :
645 :
646 :
647 102240357 : Real Elem::hmax() const
648 : {
649 102240357 : Real h_max=0;
650 :
651 : // Avoid calling a virtual a lot of times
652 102240357 : const auto n_vertices = this->n_vertices();
653 :
654 672215660 : for (unsigned int n_outer=0; n_outer<n_vertices; n_outer++)
655 2143668759 : for (unsigned int n_inner=n_outer+1; n_inner<n_vertices; n_inner++)
656 : {
657 104361216 : const auto diff = (this->point(n_outer) - this->point(n_inner));
658 :
659 1808296610 : h_max = std::max(h_max, diff.norm_sq());
660 : }
661 :
662 108826852 : return std::sqrt(h_max);
663 : }
664 :
665 :
666 :
667 261488 : Real Elem::length(const unsigned int n1,
668 : const unsigned int n2) const
669 : {
670 21728 : libmesh_assert_less ( n1, this->n_vertices() );
671 21728 : libmesh_assert_less ( n2, this->n_vertices() );
672 :
673 283216 : return (this->point(n1) - this->point(n2)).norm();
674 : }
675 :
676 :
677 :
678 0 : dof_id_type Elem::key () const
679 : {
680 0 : const unsigned short n_n = this->n_nodes();
681 :
682 : std::array<dof_id_type, Elem::max_n_nodes> node_ids;
683 :
684 0 : for (unsigned short n=0; n != n_n; ++n)
685 0 : node_ids[n] = this->node_id(n);
686 :
687 : // Always sort, so that different local node numberings hash to the
688 : // same value.
689 0 : std::sort (node_ids.begin(), node_ids.begin()+n_n);
690 :
691 0 : return Utility::hashword(node_ids.data(), n_n);
692 : }
693 :
694 :
695 :
696 130998083 : bool Elem::operator == (const Elem & rhs) const
697 : {
698 : // If the elements aren't the same type, they aren't equal
699 130998083 : if (this->type() != rhs.type())
700 4 : return false;
701 :
702 130997941 : const unsigned short n_n = this->n_nodes();
703 2676504 : libmesh_assert_equal_to(n_n, rhs.n_nodes());
704 :
705 : // Make two sorted arrays of global node ids and compare them for
706 : // equality.
707 : std::array<dof_id_type, Elem::max_n_nodes> this_ids, rhs_ids;
708 :
709 450874701 : for (unsigned short n = 0; n != n_n; n++)
710 : {
711 319876760 : this_ids[n] = this->node_id(n);
712 326955054 : rhs_ids[n] = rhs.node_id(n);
713 : }
714 :
715 : // Sort the vectors to rule out different local node numberings.
716 130997941 : std::sort(this_ids.begin(), this_ids.begin()+n_n);
717 130997941 : std::sort(rhs_ids.begin(), rhs_ids.begin()+n_n);
718 :
719 : // If the node ids match, the elements are equal!
720 450813060 : for (unsigned short n = 0; n != n_n; ++n)
721 319832115 : if (this_ids[n] != rhs_ids[n])
722 0 : return false;
723 2676504 : return true;
724 : }
725 :
726 :
727 :
728 789345 : bool Elem::topologically_equal (const Elem & rhs) const
729 : {
730 : // If the elements aren't the same type, they aren't equal
731 789345 : if (this->type() != rhs.type())
732 0 : return false;
733 :
734 338536 : libmesh_assert_equal_to(this->n_nodes(), rhs.n_nodes());
735 :
736 6682215 : for (auto n : make_range(this->n_nodes()))
737 6159614 : if (this->node_id(n) != rhs.node_id(n))
738 0 : return false;
739 :
740 3953178 : for (auto neigh : make_range(this->n_neighbors()))
741 : {
742 3221540 : if (!this->neighbor_ptr(neigh))
743 : {
744 420234 : if (rhs.neighbor_ptr(neigh))
745 6 : return false;
746 400258 : continue;
747 : }
748 4050500 : if (!rhs.neighbor_ptr(neigh) ||
749 1296766 : this->neighbor_ptr(neigh)->id() !=
750 1296766 : rhs.neighbor_ptr(neigh)->id())
751 6 : return false;
752 : }
753 :
754 803998 : if (this->parent())
755 : {
756 32076 : if (!rhs.parent())
757 0 : return false;
758 31776 : if (this->parent()->id() != rhs.parent()->id())
759 0 : return false;
760 : }
761 771922 : else if (rhs.parent())
762 0 : return false;
763 :
764 789204 : if (this->interior_parent())
765 : {
766 0 : if (!rhs.interior_parent())
767 0 : return false;
768 0 : if (this->interior_parent()->id() !=
769 0 : rhs.interior_parent()->id())
770 0 : return false;
771 : }
772 789204 : else if (rhs.interior_parent())
773 0 : return false;
774 :
775 338530 : return true;
776 : }
777 :
778 :
779 :
780 0 : bool Elem::is_semilocal(const processor_id_type my_pid) const
781 : {
782 0 : std::set<const Elem *> point_neighbors;
783 :
784 0 : this->find_point_neighbors(point_neighbors);
785 :
786 0 : for (const auto & elem : point_neighbors)
787 0 : if (elem->processor_id() == my_pid)
788 0 : return true;
789 :
790 0 : return false;
791 : }
792 :
793 :
794 :
795 7709 : unsigned int Elem::which_side_am_i (const Elem * e) const
796 : {
797 708 : libmesh_assert(e);
798 :
799 7709 : const unsigned int ns = this->n_sides();
800 7709 : const unsigned int nn = this->n_nodes();
801 :
802 7709 : const unsigned int en = e->n_nodes();
803 :
804 : // e might be on any side until proven otherwise
805 8417 : std::vector<bool> might_be_side(ns, true);
806 :
807 30670 : for (unsigned int i=0; i != en; ++i)
808 : {
809 25069 : Point side_point = e->point(i);
810 2108 : unsigned int local_node_id = libMesh::invalid_uint;
811 :
812 : // Look for a node of this that's contiguous with node i of
813 : // e. Note that the exact floating point comparison of Point
814 : // positions is intentional, see the class documentation for
815 : // this function.
816 229112 : for (unsigned int j=0; j != nn; ++j)
817 37848 : if (this->point(j) == side_point)
818 2108 : local_node_id = j;
819 :
820 : // If a node of e isn't contiguous with some node of this, then
821 : // e isn't a side of this.
822 22961 : if (local_node_id == libMesh::invalid_uint)
823 0 : return libMesh::invalid_uint;
824 :
825 : // If a node of e isn't contiguous with some node on side s of
826 : // this, then e isn't on side s.
827 114639 : for (unsigned int s=0; s != ns; ++s)
828 91678 : if (!this->is_node_on_side(local_node_id, s))
829 9816 : might_be_side[s] = false;
830 : }
831 :
832 21383 : for (unsigned int s=0; s != ns; ++s)
833 23314 : if (might_be_side[s])
834 : {
835 : #ifdef DEBUG
836 1593 : for (unsigned int s2=s+1; s2 < ns; ++s2)
837 885 : libmesh_assert (!might_be_side[s2]);
838 : #endif
839 7709 : return s;
840 : }
841 :
842 : // Didn't find any matching side
843 0 : return libMesh::invalid_uint;
844 : }
845 :
846 :
847 :
848 1376921128 : bool Elem::contains_vertex_of(const Elem * e, bool mesh_connection) const
849 : {
850 : // Our vertices are the first numbered nodes
851 1376921128 : const unsigned int nv = e->n_vertices();
852 1376921128 : const unsigned int my_nv = this->n_vertices();
853 :
854 : // Check for vertex-to-vertex containment first; contains_point() is
855 : // expensive
856 4752946196 : for (auto n : make_range(nv))
857 : {
858 17155863 : const Node * vertex = e->node_ptr(n);
859 22780888103 : for (auto my_n : make_range(my_nv))
860 19404863035 : if (&this->node_ref(my_n) == vertex)
861 1312517 : return true;
862 : }
863 :
864 : // If e is in our mesh, then we might be done testing
865 460109670 : if (mesh_connection)
866 : {
867 460109670 : const unsigned int l = this->level();
868 460109670 : const unsigned int el = e->level();
869 :
870 460109670 : if (l >= el)
871 805530 : return false;
872 :
873 : // We could also return false for l==el-1 iff we knew we had no
874 : // triangular faces, but we don't have an API to check that.
875 : }
876 :
877 : // Our vertices are the first numbered nodes
878 48578624 : for (auto n : make_range(nv))
879 42535333 : if (this->contains_point(e->point(n)))
880 0 : return true;
881 1632 : return false;
882 : }
883 :
884 :
885 :
886 0 : bool Elem::contains_edge_of(const Elem * e) const
887 : {
888 0 : unsigned int num_contained_edges = 0;
889 :
890 : // Our vertices are the first numbered nodes
891 0 : for (auto n : make_range(e->n_vertices()))
892 : {
893 0 : if (this->contains_point(e->point(n)))
894 : {
895 0 : num_contained_edges++;
896 0 : if (num_contained_edges>=2)
897 : {
898 0 : return true;
899 : }
900 : }
901 : }
902 0 : return false;
903 : }
904 :
905 :
906 :
907 575865 : void Elem::find_point_neighbors(const Point & p,
908 : std::set<const Elem *> & neighbor_set) const
909 : {
910 51170 : libmesh_assert(this->contains_point(p));
911 51170 : libmesh_assert(this->active());
912 :
913 51170 : neighbor_set.clear();
914 575865 : neighbor_set.insert(this);
915 :
916 102340 : std::set<const Elem *> untested_set, next_untested_set;
917 575865 : untested_set.insert(this);
918 :
919 : #ifdef LIBMESH_ENABLE_AMR
920 102340 : std::vector<const Elem *> active_neighbor_children;
921 : #endif // #ifdef LIBMESH_ENABLE_AMR
922 :
923 1738094 : while (!untested_set.empty())
924 : {
925 : // Loop over all the elements in the patch that haven't already
926 : // been tested
927 2473563 : for (const auto & elem : untested_set)
928 6921039 : for (auto current_neighbor : elem->neighbor_ptr_range())
929 : {
930 5922858 : if (current_neighbor &&
931 5044197 : current_neighbor != remote_elem) // we have a real neighbor on this side
932 : {
933 426766 : if (current_neighbor->active()) // ... if it is active
934 : {
935 422224 : auto it = neighbor_set.lower_bound(current_neighbor);
936 8811853 : if ((it == neighbor_set.end() || *it != current_neighbor) &&
937 3889688 : current_neighbor->contains_point(p)) // ... don't have and touches p
938 : {
939 : // Add it and test it
940 671790 : next_untested_set.insert(current_neighbor);
941 733691 : neighbor_set.emplace_hint(it, current_neighbor);
942 : }
943 : }
944 : #ifdef LIBMESH_ENABLE_AMR
945 : else // ... the neighbor is *not* active,
946 : { // ... so add *all* neighboring
947 : // active children that touch p
948 4542 : active_neighbor_children.clear();
949 : current_neighbor->active_family_tree_by_neighbor
950 14099 : (active_neighbor_children, elem);
951 :
952 42611 : for (const auto & current_child : active_neighbor_children)
953 : {
954 9084 : auto it = neighbor_set.lower_bound(current_child);
955 55778 : if ((it == neighbor_set.end() || *it != current_child) &&
956 27266 : current_child->contains_point(p))
957 : {
958 : // Add it and test it
959 1236 : next_untested_set.insert(current_child);
960 1778 : neighbor_set.emplace_hint(it, current_child);
961 : }
962 : }
963 : }
964 : #endif // #ifdef LIBMESH_ENABLE_AMR
965 : }
966 : }
967 101794 : untested_set.swap(next_untested_set);
968 101794 : next_untested_set.clear();
969 : }
970 575865 : }
971 :
972 :
973 :
974 10187738 : void Elem::find_point_neighbors(std::set<const Elem *> & neighbor_set) const
975 : {
976 10187738 : this->find_point_neighbors(neighbor_set, this);
977 10187738 : }
978 :
979 :
980 :
981 10187738 : void Elem::find_point_neighbors(std::set<const Elem *> & neighbor_set,
982 : const Elem * start_elem) const
983 : {
984 10187738 : ElemInternal::find_point_neighbors(this, neighbor_set, start_elem);
985 10187738 : }
986 :
987 :
988 :
989 0 : void Elem::find_point_neighbors(std::set<Elem *> & neighbor_set,
990 : Elem * start_elem)
991 : {
992 0 : ElemInternal::find_point_neighbors(this, neighbor_set, start_elem);
993 0 : }
994 :
995 :
996 :
997 480 : void Elem::find_edge_neighbors(const Point & p1,
998 : const Point & p2,
999 : std::set<const Elem *> & neighbor_set) const
1000 : {
1001 : // Simple but perhaps suboptimal code: find elements containing the
1002 : // first point, then winnow this set down by removing elements which
1003 : // don't also contain the second point
1004 :
1005 40 : libmesh_assert(this->contains_point(p2));
1006 480 : this->find_point_neighbors(p1, neighbor_set);
1007 :
1008 40 : std::set<const Elem *>::iterator it = neighbor_set.begin();
1009 40 : const std::set<const Elem *>::iterator end = neighbor_set.end();
1010 :
1011 1812 : while (it != end)
1012 : {
1013 : // As of C++11, set::erase returns an iterator to the element
1014 : // following the erased element, or end.
1015 1332 : if (!(*it)->contains_point(p2))
1016 781 : it = neighbor_set.erase(it);
1017 : else
1018 40 : ++it;
1019 : }
1020 480 : }
1021 :
1022 :
1023 :
1024 0 : void Elem::find_edge_neighbors(std::set<const Elem *> & neighbor_set) const
1025 : {
1026 0 : neighbor_set.clear();
1027 0 : neighbor_set.insert(this);
1028 :
1029 0 : std::set<const Elem *> untested_set, next_untested_set;
1030 0 : untested_set.insert(this);
1031 :
1032 0 : while (!untested_set.empty())
1033 : {
1034 : // Loop over all the elements in the patch that haven't already
1035 : // been tested
1036 0 : for (const auto & elem : untested_set)
1037 : {
1038 0 : for (auto current_neighbor : elem->neighbor_ptr_range())
1039 : {
1040 0 : if (current_neighbor &&
1041 0 : current_neighbor != remote_elem) // we have a real neighbor on this side
1042 : {
1043 0 : if (current_neighbor->active()) // ... if it is active
1044 : {
1045 0 : if (this->contains_edge_of(current_neighbor) // ... and touches us
1046 0 : || current_neighbor->contains_edge_of(this))
1047 : {
1048 : // Make sure we'll test it
1049 0 : if (!neighbor_set.count(current_neighbor))
1050 0 : next_untested_set.insert (current_neighbor);
1051 :
1052 : // And add it
1053 0 : neighbor_set.insert (current_neighbor);
1054 : }
1055 : }
1056 : #ifdef LIBMESH_ENABLE_AMR
1057 : else // ... the neighbor is *not* active,
1058 : { // ... so add *all* neighboring
1059 : // active children
1060 0 : std::vector<const Elem *> active_neighbor_children;
1061 :
1062 : current_neighbor->active_family_tree_by_neighbor
1063 0 : (active_neighbor_children, elem);
1064 :
1065 0 : for (const auto & current_child : active_neighbor_children)
1066 0 : if (this->contains_edge_of(current_child) || current_child->contains_edge_of(this))
1067 : {
1068 : // Make sure we'll test it
1069 0 : if (!neighbor_set.count(current_child))
1070 0 : next_untested_set.insert (current_child);
1071 :
1072 0 : neighbor_set.insert (current_child);
1073 : }
1074 : }
1075 : #endif // #ifdef LIBMESH_ENABLE_AMR
1076 : }
1077 : }
1078 : }
1079 0 : untested_set.swap(next_untested_set);
1080 0 : next_untested_set.clear();
1081 : }
1082 0 : }
1083 :
1084 :
1085 :
1086 5376 : void Elem::find_interior_neighbors(std::set<const Elem *> & neighbor_set) const
1087 : {
1088 5376 : ElemInternal::find_interior_neighbors(this, neighbor_set);
1089 5376 : }
1090 :
1091 :
1092 :
1093 6253 : void Elem::find_interior_neighbors(std::set<Elem *> & neighbor_set)
1094 : {
1095 6253 : ElemInternal::find_interior_neighbors(this, neighbor_set);
1096 6253 : }
1097 :
1098 :
1099 :
1100 86870497 : const Elem * Elem::interior_parent () const
1101 : {
1102 : // interior parents make no sense for full-dimensional elements.
1103 86870497 : if (this->dim() >= LIBMESH_DIM)
1104 837377 : return nullptr;
1105 :
1106 : // they USED TO BE only good for level-0 elements, but we now
1107 : // support keeping interior_parent() valid on refined boundary
1108 : // elements.
1109 : // if (this->level() != 0)
1110 : // return this->parent()->interior_parent();
1111 :
1112 : // We store the interior_parent pointer after both the parent
1113 : // neighbor and neighbor pointers
1114 47107132 : Elem * interior_p = _elemlinks[1+this->n_sides()];
1115 :
1116 : // If we have an interior_parent, we USED TO assume it was a
1117 : // one-higher-dimensional interior element, but we now allow e.g.
1118 : // edge elements to have a 3D interior_parent with no
1119 : // intermediate 2D element.
1120 : // libmesh_assert (!interior_p ||
1121 : // interior_p->dim() == (this->dim()+1));
1122 2466729 : libmesh_assert (!interior_p ||
1123 : (interior_p == remote_elem) ||
1124 : (interior_p->dim() > this->dim()));
1125 :
1126 : // If an element in a multi-dimensional mesh has an interior_parent
1127 : // link, it should be at our level or coarser, just like a neighbor
1128 : // link. Our collect_families() code relies on this, but it might
1129 : // be tempting for users to manually assign something that breaks
1130 : // it.
1131 : //
1132 : // However, we *also* create temporary side elements, and we don't
1133 : // bother with creating ancestors for those, so they can be at level
1134 : // 0 even when they're sides of non-level-0 elements.
1135 2466729 : libmesh_assert (!interior_p ||
1136 : (interior_p->level() <= this->level()) ||
1137 : (this->level() == 0 &&
1138 : this->id() == DofObject::invalid_id));
1139 :
1140 47107132 : return interior_p;
1141 : }
1142 :
1143 :
1144 :
1145 101288689 : Elem * Elem::interior_parent ()
1146 : {
1147 : // See the const version for comments
1148 101288689 : if (this->dim() >= LIBMESH_DIM)
1149 520192 : return nullptr;
1150 :
1151 77623403 : Elem * interior_p = _elemlinks[1+this->n_sides()];
1152 :
1153 3179363 : libmesh_assert (!interior_p ||
1154 : (interior_p == remote_elem) ||
1155 : (interior_p->dim() > this->dim()));
1156 :
1157 77623403 : return interior_p;
1158 : }
1159 :
1160 :
1161 :
1162 4480174734 : void Elem::set_interior_parent (Elem * p)
1163 : {
1164 : // interior parents make no sense for full-dimensional elements.
1165 428367584 : libmesh_assert (!p ||
1166 : this->dim() < LIBMESH_DIM);
1167 :
1168 : // If we have an interior_parent, we USED TO assume it was a
1169 : // one-higher-dimensional interior element, but we now allow e.g.
1170 : // edge elements to have a 3D interior_parent with no
1171 : // intermediate 2D element.
1172 : // libmesh_assert (!p ||
1173 : // p->dim() == (this->dim()+1));
1174 428367584 : libmesh_assert (!p ||
1175 : (p == remote_elem) ||
1176 : (p->dim() > this->dim()));
1177 :
1178 4480174734 : _elemlinks[1+this->n_sides()] = p;
1179 4480174734 : }
1180 :
1181 :
1182 :
1183 : #ifdef LIBMESH_ENABLE_PERIODIC
1184 :
1185 786435 : Elem * Elem::topological_neighbor (const unsigned int i,
1186 : MeshBase & mesh,
1187 : const PointLocatorBase & point_locator,
1188 : const PeriodicBoundaries * pb)
1189 : {
1190 488472 : libmesh_assert_less (i, this->n_neighbors());
1191 :
1192 687080 : Elem * neighbor_i = this->neighbor_ptr(i);
1193 786435 : if (neighbor_i != nullptr)
1194 476930 : return neighbor_i;
1195 :
1196 22252 : if (pb)
1197 : {
1198 : // Since the neighbor is nullptr it must be on a boundary. We need
1199 : // see if this is a periodic boundary in which case it will have a
1200 : // topological neighbor
1201 11542 : std::vector<boundary_id_type> bc_ids;
1202 22252 : mesh.get_boundary_info().boundary_ids(this, cast_int<unsigned short>(i), bc_ids);
1203 22252 : for (const auto & id : bc_ids)
1204 22252 : if (pb->boundary(id))
1205 : {
1206 : // Since the point locator inside of periodic boundaries
1207 : // returns a const pointer we will retrieve the proper
1208 : // pointer directly from the mesh object.
1209 22252 : const Elem * const cn = pb->neighbor(id, point_locator, this, i);
1210 11542 : neighbor_i = const_cast<Elem *>(cn);
1211 :
1212 : // Since coarse elements do not have more refined
1213 : // neighbors we need to make sure that we don't return one
1214 : // of these types of neighbors.
1215 22252 : if (neighbor_i)
1216 24573 : while (level() < neighbor_i->level())
1217 2026 : neighbor_i = neighbor_i->parent();
1218 11542 : return neighbor_i;
1219 : }
1220 : }
1221 :
1222 0 : return nullptr;
1223 : }
1224 :
1225 :
1226 :
1227 47450 : const Elem * Elem::topological_neighbor (const unsigned int i,
1228 : const MeshBase & mesh,
1229 : const PointLocatorBase & point_locator,
1230 : const PeriodicBoundaries * pb) const
1231 : {
1232 19968 : libmesh_assert_less (i, this->n_neighbors());
1233 :
1234 31168 : const Elem * neighbor_i = this->neighbor_ptr(i);
1235 47450 : if (neighbor_i != nullptr)
1236 11098 : return neighbor_i;
1237 :
1238 31230 : if (pb)
1239 : {
1240 : // Since the neighbor is nullptr it must be on a boundary. We need
1241 : // see if this is a periodic boundary in which case it will have a
1242 : // topological neighbor
1243 8870 : std::vector<boundary_id_type> bc_ids;
1244 31230 : mesh.get_boundary_info().boundary_ids(this, cast_int<unsigned short>(i), bc_ids);
1245 34152 : for (const auto & id : bc_ids)
1246 21754 : if (pb->boundary(id))
1247 : {
1248 18832 : neighbor_i = pb->neighbor(id, point_locator, this, i);
1249 :
1250 : // Since coarse elements do not have more refined
1251 : // neighbors we need to make sure that we don't return one
1252 : // of these types of neighbors.
1253 18832 : if (neighbor_i)
1254 20819 : while (level() < neighbor_i->level())
1255 1836 : neighbor_i = neighbor_i->parent();
1256 8134 : return neighbor_i;
1257 : }
1258 : }
1259 :
1260 736 : return nullptr;
1261 : }
1262 :
1263 :
1264 20006 : bool Elem::has_topological_neighbor (const Elem * elem,
1265 : const MeshBase & mesh,
1266 : const PointLocatorBase & point_locator,
1267 : const PeriodicBoundaries * pb) const
1268 : {
1269 : // First see if this is a normal "interior" neighbor
1270 20006 : if (has_neighbor(elem))
1271 0 : return true;
1272 :
1273 20056 : for (auto n : this->side_index_range())
1274 20056 : if (this->topological_neighbor(n, mesh, point_locator, pb))
1275 13942 : return true;
1276 :
1277 0 : return false;
1278 : }
1279 :
1280 :
1281 : #endif
1282 :
1283 : #ifndef NDEBUG
1284 :
1285 0 : void Elem::libmesh_assert_valid_node_pointers() const
1286 : {
1287 0 : libmesh_assert(this->valid_id());
1288 0 : for (auto n : this->node_index_range())
1289 : {
1290 0 : libmesh_assert(this->node_ptr(n));
1291 0 : libmesh_assert(this->node_ptr(n)->valid_id());
1292 : }
1293 0 : }
1294 :
1295 :
1296 :
1297 2867360 : void Elem::libmesh_assert_valid_neighbors() const
1298 : {
1299 14127190 : for (auto n : this->side_index_range())
1300 : {
1301 11259830 : const Elem * neigh = this->neighbor_ptr(n);
1302 :
1303 : // Any element might have a remote neighbor; checking
1304 : // to make sure that's not inaccurate is tough.
1305 11259830 : if (neigh == remote_elem)
1306 10136 : continue;
1307 :
1308 11249694 : if (neigh)
1309 : {
1310 : // Only subactive elements have subactive neighbors
1311 10621588 : libmesh_assert (this->subactive() || !neigh->subactive());
1312 :
1313 10621588 : const Elem * elem = this;
1314 :
1315 : // If we're subactive but our neighbor isn't, its
1316 : // return neighbor link will be to our first active
1317 : // ancestor OR to our inactive ancestor of the same
1318 : // level as neigh,
1319 10621588 : if (this->subactive() && !neigh->subactive())
1320 : {
1321 60504 : for (elem = this; !elem->active();
1322 30252 : elem = elem->parent())
1323 30252 : libmesh_assert(elem);
1324 : }
1325 : else
1326 : {
1327 10591336 : unsigned int rev = neigh->which_neighbor_am_i(elem);
1328 10591336 : libmesh_assert_less (rev, neigh->n_neighbors());
1329 :
1330 10591336 : if (this->subactive() && !neigh->subactive())
1331 : {
1332 0 : while (neigh->neighbor_ptr(rev) != elem)
1333 : {
1334 0 : libmesh_assert(elem->parent());
1335 0 : elem = elem->parent();
1336 : }
1337 : }
1338 : else
1339 : {
1340 10591336 : const Elem * nn = neigh->neighbor_ptr(rev);
1341 10591336 : libmesh_assert(nn);
1342 :
1343 10920376 : for (; elem != nn; elem = elem->parent())
1344 329040 : libmesh_assert(elem);
1345 : }
1346 : }
1347 : }
1348 : // If we don't have a neighbor and we're not subactive, our
1349 : // ancestors shouldn't have any neighbors in this same
1350 : // direction.
1351 628106 : else if (!this->subactive())
1352 : {
1353 616862 : const Elem * my_parent = this->parent();
1354 883468 : if (my_parent &&
1355 : // A parent with a different dimension isn't really one of
1356 : // our ancestors, it means we're on a boundary mesh and this
1357 : // is an interior mesh element for which we're on a side.
1358 : // Nothing to test for in that case.
1359 266606 : (my_parent->dim() == this->dim()))
1360 266606 : libmesh_assert (!my_parent->neighbor_ptr(n));
1361 : }
1362 : }
1363 2867360 : }
1364 :
1365 : #endif // !NDEBUG
1366 :
1367 :
1368 :
1369 230861545 : void Elem::make_links_to_me_local(unsigned int n, unsigned int nn)
1370 : {
1371 77362 : Elem * neigh = this->neighbor_ptr(n);
1372 :
1373 : // Don't bother calling this function unless it's necessary
1374 38681 : libmesh_assert(neigh);
1375 38681 : libmesh_assert(!neigh->is_remote());
1376 :
1377 230861545 : const unsigned int this_level = this->level();
1378 230861545 : const unsigned int neigh_level = neigh->level();
1379 :
1380 : // We never have neighbors more refined than us
1381 38681 : libmesh_assert_less_equal (neigh_level, this_level);
1382 :
1383 : // We never have subactive neighbors of non subactive elements
1384 38681 : libmesh_assert(!neigh->subactive() || this->subactive());
1385 :
1386 : // If we have a neighbor less refined than us then it must not
1387 : // have any more refined descendants we could have pointed to
1388 : // instead.
1389 38681 : libmesh_assert((neigh_level == this_level) ||
1390 : (neigh->active() && !this->subactive()) ||
1391 : (!neigh->has_children() && this->subactive()));
1392 :
1393 : // What side of neigh are we on? nn.
1394 : //
1395 : // We can't use the usual Elem method because we're in the middle of
1396 : // restoring topology. We can't compare side_ptr nodes because
1397 : // users want to abuse neighbor_ptr to point to
1398 : // not-technically-neighbors across mesh slits. We can't compare
1399 : // node locations because users want to move those
1400 : // not-technically-neighbors until they're
1401 : // not-even-geometrically-neighbors.
1402 :
1403 : // Find any elements that might need to point to elem
1404 77362 : std::vector<Elem *> neigh_family;
1405 : #ifdef LIBMESH_ENABLE_AMR
1406 230861545 : neigh->total_family_tree_by_side(neigh_family, nn);
1407 : #else
1408 : neigh_family.push_back(neigh);
1409 : #endif
1410 :
1411 : // Pull objects out of the loop to reduce heap operations
1412 230861545 : std::unique_ptr<Elem> my_side, neigh_side;
1413 :
1414 : // And point them to elem
1415 525821636 : for (auto & neigh_family_member : neigh_family)
1416 : {
1417 : // Only subactive elements point to other subactive elements
1418 294960091 : const bool member_subactive = neigh_family_member->subactive();
1419 294960091 : if (this->subactive() && !member_subactive)
1420 14232 : continue;
1421 :
1422 : // neigh (and possibly some of its family) might be at a lower
1423 : // level than us, with a neighbor at that lower level. We
1424 : // would exit early if neigh was at a lower level, except we do
1425 : // want to handle neighbor links from subactive descendents.
1426 294958239 : const unsigned int member_level = neigh_family_member->level();
1427 294958239 : if (member_level < this_level)
1428 1968026 : continue;
1429 :
1430 : // Ideally, the neighbor link ought to either be correct
1431 : // already or ought to be to remote_elem.
1432 : //
1433 : // However, if we're redistributing a newly created elem,
1434 : // after an AMR step but before find_neighbors has fixed up
1435 : // neighbor links, we might have an out of date neighbor
1436 : // link to elem's parent instead.
1437 : #ifdef LIBMESH_ENABLE_AMR
1438 51061 : libmesh_assert((neigh_family_member->neighbor_ptr(nn) &&
1439 : (neigh_family_member->neighbor_ptr(nn)->active() ||
1440 : this->is_ancestor_of(neigh_family_member->neighbor_ptr(nn)) ||
1441 : neigh_family_member->neighbor_ptr(nn)->is_ancestor_of(this))) ||
1442 : (neigh_family_member->neighbor_ptr(nn) == remote_elem) ||
1443 : ((this->refinement_flag() == JUST_REFINED) &&
1444 : (this->parent() != nullptr) &&
1445 : (neigh_family_member->neighbor_ptr(nn) == this->parent())));
1446 : #else
1447 : libmesh_assert((neigh_family_member->neighbor_ptr(nn) == this) ||
1448 : (neigh_family_member->neighbor_ptr(nn) == remote_elem));
1449 : #endif
1450 :
1451 : // Some of neigh's family might be at a higher (finer) level
1452 : // than us, and might need to point back to one of our children
1453 : // rather than to us.
1454 292990213 : if (member_level > this_level)
1455 : {
1456 64073502 : if (this->ancestor())
1457 63238571 : continue;
1458 :
1459 822551 : if (member_subactive && this->has_children())
1460 15036 : continue;
1461 : }
1462 :
1463 : // If neigh is at a coarser level than us, some of neigh's
1464 : // subactive family just won't line up with us!
1465 229724226 : if (neigh_level < this_level &&
1466 0 : member_level > neigh_level)
1467 : {
1468 0 : libmesh_assert(member_subactive);
1469 :
1470 25044 : this->side_ptr(my_side, n);
1471 25044 : neigh_family_member->side_ptr(neigh_side, nn);
1472 :
1473 25044 : if (*my_side != *neigh_side)
1474 16996 : continue;
1475 : }
1476 :
1477 229707230 : neigh_family_member->set_neighbor(nn, this);
1478 : }
1479 461645728 : }
1480 :
1481 :
1482 36156743 : void Elem::make_links_to_me_remote()
1483 : {
1484 1801 : libmesh_assert_not_equal_to (this, remote_elem);
1485 :
1486 : // We need to have handled any children first
1487 : #if defined(LIBMESH_ENABLE_AMR) && defined(DEBUG)
1488 1801 : if (this->has_children())
1489 268 : for (auto & child : this->child_ref_range())
1490 208 : libmesh_assert_equal_to (&child, remote_elem);
1491 : #endif
1492 :
1493 : // Remotify any neighbor links
1494 183281867 : for (auto neigh : this->neighbor_ptr_range())
1495 : {
1496 147125124 : if (neigh && neigh != remote_elem)
1497 : {
1498 : // My neighbor should never be more refined than me; my real
1499 : // neighbor would have been its parent in that case.
1500 3916 : libmesh_assert_greater_equal (this->level(), neigh->level());
1501 :
1502 141875551 : if (this->level() == neigh->level() &&
1503 70596529 : neigh->has_neighbor(this))
1504 : {
1505 : #ifdef LIBMESH_ENABLE_AMR
1506 : // My neighbor may have descendants which also consider me a
1507 : // neighbor
1508 7832 : std::vector<Elem *> family;
1509 70599818 : neigh->total_family_tree_by_neighbor (family, this);
1510 :
1511 : // FIXME - There's a lot of ugly const_casts here; we
1512 : // may want to make remote_elem non-const
1513 141207431 : for (auto & n : family)
1514 : {
1515 3916 : libmesh_assert (n);
1516 70607613 : if (n == remote_elem)
1517 0 : continue;
1518 70607613 : unsigned int my_s = n->which_neighbor_am_i(this);
1519 3916 : libmesh_assert_less (my_s, n->n_neighbors());
1520 3916 : libmesh_assert_equal_to (n->neighbor_ptr(my_s), this);
1521 70607613 : n->set_neighbor(my_s, const_cast<RemoteElem *>(remote_elem));
1522 : }
1523 : #else
1524 : unsigned int my_s = neigh->which_neighbor_am_i(this);
1525 : libmesh_assert_less (my_s, neigh->n_neighbors());
1526 : libmesh_assert_equal_to (neigh->neighbor_ptr(my_s), this);
1527 : neigh->set_neighbor(my_s, const_cast<RemoteElem *>(remote_elem));
1528 : #endif
1529 : }
1530 : #ifdef LIBMESH_ENABLE_AMR
1531 : // Even if my neighbor doesn't link back to me, it might
1532 : // have subactive descendants which do
1533 0 : else if (neigh->has_children())
1534 : {
1535 : // If my neighbor at the same level doesn't have me as a
1536 : // neighbor, I must be subactive
1537 0 : libmesh_assert(this->level() > neigh->level() ||
1538 : this->subactive());
1539 :
1540 : // My neighbor must have some ancestor of mine as a
1541 : // neighbor
1542 0 : Elem * my_ancestor = this->parent();
1543 0 : libmesh_assert(my_ancestor);
1544 1760 : while (!neigh->has_neighbor(my_ancestor))
1545 : {
1546 0 : my_ancestor = my_ancestor->parent();
1547 0 : libmesh_assert(my_ancestor);
1548 : }
1549 :
1550 : // My neighbor may have descendants which consider me a
1551 : // neighbor
1552 0 : std::vector<Elem *> family;
1553 1760 : neigh->total_family_tree_by_subneighbor (family, my_ancestor, this);
1554 :
1555 2874 : for (auto & n : family)
1556 : {
1557 0 : libmesh_assert (n);
1558 1114 : if (n->is_remote())
1559 0 : continue;
1560 1114 : unsigned int my_s = n->which_neighbor_am_i(this);
1561 0 : libmesh_assert_less (my_s, n->n_neighbors());
1562 0 : libmesh_assert_equal_to (n->neighbor_ptr(my_s), this);
1563 : // TODO: we may want to make remote_elem non-const.
1564 1114 : n->set_neighbor(my_s, const_cast<RemoteElem *>(remote_elem));
1565 : }
1566 : }
1567 : #endif
1568 : }
1569 : }
1570 :
1571 : #ifdef LIBMESH_ENABLE_AMR
1572 : // Remotify parent's child link
1573 3602 : Elem * my_parent = this->parent();
1574 26999795 : if (my_parent &&
1575 : // As long as it's not already remote
1576 63154737 : my_parent != remote_elem &&
1577 : // And it's a real parent, not an interior parent
1578 26997994 : this->dim() == my_parent->dim())
1579 : {
1580 26997994 : unsigned int me = my_parent->which_child_am_i(this);
1581 672 : libmesh_assert_equal_to (my_parent->child_ptr(me), this);
1582 26997994 : my_parent->set_child(me, const_cast<RemoteElem *>(remote_elem));
1583 : }
1584 : #endif
1585 36156743 : }
1586 :
1587 :
1588 813 : void Elem::remove_links_to_me()
1589 : {
1590 31 : libmesh_assert_not_equal_to (this, remote_elem);
1591 :
1592 : // We need to have handled any children first
1593 : #ifdef LIBMESH_ENABLE_AMR
1594 31 : libmesh_assert (!this->has_children());
1595 : #endif
1596 :
1597 : // Nullify any neighbor links
1598 4065 : for (auto neigh : this->neighbor_ptr_range())
1599 : {
1600 3252 : if (neigh && neigh != remote_elem)
1601 : {
1602 : // My neighbor should never be more refined than me; my real
1603 : // neighbor would have been its parent in that case.
1604 47 : libmesh_assert_greater_equal (this->level(), neigh->level());
1605 :
1606 2504 : if (this->level() == neigh->level() &&
1607 1205 : neigh->has_neighbor(this))
1608 : {
1609 : #ifdef LIBMESH_ENABLE_AMR
1610 : // My neighbor may have descendants which also consider me a
1611 : // neighbor
1612 94 : std::vector<Elem *> family;
1613 1252 : neigh->total_family_tree_by_neighbor (family, this);
1614 :
1615 2504 : for (auto & n : family)
1616 : {
1617 47 : libmesh_assert (n);
1618 1252 : if (n->is_remote())
1619 0 : continue;
1620 1252 : unsigned int my_s = n->which_neighbor_am_i(this);
1621 47 : libmesh_assert_less (my_s, n->n_neighbors());
1622 47 : libmesh_assert_equal_to (n->neighbor_ptr(my_s), this);
1623 1252 : n->set_neighbor(my_s, nullptr);
1624 : }
1625 : #else
1626 : unsigned int my_s = neigh->which_neighbor_am_i(this);
1627 : libmesh_assert_less (my_s, neigh->n_neighbors());
1628 : libmesh_assert_equal_to (neigh->neighbor_ptr(my_s), this);
1629 : neigh->set_neighbor(my_s, nullptr);
1630 : #endif
1631 : }
1632 : #ifdef LIBMESH_ENABLE_AMR
1633 : // Even if my neighbor doesn't link back to me, it might
1634 : // have subactive descendants which do
1635 0 : else if (neigh->has_children())
1636 : {
1637 : // If my neighbor at the same level doesn't have me as a
1638 : // neighbor, I must be subactive
1639 0 : libmesh_assert(this->level() > neigh->level() ||
1640 : this->subactive());
1641 :
1642 : // My neighbor must have some ancestor of mine as a
1643 : // neighbor
1644 0 : Elem * my_ancestor = this->parent();
1645 0 : libmesh_assert(my_ancestor);
1646 0 : while (!neigh->has_neighbor(my_ancestor))
1647 : {
1648 0 : my_ancestor = my_ancestor->parent();
1649 0 : libmesh_assert(my_ancestor);
1650 : }
1651 :
1652 : // My neighbor may have descendants which consider me a
1653 : // neighbor
1654 0 : std::vector<Elem *> family;
1655 0 : neigh->total_family_tree_by_subneighbor (family, my_ancestor, this);
1656 :
1657 0 : for (auto & n : family)
1658 : {
1659 0 : libmesh_assert (n);
1660 0 : if (n->is_remote())
1661 0 : continue;
1662 0 : unsigned int my_s = n->which_neighbor_am_i(this);
1663 0 : libmesh_assert_less (my_s, n->n_neighbors());
1664 0 : libmesh_assert_equal_to (n->neighbor_ptr(my_s), this);
1665 0 : n->set_neighbor(my_s, nullptr);
1666 : }
1667 : }
1668 : #endif
1669 : }
1670 : }
1671 :
1672 : #ifdef LIBMESH_ENABLE_AMR
1673 : // We can't currently delete a child with a parent!
1674 31 : libmesh_assert (!this->parent());
1675 : #endif
1676 813 : }
1677 :
1678 :
1679 :
1680 178350 : void Elem::write_connectivity (std::ostream & out_stream,
1681 : const IOPackage iop) const
1682 : {
1683 16310 : libmesh_assert (out_stream.good());
1684 16310 : libmesh_assert(_nodes);
1685 16310 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
1686 :
1687 178350 : switch (iop)
1688 : {
1689 16185 : case TECPLOT:
1690 : {
1691 : // This connectivity vector will be used repeatedly instead
1692 : // of being reconstructed inside the loop.
1693 32370 : std::vector<dof_id_type> conn;
1694 498180 : for (auto sc : make_range(this->n_sub_elem()))
1695 : {
1696 176850 : this->connectivity(sc, TECPLOT, conn);
1697 :
1698 176850 : std::copy(conn.begin(),
1699 : conn.end(),
1700 16185 : std::ostream_iterator<dof_id_type>(out_stream, " "));
1701 :
1702 176850 : out_stream << '\n';
1703 : }
1704 16185 : return;
1705 : }
1706 :
1707 125 : case UCD:
1708 : {
1709 13500 : for (auto i : this->node_index_range())
1710 13000 : out_stream << this->node_id(i)+1 << "\t";
1711 :
1712 1500 : out_stream << '\n';
1713 1500 : return;
1714 : }
1715 :
1716 0 : default:
1717 0 : libmesh_error_msg("Unsupported IO package " << iop);
1718 : }
1719 : }
1720 :
1721 :
1722 :
1723 1417981 : Real Elem::quality (const ElemQuality q) const
1724 : {
1725 1417981 : switch (q)
1726 : {
1727 : // Return the maximum ratio of edge lengths, or zero if that
1728 : // maximum would otherwise be infinity.
1729 33297 : case EDGE_LENGTH_RATIO:
1730 : {
1731 33297 : if (this->dim() < 2)
1732 7 : return 1;
1733 :
1734 35986 : std::vector<Real> edge_lengths(this->n_edges());
1735 274215 : for (auto e : index_range(edge_lengths))
1736 461889 : edge_lengths[e] = (this->point(this->local_edge_node(e,1)) -
1737 461889 : this->point(this->local_edge_node(e,0))).norm();
1738 :
1739 2773 : auto [min, max] =
1740 33213 : std::minmax_element(edge_lengths.begin(), edge_lengths.end());
1741 :
1742 33213 : if (*min == 0.)
1743 0 : return 0.;
1744 : else
1745 33213 : return *max / *min;
1746 : }
1747 :
1748 267652 : case MIN_ANGLE:
1749 : case MAX_ANGLE:
1750 : {
1751 : // 1D elements don't have interior angles, so just return some
1752 : // dummy value in that case.
1753 267652 : if (this->dim() < 2)
1754 14 : return 0.;
1755 :
1756 : // Initialize return values
1757 267484 : Real min_angle = std::numeric_limits<Real>::max();
1758 267484 : Real max_angle = -std::numeric_limits<Real>::max();
1759 :
1760 3031040 : for (auto n : this->node_index_range())
1761 : {
1762 : // Get list of edge ids adjacent to this node.
1763 2763556 : const auto adjacent_edge_ids = this->edges_adjacent_to_node(n);
1764 :
1765 : // Skip any nodes with fewer than 2 adjacent edges. You
1766 : // need at least two adjacent edges to form an interior
1767 : // element angle.
1768 460364 : auto N = adjacent_edge_ids.size();
1769 2763556 : if (N < 2)
1770 131946 : continue;
1771 :
1772 : // Consider all possible pairs of edges adjacent to node n
1773 3515652 : for (unsigned int first = 0; first < N-1; ++first)
1774 5850332 : for (unsigned int second = first+1; second < N; ++second)
1775 : {
1776 : // Get ids of first and second edges
1777 3808532 : auto first_edge = adjacent_edge_ids[first];
1778 3515748 : auto second_edge = adjacent_edge_ids[second];
1779 :
1780 : // Get node ids of first and second edge
1781 3515748 : auto first_edge_node_0 = this->local_edge_node(first_edge, 0);
1782 3515748 : auto first_edge_node_1 = this->local_edge_node(first_edge, 1);
1783 3515748 : auto second_edge_node_0 = this->local_edge_node(second_edge, 0);
1784 3515748 : auto second_edge_node_1 = this->local_edge_node(second_edge, 1);
1785 :
1786 : // Orient both edges so that edge node 0 == n
1787 3515748 : if (first_edge_node_0 != n)
1788 201378 : std::swap(first_edge_node_0, first_edge_node_1);
1789 3515748 : if (second_edge_node_0 != n)
1790 95984 : std::swap(second_edge_node_0, second_edge_node_1);
1791 :
1792 292784 : libmesh_assert_equal_to(first_edge_node_0, n);
1793 292784 : libmesh_assert_equal_to(second_edge_node_0, n);
1794 :
1795 : // Locally oriented edge vectors
1796 : Point
1797 585568 : first_ev = this->point(first_edge_node_1) - this->point(first_edge_node_0),
1798 292784 : second_ev = this->point(second_edge_node_1) - this->point(second_edge_node_0);
1799 :
1800 : // Angle between them in the range [0, pi]
1801 3515748 : Real theta = std::acos(first_ev.unit() * second_ev.unit());
1802 :
1803 : // Track min and max angles seen
1804 3515748 : min_angle = std::min(theta, min_angle);
1805 3515748 : max_angle = std::max(theta, max_angle);
1806 : }
1807 : }
1808 :
1809 : // Return requested extreme value (in degrees)
1810 267484 : return Real(180)/libMesh::pi * ((q == MIN_ANGLE) ? min_angle : max_angle);
1811 : }
1812 :
1813 257764 : case MIN_DIHEDRAL_ANGLE:
1814 : case MAX_DIHEDRAL_ANGLE:
1815 : {
1816 : // For 1D and 2D elements, just call this function with MIN,MAX_ANGLE instead.
1817 257764 : if (this->dim() < 3)
1818 8680 : return this->quality((q == MIN_DIHEDRAL_ANGLE) ? MIN_ANGLE : MAX_ANGLE);
1819 :
1820 : // Initialize return values
1821 251812 : Real min_angle = std::numeric_limits<Real>::max();
1822 251812 : Real max_angle = -std::numeric_limits<Real>::max();
1823 :
1824 : // Algorithm for computing dihedral angles:
1825 : // .) Loop over the edges, use the edge_sides_map to get the
1826 : // two sides adjacent to the edge.
1827 : // .) For each adjacent side, use the side_nodes_map to get the
1828 : // list of three (or more) node ids on that side.
1829 : // .) Use the node ids to compute the (approximate) inward
1830 : // normal for the side. Note: this is approximate since 3D
1831 : // elements with four-sided faces and quadratic tetrahedra
1832 : // do not necessarily have planar faces.
1833 : // .) Compute the dihedral angle for the current edge, compare
1834 : // it to the current min and max dihedral angles.
1835 1954672 : for (auto e : this->edge_index_range())
1836 : {
1837 : // Get list of edge ids adjacent to this node.
1838 1844772 : const auto adjacent_side_ids = this->sides_on_edge(e);
1839 :
1840 : // All 3D elements should have exactly two sides adjacent to each edge.
1841 141912 : libmesh_assert_equal_to(adjacent_side_ids.size(), 2);
1842 :
1843 : // Get lists of node ids on each side
1844 1844772 : const auto side_0_node_ids = this->nodes_on_side(adjacent_side_ids[0]);
1845 1702860 : const auto side_1_node_ids = this->nodes_on_side(adjacent_side_ids[1]);
1846 :
1847 : // All 3D elements have at least three nodes on each side
1848 141912 : libmesh_assert_greater_equal(side_0_node_ids.size(), 3);
1849 141912 : libmesh_assert_greater_equal(side_1_node_ids.size(), 3);
1850 :
1851 : // Construct (approximate) inward normal on each adjacent side
1852 : const auto side_0_normal =
1853 1844772 : (this->point(side_0_node_ids[2]) - this->point(side_0_node_ids[0])).cross
1854 1844772 : (this->point(side_0_node_ids[1]) - this->point(side_0_node_ids[0])).unit();
1855 : const auto side_1_normal =
1856 1844772 : (this->point(side_1_node_ids[2]) - this->point(side_1_node_ids[0])).cross
1857 1844772 : (this->point(side_1_node_ids[1]) - this->point(side_1_node_ids[0])).unit();
1858 :
1859 : // Compute dihedral angle between the planes.
1860 : // Using the absolute value ensures that we get the same result
1861 : // even if the orientation of one of the planes is flipped, i.e.
1862 : // it always gives us a value in the range [0, pi/2].
1863 : // https://en.wikipedia.org/wiki/Dihedral_angle
1864 1702860 : Real theta = std::acos(std::abs(side_0_normal * side_1_normal));
1865 :
1866 : // Track min and max angles seen
1867 1702860 : min_angle = std::min(theta, min_angle);
1868 1702860 : max_angle = std::max(theta, max_angle);
1869 : }
1870 :
1871 : // Return requested extreme value (in degrees)
1872 251812 : return Real(180)/libMesh::pi * ((q == MIN_DIHEDRAL_ANGLE) ? min_angle : max_angle);
1873 : }
1874 :
1875 231202 : case JACOBIAN:
1876 : case SCALED_JACOBIAN:
1877 : {
1878 : // 1D elements don't have interior corners, so this metric
1879 : // does not really apply to them.
1880 231202 : const auto N = this->dim();
1881 231202 : if (N < 2)
1882 14 : return 1.;
1883 :
1884 : // Initialize return value
1885 231034 : Real min_node_area = std::numeric_limits<Real>::max();
1886 :
1887 2520086 : for (auto n : this->node_index_range())
1888 : {
1889 : // Get list of edge ids adjacent to this node.
1890 2289052 : auto adjacent_edge_ids = this->edges_adjacent_to_node(n);
1891 :
1892 : // Skip any nodes that don't have dim() adjacent edges. In 2D,
1893 : // you need at least two adjacent edges to compute the cross
1894 : // product, and in 3D, you need at least three adjacent edges
1895 : // to compute the scalar triple product. The only element
1896 : // type which has an unusual topology in this regard is the
1897 : // Pyramid element in 3D, where 4 edges meet at the apex node.
1898 : // For now, we just skip this node when computing the JACOBIAN
1899 : // metric for Pyramids.
1900 2479786 : if (adjacent_edge_ids.size() != N)
1901 108744 : continue;
1902 :
1903 : // Construct oriented edges
1904 984988 : std::vector<Point> oriented_edges(N);
1905 3910064 : for (auto i : make_range(N))
1906 : {
1907 3168612 : auto node_0 = this->local_edge_node(adjacent_edge_ids[i], 0);
1908 3168612 : auto node_1 = this->local_edge_node(adjacent_edge_ids[i], 1);
1909 2925076 : if (node_0 != n)
1910 120978 : std::swap(node_0, node_1);
1911 3412148 : oriented_edges[i] = this->point(node_1) - this->point(node_0);
1912 : }
1913 :
1914 : // Compute unscaled area (2D) or volume (3D) using the
1915 : // cross product (2D) or scalar triple product (3D) of the
1916 : // oriented edges. We take the absolute value so that we
1917 : // don't have to worry about the sign of the area.
1918 984988 : Real node_area = (N == 2) ?
1919 29888 : cross_norm(oriented_edges[0], oriented_edges[1]) :
1920 955100 : std::abs(triple_product(oriented_edges[0], oriented_edges[1], oriented_edges[2]));
1921 :
1922 : // Divide by (non-zero) edge lengths if computing scaled Jacobian.
1923 : // If any length is zero, then set node_area to zero. This means that
1924 : // e.g. degenerate quadrilaterals will have a zero SCALED_JACOBIAN metric.
1925 984988 : if (q == SCALED_JACOBIAN)
1926 2018392 : for (auto i : make_range(N))
1927 : {
1928 1509770 : Real len_i = oriented_edges[i].norm();
1929 1509770 : node_area = (len_i == 0.) ? 0. : (node_area / len_i);
1930 : }
1931 :
1932 : // Update minimum
1933 984988 : min_node_area = std::min(node_area, min_node_area);
1934 : }
1935 :
1936 231034 : return min_node_area;
1937 : }
1938 :
1939 : // Relative size metric: min over the corner nodes of min(tau,
1940 : // 1/tau), where tau is the ratio of the corner's nodal Jacobian
1941 : // determinant to that of the ideal (regular) element of the same
1942 : // volume -- ReferenceElem::ideal_target(), rescaled to this
1943 : // element's volume via the volume ratio. tau = 1 at every corner
1944 : // of an element whose Jacobian is uniform (any affine element:
1945 : // parallelogram, box, or regular simplex, at any scale), so the
1946 : // metric is 1; non-uniform (tapered/sheared) elements score below
1947 : // 1, and a degenerate corner drives it to 0.
1948 108472 : case SIZE:
1949 : {
1950 : // 1D elements don't have interior corners, so this metric does
1951 : // not really apply to them.
1952 108472 : const auto N = this->dim();
1953 108472 : if (N < 2)
1954 0 : return 1.;
1955 :
1956 108472 : const Real vol = this->volume();
1957 108472 : if (vol == 0.)
1958 0 : return 0.;
1959 :
1960 : // Ideal (regular) element of the same type; its Jacobian is
1961 : // uniform. We compare nodal determinants after rescaling it to
1962 : // this element's volume, i.e. multiply by ideal_vol / this_vol.
1963 126488 : const auto ideal_pair = ReferenceElem::ideal_target(this->type());
1964 9008 : const Elem & ideal = *ideal_pair.first;
1965 108472 : const Real vol_ratio = ideal.volume() / vol;
1966 :
1967 : // Nodal Jacobian determinant at node n of element el (the same
1968 : // construction as the JACOBIAN metric above).
1969 771976 : auto nodal_det = [](const Elem & el, const unsigned int n,
1970 : const std::vector<unsigned int> & edge_ids,
1971 153568 : const unsigned int dim)
1972 : {
1973 925544 : std::vector<Point> e(dim);
1974 3673992 : for (unsigned int i = 0; i != dim; ++i)
1975 : {
1976 2976608 : auto n0 = el.local_edge_node(edge_ids[i], 0);
1977 2976608 : auto n1 = el.local_edge_node(edge_ids[i], 1);
1978 2748448 : if (n0 != n)
1979 114080 : std::swap(n0, n1);
1980 3204768 : e[i] = el.point(n1) - el.point(n0);
1981 : }
1982 925544 : return (dim == 2) ? cross_norm(e[0], e[1])
1983 1825096 : : std::abs(triple_product(e[0], e[1], e[2]));
1984 : };
1985 :
1986 108472 : Real size = 1.;
1987 9008 : bool have_corner = false;
1988 1173260 : for (auto n : this->node_index_range())
1989 : {
1990 : // Skip any nodes that don't have dim() adjacent edges (see
1991 : // the JACOBIAN metric above for the Pyramid apex caveat).
1992 1064788 : const auto adjacent_edge_ids = this->edges_adjacent_to_node(n);
1993 1153348 : if (adjacent_edge_ids.size() != N)
1994 50168 : continue;
1995 38392 : have_corner = true;
1996 :
1997 462772 : const Real a = nodal_det(*this, n, adjacent_edge_ids, N);
1998 462772 : const Real aw = nodal_det(ideal, n, adjacent_edge_ids, N);
1999 :
2000 : // Degenerate corner: worst quality.
2001 462772 : if (a == 0. || aw == 0.)
2002 0 : return 0.;
2003 :
2004 462772 : const Real tau = (a / aw) * vol_ratio;
2005 887152 : size = std::min(size, std::min(tau, Real(1) / tau));
2006 : }
2007 :
2008 108472 : return have_corner ? size : 0.;
2009 : }
2010 :
2011 : // Maximum condition number of the nodal Jacobian over the corner
2012 : // nodes, measured against the ideal (regular) element rather than
2013 : // the reference element. At each corner the physical nodal
2014 : // Jacobian A and the ideal nodal Jacobian W (taken from the same
2015 : // corner of ReferenceElem::ideal_target) give the weighted
2016 : // Jacobian A W^{-1}, whose Frobenius condition number, via the
2017 : // corner metric tensors T_A = A^T A and T_W = W^T W (which also
2018 : // handles a lower-dimensional element embedded in 3D), is
2019 : // kappa = sqrt(tr(T_A T_W^{-1}) * tr(T_W T_A^{-1})) / N.
2020 : // This is 1 for a corner similar to the ideal one -- so an
2021 : // equilateral triangle or regular tetrahedron scores 1, not just
2022 : // a right-angled corner -- and grows with distortion. A
2023 : // degenerate corner has an infinite condition number, reported as
2024 : // 0 (0 stands in for infinity, cf. EDGE_LENGTH_RATIO).
2025 108330 : case CONDITION:
2026 : {
2027 : // 1D elements don't have interior corners, so this metric does
2028 : // not really apply to them.
2029 108330 : const auto N = this->dim();
2030 108330 : if (N < 2)
2031 0 : return 1.;
2032 :
2033 : // Ideal (regular) element of the same type; W is its nodal
2034 : // Jacobian. Its scale is irrelevant here (the condition number
2035 : // is scale invariant), so the reference-volume sizing is fine.
2036 126338 : const auto ideal_pair = ReferenceElem::ideal_target(this->type());
2037 9004 : const Elem & ideal = *ideal_pair.first;
2038 :
2039 : // Corner metric tensor T = A^T A at node n of element el, padded
2040 : // with the identity in unused dimensions so that RealTensor's
2041 : // 3x3 inverse yields the correct NxN inverse.
2042 770368 : auto metric_tensor = [](const Elem & el, const unsigned int n,
2043 : const std::vector<unsigned int> & edge_ids,
2044 153472 : const unsigned int dim)
2045 : {
2046 1000576 : std::vector<Point> e(dim);
2047 3667744 : for (unsigned int i = 0; i != dim; ++i)
2048 : {
2049 2971936 : auto n0 = el.local_edge_node(edge_ids[i], 0);
2050 2971936 : auto n1 = el.local_edge_node(edge_ids[i], 1);
2051 2743904 : if (n0 != n)
2052 114016 : std::swap(n0, n1);
2053 3199968 : e[i] = el.point(n1) - el.point(n0);
2054 : }
2055 847104 : RealTensor T(1, 0, 0, 0, 1, 0, 0, 0, 1);
2056 3667744 : for (unsigned int i = 0; i != dim; ++i)
2057 10920384 : for (unsigned int j = 0; j != dim; ++j)
2058 9535968 : T(i, j) = e[i] * e[j];
2059 1000576 : return T;
2060 : };
2061 :
2062 : // kappa >= 1 for every matrix, so 1 is both the ideal value and
2063 : // a safe floor for the running maximum.
2064 108330 : Real max_cond = 1.;
2065 :
2066 1172266 : for (auto n : this->node_index_range())
2067 : {
2068 : // Skip any nodes that don't have dim() adjacent edges (see
2069 : // the JACOBIAN metric above for the Pyramid apex caveat).
2070 1063936 : const auto adjacent_edge_ids = this->edges_adjacent_to_node(n);
2071 1152472 : if (adjacent_edge_ids.size() != N)
2072 50168 : continue;
2073 :
2074 461920 : const RealTensor Ta = metric_tensor(*this, n, adjacent_edge_ids, N);
2075 :
2076 : // Degenerate corner: infinite condition number.
2077 461920 : if (Ta.det() == 0.)
2078 0 : return 0.;
2079 :
2080 500288 : const RealTensor Tw = metric_tensor(ideal, n, adjacent_edge_ids, N);
2081 500288 : const RealTensor Ta_inv = Ta.inverse();
2082 461920 : const RealTensor Tw_inv = Tw.inverse();
2083 :
2084 : // num1 = tr(T_A T_W^{-1}), num2 = tr(T_W T_A^{-1}) over the
2085 : // NxN blocks (both symmetric, so summed as elementwise dot
2086 : // products).
2087 38368 : Real num1 = 0., num2 = 0.;
2088 1833872 : for (auto i : make_range(N))
2089 5460192 : for (auto j : make_range(N))
2090 : {
2091 4088240 : num1 += Ta(i, j) * Tw_inv(i, j);
2092 4088240 : num2 += Tw(i, j) * Ta_inv(i, j);
2093 : }
2094 :
2095 592964 : max_cond = std::max(max_cond, std::sqrt(num1 * num2) / N);
2096 : }
2097 :
2098 108330 : return max_cond;
2099 : }
2100 :
2101 : // Return 1 if we made it here
2102 411264 : default:
2103 : {
2104 411471 : libmesh_do_once( libmesh_here();
2105 :
2106 : libMesh::err << "ERROR: quality metric "
2107 : << Utility::enum_to_string(q)
2108 : << " not implemented on element type "
2109 : << Utility::enum_to_string(this->type())
2110 : << std::endl
2111 : << "Returning 1."
2112 : << std::endl; );
2113 :
2114 34272 : return 1.;
2115 : }
2116 : }
2117 : }
2118 :
2119 :
2120 :
2121 138187255 : bool Elem::ancestor() const
2122 : {
2123 : #ifdef LIBMESH_ENABLE_AMR
2124 :
2125 : // Use a fast, DistributedMesh-safe definition
2126 : const bool is_ancestor =
2127 86729329 : !this->active() && !this->subactive();
2128 :
2129 : // But check for inconsistencies if we have time
2130 : #ifdef DEBUG
2131 5382328 : if (!is_ancestor && this->has_children())
2132 : {
2133 518656 : for (auto & c : this->child_ref_range())
2134 : {
2135 415616 : if (&c != remote_elem)
2136 : {
2137 415568 : libmesh_assert(!c.active());
2138 415568 : libmesh_assert(!c.ancestor());
2139 : }
2140 : }
2141 : }
2142 : #endif // DEBUG
2143 :
2144 138187255 : return is_ancestor;
2145 :
2146 : #else
2147 : return false;
2148 : #endif
2149 : }
2150 :
2151 :
2152 :
2153 : #ifdef LIBMESH_ENABLE_AMR
2154 :
2155 144494 : void Elem::add_child (Elem * elem)
2156 : {
2157 144494 : const unsigned int nc = this->n_children();
2158 :
2159 144494 : if (!_children)
2160 : {
2161 37737 : _children = std::make_unique<Elem *[]>(nc);
2162 :
2163 180795 : for (unsigned int c = 0; c != nc; c++)
2164 5724 : this->set_child(c, nullptr);
2165 : }
2166 :
2167 360525 : for (unsigned int c = 0; c != nc; c++)
2168 : {
2169 374815 : if (this->_children[c] == nullptr || this->_children[c] == remote_elem)
2170 : {
2171 5724 : libmesh_assert_equal_to (this, elem->parent());
2172 5724 : this->set_child(c, elem);
2173 144494 : return;
2174 : }
2175 : }
2176 :
2177 0 : libmesh_error_msg("Error: Tried to add a child to an element with full children array");
2178 : }
2179 :
2180 :
2181 :
2182 50245814 : void Elem::add_child (Elem * elem, unsigned int c)
2183 : {
2184 49532 : if (!this->has_children())
2185 : {
2186 6189715 : const unsigned int nc = this->n_children();
2187 6203381 : _children = std::make_unique<Elem *[]>(nc);
2188 :
2189 31152681 : for (unsigned int i = 0; i != nc; i++)
2190 49180 : this->set_child(i, nullptr);
2191 : }
2192 :
2193 49532 : libmesh_assert (this->_children[c] == nullptr || this->child_ptr(c) == remote_elem);
2194 49532 : libmesh_assert (elem == remote_elem || this == elem->parent());
2195 :
2196 49532 : this->set_child(c, elem);
2197 50245814 : }
2198 :
2199 :
2200 :
2201 231063 : void Elem::replace_child (Elem * elem, unsigned int c)
2202 : {
2203 24104 : libmesh_assert(this->has_children());
2204 :
2205 24104 : libmesh_assert(this->child_ptr(c));
2206 :
2207 24104 : this->set_child(c, elem);
2208 231063 : }
2209 :
2210 :
2211 :
2212 9501 : void Elem::family_tree (std::vector<const Elem *> & family,
2213 : bool reset) const
2214 : {
2215 9501 : ElemInternal::family_tree(this, family, reset);
2216 9501 : }
2217 :
2218 :
2219 :
2220 0 : void Elem::family_tree (std::vector<Elem *> & family,
2221 : bool reset)
2222 : {
2223 0 : ElemInternal::family_tree(this, family, reset);
2224 0 : }
2225 :
2226 :
2227 :
2228 46589 : void Elem::total_family_tree (std::vector<const Elem *> & family,
2229 : bool reset) const
2230 : {
2231 46589 : ElemInternal::total_family_tree(this, family, reset);
2232 46589 : }
2233 :
2234 :
2235 :
2236 32349676 : void Elem::total_family_tree (std::vector<Elem *> & family,
2237 : bool reset)
2238 : {
2239 32349676 : ElemInternal::total_family_tree(this, family, reset);
2240 32349676 : }
2241 :
2242 :
2243 :
2244 1032935 : void Elem::active_family_tree (std::vector<const Elem *> & active_family,
2245 : bool reset) const
2246 : {
2247 1032935 : ElemInternal::active_family_tree(this, active_family, reset);
2248 1032935 : }
2249 :
2250 :
2251 :
2252 0 : void Elem::active_family_tree (std::vector<Elem *> & active_family,
2253 : bool reset)
2254 : {
2255 0 : ElemInternal::active_family_tree(this, active_family, reset);
2256 0 : }
2257 :
2258 :
2259 :
2260 0 : void Elem::family_tree_by_side (std::vector<const Elem *> & family,
2261 : unsigned int side,
2262 : bool reset) const
2263 : {
2264 0 : ElemInternal::family_tree_by_side(this, family, side, reset);
2265 0 : }
2266 :
2267 :
2268 :
2269 0 : void Elem:: family_tree_by_side (std::vector<Elem *> & family,
2270 : unsigned int side,
2271 : bool reset)
2272 : {
2273 0 : ElemInternal::family_tree_by_side(this, family, side, reset);
2274 0 : }
2275 :
2276 :
2277 :
2278 0 : void Elem::total_family_tree_by_side (std::vector<const Elem *> & family,
2279 : unsigned int side,
2280 : bool reset) const
2281 : {
2282 0 : ElemInternal::total_family_tree_by_side(this, family, side, reset);
2283 0 : }
2284 :
2285 :
2286 :
2287 230861545 : void Elem::total_family_tree_by_side (std::vector<Elem *> & family,
2288 : unsigned int side,
2289 : bool reset)
2290 : {
2291 230861545 : ElemInternal::total_family_tree_by_side(this, family, side, reset);
2292 230861545 : }
2293 :
2294 :
2295 :
2296 591279 : void Elem::active_family_tree_by_side (std::vector<const Elem *> & family,
2297 : unsigned int side,
2298 : bool reset) const
2299 : {
2300 591279 : ElemInternal::active_family_tree_by_side(this, family, side, reset);
2301 591279 : }
2302 :
2303 :
2304 :
2305 0 : void Elem::active_family_tree_by_side (std::vector<Elem *> & family,
2306 : unsigned int side,
2307 : bool reset)
2308 : {
2309 0 : ElemInternal::active_family_tree_by_side(this, family, side, reset);
2310 0 : }
2311 :
2312 :
2313 :
2314 0 : void Elem::family_tree_by_neighbor (std::vector<const Elem *> & family,
2315 : const Elem * neighbor,
2316 : bool reset) const
2317 : {
2318 0 : ElemInternal::family_tree_by_neighbor(this, family, neighbor, reset);
2319 0 : }
2320 :
2321 :
2322 :
2323 0 : void Elem::family_tree_by_neighbor (std::vector<Elem *> & family,
2324 : Elem * neighbor,
2325 : bool reset)
2326 : {
2327 0 : ElemInternal::family_tree_by_neighbor(this, family, neighbor, reset);
2328 0 : }
2329 :
2330 :
2331 :
2332 0 : void Elem::total_family_tree_by_neighbor (std::vector<const Elem *> & family,
2333 : const Elem * neighbor,
2334 : bool reset) const
2335 : {
2336 0 : ElemInternal::total_family_tree_by_neighbor(this, family, neighbor, reset);
2337 0 : }
2338 :
2339 :
2340 :
2341 70601070 : void Elem::total_family_tree_by_neighbor (std::vector<Elem *> & family,
2342 : Elem * neighbor,
2343 : bool reset)
2344 : {
2345 70601070 : ElemInternal::total_family_tree_by_neighbor(this, family, neighbor, reset);
2346 70601070 : }
2347 :
2348 :
2349 :
2350 0 : void Elem::family_tree_by_subneighbor (std::vector<const Elem *> & family,
2351 : const Elem * neighbor,
2352 : const Elem * subneighbor,
2353 : bool reset) const
2354 : {
2355 0 : ElemInternal::family_tree_by_subneighbor(this, family, neighbor, subneighbor, reset);
2356 0 : }
2357 :
2358 :
2359 :
2360 0 : void Elem::family_tree_by_subneighbor (std::vector<Elem *> & family,
2361 : Elem * neighbor,
2362 : Elem * subneighbor,
2363 : bool reset)
2364 : {
2365 0 : ElemInternal::family_tree_by_subneighbor(this, family, neighbor, subneighbor, reset);
2366 0 : }
2367 :
2368 :
2369 :
2370 0 : void Elem::total_family_tree_by_subneighbor (std::vector<const Elem *> & family,
2371 : const Elem * neighbor,
2372 : const Elem * subneighbor,
2373 : bool reset) const
2374 : {
2375 0 : ElemInternal::total_family_tree_by_subneighbor(this, family, neighbor, subneighbor, reset);
2376 0 : }
2377 :
2378 :
2379 :
2380 2874 : void Elem::total_family_tree_by_subneighbor (std::vector<Elem *> & family,
2381 : Elem * neighbor,
2382 : Elem * subneighbor,
2383 : bool reset)
2384 : {
2385 2874 : ElemInternal::total_family_tree_by_subneighbor(this, family, neighbor, subneighbor, reset);
2386 2874 : }
2387 :
2388 :
2389 :
2390 42396543 : void Elem::active_family_tree_by_neighbor (std::vector<const Elem *> & family,
2391 : const Elem * neighbor,
2392 : bool reset) const
2393 : {
2394 42396543 : ElemInternal::active_family_tree_by_neighbor(this, family, neighbor, reset);
2395 42396543 : }
2396 :
2397 :
2398 :
2399 2736 : void Elem::active_family_tree_by_neighbor (std::vector<Elem *> & family,
2400 : Elem * neighbor,
2401 : bool reset)
2402 : {
2403 2736 : ElemInternal::active_family_tree_by_neighbor(this, family, neighbor, reset);
2404 2736 : }
2405 :
2406 :
2407 :
2408 26334 : void Elem::active_family_tree_by_topological_neighbor (std::vector<const Elem *> & family,
2409 : const Elem * neighbor,
2410 : const MeshBase & mesh,
2411 : const PointLocatorBase & point_locator,
2412 : const PeriodicBoundaries * pb,
2413 : bool reset) const
2414 : {
2415 26334 : ElemInternal::active_family_tree_by_topological_neighbor(this, family, neighbor,
2416 : mesh, point_locator, pb,
2417 : reset);
2418 26334 : }
2419 :
2420 :
2421 :
2422 0 : void Elem::active_family_tree_by_topological_neighbor (std::vector<Elem *> & family,
2423 : Elem * neighbor,
2424 : const MeshBase & mesh,
2425 : const PointLocatorBase & point_locator,
2426 : const PeriodicBoundaries * pb,
2427 : bool reset)
2428 : {
2429 0 : ElemInternal::active_family_tree_by_topological_neighbor(this, family, neighbor,
2430 : mesh, point_locator, pb,
2431 : reset);
2432 0 : }
2433 :
2434 :
2435 27663040 : bool Elem::is_child_on_edge(const unsigned int c,
2436 : const unsigned int e) const
2437 : {
2438 16936768 : libmesh_assert_less (c, this->n_children());
2439 16936768 : libmesh_assert_less (e, this->n_edges());
2440 :
2441 43551954 : std::unique_ptr<const Elem> my_edge = this->build_edge_ptr(e);
2442 26615186 : std::unique_ptr<const Elem> child_edge = this->child_ptr(c)->build_edge_ptr(e);
2443 :
2444 : // We're assuming that an overlapping child edge has the same
2445 : // number and orientation as its parent
2446 43251682 : return (child_edge->node_id(0) == my_edge->node_id(0) ||
2447 58092742 : child_edge->node_id(1) == my_edge->node_id(1));
2448 9678418 : }
2449 :
2450 :
2451 :
2452 15961514 : unsigned int Elem::min_p_level_by_neighbor(const Elem * neighbor_in,
2453 : unsigned int current_min) const
2454 : {
2455 1450032 : libmesh_assert(!this->subactive());
2456 1450032 : libmesh_assert(neighbor_in->active());
2457 :
2458 : // If we're an active element this is simple
2459 1450032 : if (this->active())
2460 17193708 : return std::min(current_min, this->p_level());
2461 :
2462 21019 : libmesh_assert(has_neighbor(neighbor_in));
2463 :
2464 : // The p_level() of an ancestor element is already the minimum
2465 : // p_level() of its children - so if that's high enough, we don't
2466 : // need to examine any children.
2467 218366 : if (current_min <= this->p_level())
2468 21019 : return current_min;
2469 :
2470 0 : unsigned int min_p_level = current_min;
2471 :
2472 0 : for (auto & c : this->child_ref_range())
2473 0 : if (&c != remote_elem && c.has_neighbor(neighbor_in))
2474 : min_p_level =
2475 0 : c.min_p_level_by_neighbor(neighbor_in, min_p_level);
2476 :
2477 0 : return min_p_level;
2478 : }
2479 :
2480 :
2481 2145470 : unsigned int Elem::min_new_p_level_by_neighbor(const Elem * neighbor_in,
2482 : unsigned int current_min) const
2483 : {
2484 95716 : libmesh_assert(!this->subactive());
2485 95716 : libmesh_assert(neighbor_in->active());
2486 :
2487 : // If we're an active element this is simple
2488 95716 : if (this->active())
2489 : {
2490 1428913 : unsigned int new_p_level = this->p_level();
2491 1493213 : if (this->p_refinement_flag() == Elem::REFINE)
2492 16 : new_p_level += 1;
2493 1428913 : if (this->p_refinement_flag() == Elem::COARSEN)
2494 : {
2495 4 : libmesh_assert_greater (new_p_level, 0);
2496 218 : new_p_level -= 1;
2497 : }
2498 1428913 : return std::min(current_min, new_p_level);
2499 : }
2500 :
2501 31416 : libmesh_assert(has_neighbor(neighbor_in));
2502 :
2503 716557 : unsigned int min_p_level = current_min;
2504 :
2505 3600199 : for (auto & c : this->child_ref_range())
2506 3012242 : if (&c != remote_elem && c.has_neighbor(neighbor_in))
2507 : min_p_level =
2508 1428913 : c.min_new_p_level_by_neighbor(neighbor_in, min_p_level);
2509 :
2510 716557 : return min_p_level;
2511 : }
2512 :
2513 :
2514 :
2515 192875749 : unsigned int Elem::as_parent_node (unsigned int child,
2516 : unsigned int child_node) const
2517 : {
2518 192875749 : const unsigned int nc = this->n_children();
2519 8507910 : libmesh_assert_less(child, nc);
2520 :
2521 : // Cached return values, indexed first by embedding_matrix version,
2522 : // then by child number, then by child node number.
2523 : std::vector<std::vector<std::vector<signed char>>> &
2524 192875749 : cached_parent_indices = this->_get_parent_indices_cache();
2525 :
2526 192875749 : unsigned int em_vers = this->embedding_matrix_version();
2527 :
2528 : // We may be updating the cache on one thread, and while that
2529 : // happens we can't safely access the cache from other threads.
2530 17015820 : Threads::spin_mutex::scoped_lock lock(parent_indices_mutex);
2531 :
2532 201383659 : if (em_vers >= cached_parent_indices.size())
2533 5898 : cached_parent_indices.resize(em_vers+1);
2534 :
2535 209891569 : if (child >= cached_parent_indices[em_vers].size())
2536 : {
2537 5902 : const signed char nn = cast_int<signed char>(this->n_nodes());
2538 :
2539 6118 : cached_parent_indices[em_vers].resize(nc);
2540 :
2541 35363 : for (unsigned int c = 0; c != nc; ++c)
2542 : {
2543 29461 : const unsigned int ncn = this->n_nodes_in_child(c);
2544 31569 : cached_parent_indices[em_vers][c].resize(ncn);
2545 300082 : for (unsigned int cn = 0; cn != ncn; ++cn)
2546 : {
2547 1288292 : for (signed char n = 0; n != nn; ++n)
2548 : {
2549 : const Real em_val = this->embedding_matrix
2550 1288292 : (c, cn, n);
2551 1288292 : if (em_val == 1)
2552 : {
2553 93530 : cached_parent_indices[em_vers][c][cn] = n;
2554 87502 : break;
2555 : }
2556 :
2557 1200790 : if (em_val != 0)
2558 : {
2559 195727 : cached_parent_indices[em_vers][c][cn] =
2560 : -1;
2561 183119 : break;
2562 : }
2563 :
2564 : // We should never see an all-zero embedding matrix
2565 : // row
2566 34690 : libmesh_assert_not_equal_to (n+1, nn);
2567 : }
2568 : }
2569 : }
2570 : }
2571 :
2572 : const signed char cache_val =
2573 209891569 : cached_parent_indices[em_vers][child][child_node];
2574 192875749 : if (cache_val == -1)
2575 5359552 : return libMesh::invalid_uint;
2576 :
2577 74871487 : return cached_parent_indices[em_vers][child][child_node];
2578 : }
2579 :
2580 :
2581 :
2582 : const std::vector<std::pair<unsigned char, unsigned char>> &
2583 111235358 : Elem::parent_bracketing_nodes(unsigned int child,
2584 : unsigned int child_node) const
2585 : {
2586 : // Indexed first by embedding matrix type, then by child id, then by
2587 : // child node, then by bracketing pair
2588 : std::vector<std::vector<std::vector<std::vector<std::pair<unsigned char, unsigned char>>>>> &
2589 111235358 : cached_bracketing_nodes = this->_get_bracketing_node_cache();
2590 :
2591 111235358 : const unsigned int em_vers = this->embedding_matrix_version();
2592 :
2593 : // We may be updating the cache on one thread, and while that
2594 : // happens we can't safely access the cache from other threads.
2595 11625076 : Threads::spin_mutex::scoped_lock lock(parent_bracketing_nodes_mutex);
2596 :
2597 117047896 : if (cached_bracketing_nodes.size() <= em_vers)
2598 5269 : cached_bracketing_nodes.resize(em_vers+1);
2599 :
2600 111235358 : const unsigned int nc = this->n_children();
2601 :
2602 : // If we haven't cached the bracketing nodes corresponding to this
2603 : // embedding matrix yet, let's do so now.
2604 122860434 : if (cached_bracketing_nodes[em_vers].size() < nc)
2605 : {
2606 : // If we're a second-order element but we're not a full-order
2607 : // element, then some of our bracketing nodes may not exist
2608 : // except on the equivalent full-order element. Let's build an
2609 : // equivalent full-order element and make a copy of its cache to
2610 : // use.
2611 8515 : if (this->default_order() != FIRST &&
2612 3242 : second_order_equivalent_type(this->type(), /*full_ordered=*/ true) != this->type())
2613 : {
2614 : // Check that we really are the non-full-order type
2615 8 : libmesh_assert_equal_to
2616 : (second_order_equivalent_type (this->type(), false),
2617 : this->type());
2618 :
2619 : // Build the full-order type
2620 : ElemType full_type =
2621 165 : second_order_equivalent_type(this->type(), /*full_ordered=*/ true);
2622 173 : std::unique_ptr<Elem> full_elem = Elem::build(full_type);
2623 :
2624 : // This won't work for elements with multiple
2625 : // embedding_matrix versions, but every such element is full
2626 : // order anyways.
2627 8 : libmesh_assert_equal_to(em_vers, 0);
2628 :
2629 : // Make sure its cache has been built. We temporarily
2630 : // release our mutex lock so that the inner call can
2631 : // re-acquire it.
2632 8 : lock.release();
2633 165 : full_elem->parent_bracketing_nodes(0,0);
2634 :
2635 : // And then we need to lock again, so that if someone *else*
2636 : // grabbed our lock before we did we don't risk accessing
2637 : // cached_bracketing_nodes while they're working on it.
2638 : // Threading is hard.
2639 8 : lock.acquire(parent_bracketing_nodes_mutex);
2640 :
2641 : // Copy its cache
2642 : cached_bracketing_nodes =
2643 165 : full_elem->_get_bracketing_node_cache();
2644 :
2645 : // Now we don't need to build the cache ourselves.
2646 181 : return cached_bracketing_nodes[em_vers][child][child_node];
2647 149 : }
2648 :
2649 5300 : cached_bracketing_nodes[em_vers].resize(nc);
2650 :
2651 5108 : const unsigned int nn = this->n_nodes();
2652 :
2653 : // We have to examine each child
2654 29097 : for (unsigned int c = 0; c != nc; ++c)
2655 : {
2656 23989 : const unsigned int ncn = this->n_nodes_in_child(c);
2657 :
2658 25777 : cached_bracketing_nodes[em_vers][c].resize(ncn);
2659 :
2660 : // We have to examine each node in that child
2661 208186 : for (unsigned int n = 0; n != ncn; ++n)
2662 : {
2663 : // If this child node isn't a vertex or an infinite
2664 : // child element's mid-infinite-edge node, then we need
2665 : // to find bracketing nodes on the child.
2666 184197 : if (!this->is_vertex_on_child(c, n)
2667 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2668 16175 : && !this->is_mid_infinite_edge_node(n)
2669 : #endif
2670 : )
2671 : {
2672 : // Use the embedding matrix to find the child node
2673 : // location in parent master element space
2674 2572 : Point bracketed_pt;
2675 :
2676 1041788 : for (unsigned int pn = 0; pn != nn; ++pn)
2677 : {
2678 : const Real em_val =
2679 968286 : this->embedding_matrix(c,n,pn);
2680 :
2681 34140 : libmesh_assert_not_equal_to (em_val, 1);
2682 968286 : if (em_val != 0.)
2683 356334 : bracketed_pt.add_scaled(this->master_point(pn), em_val);
2684 : }
2685 :
2686 : // Check each pair of nodes on the child which are
2687 : // also both parent nodes
2688 1041788 : for (unsigned int n1 = 0; n1 != ncn; ++n1)
2689 : {
2690 968286 : if (n1 == n)
2691 610786 : continue;
2692 :
2693 : unsigned int parent_n1 =
2694 894784 : this->as_parent_node(c,n1);
2695 :
2696 894784 : if (parent_n1 == libMesh::invalid_uint)
2697 518524 : continue;
2698 :
2699 357500 : Point p1 = this->master_point(parent_n1);
2700 :
2701 3773206 : for (unsigned int n2 = n1+1; n2 < nn; ++n2)
2702 : {
2703 3503924 : if (n2 == n)
2704 2808158 : continue;
2705 :
2706 : unsigned int parent_n2 =
2707 3234642 : this->as_parent_node(c,n2);
2708 :
2709 3234642 : if (parent_n2 == libMesh::invalid_uint)
2710 2450212 : continue;
2711 :
2712 695766 : Point p2 = this->master_point(parent_n2);
2713 :
2714 25444 : Point pmid = (p1 + p2)/2;
2715 :
2716 25444 : if (pmid == bracketed_pt)
2717 : {
2718 97470 : cached_bracketing_nodes[em_vers][c][n].emplace_back(parent_n1, parent_n2);
2719 88218 : break;
2720 : }
2721 : else
2722 22360 : libmesh_assert(!pmid.absolute_fuzzy_equals(bracketed_pt));
2723 : }
2724 : }
2725 : }
2726 : // If this child node is a parent node, we need to
2727 : // find bracketing nodes on the parent.
2728 : else
2729 : {
2730 110695 : unsigned int parent_node = this->as_parent_node(c,n);
2731 :
2732 4202 : Point bracketed_pt;
2733 :
2734 : // If we're not a parent node, use the embedding
2735 : // matrix to find the child node location in parent
2736 : // master element space
2737 110695 : if (parent_node == libMesh::invalid_uint)
2738 : {
2739 351898 : for (unsigned int pn = 0; pn != nn; ++pn)
2740 : {
2741 : const Real em_val =
2742 302681 : this->embedding_matrix(c,n,pn);
2743 :
2744 11828 : libmesh_assert_not_equal_to (em_val, 1);
2745 302681 : if (em_val != 0.)
2746 151178 : bracketed_pt.add_scaled(this->master_point(pn), em_val);
2747 : }
2748 : }
2749 : // If we're a parent node then we need no arithmetic
2750 : else
2751 61478 : bracketed_pt = this->master_point(parent_node);
2752 :
2753 1053222 : for (unsigned int n1 = 0; n1 != nn; ++n1)
2754 : {
2755 942527 : if (n1 == parent_node)
2756 61478 : continue;
2757 :
2758 881049 : Point p1 = this->master_point(n1);
2759 :
2760 4831850 : for (unsigned int n2 = n1+1; n2 < nn; ++n2)
2761 : {
2762 4105098 : if (n2 == parent_node)
2763 10792 : continue;
2764 :
2765 3813036 : Point pmid = (p1 + this->master_point(n2))/2;
2766 :
2767 140116 : if (pmid == bracketed_pt)
2768 : {
2769 171241 : cached_bracketing_nodes[em_vers][c][n].emplace_back(n1, n2);
2770 5648 : break;
2771 : }
2772 : else
2773 134468 : libmesh_assert(!pmid.absolute_fuzzy_equals(bracketed_pt));
2774 : }
2775 : }
2776 : }
2777 : }
2778 : }
2779 : }
2780 :
2781 128672783 : return cached_bracketing_nodes[em_vers][child][child_node];
2782 : }
2783 :
2784 :
2785 : const std::vector<std::pair<dof_id_type, dof_id_type>>
2786 120818012 : Elem::bracketing_nodes(unsigned int child,
2787 : unsigned int child_node) const
2788 : {
2789 6234618 : std::vector<std::pair<dof_id_type, dof_id_type>> returnval;
2790 :
2791 : const std::vector<std::pair<unsigned char, unsigned char>> & pbc =
2792 120818012 : this->parent_bracketing_nodes(child,child_node);
2793 :
2794 254626622 : for (const auto & pb : pbc)
2795 : {
2796 133808610 : const unsigned short n_n = this->n_nodes();
2797 133808610 : if (pb.first < n_n && pb.second < n_n)
2798 147273666 : returnval.emplace_back(this->node_id(pb.first), this->node_id(pb.second));
2799 : else
2800 : {
2801 : // We must be on a non-full-order higher order element...
2802 33008 : libmesh_assert_not_equal_to(this->default_order(), FIRST);
2803 33008 : libmesh_assert_not_equal_to
2804 : (second_order_equivalent_type (this->type(), true),
2805 : this->type());
2806 33008 : libmesh_assert_equal_to
2807 : (second_order_equivalent_type (this->type(), false),
2808 : this->type());
2809 :
2810 : // And that's a shame, because this is a nasty search:
2811 :
2812 : // Build the full-order type
2813 : ElemType full_type =
2814 908404 : second_order_equivalent_type(this->type(), /*full_ordered=*/ true);
2815 941412 : std::unique_ptr<Elem> full_elem = Elem::build(full_type);
2816 :
2817 908404 : dof_id_type pt1 = DofObject::invalid_id;
2818 908404 : dof_id_type pt2 = DofObject::invalid_id;
2819 :
2820 : // Find the bracketing nodes by figuring out what
2821 : // already-created children will have them.
2822 :
2823 : // This only doesn't break horribly because we add children
2824 : // and nodes in straightforward + hierarchical orders...
2825 4993928 : for (unsigned int c=0; c <= child; ++c)
2826 : {
2827 147320 : libmesh_assert(this->child_ptr(c));
2828 4085524 : if (this->child_ptr(c) == remote_elem)
2829 0 : continue;
2830 :
2831 73004408 : for (auto n : make_range(this->n_nodes_in_child(c)))
2832 : {
2833 66036404 : if (c == child && n == child_node)
2834 33008 : break;
2835 :
2836 67475880 : if (pb.first == full_elem->as_parent_node(c,n))
2837 : {
2838 : // We should be consistent
2839 99160 : if (pt1 != DofObject::invalid_id)
2840 66330 : libmesh_assert_equal_to(pt1, this->child_ptr(c)->node_id(n));
2841 :
2842 2846778 : pt1 = this->child_ptr(c)->node_id(n);
2843 : }
2844 :
2845 67475880 : if (pb.second == full_elem->as_parent_node(c,n))
2846 : {
2847 : // We should be consistent
2848 72088 : if (pt2 != DofObject::invalid_id)
2849 40682 : libmesh_assert_equal_to(pt2, this->child_ptr(c)->node_id(n));
2850 :
2851 2068296 : pt2 = this->child_ptr(c)->node_id(n);
2852 : }
2853 : }
2854 : }
2855 :
2856 : // We should *usually* find all bracketing nodes by the time
2857 : // we query them (again, because of the child & node add
2858 : // order)
2859 : //
2860 : // The exception is if we're a HEX20, in which case we will
2861 : // find pairs of vertex nodes and edge nodes bracketing the
2862 : // new central node but we *won't* find the pairs of face
2863 : // nodes which we would have had on a HEX27. In that case
2864 : // we'll still have enough bracketing nodes for a
2865 : // topological lookup, but we won't be able to make the
2866 : // following assertions.
2867 908404 : if (this->type() != HEX20)
2868 : {
2869 15920 : libmesh_assert_not_equal_to (pt1, DofObject::invalid_id);
2870 15920 : libmesh_assert_not_equal_to (pt2, DofObject::invalid_id);
2871 : }
2872 :
2873 908404 : if (pt1 != DofObject::invalid_id &&
2874 903549 : pt2 != DofObject::invalid_id)
2875 864709 : returnval.emplace_back(pt1, pt2);
2876 842388 : }
2877 : }
2878 :
2879 120818012 : return returnval;
2880 : }
2881 : #endif // #ifdef LIBMESH_ENABLE_AMR
2882 :
2883 :
2884 :
2885 :
2886 146656638 : bool Elem::contains_point (const Point & p, Real tol) const
2887 : {
2888 : // We currently allow the user to enlarge the bounding box by
2889 : // providing a tol > TOLERANCE (so this routine is identical to
2890 : // Elem::close_to_point()), but print a warning so that the
2891 : // user can eventually switch his code over to calling close_to_point()
2892 : // instead, which is intended to be used for this purpose.
2893 146656638 : if (tol > TOLERANCE)
2894 : {
2895 0 : libmesh_do_once(libMesh::err
2896 : << "WARNING: Resizing bounding box to match user-specified tolerance!\n"
2897 : << "In the future, calls to Elem::contains_point() with tol > TOLERANCE\n"
2898 : << "will be more optimized, but should not be used\n"
2899 : << "to search for points 'close to' elements!\n"
2900 : << "Instead, use Elem::close_to_point() for this purpose.\n"
2901 : << std::endl;);
2902 0 : return this->point_test(p, tol, tol);
2903 : }
2904 : else
2905 146656638 : return this->point_test(p, TOLERANCE, tol);
2906 : }
2907 :
2908 :
2909 :
2910 :
2911 6059349 : bool Elem::close_to_point (const Point & p, Real tol) const
2912 : {
2913 : // This test uses the user's passed-in tolerance for the
2914 : // bounding box test as well, thereby allowing the routine to
2915 : // find points which are not only "in" the element, but also
2916 : // "nearby" to within some tolerance.
2917 6059349 : return this->point_test(p, tol, tol);
2918 : }
2919 :
2920 :
2921 :
2922 :
2923 152715987 : bool Elem::point_test(const Point & p, Real box_tol, Real map_tol) const
2924 : {
2925 7579350 : libmesh_assert_greater (box_tol, 0.);
2926 7579350 : libmesh_assert_greater (map_tol, 0.);
2927 :
2928 : // This is a great optimization on first order elements, but it
2929 : // could return false negatives on higher orders
2930 152715987 : if (this->default_order() == FIRST)
2931 : {
2932 : // Check to make sure the element *could* contain this point, so we
2933 : // can avoid an expensive inverse_map call if it doesn't.
2934 : bool
2935 : #if LIBMESH_DIM > 2
2936 3293746 : point_above_min_z = false,
2937 3293746 : point_below_max_z = false,
2938 : #endif
2939 : #if LIBMESH_DIM > 1
2940 3293746 : point_above_min_y = false,
2941 3293746 : point_below_max_y = false,
2942 : #endif
2943 3293746 : point_above_min_x = false,
2944 3293746 : point_below_max_x = false;
2945 :
2946 : // For relative bounding box checks in physical space
2947 62447328 : const Real my_hmax = this->hmax();
2948 :
2949 494789637 : for (auto & n : this->node_ref_range())
2950 : {
2951 432342309 : point_above_min_x = point_above_min_x || (n(0) - my_hmax*box_tol <= p(0));
2952 432342309 : point_below_max_x = point_below_max_x || (n(0) + my_hmax*box_tol >= p(0));
2953 : #if LIBMESH_DIM > 1
2954 432342309 : point_above_min_y = point_above_min_y || (n(1) - my_hmax*box_tol <= p(1));
2955 432342309 : point_below_max_y = point_below_max_y || (n(1) + my_hmax*box_tol >= p(1));
2956 : #endif
2957 : #if LIBMESH_DIM > 2
2958 432342309 : point_above_min_z = point_above_min_z || (n(2) - my_hmax*box_tol <= p(2));
2959 432342309 : point_below_max_z = point_below_max_z || (n(2) + my_hmax*box_tol >= p(2));
2960 : #endif
2961 : }
2962 :
2963 62447328 : if (
2964 : #if LIBMESH_DIM > 2
2965 3293746 : !point_above_min_z ||
2966 3176240 : !point_below_max_z ||
2967 : #endif
2968 : #if LIBMESH_DIM > 1
2969 40278562 : !point_above_min_y ||
2970 2095221 : !point_below_max_y ||
2971 : #endif
2972 18053619 : !point_above_min_x ||
2973 885412 : !point_below_max_x)
2974 2565832 : return false;
2975 : }
2976 :
2977 : // To be on the safe side, we converge the inverse_map() iteration
2978 : // to a slightly tighter tolerance than that requested by the
2979 : // user...
2980 103738658 : const Point mapped_point = FEMap::inverse_map(this->dim(),
2981 : this,
2982 : p,
2983 : 0.1*map_tol, // <- this is |dx| tolerance, the Newton residual should be ~ |dx|^2
2984 : /*secure=*/ false,
2985 9928898 : /*extra_checks=*/ false);
2986 :
2987 : // Check that the refspace point maps back to p! This is only necessary
2988 : // for 1D and 2D elements, 3D elements always live in 3D.
2989 : //
2990 : // TODO: The contains_point() function could most likely be implemented
2991 : // more efficiently in the element sub-classes themselves, at least for
2992 : // the linear element types.
2993 98725140 : if (this->dim() < 3)
2994 : {
2995 30397000 : Point xyz = FEMap::map(this->dim(), this, mapped_point);
2996 :
2997 : // Compute the distance between the original point and the re-mapped point.
2998 : // They should be in the same place.
2999 28098632 : Real dist = (xyz - p).norm();
3000 :
3001 :
3002 : // If dist is larger than some fraction of the tolerance, then return false.
3003 : // This can happen when e.g. a 2D element is living in 3D, and
3004 : // FEMap::inverse_map() maps p onto the projection of the element,
3005 : // effectively "tricking" on_reference_element().
3006 30397000 : if (dist > this->hmax() * map_tol)
3007 8440 : return false;
3008 : }
3009 :
3010 98716700 : return this->on_reference_element(mapped_point, map_tol);
3011 : }
3012 :
3013 :
3014 :
3015 :
3016 78600 : bool Elem::has_invertible_map(Real /*tol*/) const
3017 : {
3018 85150 : QNodal qnodal {this->dim()};
3019 91700 : FEMap fe_map;
3020 6550 : auto & jac = fe_map.get_jacobian();
3021 :
3022 : // We have a separate check for is_singular_node() below, so in this
3023 : // case its "OK" to do nodal quadrature on pyramids.
3024 78600 : qnodal.allow_nodal_pyramid_quadrature = true;
3025 78600 : qnodal.init(*this);
3026 6550 : auto & qp = qnodal.get_points();
3027 6550 : libmesh_assert_equal_to(qp.size(), this->n_nodes());
3028 :
3029 85150 : std::vector<Point> one_point(1);
3030 85150 : std::vector<Real> one_weight(1,1);
3031 1086900 : for (auto i : index_range(qp))
3032 : {
3033 1008300 : if (this->is_singular_node(i))
3034 8448 : continue;
3035 :
3036 999084 : one_point[0] = qp[i];
3037 :
3038 999084 : fe_map.init_reference_to_physical_map(this->dim(), one_point, this);
3039 999084 : fe_map.compute_map(this->dim(), one_weight, this, false);
3040 :
3041 999084 : if (jac[0] <= 0)
3042 0 : return false;
3043 : }
3044 :
3045 13100 : return true;
3046 65500 : }
3047 :
3048 :
3049 :
3050 0 : void Elem::print_info (std::ostream & os) const
3051 : {
3052 0 : os << this->get_info()
3053 0 : << std::endl;
3054 0 : }
3055 :
3056 :
3057 :
3058 0 : std::string Elem::get_info () const
3059 : {
3060 0 : std::ostringstream oss;
3061 :
3062 : oss << " Elem Information" << '\n'
3063 0 : << " id()=";
3064 :
3065 0 : if (this->valid_id())
3066 0 : oss << this->id();
3067 : else
3068 0 : oss << "invalid";
3069 :
3070 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
3071 0 : oss << ", unique_id()=";
3072 0 : if (this->valid_unique_id())
3073 0 : oss << this->unique_id();
3074 : else
3075 0 : oss << "invalid";
3076 : #endif
3077 :
3078 0 : oss << ", subdomain_id()=" << this->subdomain_id();
3079 0 : oss << ", processor_id()=" << this->processor_id() << '\n';
3080 :
3081 0 : oss << " type()=" << Utility::enum_to_string(this->type()) << '\n'
3082 0 : << " dim()=" << this->dim() << '\n'
3083 0 : << " n_nodes()=" << this->n_nodes() << '\n';
3084 :
3085 0 : oss << " mapping=" << Utility::enum_to_string(this->mapping_type()) << '\n';
3086 :
3087 0 : for (auto n : this->node_index_range())
3088 : {
3089 0 : oss << " " << n << this->node_ref(n);
3090 0 : if (this->mapping_type() == RATIONAL_BERNSTEIN_MAP)
3091 : {
3092 0 : const unsigned char datum_index = this->mapping_data();
3093 0 : oss << " weight=" <<
3094 0 : this->node_ref(n).get_extra_datum<Real>(datum_index) << '\n';
3095 : }
3096 : }
3097 :
3098 0 : oss << " n_sides()=" << this->n_sides() << '\n';
3099 :
3100 0 : for (auto s : this->side_index_range())
3101 : {
3102 0 : oss << " neighbor(" << s << ")=";
3103 0 : if (this->neighbor_ptr(s))
3104 0 : oss << this->neighbor_ptr(s)->id() << '\n';
3105 : else
3106 0 : oss << "nullptr\n";
3107 : }
3108 :
3109 0 : if (!this->infinite())
3110 : {
3111 0 : oss << " hmin()=" << this->hmin()
3112 0 : << ", hmax()=" << this->hmax() << '\n'
3113 0 : << " volume()=" << this->volume() << '\n';
3114 : }
3115 0 : oss << " active()=" << this->active()
3116 0 : << ", ancestor()=" << this->ancestor()
3117 0 : << ", subactive()=" << this->subactive()
3118 0 : << ", has_children()=" << this->has_children() << '\n'
3119 0 : << " parent()=";
3120 0 : if (this->parent())
3121 0 : oss << this->parent()->id() << '\n';
3122 : else
3123 0 : oss << "nullptr\n";
3124 0 : oss << " level()=" << this->level()
3125 0 : << ", p_level()=" << this->p_level() << '\n'
3126 : #ifdef LIBMESH_ENABLE_AMR
3127 0 : << " refinement_flag()=" << Utility::enum_to_string(this->refinement_flag()) << '\n'
3128 0 : << " p_refinement_flag()=" << Utility::enum_to_string(this->p_refinement_flag()) << '\n'
3129 : #endif
3130 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
3131 0 : << " infinite()=" << this->infinite() << '\n';
3132 0 : if (this->infinite())
3133 0 : oss << " origin()=" << this->origin() << '\n'
3134 : #endif
3135 : ;
3136 :
3137 0 : oss << " DoFs=";
3138 0 : for (auto s : make_range(this->n_systems()))
3139 0 : for (auto v : make_range(this->n_vars(s)))
3140 0 : for (auto c : make_range(this->n_comp(s,v)))
3141 0 : oss << '(' << s << '/' << v << '/' << this->dof_number(s,v,c) << ") ";
3142 :
3143 :
3144 0 : return oss.str();
3145 0 : }
3146 :
3147 :
3148 :
3149 0 : void Elem::nullify_neighbors ()
3150 : {
3151 : // Tell any of my neighbors about my death...
3152 : // Looks strange, huh?
3153 0 : for (auto n : this->side_index_range())
3154 : {
3155 0 : Elem * current_neighbor = this->neighbor_ptr(n);
3156 0 : if (current_neighbor && current_neighbor != remote_elem)
3157 : {
3158 : // Note: it is possible that I see the neighbor
3159 : // (which is coarser than me)
3160 : // but they don't see me, so avoid that case.
3161 0 : if (current_neighbor->level() == this->level())
3162 : {
3163 0 : const unsigned int w_n_a_i = current_neighbor->which_neighbor_am_i(this);
3164 0 : libmesh_assert_less (w_n_a_i, current_neighbor->n_neighbors());
3165 0 : current_neighbor->set_neighbor(w_n_a_i, nullptr);
3166 0 : this->set_neighbor(n, nullptr);
3167 : }
3168 : }
3169 : }
3170 0 : }
3171 :
3172 :
3173 :
3174 :
3175 0 : unsigned int Elem::n_second_order_adjacent_vertices (const unsigned int) const
3176 : {
3177 : // for linear elements, always return 0
3178 0 : return 0;
3179 : }
3180 :
3181 :
3182 :
3183 0 : unsigned short int Elem::second_order_adjacent_vertex (const unsigned int,
3184 : const unsigned int) const
3185 : {
3186 : // for linear elements, always return 0
3187 0 : return 0;
3188 : }
3189 :
3190 :
3191 :
3192 : std::pair<unsigned short int, unsigned short int>
3193 0 : Elem::second_order_child_vertex (const unsigned int) const
3194 : {
3195 : // for linear elements, always return 0
3196 0 : return std::pair<unsigned short int, unsigned short int>(0,0);
3197 : }
3198 :
3199 :
3200 :
3201 272662 : ElemType Elem::first_order_equivalent_type (const ElemType et)
3202 : {
3203 26532 : switch (et)
3204 : {
3205 28 : case NODEELEM:
3206 28 : return NODEELEM;
3207 3118 : case EDGE2:
3208 : case EDGE3:
3209 : case EDGE4:
3210 3118 : return EDGE2;
3211 480 : case TRI3:
3212 : case TRI6:
3213 : case TRI7:
3214 480 : return TRI3;
3215 0 : case TRISHELL3:
3216 0 : return TRISHELL3;
3217 22338 : case QUAD4:
3218 : case QUAD8:
3219 : case QUAD9:
3220 22338 : return QUAD4;
3221 0 : case QUADSHELL4:
3222 : case QUADSHELL8:
3223 : case QUADSHELL9:
3224 0 : return QUADSHELL4;
3225 0 : case TET4:
3226 : case TET10:
3227 : case TET14:
3228 0 : return TET4;
3229 0 : case HEX8:
3230 : case HEX27:
3231 : case HEX20:
3232 0 : return HEX8;
3233 480 : case PRISM6:
3234 : case PRISM15:
3235 : case PRISM18:
3236 : case PRISM20:
3237 : case PRISM21:
3238 480 : return PRISM6;
3239 0 : case PYRAMID5:
3240 : case PYRAMID13:
3241 : case PYRAMID14:
3242 : case PYRAMID18:
3243 0 : return PYRAMID5;
3244 :
3245 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
3246 :
3247 16 : case INFEDGE2:
3248 16 : return INFEDGE2;
3249 72 : case INFQUAD4:
3250 : case INFQUAD6:
3251 72 : return INFQUAD4;
3252 0 : case INFHEX8:
3253 : case INFHEX16:
3254 : case INFHEX18:
3255 0 : return INFHEX8;
3256 0 : case INFPRISM6:
3257 : case INFPRISM12:
3258 0 : return INFPRISM6;
3259 :
3260 : #endif
3261 :
3262 0 : default:
3263 : // unknown element
3264 0 : return INVALID_ELEM;
3265 : }
3266 : }
3267 :
3268 :
3269 : // Most of our elements have the same topology on every side, but the
3270 : // prisms and pyramids have both triangular and quadrilateral faces, and
3271 : // a 2D or 3D infinite element's side 0 is its finite base while the rest
3272 : // are infinite. The polygons and polyhedra have one side type but no
3273 : // fixed side count, so they can answer here too while their side index
3274 : // goes unchecked; ask an actual Elem when you have one.
3275 99551 : ElemType Elem::side_type (const ElemType t,
3276 : const unsigned int s)
3277 : {
3278 50961 : libmesh_assert_less (s, type_to_n_sides_map[t]);
3279 :
3280 99551 : switch (t)
3281 : {
3282 36 : case EDGE2:
3283 : case EDGE3:
3284 : case EDGE4:
3285 36 : return NODEELEM;
3286 3710 : case TRI3:
3287 : case TRISHELL3:
3288 : case QUAD4:
3289 : case QUADSHELL4:
3290 : // A polygon's sides are all edges and a polyhedron's are all
3291 : // polygons, however many a given element turns out to have
3292 : case C0POLYGON:
3293 3710 : return EDGE2;
3294 5712 : case TRI6:
3295 : case TRI7:
3296 : case QUAD8:
3297 : case QUADSHELL8:
3298 : case QUAD9:
3299 : case QUADSHELL9:
3300 5712 : return EDGE3;
3301 6144 : case TET4:
3302 6144 : return TRI3;
3303 10752 : case TET10:
3304 10752 : return TRI6;
3305 12288 : case TET14:
3306 12288 : return TRI7;
3307 480 : case HEX8:
3308 480 : return QUAD4;
3309 864 : case HEX20:
3310 864 : return QUAD8;
3311 960 : case HEX27:
3312 960 : return QUAD9;
3313 : // Sides 0 and 4 of a prism are its triangles
3314 1536 : case PRISM6:
3315 1536 : return (s == 0 || s == 4) ? TRI3 : QUAD4;
3316 2112 : case PRISM15:
3317 2112 : return (s == 0 || s == 4) ? TRI6 : QUAD8;
3318 2208 : case PRISM18:
3319 2208 : return (s == 0 || s == 4) ? TRI6 : QUAD9;
3320 4544 : case PRISM20:
3321 : case PRISM21:
3322 4544 : return (s == 0 || s == 4) ? TRI7 : QUAD9;
3323 : // Side 4 of a pyramid is its quadrilateral base
3324 4416 : case PYRAMID5:
3325 4416 : return (s < 4) ? TRI3 : QUAD4;
3326 5952 : case PYRAMID13:
3327 5952 : return (s < 4) ? TRI6 : QUAD8;
3328 6048 : case PYRAMID14:
3329 6048 : return (s < 4) ? TRI6 : QUAD9;
3330 6432 : case PYRAMID18:
3331 6432 : return (s < 4) ? TRI7 : QUAD9;
3332 72 : case C0POLYHEDRON:
3333 72 : return C0POLYGON;
3334 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
3335 : // An InfEdge2's sides are its end nodes, like any other 1D element's
3336 0 : case INFEDGE2:
3337 0 : return NODEELEM;
3338 : // A 2D or 3D infinite element's side 0 is the finite base it was
3339 : // built from; its remaining sides run out to infinity with it
3340 9 : case INFQUAD4:
3341 9 : return (s == 0) ? EDGE2 : INFEDGE2;
3342 9 : case INFQUAD6:
3343 9 : return (s == 0) ? EDGE3 : INFEDGE2;
3344 15 : case INFHEX8:
3345 15 : return (s == 0) ? QUAD4 : INFQUAD4;
3346 15 : case INFHEX16:
3347 15 : return (s == 0) ? QUAD8 : INFQUAD6;
3348 15 : case INFHEX18:
3349 15 : return (s == 0) ? QUAD9 : INFQUAD6;
3350 12 : case INFPRISM6:
3351 12 : return (s == 0) ? TRI3 : INFQUAD4;
3352 12 : case INFPRISM12:
3353 12 : return (s == 0) ? TRI6 : INFQUAD6;
3354 : #endif
3355 0 : default:
3356 0 : libmesh_error_msg("No side type for element type " << Utility::enum_to_string(t));
3357 : }
3358 :
3359 : return INVALID_ELEM;
3360 : }
3361 :
3362 :
3363 :
3364 : // This reads the same constexpr side_nodes_map tables the element
3365 : // classes use for their virtual local_side_node(), so that code which
3366 : // knows an element type but has no Elem to call a virtual function on --
3367 : // building a side of an element that doesn't exist yet, or a device
3368 : // kernel that has only the type -- can still get at the reference
3369 : // element topology.
3370 : //
3371 : // The polygons and polyhedra have no such table to read, and the
3372 : // infinite elements' maps have no caller here yet, so both error out:
3373 : // asking for a topology we can't answer is a programming error rather
3374 : // than something to signal with a return value.
3375 : // A finite element's edges all have the same type, unlike its sides,
3376 : // which differ on the prisms and pyramids, so this takes no edge index.
3377 : // The infinite elements are the exception -- the edges of their finite
3378 : // base are not the ones running out to infinity -- and like the rest of
3379 : // the static topology lookups they are left to their virtual overrides.
3380 182252 : ElemType Elem::edge_type (const ElemType t)
3381 : {
3382 182252 : switch (t)
3383 : {
3384 : // 1D elements have no edges
3385 0 : case EDGE2:
3386 : case EDGE3:
3387 : case EDGE4:
3388 0 : return INVALID_ELEM;
3389 : // A 2D element's edges are its sides
3390 5660 : case TRI3:
3391 : case TRISHELL3:
3392 : case TRI6:
3393 : case TRI7:
3394 : case QUAD4:
3395 : case QUADSHELL4:
3396 : case QUAD8:
3397 : case QUADSHELL8:
3398 : case QUAD9:
3399 : case QUADSHELL9:
3400 : case C0POLYGON:
3401 5660 : return side_type(t, 0);
3402 : // A first-order 3D element's edges hold their two vertices
3403 49872 : case TET4:
3404 : case HEX8:
3405 : case PRISM6:
3406 : case PYRAMID5:
3407 : case C0POLYHEDRON:
3408 49872 : return EDGE2;
3409 : // and a second-order element's add the midpoint
3410 126720 : case TET10:
3411 : case TET14:
3412 : case HEX20:
3413 : case HEX27:
3414 : case PRISM15:
3415 : case PRISM18:
3416 : case PRISM20:
3417 : case PRISM21:
3418 : case PYRAMID13:
3419 : case PYRAMID14:
3420 : case PYRAMID18:
3421 126720 : return EDGE3;
3422 0 : default:
3423 0 : libmesh_error_msg("No edge type for element type " << Utility::enum_to_string(t));
3424 : }
3425 :
3426 : return INVALID_ELEM;
3427 : }
3428 :
3429 :
3430 :
3431 546336 : unsigned int Elem::local_side_node (const ElemType t,
3432 : const unsigned int side,
3433 : const unsigned int side_node)
3434 : {
3435 : // The nodes on a side are the nodes of the side's own element type, so
3436 : // we compose the two lookups here rather than tabulating side node
3437 : // counts a second time. That also keeps us inside the meaningful part
3438 : // of a row: the prisms' and pyramids' triangular sides carry fewer
3439 : // nodes than their quadrilateral ones, and their rows are padded out
3440 : // to the longer length
3441 45528 : libmesh_assert_less (side, type_to_n_sides_map[t]);
3442 45528 : libmesh_assert_less (side_node, type_to_n_nodes_map[side_type(t, side)]);
3443 :
3444 546336 : switch (t)
3445 : {
3446 : // A 1D element's sides are its end nodes
3447 24 : case EDGE2:
3448 : case EDGE3:
3449 : case EDGE4:
3450 24 : return side;
3451 : // The shell elements are numbered like the elements they shadow
3452 4608 : case TRI3:
3453 : case TRISHELL3:
3454 4608 : return Tri3::side_nodes_map[side][side_node];
3455 3456 : case TRI6:
3456 3456 : return Tri6::side_nodes_map[side][side_node];
3457 3456 : case TRI7:
3458 3456 : return Tri7::side_nodes_map[side][side_node];
3459 3072 : case QUAD4:
3460 : case QUADSHELL4:
3461 3072 : return Quad4::side_nodes_map[side][side_node];
3462 4608 : case QUAD8:
3463 : case QUADSHELL8:
3464 4608 : return Quad8::side_nodes_map[side][side_node];
3465 4608 : case QUAD9:
3466 : case QUADSHELL9:
3467 4608 : return Quad9::side_nodes_map[side][side_node];
3468 55296 : case TET4:
3469 55296 : return Tet4::side_nodes_map[side][side_node];
3470 110592 : case TET10:
3471 110592 : return Tet10::side_nodes_map[side][side_node];
3472 129024 : case TET14:
3473 129024 : return Tet14::side_nodes_map[side][side_node];
3474 4608 : case HEX8:
3475 4608 : return Hex8::side_nodes_map[side][side_node];
3476 9216 : case HEX20:
3477 9216 : return Hex20::side_nodes_map[side][side_node];
3478 10368 : case HEX27:
3479 10368 : return Hex27::side_nodes_map[side][side_node];
3480 6912 : case PRISM6:
3481 6912 : return Prism6::side_nodes_map[side][side_node];
3482 13824 : case PRISM15:
3483 13824 : return Prism15::side_nodes_map[side][side_node];
3484 14976 : case PRISM18:
3485 14976 : return Prism18::side_nodes_map[side][side_node];
3486 15744 : case PRISM20:
3487 15744 : return Prism20::side_nodes_map[side][side_node];
3488 15744 : case PRISM21:
3489 15744 : return Prism21::side_nodes_map[side][side_node];
3490 18432 : case PYRAMID5:
3491 18432 : return Pyramid5::side_nodes_map[side][side_node];
3492 36864 : case PYRAMID13:
3493 36864 : return Pyramid13::side_nodes_map[side][side_node];
3494 38016 : case PYRAMID14:
3495 38016 : return Pyramid14::side_nodes_map[side][side_node];
3496 42624 : case PYRAMID18:
3497 42624 : return Pyramid18::side_nodes_map[side][side_node];
3498 0 : default:
3499 0 : libmesh_error_msg("No static side node map for element type " << Utility::enum_to_string(t));
3500 : }
3501 :
3502 : return invalid_uint;
3503 : }
3504 :
3505 :
3506 :
3507 : // A 2D element's edges are its sides -- Face::local_edge_node() defines
3508 : // the two to be the same thing -- so we answer those from the side map
3509 : // and keep one copy of the 2D numbering. 1D elements have no edges to
3510 : // ask about, and like local_side_node() above we error rather than
3511 : // return a flag for a type we have no table for.
3512 401280 : unsigned int Elem::local_edge_node (const ElemType t,
3513 : const unsigned int edge,
3514 : const unsigned int edge_node)
3515 : {
3516 33440 : libmesh_assert_less (edge, type_to_n_edges_map[t]);
3517 33440 : libmesh_assert_less (edge_node, type_to_n_nodes_map[edge_type(t)]);
3518 :
3519 401280 : switch (t)
3520 : {
3521 11904 : case TRI3:
3522 : case TRISHELL3:
3523 : case TRI6:
3524 : case TRI7:
3525 : case QUAD4:
3526 : case QUADSHELL4:
3527 : case QUAD8:
3528 : case QUADSHELL8:
3529 : case QUAD9:
3530 : case QUADSHELL9:
3531 11904 : return local_side_node(t, edge, edge_node);
3532 55296 : case TET4:
3533 55296 : return Tet4::edge_nodes_map[edge][edge_node];
3534 82944 : case TET10:
3535 82944 : return Tet10::edge_nodes_map[edge][edge_node];
3536 82944 : case TET14:
3537 82944 : return Tet14::edge_nodes_map[edge][edge_node];
3538 4608 : case HEX8:
3539 4608 : return Hex8::edge_nodes_map[edge][edge_node];
3540 6912 : case HEX20:
3541 6912 : return Hex20::edge_nodes_map[edge][edge_node];
3542 6912 : case HEX27:
3543 6912 : return Hex27::edge_nodes_map[edge][edge_node];
3544 6912 : case PRISM6:
3545 6912 : return Prism6::edge_nodes_map[edge][edge_node];
3546 10368 : case PRISM15:
3547 10368 : return Prism15::edge_nodes_map[edge][edge_node];
3548 10368 : case PRISM18:
3549 10368 : return Prism18::edge_nodes_map[edge][edge_node];
3550 10368 : case PRISM20:
3551 10368 : return Prism20::edge_nodes_map[edge][edge_node];
3552 10368 : case PRISM21:
3553 10368 : return Prism21::edge_nodes_map[edge][edge_node];
3554 18432 : case PYRAMID5:
3555 18432 : return Pyramid5::edge_nodes_map[edge][edge_node];
3556 27648 : case PYRAMID13:
3557 27648 : return Pyramid13::edge_nodes_map[edge][edge_node];
3558 27648 : case PYRAMID14:
3559 27648 : return Pyramid14::edge_nodes_map[edge][edge_node];
3560 27648 : case PYRAMID18:
3561 27648 : return Pyramid18::edge_nodes_map[edge][edge_node];
3562 0 : default:
3563 0 : libmesh_error_msg("No static edge node map for element type " << Utility::enum_to_string(t));
3564 : }
3565 :
3566 : return invalid_uint;
3567 : }
3568 :
3569 :
3570 :
3571 3516844 : ElemType Elem::second_order_equivalent_type (const ElemType et,
3572 : const bool full_ordered)
3573 : {
3574 3516844 : switch (et)
3575 : {
3576 2 : case NODEELEM:
3577 2 : return NODEELEM;
3578 3806 : case EDGE2:
3579 : case EDGE3:
3580 : {
3581 : // full_ordered not relevant
3582 3806 : return EDGE3;
3583 : }
3584 :
3585 0 : case EDGE4:
3586 : {
3587 : // full_ordered not relevant
3588 0 : return EDGE4;
3589 : }
3590 :
3591 21150 : case TRI3:
3592 : case TRI6:
3593 : {
3594 : // full_ordered not relevant
3595 21150 : return TRI6;
3596 : }
3597 :
3598 0 : case TRI7:
3599 0 : return TRI7;
3600 :
3601 : // Currently there is no TRISHELL6, so similarly to other types
3602 : // where this is the case, we just return the input.
3603 0 : case TRISHELL3:
3604 0 : return TRISHELL3;
3605 :
3606 33454 : case QUAD4:
3607 : case QUAD8:
3608 : {
3609 33454 : if (full_ordered)
3610 2590 : return QUAD9;
3611 : else
3612 466 : return QUAD8;
3613 : }
3614 :
3615 8802 : case QUADSHELL4:
3616 : case QUADSHELL8:
3617 : {
3618 8802 : if (full_ordered)
3619 896 : return QUADSHELL9;
3620 : else
3621 448 : return QUADSHELL8;
3622 : }
3623 :
3624 258 : case QUAD9:
3625 : {
3626 : // full_ordered not relevant
3627 258 : return QUAD9;
3628 : }
3629 :
3630 2 : case QUADSHELL9:
3631 : {
3632 : // full_ordered not relevant
3633 2 : return QUADSHELL9;
3634 : }
3635 :
3636 1802795 : case TET4:
3637 : case TET10:
3638 : {
3639 : // full_ordered not relevant
3640 1802795 : return TET10;
3641 : }
3642 :
3643 0 : case TET14:
3644 0 : return TET14;
3645 :
3646 662376 : case HEX8:
3647 : case HEX20:
3648 : {
3649 : // see below how this correlates with INFHEX8
3650 662376 : if (full_ordered)
3651 43424 : return HEX27;
3652 : else
3653 17090 : return HEX20;
3654 : }
3655 :
3656 460 : case HEX27:
3657 : {
3658 : // full_ordered not relevant
3659 460 : return HEX27;
3660 : }
3661 :
3662 456152 : case PRISM6:
3663 : case PRISM15:
3664 : {
3665 456152 : if (full_ordered)
3666 29598 : return PRISM18;
3667 : else
3668 14786 : return PRISM15;
3669 : }
3670 :
3671 : // full_ordered not relevant, already fully second order
3672 12 : case PRISM18:
3673 12 : return PRISM18;
3674 0 : case PRISM20:
3675 0 : return PRISM20;
3676 0 : case PRISM21:
3677 0 : return PRISM21;
3678 0 : case PYRAMID18:
3679 0 : return PYRAMID18;
3680 :
3681 518939 : case PYRAMID5:
3682 : case PYRAMID13:
3683 : {
3684 518939 : if (full_ordered)
3685 7322 : return PYRAMID14;
3686 : else
3687 14592 : return PYRAMID13;
3688 : }
3689 :
3690 0 : case PYRAMID14:
3691 : {
3692 : // full_ordered not relevant
3693 0 : return PYRAMID14;
3694 : }
3695 :
3696 :
3697 :
3698 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
3699 :
3700 : // infinite elements
3701 0 : case INFEDGE2:
3702 : {
3703 0 : return INFEDGE2;
3704 : }
3705 :
3706 5 : case INFQUAD4:
3707 : case INFQUAD6:
3708 : {
3709 : // full_ordered not relevant
3710 5 : return INFQUAD6;
3711 : }
3712 :
3713 1032 : case INFHEX8:
3714 : case INFHEX16:
3715 : {
3716 : /*
3717 : * Note that this matches with \p Hex8:
3718 : * For full-ordered, \p InfHex18 and \p Hex27
3719 : * belong together, and for not full-ordered,
3720 : * \p InfHex16 and \p Hex20 belong together.
3721 : */
3722 1032 : if (full_ordered)
3723 452 : return INFHEX18;
3724 : else
3725 226 : return INFHEX16;
3726 : }
3727 :
3728 4 : case INFHEX18:
3729 : {
3730 : // full_ordered not relevant
3731 4 : return INFHEX18;
3732 : }
3733 :
3734 5 : case INFPRISM6:
3735 : case INFPRISM12:
3736 : {
3737 : // full_ordered not relevant
3738 5 : return INFPRISM12;
3739 : }
3740 :
3741 : #endif
3742 :
3743 :
3744 0 : default:
3745 : {
3746 : // what did we miss?
3747 0 : libmesh_error_msg("No second order equivalent element type for et = "
3748 : << Utility::enum_to_string(et));
3749 : }
3750 : }
3751 : }
3752 :
3753 :
3754 :
3755 5741251 : ElemType Elem::complete_order_equivalent_type (const ElemType et)
3756 : {
3757 5741251 : switch (et)
3758 : {
3759 0 : case NODEELEM:
3760 0 : return NODEELEM;
3761 1278 : case EDGE2:
3762 : case EDGE3:
3763 1278 : return EDGE3;
3764 :
3765 0 : case EDGE4:
3766 0 : return EDGE4;
3767 :
3768 7086 : case TRI3:
3769 : case TRI6:
3770 : case TRI7:
3771 7086 : return TRI7;
3772 :
3773 : // Currently there is no TRISHELL6, so similarly to other types
3774 : // where this is the case, we just return the input.
3775 0 : case TRISHELL3:
3776 0 : return TRISHELL3;
3777 :
3778 900 : case QUAD4:
3779 : case QUAD8:
3780 : case QUAD9:
3781 900 : return QUAD9;
3782 :
3783 0 : case QUADSHELL4:
3784 : case QUADSHELL8:
3785 : case QUADSHELL9:
3786 0 : return QUADSHELL9;
3787 :
3788 5477767 : case TET4:
3789 : case TET10:
3790 : case TET14:
3791 5477767 : return TET14;
3792 :
3793 316 : case HEX8:
3794 : case HEX20:
3795 : case HEX27:
3796 316 : return HEX27;
3797 :
3798 : // We don't strictly need the 21st node for DoFs, but some code
3799 : // depends on it for e.g. refinement patterns
3800 33449 : case PRISM6:
3801 : case PRISM15:
3802 : case PRISM18:
3803 : case PRISM20:
3804 : case PRISM21:
3805 33449 : return PRISM21;
3806 :
3807 220455 : case PYRAMID5:
3808 : case PYRAMID13:
3809 : case PYRAMID14:
3810 : case PYRAMID18:
3811 220455 : return PYRAMID18;
3812 :
3813 :
3814 :
3815 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
3816 0 : case INFEDGE2:
3817 0 : return INFEDGE2;
3818 :
3819 0 : case INFQUAD4:
3820 : case INFQUAD6:
3821 0 : return INFQUAD6;
3822 :
3823 0 : case INFHEX8:
3824 : case INFHEX16:
3825 : case INFHEX18:
3826 0 : return INFHEX18;
3827 :
3828 : // Probably ought to implement INFPRISM13 and/or 14; until then
3829 : // we'll just return the not-complete-order ElemType, in
3830 : // accordance which what seems like bad precedent...
3831 0 : case INFPRISM6:
3832 : case INFPRISM12:
3833 0 : return INFPRISM12;
3834 : #endif // LIBMESH_ENABLE_INFINITE_ELEMENTS
3835 :
3836 0 : default:
3837 : {
3838 : // what did we miss?
3839 0 : libmesh_error_msg("No complete order equivalent element type for et = "
3840 : << Utility::enum_to_string(et));
3841 : }
3842 : }
3843 : }
3844 :
3845 :
3846 :
3847 0 : Elem::side_iterator Elem::boundary_sides_begin()
3848 : {
3849 0 : Predicates::BoundarySide<SideIter> bsp;
3850 0 : return side_iterator(this->_first_side(), this->_last_side(), bsp);
3851 : }
3852 :
3853 :
3854 :
3855 :
3856 0 : Elem::side_iterator Elem::boundary_sides_end()
3857 : {
3858 0 : Predicates::BoundarySide<SideIter> bsp;
3859 0 : return side_iterator(this->_last_side(), this->_last_side(), bsp);
3860 : }
3861 :
3862 :
3863 :
3864 :
3865 157903 : Real Elem::volume () const
3866 : {
3867 : // The default implementation builds a finite element of the correct
3868 : // order and sums up the JxW contributions. This can be expensive,
3869 : // so the various element types can overload this method and compute
3870 : // the volume more efficiently.
3871 157903 : const FEFamily mapping_family = FEMap::map_fe_type(*this);
3872 157903 : const FEType fe_type(this->default_order(), mapping_family);
3873 :
3874 157903 : std::unique_ptr<FEBase> fe (FEBase::build(this->dim(),
3875 175588 : fe_type));
3876 :
3877 : // Use a map with a negative Jacobian tolerance, in case we're asked
3878 : // for the net volume of a tangled element
3879 17685 : fe->get_fe_map().set_jacobian_tolerance(std::numeric_limits<Real>::lowest());
3880 :
3881 41928 : const std::vector<Real> & JxW = fe->get_JxW();
3882 :
3883 : // The default quadrature rule should integrate the mass matrix,
3884 : // thus it should be plenty to compute the area
3885 157903 : QGauss qrule (this->dim(), fe_type.default_quadrature_order());
3886 :
3887 157903 : fe->attach_quadrature_rule(&qrule);
3888 :
3889 157903 : fe->reinit(this);
3890 :
3891 17685 : Real vol=0.;
3892 6519764 : for (auto jxw : JxW)
3893 6361861 : vol += jxw;
3894 :
3895 175588 : return vol;
3896 :
3897 128003 : }
3898 :
3899 :
3900 :
3901 40845840 : BoundingBox Elem::loose_bounding_box () const
3902 : {
3903 40845840 : Point pmin = this->point(0);
3904 40845840 : Point pmax = pmin;
3905 :
3906 40845840 : unsigned int n_points = this->n_nodes();
3907 215185789 : for (unsigned int p=0; p != n_points; ++p)
3908 697359796 : for (unsigned d=0; d<LIBMESH_DIM; ++d)
3909 : {
3910 76107846 : const Point & pt = this->point(p);
3911 523019847 : if (pmin(d) > pt(d))
3912 34581741 : pmin(d) = pt(d);
3913 :
3914 523019847 : if (pmax(d) < pt(d))
3915 83026046 : pmax(d) = pt(d);
3916 : }
3917 :
3918 43950588 : return BoundingBox(pmin, pmax);
3919 : }
3920 :
3921 : Point
3922 9399 : Elem::side_vertex_average_normal(const unsigned int s) const
3923 : {
3924 9399 : unsigned int dim = this->dim();
3925 9399 : const std::unique_ptr<const Elem> face = this->build_side_ptr(s);
3926 : std::unique_ptr<libMesh::FEBase> fe(
3927 9717 : libMesh::FEBase::build(dim, libMesh::FEType(this->default_order())));
3928 318 : fe->add_p_level_in_reinit(false); // Yes, really use the default order
3929 795 : const std::vector<Point> & normals = fe->get_normals();
3930 9399 : std::vector<Point> ref_side_vertex_average_v = {face->reference_elem()->vertex_average()};
3931 9399 : fe->reinit(this, s, TOLERANCE, &ref_side_vertex_average_v);
3932 18798 : return normals[0];
3933 8763 : }
3934 :
3935 3431002 : bool Elem::is_vertex_on_parent(unsigned int c,
3936 : unsigned int n) const
3937 : {
3938 : #ifdef LIBMESH_ENABLE_AMR
3939 :
3940 3431002 : unsigned int my_n_vertices = this->n_vertices();
3941 17136072 : for (unsigned int n_parent = 0; n_parent != my_n_vertices;
3942 : ++n_parent)
3943 16890944 : if (this->node_ptr(n_parent) == this->child_ptr(c)->node_ptr(n))
3944 54947 : return true;
3945 256846 : return false;
3946 :
3947 : #else
3948 :
3949 : // No AMR?
3950 : libmesh_ignore(c,n);
3951 : libmesh_error_msg("ERROR: AMR disabled, how did we get here?");
3952 : return true;
3953 :
3954 : #endif
3955 : }
3956 :
3957 :
3958 :
3959 0 : unsigned int Elem::opposite_side(const unsigned int /*s*/) const
3960 : {
3961 : // If the subclass didn't rederive this, using it is an error
3962 0 : libmesh_not_implemented();
3963 : }
3964 :
3965 :
3966 :
3967 0 : unsigned int Elem::opposite_node(const unsigned int /*n*/,
3968 : const unsigned int /*s*/) const
3969 : {
3970 : // If the subclass didn't rederive this, using it is an error
3971 0 : libmesh_not_implemented();
3972 : }
3973 :
3974 :
3975 37659 : unsigned int Elem::center_node_on_side(const unsigned short libmesh_dbg_var(side)) const
3976 : {
3977 3160 : libmesh_assert_less (side, this->n_sides());
3978 37659 : return invalid_uint;
3979 : }
3980 :
3981 :
3982 188382 : void Elem::swap2boundarysides(unsigned short s1,
3983 : unsigned short s2,
3984 : BoundaryInfo * boundary_info) const
3985 : {
3986 18964 : std::vector<boundary_id_type> ids1, ids2;
3987 188382 : boundary_info->boundary_ids(this, s1, ids1);
3988 188382 : boundary_info->boundary_ids(this, s2, ids2);
3989 188382 : boundary_info->remove_side(this, s1);
3990 188382 : boundary_info->remove_side(this, s2);
3991 188382 : if (!ids1.empty())
3992 7713 : boundary_info->add_side(this, s2, ids1);
3993 188382 : if (!ids2.empty())
3994 4392 : boundary_info->add_side(this, s1, ids2);
3995 188382 : }
3996 :
3997 :
3998 370518 : void Elem::swap2boundaryedges(unsigned short e1,
3999 : unsigned short e2,
4000 : BoundaryInfo * boundary_info) const
4001 : {
4002 38012 : std::vector<boundary_id_type> ids1, ids2;
4003 370518 : boundary_info->edge_boundary_ids(this, e1, ids1);
4004 370518 : boundary_info->edge_boundary_ids(this, e2, ids2);
4005 370518 : boundary_info->remove_edge(this, e1);
4006 370518 : boundary_info->remove_edge(this, e2);
4007 370518 : if (!ids1.empty())
4008 0 : boundary_info->add_edge(this, e2, ids1);
4009 370518 : if (!ids2.empty())
4010 0 : boundary_info->add_edge(this, e1, ids2);
4011 370518 : }
4012 :
4013 : bool
4014 8106214 : Elem::is_internal(const unsigned int i) const
4015 : {
4016 8106214 : switch (this->dim())
4017 : {
4018 1 : case 0:
4019 1 : return false;
4020 :
4021 216 : case 1:
4022 216 : return !this->is_vertex(i);
4023 :
4024 7990185 : case 2:
4025 7990185 : return !this->is_vertex(i) && !this->is_edge(i);
4026 :
4027 115801 : case 3:
4028 115801 : return !this->is_vertex(i) && !this->is_edge(i) && !this->is_face(i);
4029 :
4030 0 : default:
4031 0 : libmesh_error_msg("impossible element dimension " << std::to_string(this->dim()));
4032 : return 0;
4033 : }
4034 : }
4035 :
4036 : bool
4037 830337916 : Elem::positive_edge_orientation(const unsigned int i) const
4038 : {
4039 65776489 : libmesh_assert_less (i, this->n_edges());
4040 :
4041 830337916 : return this->point(this->local_edge_node(i, 0)) >
4042 896114405 : this->point(this->local_edge_node(i, 1));
4043 : }
4044 :
4045 : bool
4046 348944750 : Elem::positive_face_orientation(const unsigned int i) const
4047 : {
4048 348944750 : return this->face_orientation(i) % 2;
4049 : }
4050 :
4051 :
4052 : unsigned int
4053 2089125692 : Elem::face_orientation(const unsigned int i) const
4054 : {
4055 169079362 : libmesh_assert_less (i, this->n_faces());
4056 :
4057 : // Get the number of vertices N of face i. Note that for 3d elements, i.e.
4058 : // elements for which this->n_faces() > 0, the number of vertices on any of
4059 : // its sides (or faces) is just the number of that face's sides (or edges).
4060 2089125692 : const unsigned int N = Elem::type_to_n_sides_map[this->side_type(i)];
4061 :
4062 2089125692 : const std::vector<unsigned int> nodes = this->nodes_on_side(i);
4063 :
4064 6778589630 : auto cmp = [&](const unsigned int & m, const unsigned int & n) -> bool
4065 10153343812 : { return this->point(m) < this->point(n); };
4066 :
4067 2089125692 : const unsigned int v = std::distance(nodes.begin(),
4068 2089125692 : std::min_element(nodes.begin(), nodes.begin() + N, cmp));
4069 :
4070 : // v is the position of the least vertex. Compare its neighbors just
4071 : // before and after it, cyclically: the face is positively oriented iff
4072 : // the lesser neighbor precedes v, i.e. the vertices run backward from
4073 : // the least vertex through the face's node map.
4074 2765443140 : return 2 * v + cmp(nodes[(v - 1 + N) % N], nodes[(v + 1) % N]);
4075 : }
4076 :
4077 : std::vector<unsigned int>
4078 1711184772 : Elem::oriented_face_vertices(const unsigned int i) const
4079 : {
4080 1711184772 : const unsigned int N = Elem::type_to_n_sides_map[this->side_type(i)];
4081 1711184772 : const unsigned int orientation = this->face_orientation(i);
4082 1711184772 : const unsigned int least = orientation / 2;
4083 :
4084 : // A positively oriented face has the lesser neighbor before the least
4085 : // vertex in its side node map
4086 1711184772 : const bool backward = orientation % 2;
4087 :
4088 1711184772 : std::vector<unsigned int> vertices(N);
4089 7674821190 : for (const auto k : make_range(N))
4090 5963636418 : vertices[k] = this->local_side_node(i, (least + (backward ? N - k : k)) % N);
4091 :
4092 1711184772 : return vertices;
4093 : }
4094 :
4095 : bool
4096 143208 : Elem::relative_edge_face_order(const unsigned int e, const unsigned int s) const
4097 : {
4098 0 : libmesh_assert_less (e, this->n_edges());
4099 0 : libmesh_assert_less (s, this->n_faces());
4100 0 : libmesh_assert (is_edge_on_side(e, s));
4101 :
4102 143208 : const unsigned int N = Elem::type_to_n_sides_map[this->side_type(s)];
4103 :
4104 : unsigned int v;
4105 255312 : for (v = 0; v < N; v++)
4106 255312 : if (this->local_side_node(s, v) == this->local_edge_node(e, 0))
4107 0 : break;
4108 :
4109 0 : libmesh_assert (this->local_side_node(s, (v + 1) % N) ==
4110 : this->local_edge_node(e, 1) ||
4111 : this->local_side_node(s, (v - 1 + N) % N) ==
4112 : this->local_edge_node(e, 1));
4113 :
4114 143208 : return this->local_side_node(s, (v + 1) % N) ==
4115 143208 : this->local_edge_node(e, 1);
4116 : }
4117 :
4118 : } // namespace libMesh
|
