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