Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
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 10050099 : std::unique_ptr<Elem> Elem::disconnected_clone() const
411 : {
412 10050099 : std::unique_ptr<Elem> returnval;
413 :
414 10050099 : switch (this->type())
415 : {
416 1887 : case C0POLYGON:
417 1887 : returnval = std::make_unique<C0Polygon>(this->n_sides());
418 1887 : break;
419 :
420 10048212 : default:
421 19776778 : returnval = Elem::build(this->type());
422 : }
423 :
424 10050099 : returnval->set_id() = this->id();
425 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
426 10050099 : if (this->valid_unique_id())
427 319648 : returnval->set_unique_id(this->unique_id());
428 : #endif
429 :
430 10050099 : const auto n_elem_ints = this->n_extra_integers();
431 10050099 : returnval->add_extra_integers(n_elem_ints);
432 10067352 : for (unsigned int i = 0; i != n_elem_ints; ++i)
433 17253 : returnval->set_extra_integer(i, this->get_extra_integer(i));
434 :
435 9741217 : returnval->inherit_data_from(*this);
436 :
437 10050099 : return returnval;
438 0 : }
439 :
440 :
441 :
442 64075550 : std::unique_ptr<Elem> Elem::build(const ElemType type,
443 : Elem * p)
444 : {
445 64075550 : switch (type)
446 : {
447 : // 0D elements
448 212161 : case NODEELEM:
449 212161 : return std::make_unique<NodeElem>(p);
450 :
451 : // 1D elements
452 496055 : case EDGE2:
453 496055 : return std::make_unique<Edge2>(p);
454 227122 : case EDGE3:
455 227122 : return std::make_unique<Edge3>(p);
456 45552 : case EDGE4:
457 45552 : return std::make_unique<Edge4>(p);
458 :
459 : // 2D elements
460 4338395 : case TRI3:
461 4338395 : return std::make_unique<Tri3>(p);
462 66496 : case TRISHELL3:
463 66496 : return std::make_unique<TriShell3>(p);
464 60799 : case TRI3SUBDIVISION:
465 60799 : return std::make_unique<Tri3Subdivision>(p);
466 866808 : case TRI6:
467 866808 : return std::make_unique<Tri6>(p);
468 569770 : case TRI7:
469 569770 : return std::make_unique<Tri7>(p);
470 25358150 : case QUAD4:
471 25358150 : return std::make_unique<Quad4>(p);
472 181374 : case QUADSHELL4:
473 181374 : return std::make_unique<QuadShell4>(p);
474 258021 : case QUAD8:
475 258021 : return std::make_unique<Quad8>(p);
476 176074 : case QUADSHELL8:
477 176074 : return std::make_unique<QuadShell8>(p);
478 3097338 : case QUAD9:
479 3097338 : return std::make_unique<Quad9>(p);
480 61081 : case QUADSHELL9:
481 61081 : return std::make_unique<QuadShell9>(p);
482 :
483 : // Well, a hexagon is *a* polygon...
484 1207 : case C0POLYGON:
485 1207 : 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 15532691 : case TET4:
495 15532691 : return std::make_unique<Tet4>(p);
496 1562682 : case TET10:
497 1562682 : return std::make_unique<Tet10>(p);
498 3332168 : case TET14:
499 3332168 : return std::make_unique<Tet14>(p);
500 3229492 : case HEX8:
501 3229492 : return std::make_unique<Hex8>(p);
502 127988 : case HEX20:
503 127988 : return std::make_unique<Hex20>(p);
504 1859382 : case HEX27:
505 1859382 : return std::make_unique<Hex27>(p);
506 361097 : case PRISM6:
507 361097 : return std::make_unique<Prism6>(p);
508 166622 : case PRISM15:
509 166622 : return std::make_unique<Prism15>(p);
510 612956 : case PRISM18:
511 612956 : return std::make_unique<Prism18>(p);
512 111507 : case PRISM20:
513 111507 : return std::make_unique<Prism20>(p);
514 245769 : case PRISM21:
515 245769 : return std::make_unique<Prism21>(p);
516 468517 : case PYRAMID5:
517 468517 : return std::make_unique<Pyramid5>(p);
518 153864 : case PYRAMID13:
519 153864 : return std::make_unique<Pyramid13>(p);
520 155928 : case PYRAMID14:
521 155928 : return std::make_unique<Pyramid14>(p);
522 131071 : case PYRAMID18:
523 131071 : 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 140 : case INFQUAD4:
532 140 : return std::make_unique<InfQuad4>(p);
533 140 : case INFQUAD6:
534 140 : return std::make_unique<InfQuad6>(p);
535 :
536 : // 3D infinite elements
537 2986 : case INFHEX8:
538 2986 : return std::make_unique<InfHex8>(p);
539 220 : case INFHEX16:
540 220 : return std::make_unique<InfHex16>(p);
541 1667 : case INFHEX18:
542 1667 : return std::make_unique<InfHex18>(p);
543 220 : case INFPRISM6:
544 220 : return std::make_unique<InfPrism6>(p);
545 2040 : case INFPRISM12:
546 2040 : 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 6829510 : 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 6829510 : auto temp = Elem::build(type, nullptr);
562 192767 : temp->set_id(id);
563 6829510 : return temp;
564 : }
565 :
566 :
567 :
568 139184 : const Elem * Elem::reference_elem () const
569 : {
570 139184 : 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 499304125 : Point Elem::vertex_average() const
670 : {
671 39654228 : Point cp;
672 :
673 499304125 : const auto n_vertices = this->n_vertices();
674 :
675 2923921456 : for (unsigned int n=0; n<n_vertices; n++)
676 388132734 : cp.add (this->point(n));
677 :
678 499304125 : return (cp /= static_cast<Real>(n_vertices));
679 : }
680 :
681 :
682 :
683 1840880 : Real Elem::hmin() const
684 : {
685 1840880 : Real h_min=std::numeric_limits<Real>::max();
686 :
687 : // Avoid calling a virtual a lot of times
688 1840880 : const auto n_vertices = this->n_vertices();
689 :
690 8774774 : for (unsigned int n_outer=0; n_outer<n_vertices; n_outer++)
691 17292877 : for (unsigned int n_inner=n_outer+1; n_inner<n_vertices; n_inner++)
692 : {
693 2015176 : const auto diff = (this->point(n_outer) - this->point(n_inner));
694 :
695 13500684 : h_min = std::min(h_min, diff.norm_sq());
696 : }
697 :
698 2021034 : return std::sqrt(h_min);
699 : }
700 :
701 :
702 :
703 93579744 : Real Elem::hmax() const
704 : {
705 93579744 : Real h_max=0;
706 :
707 : // Avoid calling a virtual a lot of times
708 93579744 : const auto n_vertices = this->n_vertices();
709 :
710 631483510 : for (unsigned int n_outer=0; n_outer<n_vertices; n_outer++)
711 2064450326 : for (unsigned int n_inner=n_outer+1; n_inner<n_vertices; n_inner++)
712 : {
713 96763450 : const auto diff = (this->point(n_outer) - this->point(n_inner));
714 :
715 1746473785 : h_max = std::max(h_max, diff.norm_sq());
716 : }
717 :
718 99422403 : return std::sqrt(h_max);
719 : }
720 :
721 :
722 :
723 0 : Real Elem::length(const unsigned int n1,
724 : const unsigned int n2) const
725 : {
726 0 : libmesh_assert_less ( n1, this->n_vertices() );
727 0 : libmesh_assert_less ( n2, this->n_vertices() );
728 :
729 0 : return (this->point(n1) - this->point(n2)).norm();
730 : }
731 :
732 :
733 :
734 0 : dof_id_type Elem::key () const
735 : {
736 0 : const unsigned short n_n = this->n_nodes();
737 :
738 : std::array<dof_id_type, Elem::max_n_nodes> node_ids;
739 :
740 0 : for (unsigned short n=0; n != n_n; ++n)
741 0 : node_ids[n] = this->node_id(n);
742 :
743 : // Always sort, so that different local node numberings hash to the
744 : // same value.
745 0 : std::sort (node_ids.begin(), node_ids.begin()+n_n);
746 :
747 0 : return Utility::hashword(node_ids.data(), n_n);
748 : }
749 :
750 :
751 :
752 122994632 : bool Elem::operator == (const Elem & rhs) const
753 : {
754 : // If the elements aren't the same type, they aren't equal
755 122994632 : if (this->type() != rhs.type())
756 0 : return false;
757 :
758 122994632 : const unsigned short n_n = this->n_nodes();
759 2557176 : libmesh_assert_equal_to(n_n, rhs.n_nodes());
760 :
761 : // Make two sorted arrays of global node ids and compare them for
762 : // equality.
763 : std::array<dof_id_type, Elem::max_n_nodes> this_ids, rhs_ids;
764 :
765 421368881 : for (unsigned short n = 0; n != n_n; n++)
766 : {
767 298374249 : this_ids[n] = this->node_id(n);
768 305171907 : rhs_ids[n] = rhs.node_id(n);
769 : }
770 :
771 : // Sort the vectors to rule out different local node numberings.
772 122994632 : std::sort(this_ids.begin(), this_ids.begin()+n_n);
773 122994632 : std::sort(rhs_ids.begin(), rhs_ids.begin()+n_n);
774 :
775 : // If the node ids match, the elements are equal!
776 421368881 : for (unsigned short n = 0; n != n_n; ++n)
777 298374249 : if (this_ids[n] != rhs_ids[n])
778 0 : return false;
779 2557176 : return true;
780 : }
781 :
782 :
783 :
784 523191 : bool Elem::topologically_equal (const Elem & rhs) const
785 : {
786 : // If the elements aren't the same type, they aren't equal
787 523191 : if (this->type() != rhs.type())
788 0 : return false;
789 :
790 330964 : libmesh_assert_equal_to(this->n_nodes(), rhs.n_nodes());
791 :
792 3615907 : for (auto n : make_range(this->n_nodes()))
793 3201436 : if (this->node_id(n) != rhs.node_id(n))
794 0 : return false;
795 :
796 2544152 : for (auto neigh : make_range(this->n_neighbors()))
797 : {
798 2046412 : if (!this->neighbor_ptr(neigh))
799 : {
800 212960 : if (rhs.neighbor_ptr(neigh))
801 6 : return false;
802 204408 : continue;
803 : }
804 3056234 : if (!rhs.neighbor_ptr(neigh) ||
805 1243810 : this->neighbor_ptr(neigh)->id() !=
806 1243810 : rhs.neighbor_ptr(neigh)->id())
807 6 : return false;
808 : }
809 :
810 530332 : if (this->parent())
811 : {
812 32220 : if (!rhs.parent())
813 0 : return false;
814 31920 : if (this->parent()->id() != rhs.parent()->id())
815 0 : return false;
816 : }
817 498112 : else if (rhs.parent())
818 0 : return false;
819 :
820 523050 : if (this->interior_parent())
821 : {
822 0 : if (!rhs.interior_parent())
823 0 : return false;
824 0 : if (this->interior_parent()->id() !=
825 0 : rhs.interior_parent()->id())
826 0 : return false;
827 : }
828 523050 : else if (rhs.interior_parent())
829 0 : return false;
830 :
831 330958 : return true;
832 : }
833 :
834 :
835 :
836 0 : bool Elem::is_semilocal(const processor_id_type my_pid) const
837 : {
838 0 : std::set<const Elem *> point_neighbors;
839 :
840 0 : this->find_point_neighbors(point_neighbors);
841 :
842 0 : for (const auto & elem : point_neighbors)
843 0 : if (elem->processor_id() == my_pid)
844 0 : return true;
845 :
846 0 : return false;
847 : }
848 :
849 :
850 :
851 7709 : unsigned int Elem::which_side_am_i (const Elem * e) const
852 : {
853 708 : libmesh_assert(e);
854 :
855 7709 : const unsigned int ns = this->n_sides();
856 7709 : const unsigned int nn = this->n_nodes();
857 :
858 7709 : const unsigned int en = e->n_nodes();
859 :
860 : // e might be on any side until proven otherwise
861 8417 : std::vector<bool> might_be_side(ns, true);
862 :
863 30670 : for (unsigned int i=0; i != en; ++i)
864 : {
865 25069 : Point side_point = e->point(i);
866 2108 : unsigned int local_node_id = libMesh::invalid_uint;
867 :
868 : // Look for a node of this that's contiguous with node i of
869 : // e. Note that the exact floating point comparison of Point
870 : // positions is intentional, see the class documentation for
871 : // this function.
872 229112 : for (unsigned int j=0; j != nn; ++j)
873 37848 : if (this->point(j) == side_point)
874 2108 : local_node_id = j;
875 :
876 : // If a node of e isn't contiguous with some node of this, then
877 : // e isn't a side of this.
878 22961 : if (local_node_id == libMesh::invalid_uint)
879 0 : return libMesh::invalid_uint;
880 :
881 : // If a node of e isn't contiguous with some node on side s of
882 : // this, then e isn't on side s.
883 114639 : for (unsigned int s=0; s != ns; ++s)
884 91678 : if (!this->is_node_on_side(local_node_id, s))
885 9816 : might_be_side[s] = false;
886 : }
887 :
888 21383 : for (unsigned int s=0; s != ns; ++s)
889 23314 : if (might_be_side[s])
890 : {
891 : #ifdef DEBUG
892 1593 : for (unsigned int s2=s+1; s2 < ns; ++s2)
893 885 : libmesh_assert (!might_be_side[s2]);
894 : #endif
895 7709 : return s;
896 : }
897 :
898 : // Didn't find any matching side
899 0 : return libMesh::invalid_uint;
900 : }
901 :
902 :
903 :
904 1258090086 : bool Elem::contains_vertex_of(const Elem * e, bool mesh_connection) const
905 : {
906 : // Our vertices are the first numbered nodes
907 1258090086 : const unsigned int nv = e->n_vertices();
908 1258090086 : const unsigned int my_nv = this->n_vertices();
909 :
910 : // Check for vertex-to-vertex containment first; contains_point() is
911 : // expensive
912 4335354711 : for (auto n : make_range(nv))
913 : {
914 14948611 : const Node * vertex = e->node_ptr(n);
915 20624886123 : for (auto my_n : make_range(my_nv))
916 17547621498 : if (&this->node_ref(my_n) == vertex)
917 1199291 : return true;
918 : }
919 :
920 : // If e is in our mesh, then we might be done testing
921 425601024 : if (mesh_connection)
922 : {
923 425601024 : const unsigned int l = this->level();
924 425601024 : const unsigned int el = e->level();
925 :
926 425601024 : if (l >= el)
927 706150 : return false;
928 :
929 : // We could also return false for l==el-1 iff we knew we had no
930 : // triangular faces, but we don't have an API to check that.
931 : }
932 :
933 : // Our vertices are the first numbered nodes
934 48834830 : for (auto n : make_range(nv))
935 42748900 : if (this->contains_point(e->point(n)))
936 0 : return true;
937 1632 : return false;
938 : }
939 :
940 :
941 :
942 0 : bool Elem::contains_edge_of(const Elem * e) const
943 : {
944 0 : unsigned int num_contained_edges = 0;
945 :
946 : // Our vertices are the first numbered nodes
947 0 : for (auto n : make_range(e->n_vertices()))
948 : {
949 0 : if (this->contains_point(e->point(n)))
950 : {
951 0 : num_contained_edges++;
952 0 : if (num_contained_edges>=2)
953 : {
954 0 : return true;
955 : }
956 : }
957 : }
958 0 : return false;
959 : }
960 :
961 :
962 :
963 575791 : void Elem::find_point_neighbors(const Point & p,
964 : std::set<const Elem *> & neighbor_set) const
965 : {
966 51170 : libmesh_assert(this->contains_point(p));
967 51170 : libmesh_assert(this->active());
968 :
969 51170 : neighbor_set.clear();
970 575791 : neighbor_set.insert(this);
971 :
972 102340 : std::set<const Elem *> untested_set, next_untested_set;
973 575791 : untested_set.insert(this);
974 :
975 : #ifdef LIBMESH_ENABLE_AMR
976 102340 : std::vector<const Elem *> active_neighbor_children;
977 : #endif // #ifdef LIBMESH_ENABLE_AMR
978 :
979 1737947 : while (!untested_set.empty())
980 : {
981 : // Loop over all the elements in the patch that haven't already
982 : // been tested
983 2473416 : for (const auto & elem : untested_set)
984 6920669 : for (auto current_neighbor : elem->neighbor_ptr_range())
985 : {
986 5922562 : if (current_neighbor &&
987 5043901 : current_neighbor != remote_elem) // we have a real neighbor on this side
988 : {
989 426766 : if (current_neighbor->active()) // ... if it is active
990 : {
991 422224 : auto it = neighbor_set.lower_bound(current_neighbor);
992 8811239 : if ((it == neighbor_set.end() || *it != current_neighbor) &&
993 3889370 : current_neighbor->contains_point(p)) // ... don't have and touches p
994 : {
995 : // Add it and test it
996 671768 : next_untested_set.insert(current_neighbor);
997 733669 : neighbor_set.emplace_hint(it, current_neighbor);
998 : }
999 : }
1000 : #ifdef LIBMESH_ENABLE_AMR
1001 : else // ... the neighbor is *not* active,
1002 : { // ... so add *all* neighboring
1003 : // active children that touch p
1004 4542 : active_neighbor_children.clear();
1005 : current_neighbor->active_family_tree_by_neighbor
1006 14099 : (active_neighbor_children, elem);
1007 :
1008 42611 : for (const auto & current_child : active_neighbor_children)
1009 : {
1010 9084 : auto it = neighbor_set.lower_bound(current_child);
1011 55800 : if ((it == neighbor_set.end() || *it != current_child) &&
1012 27288 : current_child->contains_point(p))
1013 : {
1014 : // Add it and test it
1015 1258 : next_untested_set.insert(current_child);
1016 1800 : neighbor_set.emplace_hint(it, current_child);
1017 : }
1018 : }
1019 : }
1020 : #endif // #ifdef LIBMESH_ENABLE_AMR
1021 : }
1022 : }
1023 101794 : untested_set.swap(next_untested_set);
1024 101794 : next_untested_set.clear();
1025 : }
1026 575791 : }
1027 :
1028 :
1029 :
1030 9407683 : void Elem::find_point_neighbors(std::set<const Elem *> & neighbor_set) const
1031 : {
1032 9407683 : this->find_point_neighbors(neighbor_set, this);
1033 9407683 : }
1034 :
1035 :
1036 :
1037 9407683 : void Elem::find_point_neighbors(std::set<const Elem *> & neighbor_set,
1038 : const Elem * start_elem) const
1039 : {
1040 9407683 : ElemInternal::find_point_neighbors(this, neighbor_set, start_elem);
1041 9407683 : }
1042 :
1043 :
1044 :
1045 0 : void Elem::find_point_neighbors(std::set<Elem *> & neighbor_set,
1046 : Elem * start_elem)
1047 : {
1048 0 : ElemInternal::find_point_neighbors(this, neighbor_set, start_elem);
1049 0 : }
1050 :
1051 :
1052 :
1053 480 : void Elem::find_edge_neighbors(const Point & p1,
1054 : const Point & p2,
1055 : std::set<const Elem *> & neighbor_set) const
1056 : {
1057 : // Simple but perhaps suboptimal code: find elements containing the
1058 : // first point, then winnow this set down by removing elements which
1059 : // don't also contain the second point
1060 :
1061 40 : libmesh_assert(this->contains_point(p2));
1062 480 : this->find_point_neighbors(p1, neighbor_set);
1063 :
1064 40 : std::set<const Elem *>::iterator it = neighbor_set.begin();
1065 40 : const std::set<const Elem *>::iterator end = neighbor_set.end();
1066 :
1067 1812 : while (it != end)
1068 : {
1069 : // As of C++11, set::erase returns an iterator to the element
1070 : // following the erased element, or end.
1071 1332 : if (!(*it)->contains_point(p2))
1072 781 : it = neighbor_set.erase(it);
1073 : else
1074 40 : ++it;
1075 : }
1076 480 : }
1077 :
1078 :
1079 :
1080 0 : void Elem::find_edge_neighbors(std::set<const Elem *> & neighbor_set) const
1081 : {
1082 0 : neighbor_set.clear();
1083 0 : neighbor_set.insert(this);
1084 :
1085 0 : std::set<const Elem *> untested_set, next_untested_set;
1086 0 : untested_set.insert(this);
1087 :
1088 0 : while (!untested_set.empty())
1089 : {
1090 : // Loop over all the elements in the patch that haven't already
1091 : // been tested
1092 0 : for (const auto & elem : untested_set)
1093 : {
1094 0 : for (auto current_neighbor : elem->neighbor_ptr_range())
1095 : {
1096 0 : if (current_neighbor &&
1097 0 : current_neighbor != remote_elem) // we have a real neighbor on this side
1098 : {
1099 0 : if (current_neighbor->active()) // ... if it is active
1100 : {
1101 0 : if (this->contains_edge_of(current_neighbor) // ... and touches us
1102 0 : || current_neighbor->contains_edge_of(this))
1103 : {
1104 : // Make sure we'll test it
1105 0 : if (!neighbor_set.count(current_neighbor))
1106 0 : next_untested_set.insert (current_neighbor);
1107 :
1108 : // And add it
1109 0 : neighbor_set.insert (current_neighbor);
1110 : }
1111 : }
1112 : #ifdef LIBMESH_ENABLE_AMR
1113 : else // ... the neighbor is *not* active,
1114 : { // ... so add *all* neighboring
1115 : // active children
1116 0 : std::vector<const Elem *> active_neighbor_children;
1117 :
1118 : current_neighbor->active_family_tree_by_neighbor
1119 0 : (active_neighbor_children, elem);
1120 :
1121 0 : for (const auto & current_child : active_neighbor_children)
1122 0 : if (this->contains_edge_of(current_child) || current_child->contains_edge_of(this))
1123 : {
1124 : // Make sure we'll test it
1125 0 : if (!neighbor_set.count(current_child))
1126 0 : next_untested_set.insert (current_child);
1127 :
1128 0 : neighbor_set.insert (current_child);
1129 : }
1130 : }
1131 : #endif // #ifdef LIBMESH_ENABLE_AMR
1132 : }
1133 : }
1134 : }
1135 0 : untested_set.swap(next_untested_set);
1136 0 : next_untested_set.clear();
1137 : }
1138 0 : }
1139 :
1140 :
1141 :
1142 5376 : void Elem::find_interior_neighbors(std::set<const Elem *> & neighbor_set) const
1143 : {
1144 5376 : ElemInternal::find_interior_neighbors(this, neighbor_set);
1145 5376 : }
1146 :
1147 :
1148 :
1149 6253 : void Elem::find_interior_neighbors(std::set<Elem *> & neighbor_set)
1150 : {
1151 6253 : ElemInternal::find_interior_neighbors(this, neighbor_set);
1152 6253 : }
1153 :
1154 :
1155 :
1156 75997965 : const Elem * Elem::interior_parent () const
1157 : {
1158 : // interior parents make no sense for full-dimensional elements.
1159 75997965 : if (this->dim() >= LIBMESH_DIM)
1160 819883 : return nullptr;
1161 :
1162 : // they USED TO BE only good for level-0 elements, but we now
1163 : // support keeping interior_parent() valid on refined boundary
1164 : // elements.
1165 : // if (this->level() != 0)
1166 : // return this->parent()->interior_parent();
1167 :
1168 : // We store the interior_parent pointer after both the parent
1169 : // neighbor and neighbor pointers
1170 43040477 : Elem * interior_p = _elemlinks[1+this->n_sides()];
1171 :
1172 : // If we have an interior_parent, we USED TO assume it was a
1173 : // one-higher-dimensional interior element, but we now allow e.g.
1174 : // edge elements to have a 3D interior_parent with no
1175 : // intermediate 2D element.
1176 : // libmesh_assert (!interior_p ||
1177 : // interior_p->dim() == (this->dim()+1));
1178 2162946 : libmesh_assert (!interior_p ||
1179 : (interior_p == remote_elem) ||
1180 : (interior_p->dim() > this->dim()));
1181 :
1182 : // If an element in a multi-dimensional mesh has an interior_parent
1183 : // link, it should be at our level or coarser, just like a neighbor
1184 : // link. Our collect_families() code relies on this, but it might
1185 : // be tempting for users to manually assign something that breaks
1186 : // it.
1187 : //
1188 : // However, we *also* create temporary side elements, and we don't
1189 : // bother with creating ancestors for those, so they can be at level
1190 : // 0 even when they're sides of non-level-0 elements.
1191 2162946 : libmesh_assert (!interior_p ||
1192 : (interior_p->level() <= this->level()) ||
1193 : (this->level() == 0 &&
1194 : this->id() == DofObject::invalid_id));
1195 :
1196 43040477 : return interior_p;
1197 : }
1198 :
1199 :
1200 :
1201 81933461 : Elem * Elem::interior_parent ()
1202 : {
1203 : // See the const version for comments
1204 81933461 : if (this->dim() >= LIBMESH_DIM)
1205 208206 : return nullptr;
1206 :
1207 70177943 : Elem * interior_p = _elemlinks[1+this->n_sides()];
1208 :
1209 2935316 : libmesh_assert (!interior_p ||
1210 : (interior_p == remote_elem) ||
1211 : (interior_p->dim() > this->dim()));
1212 :
1213 70177943 : return interior_p;
1214 : }
1215 :
1216 :
1217 :
1218 5248471697 : void Elem::set_interior_parent (Elem * p)
1219 : {
1220 : // interior parents make no sense for full-dimensional elements.
1221 490902960 : libmesh_assert (!p ||
1222 : this->dim() < LIBMESH_DIM);
1223 :
1224 : // If we have an interior_parent, we USED TO assume it was a
1225 : // one-higher-dimensional interior element, but we now allow e.g.
1226 : // edge elements to have a 3D interior_parent with no
1227 : // intermediate 2D element.
1228 : // libmesh_assert (!p ||
1229 : // p->dim() == (this->dim()+1));
1230 490902960 : libmesh_assert (!p ||
1231 : (p == remote_elem) ||
1232 : (p->dim() > this->dim()));
1233 :
1234 5248471697 : _elemlinks[1+this->n_sides()] = p;
1235 5248471697 : }
1236 :
1237 :
1238 :
1239 : #ifdef LIBMESH_ENABLE_PERIODIC
1240 :
1241 786435 : Elem * Elem::topological_neighbor (const unsigned int i,
1242 : MeshBase & mesh,
1243 : const PointLocatorBase & point_locator,
1244 : const PeriodicBoundaries * pb)
1245 : {
1246 488472 : libmesh_assert_less (i, this->n_neighbors());
1247 :
1248 687080 : Elem * neighbor_i = this->neighbor_ptr(i);
1249 786435 : if (neighbor_i != nullptr)
1250 476930 : return neighbor_i;
1251 :
1252 22252 : if (pb)
1253 : {
1254 : // Since the neighbor is nullptr it must be on a boundary. We need
1255 : // see if this is a periodic boundary in which case it will have a
1256 : // topological neighbor
1257 11542 : std::vector<boundary_id_type> bc_ids;
1258 22252 : mesh.get_boundary_info().boundary_ids(this, cast_int<unsigned short>(i), bc_ids);
1259 22252 : for (const auto & id : bc_ids)
1260 22252 : if (pb->boundary(id))
1261 : {
1262 : // Since the point locator inside of periodic boundaries
1263 : // returns a const pointer we will retrieve the proper
1264 : // pointer directly from the mesh object.
1265 22252 : const Elem * const cn = pb->neighbor(id, point_locator, this, i);
1266 11542 : neighbor_i = const_cast<Elem *>(cn);
1267 :
1268 : // Since coarse elements do not have more refined
1269 : // neighbors we need to make sure that we don't return one
1270 : // of these types of neighbors.
1271 22252 : if (neighbor_i)
1272 24573 : while (level() < neighbor_i->level())
1273 2026 : neighbor_i = neighbor_i->parent();
1274 11542 : return neighbor_i;
1275 : }
1276 : }
1277 :
1278 0 : return nullptr;
1279 : }
1280 :
1281 :
1282 :
1283 47450 : const Elem * Elem::topological_neighbor (const unsigned int i,
1284 : const MeshBase & mesh,
1285 : const PointLocatorBase & point_locator,
1286 : const PeriodicBoundaries * pb) const
1287 : {
1288 19968 : libmesh_assert_less (i, this->n_neighbors());
1289 :
1290 31168 : const Elem * neighbor_i = this->neighbor_ptr(i);
1291 47450 : if (neighbor_i != nullptr)
1292 11098 : return neighbor_i;
1293 :
1294 31230 : if (pb)
1295 : {
1296 : // Since the neighbor is nullptr it must be on a boundary. We need
1297 : // see if this is a periodic boundary in which case it will have a
1298 : // topological neighbor
1299 8870 : std::vector<boundary_id_type> bc_ids;
1300 31230 : mesh.get_boundary_info().boundary_ids(this, cast_int<unsigned short>(i), bc_ids);
1301 34152 : for (const auto & id : bc_ids)
1302 21754 : if (pb->boundary(id))
1303 : {
1304 18832 : neighbor_i = pb->neighbor(id, point_locator, this, i);
1305 :
1306 : // Since coarse elements do not have more refined
1307 : // neighbors we need to make sure that we don't return one
1308 : // of these types of neighbors.
1309 18832 : if (neighbor_i)
1310 20819 : while (level() < neighbor_i->level())
1311 1836 : neighbor_i = neighbor_i->parent();
1312 8134 : return neighbor_i;
1313 : }
1314 : }
1315 :
1316 736 : return nullptr;
1317 : }
1318 :
1319 :
1320 20006 : bool Elem::has_topological_neighbor (const Elem * elem,
1321 : const MeshBase & mesh,
1322 : const PointLocatorBase & point_locator,
1323 : const PeriodicBoundaries * pb) const
1324 : {
1325 : // First see if this is a normal "interior" neighbor
1326 20006 : if (has_neighbor(elem))
1327 0 : return true;
1328 :
1329 20056 : for (auto n : this->side_index_range())
1330 20056 : if (this->topological_neighbor(n, mesh, point_locator, pb))
1331 13942 : return true;
1332 :
1333 0 : return false;
1334 : }
1335 :
1336 :
1337 : #endif
1338 :
1339 : #ifndef NDEBUG
1340 :
1341 0 : void Elem::libmesh_assert_valid_node_pointers() const
1342 : {
1343 0 : libmesh_assert(this->valid_id());
1344 0 : for (auto n : this->node_index_range())
1345 : {
1346 0 : libmesh_assert(this->node_ptr(n));
1347 0 : libmesh_assert(this->node_ptr(n)->valid_id());
1348 : }
1349 0 : }
1350 :
1351 :
1352 :
1353 2662706 : void Elem::libmesh_assert_valid_neighbors() const
1354 : {
1355 13163086 : for (auto n : this->side_index_range())
1356 : {
1357 10500380 : const Elem * neigh = this->neighbor_ptr(n);
1358 :
1359 : // Any element might have a remote neighbor; checking
1360 : // to make sure that's not inaccurate is tough.
1361 10500380 : if (neigh == remote_elem)
1362 8524 : continue;
1363 :
1364 10491856 : if (neigh)
1365 : {
1366 : // Only subactive elements have subactive neighbors
1367 9928312 : libmesh_assert (this->subactive() || !neigh->subactive());
1368 :
1369 9928312 : const Elem * elem = this;
1370 :
1371 : // If we're subactive but our neighbor isn't, its
1372 : // return neighbor link will be to our first active
1373 : // ancestor OR to our inactive ancestor of the same
1374 : // level as neigh,
1375 9928312 : if (this->subactive() && !neigh->subactive())
1376 : {
1377 40808 : for (elem = this; !elem->active();
1378 20404 : elem = elem->parent())
1379 20404 : libmesh_assert(elem);
1380 : }
1381 : else
1382 : {
1383 9907908 : unsigned int rev = neigh->which_neighbor_am_i(elem);
1384 9907908 : libmesh_assert_less (rev, neigh->n_neighbors());
1385 :
1386 9907908 : if (this->subactive() && !neigh->subactive())
1387 : {
1388 0 : while (neigh->neighbor_ptr(rev) != elem)
1389 : {
1390 0 : libmesh_assert(elem->parent());
1391 0 : elem = elem->parent();
1392 : }
1393 : }
1394 : else
1395 : {
1396 9907908 : const Elem * nn = neigh->neighbor_ptr(rev);
1397 9907908 : libmesh_assert(nn);
1398 :
1399 10177160 : for (; elem != nn; elem = elem->parent())
1400 269252 : libmesh_assert(elem);
1401 : }
1402 : }
1403 : }
1404 : // If we don't have a neighbor and we're not subactive, our
1405 : // ancestors shouldn't have any neighbors in this same
1406 : // direction.
1407 563544 : else if (!this->subactive())
1408 : {
1409 552828 : const Elem * my_parent = this->parent();
1410 804214 : if (my_parent &&
1411 : // A parent with a different dimension isn't really one of
1412 : // our ancestors, it means we're on a boundary mesh and this
1413 : // is an interior mesh element for which we're on a side.
1414 : // Nothing to test for in that case.
1415 251386 : (my_parent->dim() == this->dim()))
1416 251386 : libmesh_assert (!my_parent->neighbor_ptr(n));
1417 : }
1418 : }
1419 2662706 : }
1420 :
1421 : #endif // !NDEBUG
1422 :
1423 :
1424 :
1425 54001105 : void Elem::make_links_to_me_local(unsigned int n, unsigned int nn)
1426 : {
1427 54001105 : Elem * neigh = this->neighbor_ptr(n);
1428 :
1429 : // Don't bother calling this function unless it's necessary
1430 787 : libmesh_assert(neigh);
1431 787 : libmesh_assert(!neigh->is_remote());
1432 :
1433 : // We never have neighbors more refined than us
1434 787 : libmesh_assert_less_equal (neigh->level(), this->level());
1435 :
1436 : // We never have subactive neighbors of non subactive elements
1437 787 : libmesh_assert(!neigh->subactive() || this->subactive());
1438 :
1439 : // If we have a neighbor less refined than us then it must not
1440 : // have any more refined descendants we could have pointed to
1441 : // instead.
1442 787 : libmesh_assert((neigh->level() == this->level()) ||
1443 : (neigh->active() && !this->subactive()) ||
1444 : (!neigh->has_children() && this->subactive()));
1445 :
1446 : // If neigh is at our level, then its family might have
1447 : // remote_elem neighbor links which need to point to us
1448 : // instead, but if not, then we're done.
1449 54001105 : if (neigh->level() != this->level())
1450 686296 : return;
1451 :
1452 : // What side of neigh are we on? nn.
1453 : //
1454 : // We can't use the usual Elem method because we're in the middle of
1455 : // restoring topology. We can't compare side_ptr nodes because
1456 : // users want to abuse neighbor_ptr to point to
1457 : // not-technically-neighbors across mesh slits. We can't compare
1458 : // node locations because users want to move those
1459 : // not-technically-neighbors until they're
1460 : // not-even-geometrically-neighbors.
1461 :
1462 : // Find any elements that ought to point to elem
1463 1574 : std::vector<Elem *> neigh_family;
1464 : #ifdef LIBMESH_ENABLE_AMR
1465 787 : if (this->active())
1466 41698445 : neigh->family_tree_by_side(neigh_family, nn);
1467 : else
1468 : #endif
1469 11616364 : neigh_family.push_back(neigh);
1470 :
1471 : // And point them to elem
1472 106642017 : for (auto & neigh_family_member : neigh_family)
1473 : {
1474 : // Only subactive elements point to other subactive elements
1475 53327208 : if (this->subactive() && !neigh_family_member->subactive())
1476 1152 : continue;
1477 :
1478 : // Ideally, the neighbor link ought to either be correct
1479 : // already or ought to be to remote_elem.
1480 : //
1481 : // However, if we're redistributing a newly created elem,
1482 : // after an AMR step but before find_neighbors has fixed up
1483 : // neighbor links, we might have an out of date neighbor
1484 : // link to elem's parent instead.
1485 : #ifdef LIBMESH_ENABLE_AMR
1486 787 : libmesh_assert((neigh_family_member->neighbor_ptr(nn) &&
1487 : (neigh_family_member->neighbor_ptr(nn)->active() ||
1488 : neigh_family_member->neighbor_ptr(nn)->is_ancestor_of(this))) ||
1489 : (neigh_family_member->neighbor_ptr(nn) == remote_elem) ||
1490 : ((this->refinement_flag() == JUST_REFINED) &&
1491 : (this->parent() != nullptr) &&
1492 : (neigh_family_member->neighbor_ptr(nn) == this->parent())));
1493 : #else
1494 : libmesh_assert((neigh_family_member->neighbor_ptr(nn) == this) ||
1495 : (neigh_family_member->neighbor_ptr(nn) == remote_elem));
1496 : #endif
1497 :
1498 53326056 : neigh_family_member->set_neighbor(nn, this);
1499 : }
1500 : }
1501 :
1502 :
1503 35019046 : void Elem::make_links_to_me_remote()
1504 : {
1505 1633 : libmesh_assert_not_equal_to (this, remote_elem);
1506 :
1507 : // We need to have handled any children first
1508 : #if defined(LIBMESH_ENABLE_AMR) && defined(DEBUG)
1509 1633 : if (this->has_children())
1510 268 : for (auto & child : this->child_ref_range())
1511 208 : libmesh_assert_equal_to (&child, remote_elem);
1512 : #endif
1513 :
1514 : // Remotify any neighbor links
1515 177333274 : for (auto neigh : this->neighbor_ptr_range())
1516 : {
1517 142314228 : if (neigh && neigh != remote_elem)
1518 : {
1519 : // My neighbor should never be more refined than me; my real
1520 : // neighbor would have been its parent in that case.
1521 3572 : libmesh_assert_greater_equal (this->level(), neigh->level());
1522 :
1523 137671361 : if (this->level() == neigh->level() &&
1524 68494755 : neigh->has_neighbor(this))
1525 : {
1526 : #ifdef LIBMESH_ENABLE_AMR
1527 : // My neighbor may have descendants which also consider me a
1528 : // neighbor
1529 7144 : std::vector<Elem *> family;
1530 68497700 : neigh->total_family_tree_by_neighbor (family, this);
1531 :
1532 : // FIXME - There's a lot of ugly const_casts here; we
1533 : // may want to make remote_elem non-const
1534 137003188 : for (auto & n : family)
1535 : {
1536 3572 : libmesh_assert (n);
1537 68505488 : if (n == remote_elem)
1538 0 : continue;
1539 68505488 : unsigned int my_s = n->which_neighbor_am_i(this);
1540 3572 : libmesh_assert_less (my_s, n->n_neighbors());
1541 3572 : libmesh_assert_equal_to (n->neighbor_ptr(my_s), this);
1542 68505488 : n->set_neighbor(my_s, const_cast<RemoteElem *>(remote_elem));
1543 : }
1544 : #else
1545 : unsigned int my_s = neigh->which_neighbor_am_i(this);
1546 : libmesh_assert_less (my_s, neigh->n_neighbors());
1547 : libmesh_assert_equal_to (neigh->neighbor_ptr(my_s), this);
1548 : neigh->set_neighbor(my_s, const_cast<RemoteElem *>(remote_elem));
1549 : #endif
1550 : }
1551 : #ifdef LIBMESH_ENABLE_AMR
1552 : // Even if my neighbor doesn't link back to me, it might
1553 : // have subactive descendants which do
1554 0 : else if (neigh->has_children())
1555 : {
1556 : // If my neighbor at the same level doesn't have me as a
1557 : // neighbor, I must be subactive
1558 0 : libmesh_assert(this->level() > neigh->level() ||
1559 : this->subactive());
1560 :
1561 : // My neighbor must have some ancestor of mine as a
1562 : // neighbor
1563 0 : Elem * my_ancestor = this->parent();
1564 0 : libmesh_assert(my_ancestor);
1565 1760 : while (!neigh->has_neighbor(my_ancestor))
1566 : {
1567 0 : my_ancestor = my_ancestor->parent();
1568 0 : libmesh_assert(my_ancestor);
1569 : }
1570 :
1571 : // My neighbor may have descendants which consider me a
1572 : // neighbor
1573 0 : std::vector<Elem *> family;
1574 1760 : neigh->total_family_tree_by_subneighbor (family, my_ancestor, this);
1575 :
1576 2874 : for (auto & n : family)
1577 : {
1578 0 : libmesh_assert (n);
1579 1114 : if (n->is_remote())
1580 0 : continue;
1581 1114 : unsigned int my_s = n->which_neighbor_am_i(this);
1582 0 : libmesh_assert_less (my_s, n->n_neighbors());
1583 0 : libmesh_assert_equal_to (n->neighbor_ptr(my_s), this);
1584 : // TODO: we may want to make remote_elem non-const.
1585 1114 : n->set_neighbor(my_s, const_cast<RemoteElem *>(remote_elem));
1586 : }
1587 : }
1588 : #endif
1589 : }
1590 : }
1591 :
1592 : #ifdef LIBMESH_ENABLE_AMR
1593 : // Remotify parent's child link
1594 3266 : Elem * my_parent = this->parent();
1595 26948262 : if (my_parent &&
1596 : // As long as it's not already remote
1597 61965675 : my_parent != remote_elem &&
1598 : // And it's a real parent, not an interior parent
1599 26946629 : this->dim() == my_parent->dim())
1600 : {
1601 26946629 : unsigned int me = my_parent->which_child_am_i(this);
1602 672 : libmesh_assert_equal_to (my_parent->child_ptr(me), this);
1603 26946629 : my_parent->set_child(me, const_cast<RemoteElem *>(remote_elem));
1604 : }
1605 : #endif
1606 35019046 : }
1607 :
1608 :
1609 813 : void Elem::remove_links_to_me()
1610 : {
1611 31 : libmesh_assert_not_equal_to (this, remote_elem);
1612 :
1613 : // We need to have handled any children first
1614 : #ifdef LIBMESH_ENABLE_AMR
1615 31 : libmesh_assert (!this->has_children());
1616 : #endif
1617 :
1618 : // Nullify any neighbor links
1619 4065 : for (auto neigh : this->neighbor_ptr_range())
1620 : {
1621 3252 : if (neigh && neigh != remote_elem)
1622 : {
1623 : // My neighbor should never be more refined than me; my real
1624 : // neighbor would have been its parent in that case.
1625 47 : libmesh_assert_greater_equal (this->level(), neigh->level());
1626 :
1627 2504 : if (this->level() == neigh->level() &&
1628 1205 : neigh->has_neighbor(this))
1629 : {
1630 : #ifdef LIBMESH_ENABLE_AMR
1631 : // My neighbor may have descendants which also consider me a
1632 : // neighbor
1633 94 : std::vector<Elem *> family;
1634 1252 : neigh->total_family_tree_by_neighbor (family, this);
1635 :
1636 2504 : for (auto & n : family)
1637 : {
1638 47 : libmesh_assert (n);
1639 1252 : if (n->is_remote())
1640 0 : continue;
1641 1252 : unsigned int my_s = n->which_neighbor_am_i(this);
1642 47 : libmesh_assert_less (my_s, n->n_neighbors());
1643 47 : libmesh_assert_equal_to (n->neighbor_ptr(my_s), this);
1644 1252 : n->set_neighbor(my_s, nullptr);
1645 : }
1646 : #else
1647 : unsigned int my_s = neigh->which_neighbor_am_i(this);
1648 : libmesh_assert_less (my_s, neigh->n_neighbors());
1649 : libmesh_assert_equal_to (neigh->neighbor_ptr(my_s), this);
1650 : neigh->set_neighbor(my_s, nullptr);
1651 : #endif
1652 : }
1653 : #ifdef LIBMESH_ENABLE_AMR
1654 : // Even if my neighbor doesn't link back to me, it might
1655 : // have subactive descendants which do
1656 0 : else if (neigh->has_children())
1657 : {
1658 : // If my neighbor at the same level doesn't have me as a
1659 : // neighbor, I must be subactive
1660 0 : libmesh_assert(this->level() > neigh->level() ||
1661 : this->subactive());
1662 :
1663 : // My neighbor must have some ancestor of mine as a
1664 : // neighbor
1665 0 : Elem * my_ancestor = this->parent();
1666 0 : libmesh_assert(my_ancestor);
1667 0 : while (!neigh->has_neighbor(my_ancestor))
1668 : {
1669 0 : my_ancestor = my_ancestor->parent();
1670 0 : libmesh_assert(my_ancestor);
1671 : }
1672 :
1673 : // My neighbor may have descendants which consider me a
1674 : // neighbor
1675 0 : std::vector<Elem *> family;
1676 0 : neigh->total_family_tree_by_subneighbor (family, my_ancestor, this);
1677 :
1678 0 : for (auto & n : family)
1679 : {
1680 0 : libmesh_assert (n);
1681 0 : if (n->is_remote())
1682 0 : continue;
1683 0 : unsigned int my_s = n->which_neighbor_am_i(this);
1684 0 : libmesh_assert_less (my_s, n->n_neighbors());
1685 0 : libmesh_assert_equal_to (n->neighbor_ptr(my_s), this);
1686 0 : n->set_neighbor(my_s, nullptr);
1687 : }
1688 : }
1689 : #endif
1690 : }
1691 : }
1692 :
1693 : #ifdef LIBMESH_ENABLE_AMR
1694 : // We can't currently delete a child with a parent!
1695 31 : libmesh_assert (!this->parent());
1696 : #endif
1697 813 : }
1698 :
1699 :
1700 :
1701 178350 : void Elem::write_connectivity (std::ostream & out_stream,
1702 : const IOPackage iop) const
1703 : {
1704 16310 : libmesh_assert (out_stream.good());
1705 16310 : libmesh_assert(_nodes);
1706 16310 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
1707 :
1708 178350 : switch (iop)
1709 : {
1710 16185 : case TECPLOT:
1711 : {
1712 : // This connectivity vector will be used repeatedly instead
1713 : // of being reconstructed inside the loop.
1714 32370 : std::vector<dof_id_type> conn;
1715 498180 : for (auto sc : make_range(this->n_sub_elem()))
1716 : {
1717 176850 : this->connectivity(sc, TECPLOT, conn);
1718 :
1719 176850 : std::copy(conn.begin(),
1720 : conn.end(),
1721 16185 : std::ostream_iterator<dof_id_type>(out_stream, " "));
1722 :
1723 176850 : out_stream << '\n';
1724 : }
1725 16185 : return;
1726 : }
1727 :
1728 125 : case UCD:
1729 : {
1730 13500 : for (auto i : this->node_index_range())
1731 13000 : out_stream << this->node_id(i)+1 << "\t";
1732 :
1733 1500 : out_stream << '\n';
1734 1500 : return;
1735 : }
1736 :
1737 0 : default:
1738 0 : libmesh_error_msg("Unsupported IO package " << iop);
1739 : }
1740 : }
1741 :
1742 :
1743 :
1744 79011 : Real Elem::quality (const ElemQuality q) const
1745 : {
1746 79011 : switch (q)
1747 : {
1748 : // Return the maximum ratio of edge lengths, or zero if that
1749 : // maximum would otherwise be infinity.
1750 11265 : case EDGE_LENGTH_RATIO:
1751 : {
1752 11265 : if (this->dim() < 2)
1753 7 : return 1;
1754 :
1755 12118 : std::vector<Real> edge_lengths(this->n_edges());
1756 85815 : for (auto e : index_range(edge_lengths))
1757 143017 : edge_lengths[e] = (this->point(this->local_edge_node(e,1)) -
1758 143017 : this->point(this->local_edge_node(e,0))).norm();
1759 :
1760 937 : auto [min, max] =
1761 11181 : std::minmax_element(edge_lengths.begin(), edge_lengths.end());
1762 :
1763 11181 : if (*min == 0.)
1764 0 : return 0.;
1765 : else
1766 11181 : return *max / *min;
1767 : }
1768 :
1769 23524 : case MIN_ANGLE:
1770 : case MAX_ANGLE:
1771 : {
1772 : // 1D elements don't have interior angles, so just return some
1773 : // dummy value in that case.
1774 23524 : if (this->dim() < 2)
1775 14 : return 0.;
1776 :
1777 : // Initialize return values
1778 23356 : Real min_angle = std::numeric_limits<Real>::max();
1779 23356 : Real max_angle = -std::numeric_limits<Real>::max();
1780 :
1781 263648 : for (auto n : this->node_index_range())
1782 : {
1783 : // Get list of edge ids adjacent to this node.
1784 240292 : const auto adjacent_edge_ids = this->edges_adjacent_to_node(n);
1785 :
1786 : // Skip any nodes with fewer than 2 adjacent edges. You
1787 : // need at least two adjacent edges to form an interior
1788 : // element angle.
1789 39820 : auto N = adjacent_edge_ids.size();
1790 240292 : if (N < 2)
1791 11472 : continue;
1792 :
1793 : // Consider all possible pairs of edges adjacent to node n
1794 306924 : for (unsigned int first = 0; first < N-1; ++first)
1795 511412 : for (unsigned int second = first+1; second < N; ++second)
1796 : {
1797 : // Get ids of first and second edges
1798 333450 : auto first_edge = adjacent_edge_ids[first];
1799 307980 : auto second_edge = adjacent_edge_ids[second];
1800 :
1801 : // Get node ids of first and second edge
1802 307980 : auto first_edge_node_0 = this->local_edge_node(first_edge, 0);
1803 307980 : auto first_edge_node_1 = this->local_edge_node(first_edge, 1);
1804 307980 : auto second_edge_node_0 = this->local_edge_node(second_edge, 0);
1805 307980 : auto second_edge_node_1 = this->local_edge_node(second_edge, 1);
1806 :
1807 : // Orient both edges so that edge node 0 == n
1808 307980 : if (first_edge_node_0 != n)
1809 17624 : std::swap(first_edge_node_0, first_edge_node_1);
1810 307980 : if (second_edge_node_0 != n)
1811 8574 : std::swap(second_edge_node_0, second_edge_node_1);
1812 :
1813 25470 : libmesh_assert_equal_to(first_edge_node_0, n);
1814 25470 : libmesh_assert_equal_to(second_edge_node_0, n);
1815 :
1816 : // Locally oriented edge vectors
1817 : Point
1818 50940 : first_ev = this->point(first_edge_node_1) - this->point(first_edge_node_0),
1819 25470 : second_ev = this->point(second_edge_node_1) - this->point(second_edge_node_0);
1820 :
1821 : // Angle between them in the range [0, pi]
1822 307980 : Real theta = std::acos(first_ev.unit() * second_ev.unit());
1823 :
1824 : // Track min and max angles seen
1825 307980 : min_angle = std::min(theta, min_angle);
1826 307980 : max_angle = std::max(theta, max_angle);
1827 : }
1828 : }
1829 :
1830 : // Return requested extreme value (in degrees)
1831 23356 : return Real(180)/libMesh::pi * ((q == MIN_ANGLE) ? min_angle : max_angle);
1832 : }
1833 :
1834 21124 : case MIN_DIHEDRAL_ANGLE:
1835 : case MAX_DIHEDRAL_ANGLE:
1836 : {
1837 : // For 1D and 2D elements, just call this function with MIN,MAX_ANGLE instead.
1838 21124 : if (this->dim() < 3)
1839 0 : return this->quality((q == MIN_DIHEDRAL_ANGLE) ? MIN_ANGLE : MAX_ANGLE);
1840 :
1841 : // Initialize return values
1842 21124 : Real min_angle = std::numeric_limits<Real>::max();
1843 21124 : Real max_angle = -std::numeric_limits<Real>::max();
1844 :
1845 : // Algorithm for computing dihedral angles:
1846 : // .) Loop over the edges, use the edge_sides_map to get the
1847 : // two sides adjacent to the edge.
1848 : // .) For each adjacent side, use the side_nodes_map to get the
1849 : // list of three (or more) node ids on that side.
1850 : // .) Use the node ids to compute the (approximate) inward
1851 : // normal for the side. Note: this is approximate since 3D
1852 : // elements with four-sided faces and quadratic tetrahedra
1853 : // do not necessarily have planar faces.
1854 : // .) Compute the dihedral angle for the current edge, compare
1855 : // it to the current min and max dihedral angles.
1856 166480 : for (auto e : this->edge_index_range())
1857 : {
1858 : // Get list of edge ids adjacent to this node.
1859 157476 : const auto adjacent_side_ids = this->sides_on_edge(e);
1860 :
1861 : // All 3D elements should have exactly two sides adjacent to each edge.
1862 12120 : libmesh_assert_equal_to(adjacent_side_ids.size(), 2);
1863 :
1864 : // Get lists of node ids on each side
1865 157476 : const auto side_0_node_ids = this->nodes_on_side(adjacent_side_ids[0]);
1866 145356 : const auto side_1_node_ids = this->nodes_on_side(adjacent_side_ids[1]);
1867 :
1868 : // All 3D elements have at least three nodes on each side
1869 12120 : libmesh_assert_greater_equal(side_0_node_ids.size(), 3);
1870 12120 : libmesh_assert_greater_equal(side_1_node_ids.size(), 3);
1871 :
1872 : // Construct (approximate) inward normal on each adjacent side
1873 : const auto side_0_normal =
1874 157476 : (this->point(side_0_node_ids[2]) - this->point(side_0_node_ids[0])).cross
1875 157476 : (this->point(side_0_node_ids[1]) - this->point(side_0_node_ids[0])).unit();
1876 : const auto side_1_normal =
1877 157476 : (this->point(side_1_node_ids[2]) - this->point(side_1_node_ids[0])).cross
1878 157476 : (this->point(side_1_node_ids[1]) - this->point(side_1_node_ids[0])).unit();
1879 :
1880 : // Compute dihedral angle between the planes.
1881 : // Using the absolute value ensures that we get the same result
1882 : // even if the orientation of one of the planes is flipped, i.e.
1883 : // it always gives us a value in the range [0, pi/2].
1884 : // https://en.wikipedia.org/wiki/Dihedral_angle
1885 145356 : Real theta = std::acos(std::abs(side_0_normal * side_1_normal));
1886 :
1887 : // Track min and max angles seen
1888 145356 : min_angle = std::min(theta, min_angle);
1889 145356 : max_angle = std::max(theta, max_angle);
1890 : }
1891 :
1892 : // Return requested extreme value (in degrees)
1893 21124 : return Real(180)/libMesh::pi * ((q == MIN_DIHEDRAL_ANGLE) ? min_angle : max_angle);
1894 : }
1895 :
1896 23098 : case JACOBIAN:
1897 : case SCALED_JACOBIAN:
1898 : {
1899 : // 1D elements don't have interior corners, so this metric
1900 : // does not really apply to them.
1901 23098 : const auto N = this->dim();
1902 23098 : if (N < 2)
1903 14 : return 1.;
1904 :
1905 : // Initialize return value
1906 22930 : Real min_node_area = std::numeric_limits<Real>::max();
1907 :
1908 261518 : for (auto n : this->node_index_range())
1909 : {
1910 : // Get list of edge ids adjacent to this node.
1911 238588 : auto adjacent_edge_ids = this->edges_adjacent_to_node(n);
1912 :
1913 : // Skip any nodes that don't have dim() adjacent edges. In 2D,
1914 : // you need at least two adjacent edges to compute the cross
1915 : // product, and in 3D, you need at least three adjacent edges
1916 : // to compute the scalar triple product. The only element
1917 : // type which has an unusual topology in this regard is the
1918 : // Pyramid element in 3D, where 4 edges meet at the apex node.
1919 : // For now, we just skip this node when computing the JACOBIAN
1920 : // metric for Pyramids.
1921 258450 : if (adjacent_edge_ids.size() != N)
1922 11856 : continue;
1923 :
1924 : // Construct oriented edges
1925 97180 : std::vector<Point> oriented_edges(N);
1926 382832 : for (auto i : make_range(N))
1927 : {
1928 309236 : auto node_0 = this->local_edge_node(adjacent_edge_ids[i], 0);
1929 309236 : auto node_1 = this->local_edge_node(adjacent_edge_ids[i], 1);
1930 285652 : if (node_0 != n)
1931 11002 : std::swap(node_0, node_1);
1932 332820 : oriented_edges[i] = this->point(node_1) - this->point(node_0);
1933 : }
1934 :
1935 : // Compute unscaled area (2D) or volume (3D) using the
1936 : // cross product (2D) or scalar triple product (3D) of the
1937 : // oriented edges. We take the absolute value so that we
1938 : // don't have to worry about the sign of the area.
1939 97180 : Real node_area = (N == 2) ?
1940 5888 : cross_norm(oriented_edges[0], oriented_edges[1]) :
1941 91292 : std::abs(triple_product(oriented_edges[0], oriented_edges[1], oriented_edges[2]));
1942 :
1943 : // Divide by (non-zero) edge lengths if computing scaled Jacobian.
1944 : // If any length is zero, then set node_area to zero. This means that
1945 : // e.g. degenerate quadrilaterals will have a zero SCALED_JACOBIAN metric.
1946 97180 : if (q == SCALED_JACOBIAN)
1947 191416 : for (auto i : make_range(N))
1948 : {
1949 142826 : Real len_i = oriented_edges[i].norm();
1950 142826 : node_area = (len_i == 0.) ? 0. : (node_area / len_i);
1951 : }
1952 :
1953 : // Update minimum
1954 97180 : min_node_area = std::min(node_area, min_node_area);
1955 : }
1956 :
1957 22930 : return min_node_area;
1958 : }
1959 :
1960 : // Return 1 if we made it here
1961 0 : default:
1962 : {
1963 0 : libmesh_do_once( libmesh_here();
1964 :
1965 : libMesh::err << "ERROR: quality metric "
1966 : << Utility::enum_to_string(q)
1967 : << " not implemented on element type "
1968 : << Utility::enum_to_string(this->type())
1969 : << std::endl
1970 : << "Returning 1."
1971 : << std::endl; );
1972 :
1973 0 : return 1.;
1974 : }
1975 : }
1976 : }
1977 :
1978 :
1979 :
1980 62018359 : bool Elem::ancestor() const
1981 : {
1982 : #ifdef LIBMESH_ENABLE_AMR
1983 :
1984 : // Use a fast, DistributedMesh-safe definition
1985 : const bool is_ancestor =
1986 20062853 : !this->active() && !this->subactive();
1987 :
1988 : // But check for inconsistencies if we have time
1989 : #ifdef DEBUG
1990 4868824 : if (!is_ancestor && this->has_children())
1991 : {
1992 467716 : for (auto & c : this->child_ref_range())
1993 : {
1994 374864 : if (&c != remote_elem)
1995 : {
1996 374816 : libmesh_assert(!c.active());
1997 374816 : libmesh_assert(!c.ancestor());
1998 : }
1999 : }
2000 : }
2001 : #endif // DEBUG
2002 :
2003 62018359 : return is_ancestor;
2004 :
2005 : #else
2006 : return false;
2007 : #endif
2008 : }
2009 :
2010 :
2011 :
2012 : #ifdef LIBMESH_ENABLE_AMR
2013 :
2014 105654 : void Elem::add_child (Elem * elem)
2015 : {
2016 105654 : const unsigned int nc = this->n_children();
2017 :
2018 105654 : if (!_children)
2019 : {
2020 27749 : _children = std::make_unique<Elem *[]>(nc);
2021 :
2022 132245 : for (unsigned int c = 0; c != nc; c++)
2023 4612 : this->set_child(c, nullptr);
2024 : }
2025 :
2026 263425 : for (unsigned int c = 0; c != nc; c++)
2027 : {
2028 274935 : if (this->_children[c] == nullptr || this->_children[c] == remote_elem)
2029 : {
2030 4612 : libmesh_assert_equal_to (this, elem->parent());
2031 4612 : this->set_child(c, elem);
2032 105654 : return;
2033 : }
2034 : }
2035 :
2036 0 : libmesh_error_msg("Error: Tried to add a child to an element with full children array");
2037 : }
2038 :
2039 :
2040 :
2041 50179838 : void Elem::add_child (Elem * elem, unsigned int c)
2042 : {
2043 49820 : if (!this->has_children())
2044 : {
2045 6182503 : const unsigned int nc = this->n_children();
2046 6196225 : _children = std::make_unique<Elem *[]>(nc);
2047 :
2048 31112617 : for (unsigned int i = 0; i != nc; i++)
2049 49468 : this->set_child(i, nullptr);
2050 : }
2051 :
2052 49820 : libmesh_assert (this->_children[c] == nullptr || this->child_ptr(c) == remote_elem);
2053 49820 : libmesh_assert (elem == remote_elem || this == elem->parent());
2054 :
2055 49820 : this->set_child(c, elem);
2056 50179838 : }
2057 :
2058 :
2059 :
2060 231463 : void Elem::replace_child (Elem * elem, unsigned int c)
2061 : {
2062 24248 : libmesh_assert(this->has_children());
2063 :
2064 24248 : libmesh_assert(this->child_ptr(c));
2065 :
2066 24248 : this->set_child(c, elem);
2067 231463 : }
2068 :
2069 :
2070 :
2071 9501 : void Elem::family_tree (std::vector<const Elem *> & family,
2072 : bool reset) const
2073 : {
2074 9501 : ElemInternal::family_tree(this, family, reset);
2075 9501 : }
2076 :
2077 :
2078 :
2079 0 : void Elem::family_tree (std::vector<Elem *> & family,
2080 : bool reset)
2081 : {
2082 0 : ElemInternal::family_tree(this, family, reset);
2083 0 : }
2084 :
2085 :
2086 :
2087 46312 : void Elem::total_family_tree (std::vector<const Elem *> & family,
2088 : bool reset) const
2089 : {
2090 46312 : ElemInternal::total_family_tree(this, family, reset);
2091 46312 : }
2092 :
2093 :
2094 :
2095 28131573 : void Elem::total_family_tree (std::vector<Elem *> & family,
2096 : bool reset)
2097 : {
2098 28131573 : ElemInternal::total_family_tree(this, family, reset);
2099 28131573 : }
2100 :
2101 :
2102 :
2103 927070 : void Elem::active_family_tree (std::vector<const Elem *> & active_family,
2104 : bool reset) const
2105 : {
2106 927070 : ElemInternal::active_family_tree(this, active_family, reset);
2107 927070 : }
2108 :
2109 :
2110 :
2111 0 : void Elem::active_family_tree (std::vector<Elem *> & active_family,
2112 : bool reset)
2113 : {
2114 0 : ElemInternal::active_family_tree(this, active_family, reset);
2115 0 : }
2116 :
2117 :
2118 :
2119 0 : void Elem::family_tree_by_side (std::vector<const Elem *> & family,
2120 : unsigned int side,
2121 : bool reset) const
2122 : {
2123 0 : ElemInternal::family_tree_by_side(this, family, side, reset);
2124 0 : }
2125 :
2126 :
2127 :
2128 41698445 : void Elem:: family_tree_by_side (std::vector<Elem *> & family,
2129 : unsigned int side,
2130 : bool reset)
2131 : {
2132 41698445 : ElemInternal::family_tree_by_side(this, family, side, reset);
2133 41698445 : }
2134 :
2135 :
2136 :
2137 596289 : void Elem::active_family_tree_by_side (std::vector<const Elem *> & family,
2138 : unsigned int side,
2139 : bool reset) const
2140 : {
2141 596289 : ElemInternal::active_family_tree_by_side(this, family, side, reset);
2142 596289 : }
2143 :
2144 :
2145 :
2146 0 : void Elem::active_family_tree_by_side (std::vector<Elem *> & family,
2147 : unsigned int side,
2148 : bool reset)
2149 : {
2150 0 : ElemInternal::active_family_tree_by_side(this, family, side, reset);
2151 0 : }
2152 :
2153 :
2154 :
2155 0 : void Elem::family_tree_by_neighbor (std::vector<const Elem *> & family,
2156 : const Elem * neighbor,
2157 : bool reset) const
2158 : {
2159 0 : ElemInternal::family_tree_by_neighbor(this, family, neighbor, reset);
2160 0 : }
2161 :
2162 :
2163 :
2164 0 : void Elem::family_tree_by_neighbor (std::vector<Elem *> & family,
2165 : Elem * neighbor,
2166 : bool reset)
2167 : {
2168 0 : ElemInternal::family_tree_by_neighbor(this, family, neighbor, reset);
2169 0 : }
2170 :
2171 :
2172 :
2173 0 : void Elem::total_family_tree_by_neighbor (std::vector<const Elem *> & family,
2174 : const Elem * neighbor,
2175 : bool reset) const
2176 : {
2177 0 : ElemInternal::total_family_tree_by_neighbor(this, family, neighbor, reset);
2178 0 : }
2179 :
2180 :
2181 :
2182 68498952 : void Elem::total_family_tree_by_neighbor (std::vector<Elem *> & family,
2183 : Elem * neighbor,
2184 : bool reset)
2185 : {
2186 68498952 : ElemInternal::total_family_tree_by_neighbor(this, family, neighbor, reset);
2187 68498952 : }
2188 :
2189 :
2190 :
2191 0 : void Elem::family_tree_by_subneighbor (std::vector<const Elem *> & family,
2192 : const Elem * neighbor,
2193 : const Elem * subneighbor,
2194 : bool reset) const
2195 : {
2196 0 : ElemInternal::family_tree_by_subneighbor(this, family, neighbor, subneighbor, reset);
2197 0 : }
2198 :
2199 :
2200 :
2201 0 : void Elem::family_tree_by_subneighbor (std::vector<Elem *> & family,
2202 : Elem * neighbor,
2203 : Elem * subneighbor,
2204 : bool reset)
2205 : {
2206 0 : ElemInternal::family_tree_by_subneighbor(this, family, neighbor, subneighbor, reset);
2207 0 : }
2208 :
2209 :
2210 :
2211 0 : void Elem::total_family_tree_by_subneighbor (std::vector<const Elem *> & family,
2212 : const Elem * neighbor,
2213 : const Elem * subneighbor,
2214 : bool reset) const
2215 : {
2216 0 : ElemInternal::total_family_tree_by_subneighbor(this, family, neighbor, subneighbor, reset);
2217 0 : }
2218 :
2219 :
2220 :
2221 2874 : void Elem::total_family_tree_by_subneighbor (std::vector<Elem *> & family,
2222 : Elem * neighbor,
2223 : Elem * subneighbor,
2224 : bool reset)
2225 : {
2226 2874 : ElemInternal::total_family_tree_by_subneighbor(this, family, neighbor, subneighbor, reset);
2227 2874 : }
2228 :
2229 :
2230 :
2231 39175978 : void Elem::active_family_tree_by_neighbor (std::vector<const Elem *> & family,
2232 : const Elem * neighbor,
2233 : bool reset) const
2234 : {
2235 39175978 : ElemInternal::active_family_tree_by_neighbor(this, family, neighbor, reset);
2236 39175978 : }
2237 :
2238 :
2239 :
2240 2736 : void Elem::active_family_tree_by_neighbor (std::vector<Elem *> & family,
2241 : Elem * neighbor,
2242 : bool reset)
2243 : {
2244 2736 : ElemInternal::active_family_tree_by_neighbor(this, family, neighbor, reset);
2245 2736 : }
2246 :
2247 :
2248 :
2249 26334 : void Elem::active_family_tree_by_topological_neighbor (std::vector<const Elem *> & family,
2250 : const Elem * neighbor,
2251 : const MeshBase & mesh,
2252 : const PointLocatorBase & point_locator,
2253 : const PeriodicBoundaries * pb,
2254 : bool reset) const
2255 : {
2256 26334 : ElemInternal::active_family_tree_by_topological_neighbor(this, family, neighbor,
2257 : mesh, point_locator, pb,
2258 : reset);
2259 26334 : }
2260 :
2261 :
2262 :
2263 0 : void Elem::active_family_tree_by_topological_neighbor (std::vector<Elem *> & family,
2264 : Elem * neighbor,
2265 : const MeshBase & mesh,
2266 : const PointLocatorBase & point_locator,
2267 : const PeriodicBoundaries * pb,
2268 : bool reset)
2269 : {
2270 0 : ElemInternal::active_family_tree_by_topological_neighbor(this, family, neighbor,
2271 : mesh, point_locator, pb,
2272 : reset);
2273 0 : }
2274 :
2275 :
2276 26856167 : bool Elem::is_child_on_edge(const unsigned int c,
2277 : const unsigned int e) const
2278 : {
2279 16129880 : libmesh_assert_less (c, this->n_children());
2280 16129880 : libmesh_assert_less (e, this->n_edges());
2281 :
2282 41938193 : std::unique_ptr<const Elem> my_edge = this->build_edge_ptr(e);
2283 25808313 : std::unique_ptr<const Elem> child_edge = this->child_ptr(c)->build_edge_ptr(e);
2284 :
2285 : // We're assuming that an overlapping child edge has the same
2286 : // number and orientation as its parent
2287 41841863 : return (child_edge->node_id(0) == my_edge->node_id(0) ||
2288 55876035 : child_edge->node_id(1) == my_edge->node_id(1));
2289 9678433 : }
2290 :
2291 :
2292 :
2293 14724254 : unsigned int Elem::min_p_level_by_neighbor(const Elem * neighbor_in,
2294 : unsigned int current_min) const
2295 : {
2296 1370299 : libmesh_assert(!this->subactive());
2297 1370299 : libmesh_assert(neighbor_in->active());
2298 :
2299 : // If we're an active element this is simple
2300 1370299 : if (this->active())
2301 15927436 : return std::min(current_min, this->p_level());
2302 :
2303 17078 : libmesh_assert(has_neighbor(neighbor_in));
2304 :
2305 : // The p_level() of an ancestor element is already the minimum
2306 : // p_level() of its children - so if that's high enough, we don't
2307 : // need to examine any children.
2308 167583 : if (current_min <= this->p_level())
2309 17078 : return current_min;
2310 :
2311 0 : unsigned int min_p_level = current_min;
2312 :
2313 0 : for (auto & c : this->child_ref_range())
2314 0 : if (&c != remote_elem && c.has_neighbor(neighbor_in))
2315 : min_p_level =
2316 0 : c.min_p_level_by_neighbor(neighbor_in, min_p_level);
2317 :
2318 0 : return min_p_level;
2319 : }
2320 :
2321 :
2322 1356921 : unsigned int Elem::min_new_p_level_by_neighbor(const Elem * neighbor_in,
2323 : unsigned int current_min) const
2324 : {
2325 76336 : libmesh_assert(!this->subactive());
2326 76336 : libmesh_assert(neighbor_in->active());
2327 :
2328 : // If we're an active element this is simple
2329 76336 : if (this->active())
2330 : {
2331 903205 : unsigned int new_p_level = this->p_level();
2332 954513 : if (this->p_refinement_flag() == Elem::REFINE)
2333 16 : new_p_level += 1;
2334 903205 : if (this->p_refinement_flag() == Elem::COARSEN)
2335 : {
2336 4 : libmesh_assert_greater (new_p_level, 0);
2337 218 : new_p_level -= 1;
2338 : }
2339 903205 : return std::min(current_min, new_p_level);
2340 : }
2341 :
2342 24958 : libmesh_assert(has_neighbor(neighbor_in));
2343 :
2344 453716 : unsigned int min_p_level = current_min;
2345 :
2346 2285934 : for (auto & c : this->child_ref_range())
2347 1934834 : if (&c != remote_elem && c.has_neighbor(neighbor_in))
2348 : min_p_level =
2349 903205 : c.min_new_p_level_by_neighbor(neighbor_in, min_p_level);
2350 :
2351 453716 : return min_p_level;
2352 : }
2353 :
2354 :
2355 :
2356 190559429 : unsigned int Elem::as_parent_node (unsigned int child,
2357 : unsigned int child_node) const
2358 : {
2359 190559429 : const unsigned int nc = this->n_children();
2360 8438418 : libmesh_assert_less(child, nc);
2361 :
2362 : // Cached return values, indexed first by embedding_matrix version,
2363 : // then by child number, then by child node number.
2364 : std::vector<std::vector<std::vector<signed char>>> &
2365 190559429 : cached_parent_indices = this->_get_parent_indices_cache();
2366 :
2367 190559429 : unsigned int em_vers = this->embedding_matrix_version();
2368 :
2369 : // We may be updating the cache on one thread, and while that
2370 : // happens we can't safely access the cache from other threads.
2371 16876836 : Threads::spin_mutex::scoped_lock lock(parent_indices_mutex);
2372 :
2373 198997631 : if (em_vers >= cached_parent_indices.size())
2374 5618 : cached_parent_indices.resize(em_vers+1);
2375 :
2376 207435833 : if (child >= cached_parent_indices[em_vers].size())
2377 : {
2378 5622 : const signed char nn = cast_int<signed char>(this->n_nodes());
2379 :
2380 5830 : cached_parent_indices[em_vers].resize(nc);
2381 :
2382 33963 : for (unsigned int c = 0; c != nc; ++c)
2383 : {
2384 28341 : const unsigned int ncn = this->n_nodes_in_child(c);
2385 30385 : cached_parent_indices[em_vers][c].resize(ncn);
2386 290562 : for (unsigned int cn = 0; cn != ncn; ++cn)
2387 : {
2388 1256512 : for (signed char n = 0; n != nn; ++n)
2389 : {
2390 : const Real em_val = this->embedding_matrix
2391 1256512 : (c, cn, n);
2392 1256512 : if (em_val == 1)
2393 : {
2394 89386 : cached_parent_indices[em_vers][c][cn] = n;
2395 83582 : break;
2396 : }
2397 :
2398 1172930 : if (em_val != 0)
2399 : {
2400 190991 : cached_parent_indices[em_vers][c][cn] =
2401 : -1;
2402 178639 : break;
2403 : }
2404 :
2405 : // We should never see an all-zero embedding matrix
2406 : // row
2407 34022 : libmesh_assert_not_equal_to (n+1, nn);
2408 : }
2409 : }
2410 : }
2411 : }
2412 :
2413 : const signed char cache_val =
2414 207435833 : cached_parent_indices[em_vers][child][child_node];
2415 190559429 : if (cache_val == -1)
2416 5315908 : return libMesh::invalid_uint;
2417 :
2418 73936611 : return cached_parent_indices[em_vers][child][child_node];
2419 : }
2420 :
2421 :
2422 :
2423 : const std::vector<std::pair<unsigned char, unsigned char>> &
2424 86672012 : Elem::parent_bracketing_nodes(unsigned int child,
2425 : unsigned int child_node) const
2426 : {
2427 : // Indexed first by embedding matrix type, then by child id, then by
2428 : // child node, then by bracketing pair
2429 : std::vector<std::vector<std::vector<std::vector<std::pair<unsigned char, unsigned char>>>>> &
2430 86672012 : cached_bracketing_nodes = this->_get_bracketing_node_cache();
2431 :
2432 86672012 : const unsigned int em_vers = this->embedding_matrix_version();
2433 :
2434 : // We may be updating the cache on one thread, and while that
2435 : // happens we can't safely access the cache from other threads.
2436 10362868 : Threads::spin_mutex::scoped_lock lock(parent_bracketing_nodes_mutex);
2437 :
2438 91852822 : if (cached_bracketing_nodes.size() <= em_vers)
2439 4989 : cached_bracketing_nodes.resize(em_vers+1);
2440 :
2441 86672012 : const unsigned int nc = this->n_children();
2442 :
2443 : // If we haven't cached the bracketing nodes corresponding to this
2444 : // embedding matrix yet, let's do so now.
2445 97033632 : if (cached_bracketing_nodes[em_vers].size() < nc)
2446 : {
2447 : // If we're a second-order element but we're not a full-order
2448 : // element, then some of our bracketing nodes may not exist
2449 : // except on the equivalent full-order element. Let's build an
2450 : // equivalent full-order element and make a copy of its cache to
2451 : // use.
2452 7955 : if (this->default_order() != FIRST &&
2453 2962 : second_order_equivalent_type(this->type(), /*full_ordered=*/ true) != this->type())
2454 : {
2455 : // Check that we really are the non-full-order type
2456 8 : libmesh_assert_equal_to
2457 : (second_order_equivalent_type (this->type(), false),
2458 : this->type());
2459 :
2460 : // Build the full-order type
2461 : ElemType full_type =
2462 165 : second_order_equivalent_type(this->type(), /*full_ordered=*/ true);
2463 173 : std::unique_ptr<Elem> full_elem = Elem::build(full_type);
2464 :
2465 : // This won't work for elements with multiple
2466 : // embedding_matrix versions, but every such element is full
2467 : // order anyways.
2468 8 : libmesh_assert_equal_to(em_vers, 0);
2469 :
2470 : // Make sure its cache has been built. We temporarily
2471 : // release our mutex lock so that the inner call can
2472 : // re-acquire it.
2473 8 : lock.release();
2474 165 : full_elem->parent_bracketing_nodes(0,0);
2475 :
2476 : // And then we need to lock again, so that if someone *else*
2477 : // grabbed our lock before we did we don't risk accessing
2478 : // cached_bracketing_nodes while they're working on it.
2479 : // Threading is hard.
2480 8 : lock.acquire(parent_bracketing_nodes_mutex);
2481 :
2482 : // Copy its cache
2483 : cached_bracketing_nodes =
2484 165 : full_elem->_get_bracketing_node_cache();
2485 :
2486 : // Now we don't need to build the cache ourselves.
2487 181 : return cached_bracketing_nodes[em_vers][child][child_node];
2488 149 : }
2489 :
2490 5012 : cached_bracketing_nodes[em_vers].resize(nc);
2491 :
2492 4828 : const unsigned int nn = this->n_nodes();
2493 :
2494 : // We have to examine each child
2495 27697 : for (unsigned int c = 0; c != nc; ++c)
2496 : {
2497 22869 : const unsigned int ncn = this->n_nodes_in_child(c);
2498 :
2499 24593 : cached_bracketing_nodes[em_vers][c].resize(ncn);
2500 :
2501 : // We have to examine each node in that child
2502 198666 : for (unsigned int n = 0; n != ncn; ++n)
2503 : {
2504 : // If this child node isn't a vertex or an infinite
2505 : // child element's mid-infinite-edge node, then we need
2506 : // to find bracketing nodes on the child.
2507 175797 : if (!this->is_vertex_on_child(c, n)
2508 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2509 15695 : && !this->is_mid_infinite_edge_node(n)
2510 : #endif
2511 : )
2512 : {
2513 : // Use the embedding matrix to find the child node
2514 : // location in parent master element space
2515 2444 : Point bracketed_pt;
2516 :
2517 1002028 : for (unsigned int pn = 0; pn != nn; ++pn)
2518 : {
2519 : const Real em_val =
2520 933006 : this->embedding_matrix(c,n,pn);
2521 :
2522 33132 : libmesh_assert_not_equal_to (em_val, 1);
2523 933006 : if (em_val != 0.)
2524 337854 : bracketed_pt.add_scaled(this->master_point(pn), em_val);
2525 : }
2526 :
2527 : // Check each pair of nodes on the child which are
2528 : // also both parent nodes
2529 1002028 : for (unsigned int n1 = 0; n1 != ncn; ++n1)
2530 : {
2531 933006 : if (n1 == n)
2532 591746 : continue;
2533 :
2534 : unsigned int parent_n1 =
2535 863984 : this->as_parent_node(c,n1);
2536 :
2537 863984 : if (parent_n1 == libMesh::invalid_uint)
2538 504380 : continue;
2539 :
2540 341260 : Point p1 = this->master_point(parent_n1);
2541 :
2542 3693686 : for (unsigned int n2 = n1+1; n2 < nn; ++n2)
2543 : {
2544 3435604 : if (n2 == n)
2545 2758878 : continue;
2546 :
2547 : unsigned int parent_n2 =
2548 3177522 : this->as_parent_node(c,n2);
2549 :
2550 3177522 : if (parent_n2 == libMesh::invalid_uint)
2551 2413220 : continue;
2552 :
2553 676726 : Point p2 = this->master_point(parent_n2);
2554 :
2555 24900 : Point pmid = (p1 + p2)/2;
2556 :
2557 24900 : if (pmid == bracketed_pt)
2558 : {
2559 91998 : cached_bracketing_nodes[em_vers][c][n].emplace_back(parent_n1, parent_n2);
2560 83178 : break;
2561 : }
2562 : else
2563 21960 : libmesh_assert(!pmid.absolute_fuzzy_equals(bracketed_pt));
2564 : }
2565 : }
2566 : }
2567 : // If this child node is a parent node, we need to
2568 : // find bracketing nodes on the parent.
2569 : else
2570 : {
2571 106775 : unsigned int parent_node = this->as_parent_node(c,n);
2572 :
2573 4090 : Point bracketed_pt;
2574 :
2575 : // If we're not a parent node, use the embedding
2576 : // matrix to find the child node location in parent
2577 : // master element space
2578 106775 : if (parent_node == libMesh::invalid_uint)
2579 : {
2580 351898 : for (unsigned int pn = 0; pn != nn; ++pn)
2581 : {
2582 : const Real em_val =
2583 302681 : this->embedding_matrix(c,n,pn);
2584 :
2585 11828 : libmesh_assert_not_equal_to (em_val, 1);
2586 302681 : if (em_val != 0.)
2587 151178 : bracketed_pt.add_scaled(this->master_point(pn), em_val);
2588 : }
2589 : }
2590 : // If we're a parent node then we need no arithmetic
2591 : else
2592 57558 : bracketed_pt = this->master_point(parent_node);
2593 :
2594 1019062 : for (unsigned int n1 = 0; n1 != nn; ++n1)
2595 : {
2596 912287 : if (n1 == parent_node)
2597 57558 : continue;
2598 :
2599 854729 : Point p1 = this->master_point(n1);
2600 :
2601 4732170 : for (unsigned int n2 = n1+1; n2 < nn; ++n2)
2602 : {
2603 4027118 : if (n2 == parent_node)
2604 10440 : continue;
2605 :
2606 3747376 : Point pmid = (p1 + this->master_point(n2))/2;
2607 :
2608 138240 : if (pmid == bracketed_pt)
2609 : {
2610 166225 : cached_bracketing_nodes[em_vers][c][n].emplace_back(n1, n2);
2611 5516 : break;
2612 : }
2613 : else
2614 132724 : libmesh_assert(!pmid.absolute_fuzzy_equals(bracketed_pt));
2615 : }
2616 : }
2617 : }
2618 : }
2619 : }
2620 : }
2621 :
2622 102214253 : return cached_bracketing_nodes[em_vers][child][child_node];
2623 : }
2624 :
2625 :
2626 : const std::vector<std::pair<dof_id_type, dof_id_type>>
2627 96254666 : Elem::bracketing_nodes(unsigned int child,
2628 : unsigned int child_node) const
2629 : {
2630 5603514 : std::vector<std::pair<dof_id_type, dof_id_type>> returnval;
2631 :
2632 : const std::vector<std::pair<unsigned char, unsigned char>> & pbc =
2633 96254666 : this->parent_bracketing_nodes(child,child_node);
2634 :
2635 200918167 : for (const auto & pb : pbc)
2636 : {
2637 104663501 : const unsigned short n_n = this->n_nodes();
2638 104663501 : if (pb.first < n_n && pb.second < n_n)
2639 116605317 : returnval.emplace_back(this->node_id(pb.first), this->node_id(pb.second));
2640 : else
2641 : {
2642 : // We must be on a non-full-order higher order element...
2643 33008 : libmesh_assert_not_equal_to(this->default_order(), FIRST);
2644 33008 : libmesh_assert_not_equal_to
2645 : (second_order_equivalent_type (this->type(), true),
2646 : this->type());
2647 33008 : libmesh_assert_equal_to
2648 : (second_order_equivalent_type (this->type(), false),
2649 : this->type());
2650 :
2651 : // And that's a shame, because this is a nasty search:
2652 :
2653 : // Build the full-order type
2654 : ElemType full_type =
2655 908404 : second_order_equivalent_type(this->type(), /*full_ordered=*/ true);
2656 941412 : std::unique_ptr<Elem> full_elem = Elem::build(full_type);
2657 :
2658 908404 : dof_id_type pt1 = DofObject::invalid_id;
2659 908404 : dof_id_type pt2 = DofObject::invalid_id;
2660 :
2661 : // Find the bracketing nodes by figuring out what
2662 : // already-created children will have them.
2663 :
2664 : // This only doesn't break horribly because we add children
2665 : // and nodes in straightforward + hierarchical orders...
2666 4993928 : for (unsigned int c=0; c <= child; ++c)
2667 : {
2668 147320 : libmesh_assert(this->child_ptr(c));
2669 4085524 : if (this->child_ptr(c) == remote_elem)
2670 0 : continue;
2671 :
2672 73004408 : for (auto n : make_range(this->n_nodes_in_child(c)))
2673 : {
2674 66036404 : if (c == child && n == child_node)
2675 33008 : break;
2676 :
2677 67475880 : if (pb.first == full_elem->as_parent_node(c,n))
2678 : {
2679 : // We should be consistent
2680 99160 : if (pt1 != DofObject::invalid_id)
2681 66330 : libmesh_assert_equal_to(pt1, this->child_ptr(c)->node_id(n));
2682 :
2683 2846778 : pt1 = this->child_ptr(c)->node_id(n);
2684 : }
2685 :
2686 67475880 : if (pb.second == full_elem->as_parent_node(c,n))
2687 : {
2688 : // We should be consistent
2689 72088 : if (pt2 != DofObject::invalid_id)
2690 40682 : libmesh_assert_equal_to(pt2, this->child_ptr(c)->node_id(n));
2691 :
2692 2068296 : pt2 = this->child_ptr(c)->node_id(n);
2693 : }
2694 : }
2695 : }
2696 :
2697 : // We should *usually* find all bracketing nodes by the time
2698 : // we query them (again, because of the child & node add
2699 : // order)
2700 : //
2701 : // The exception is if we're a HEX20, in which case we will
2702 : // find pairs of vertex nodes and edge nodes bracketing the
2703 : // new central node but we *won't* find the pairs of face
2704 : // nodes which we would have had on a HEX27. In that case
2705 : // we'll still have enough bracketing nodes for a
2706 : // topological lookup, but we won't be able to make the
2707 : // following assertions.
2708 908404 : if (this->type() != HEX20)
2709 : {
2710 15920 : libmesh_assert_not_equal_to (pt1, DofObject::invalid_id);
2711 15920 : libmesh_assert_not_equal_to (pt2, DofObject::invalid_id);
2712 : }
2713 :
2714 908404 : if (pt1 != DofObject::invalid_id &&
2715 903549 : pt2 != DofObject::invalid_id)
2716 864709 : returnval.emplace_back(pt1, pt2);
2717 842388 : }
2718 : }
2719 :
2720 96254666 : return returnval;
2721 : }
2722 : #endif // #ifdef LIBMESH_ENABLE_AMR
2723 :
2724 :
2725 :
2726 :
2727 119184139 : bool Elem::contains_point (const Point & p, Real tol) const
2728 : {
2729 : // We currently allow the user to enlarge the bounding box by
2730 : // providing a tol > TOLERANCE (so this routine is identical to
2731 : // Elem::close_to_point()), but print a warning so that the
2732 : // user can eventually switch his code over to calling close_to_point()
2733 : // instead, which is intended to be used for this purpose.
2734 119184139 : if (tol > TOLERANCE)
2735 : {
2736 0 : libmesh_do_once(libMesh::err
2737 : << "WARNING: Resizing bounding box to match user-specified tolerance!\n"
2738 : << "In the future, calls to Elem::contains_point() with tol > TOLERANCE\n"
2739 : << "will be more optimized, but should not be used\n"
2740 : << "to search for points 'close to' elements!\n"
2741 : << "Instead, use Elem::close_to_point() for this purpose.\n"
2742 : << std::endl;);
2743 0 : return this->point_test(p, tol, tol);
2744 : }
2745 : else
2746 119184139 : return this->point_test(p, TOLERANCE, tol);
2747 : }
2748 :
2749 :
2750 :
2751 :
2752 2231639 : bool Elem::close_to_point (const Point & p, Real tol) const
2753 : {
2754 : // This test uses the user's passed-in tolerance for the
2755 : // bounding box test as well, thereby allowing the routine to
2756 : // find points which are not only "in" the element, but also
2757 : // "nearby" to within some tolerance.
2758 2231639 : return this->point_test(p, tol, tol);
2759 : }
2760 :
2761 :
2762 :
2763 :
2764 121415778 : bool Elem::point_test(const Point & p, Real box_tol, Real map_tol) const
2765 : {
2766 6554687 : libmesh_assert_greater (box_tol, 0.);
2767 6554687 : libmesh_assert_greater (map_tol, 0.);
2768 :
2769 : // This is a great optimization on first order elements, but it
2770 : // could return false negatives on higher orders
2771 121415778 : if (this->default_order() == FIRST)
2772 : {
2773 : // Check to make sure the element *could* contain this point, so we
2774 : // can avoid an expensive inverse_map call if it doesn't.
2775 : bool
2776 : #if LIBMESH_DIM > 2
2777 3253530 : point_above_min_z = false,
2778 3253530 : point_below_max_z = false,
2779 : #endif
2780 : #if LIBMESH_DIM > 1
2781 3253530 : point_above_min_y = false,
2782 3253530 : point_below_max_y = false,
2783 : #endif
2784 3253530 : point_above_min_x = false,
2785 3253530 : point_below_max_x = false;
2786 :
2787 : // For relative bounding box checks in physical space
2788 61725134 : const Real my_hmax = this->hmax();
2789 :
2790 489535647 : for (auto & n : this->node_ref_range())
2791 : {
2792 427810513 : point_above_min_x = point_above_min_x || (n(0) - my_hmax*box_tol <= p(0));
2793 427810513 : point_below_max_x = point_below_max_x || (n(0) + my_hmax*box_tol >= p(0));
2794 : #if LIBMESH_DIM > 1
2795 427810513 : point_above_min_y = point_above_min_y || (n(1) - my_hmax*box_tol <= p(1));
2796 427810513 : point_below_max_y = point_below_max_y || (n(1) + my_hmax*box_tol >= p(1));
2797 : #endif
2798 : #if LIBMESH_DIM > 2
2799 427810513 : point_above_min_z = point_above_min_z || (n(2) - my_hmax*box_tol <= p(2));
2800 427810513 : point_below_max_z = point_below_max_z || (n(2) + my_hmax*box_tol >= p(2));
2801 : #endif
2802 : }
2803 :
2804 61725134 : if (
2805 : #if LIBMESH_DIM > 2
2806 3253530 : !point_above_min_z ||
2807 3136352 : !point_below_max_z ||
2808 : #endif
2809 : #if LIBMESH_DIM > 1
2810 39831121 : !point_above_min_y ||
2811 2064481 : !point_below_max_y ||
2812 : #endif
2813 17707698 : !point_above_min_x ||
2814 860216 : !point_below_max_x)
2815 2547843 : return false;
2816 : }
2817 :
2818 : // To be on the safe side, we converge the inverse_map() iteration
2819 : // to a slightly tighter tolerance than that requested by the
2820 : // user...
2821 71861284 : const Point mapped_point = FEMap::inverse_map(this->dim(),
2822 : this,
2823 : p,
2824 : 0.1*map_tol, // <- this is |dx| tolerance, the Newton residual should be ~ |dx|^2
2825 : /*secure=*/ false,
2826 7935342 : /*extra_checks=*/ false);
2827 :
2828 : // Check that the refspace point maps back to p! This is only necessary
2829 : // for 1D and 2D elements, 3D elements always live in 3D.
2830 : //
2831 : // TODO: The contains_point() function could most likely be implemented
2832 : // more efficiently in the element sub-classes themselves, at least for
2833 : // the linear element types.
2834 67854440 : if (this->dim() < 3)
2835 : {
2836 26499849 : Point xyz = FEMap::map(this->dim(), this, mapped_point);
2837 :
2838 : // Compute the distance between the original point and the re-mapped point.
2839 : // They should be in the same place.
2840 24536919 : Real dist = (xyz - p).norm();
2841 :
2842 :
2843 : // If dist is larger than some fraction of the tolerance, then return false.
2844 : // This can happen when e.g. a 2D element is living in 3D, and
2845 : // FEMap::inverse_map() maps p onto the projection of the element,
2846 : // effectively "tricking" on_reference_element().
2847 26499849 : if (dist > this->hmax() * map_tol)
2848 8449 : return false;
2849 : }
2850 :
2851 67845991 : return this->on_reference_element(mapped_point, map_tol);
2852 : }
2853 :
2854 :
2855 :
2856 :
2857 78600 : bool Elem::has_invertible_map(Real /*tol*/) const
2858 : {
2859 85150 : QNodal qnodal {this->dim()};
2860 91700 : FEMap fe_map;
2861 6550 : auto & jac = fe_map.get_jacobian();
2862 :
2863 : // We have a separate check for is_singular_node() below, so in this
2864 : // case its "OK" to do nodal quadrature on pyramids.
2865 78600 : qnodal.allow_nodal_pyramid_quadrature = true;
2866 78600 : qnodal.init(*this);
2867 6550 : auto & qp = qnodal.get_points();
2868 6550 : libmesh_assert_equal_to(qp.size(), this->n_nodes());
2869 :
2870 85150 : std::vector<Point> one_point(1);
2871 85150 : std::vector<Real> one_weight(1,1);
2872 1086900 : for (auto i : index_range(qp))
2873 : {
2874 1008300 : if (this->is_singular_node(i))
2875 8448 : continue;
2876 :
2877 999084 : one_point[0] = qp[i];
2878 :
2879 999084 : fe_map.init_reference_to_physical_map(this->dim(), one_point, this);
2880 999084 : fe_map.compute_map(this->dim(), one_weight, this, false);
2881 :
2882 999084 : if (jac[0] <= 0)
2883 0 : return false;
2884 : }
2885 :
2886 13100 : return true;
2887 65500 : }
2888 :
2889 :
2890 :
2891 0 : void Elem::print_info (std::ostream & os) const
2892 : {
2893 0 : os << this->get_info()
2894 0 : << std::endl;
2895 0 : }
2896 :
2897 :
2898 :
2899 0 : std::string Elem::get_info () const
2900 : {
2901 0 : std::ostringstream oss;
2902 :
2903 : oss << " Elem Information" << '\n'
2904 0 : << " id()=";
2905 :
2906 0 : if (this->valid_id())
2907 0 : oss << this->id();
2908 : else
2909 0 : oss << "invalid";
2910 :
2911 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
2912 0 : oss << ", unique_id()=";
2913 0 : if (this->valid_unique_id())
2914 0 : oss << this->unique_id();
2915 : else
2916 0 : oss << "invalid";
2917 : #endif
2918 :
2919 0 : oss << ", subdomain_id()=" << this->subdomain_id();
2920 0 : oss << ", processor_id()=" << this->processor_id() << '\n';
2921 :
2922 0 : oss << " type()=" << Utility::enum_to_string(this->type()) << '\n'
2923 0 : << " dim()=" << this->dim() << '\n'
2924 0 : << " n_nodes()=" << this->n_nodes() << '\n';
2925 :
2926 0 : oss << " mapping=" << Utility::enum_to_string(this->mapping_type()) << '\n';
2927 :
2928 0 : for (auto n : this->node_index_range())
2929 : {
2930 0 : oss << " " << n << this->node_ref(n);
2931 0 : if (this->mapping_type() == RATIONAL_BERNSTEIN_MAP)
2932 : {
2933 0 : const unsigned char datum_index = this->mapping_data();
2934 0 : oss << " weight=" <<
2935 0 : this->node_ref(n).get_extra_datum<Real>(datum_index) << '\n';
2936 : }
2937 : }
2938 :
2939 0 : oss << " n_sides()=" << this->n_sides() << '\n';
2940 :
2941 0 : for (auto s : this->side_index_range())
2942 : {
2943 0 : oss << " neighbor(" << s << ")=";
2944 0 : if (this->neighbor_ptr(s))
2945 0 : oss << this->neighbor_ptr(s)->id() << '\n';
2946 : else
2947 0 : oss << "nullptr\n";
2948 : }
2949 :
2950 0 : if (!this->infinite())
2951 : {
2952 0 : oss << " hmin()=" << this->hmin()
2953 0 : << ", hmax()=" << this->hmax() << '\n'
2954 0 : << " volume()=" << this->volume() << '\n';
2955 : }
2956 0 : oss << " active()=" << this->active()
2957 0 : << ", ancestor()=" << this->ancestor()
2958 0 : << ", subactive()=" << this->subactive()
2959 0 : << ", has_children()=" << this->has_children() << '\n'
2960 0 : << " parent()=";
2961 0 : if (this->parent())
2962 0 : oss << this->parent()->id() << '\n';
2963 : else
2964 0 : oss << "nullptr\n";
2965 0 : oss << " level()=" << this->level()
2966 0 : << ", p_level()=" << this->p_level() << '\n'
2967 : #ifdef LIBMESH_ENABLE_AMR
2968 0 : << " refinement_flag()=" << Utility::enum_to_string(this->refinement_flag()) << '\n'
2969 0 : << " p_refinement_flag()=" << Utility::enum_to_string(this->p_refinement_flag()) << '\n'
2970 : #endif
2971 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2972 0 : << " infinite()=" << this->infinite() << '\n';
2973 0 : if (this->infinite())
2974 0 : oss << " origin()=" << this->origin() << '\n'
2975 : #endif
2976 : ;
2977 :
2978 0 : oss << " DoFs=";
2979 0 : for (auto s : make_range(this->n_systems()))
2980 0 : for (auto v : make_range(this->n_vars(s)))
2981 0 : for (auto c : make_range(this->n_comp(s,v)))
2982 0 : oss << '(' << s << '/' << v << '/' << this->dof_number(s,v,c) << ") ";
2983 :
2984 :
2985 0 : return oss.str();
2986 0 : }
2987 :
2988 :
2989 :
2990 0 : void Elem::nullify_neighbors ()
2991 : {
2992 : // Tell any of my neighbors about my death...
2993 : // Looks strange, huh?
2994 0 : for (auto n : this->side_index_range())
2995 : {
2996 0 : Elem * current_neighbor = this->neighbor_ptr(n);
2997 0 : if (current_neighbor && current_neighbor != remote_elem)
2998 : {
2999 : // Note: it is possible that I see the neighbor
3000 : // (which is coarser than me)
3001 : // but they don't see me, so avoid that case.
3002 0 : if (current_neighbor->level() == this->level())
3003 : {
3004 0 : const unsigned int w_n_a_i = current_neighbor->which_neighbor_am_i(this);
3005 0 : libmesh_assert_less (w_n_a_i, current_neighbor->n_neighbors());
3006 0 : current_neighbor->set_neighbor(w_n_a_i, nullptr);
3007 0 : this->set_neighbor(n, nullptr);
3008 : }
3009 : }
3010 : }
3011 0 : }
3012 :
3013 :
3014 :
3015 :
3016 0 : unsigned int Elem::n_second_order_adjacent_vertices (const unsigned int) const
3017 : {
3018 : // for linear elements, always return 0
3019 0 : return 0;
3020 : }
3021 :
3022 :
3023 :
3024 0 : unsigned short int Elem::second_order_adjacent_vertex (const unsigned int,
3025 : const unsigned int) const
3026 : {
3027 : // for linear elements, always return 0
3028 0 : return 0;
3029 : }
3030 :
3031 :
3032 :
3033 : std::pair<unsigned short int, unsigned short int>
3034 0 : Elem::second_order_child_vertex (const unsigned int) const
3035 : {
3036 : // for linear elements, always return 0
3037 0 : return std::pair<unsigned short int, unsigned short int>(0,0);
3038 : }
3039 :
3040 :
3041 :
3042 271926 : ElemType Elem::first_order_equivalent_type (const ElemType et)
3043 : {
3044 26644 : switch (et)
3045 : {
3046 28 : case NODEELEM:
3047 28 : return NODEELEM;
3048 3118 : case EDGE2:
3049 : case EDGE3:
3050 : case EDGE4:
3051 3118 : return EDGE2;
3052 448 : case TRI3:
3053 : case TRI6:
3054 : case TRI7:
3055 448 : return TRI3;
3056 0 : case TRISHELL3:
3057 0 : return TRISHELL3;
3058 22482 : case QUAD4:
3059 : case QUAD8:
3060 : case QUAD9:
3061 22482 : return QUAD4;
3062 0 : case QUADSHELL4:
3063 : case QUADSHELL8:
3064 : case QUADSHELL9:
3065 0 : return QUADSHELL4;
3066 0 : case TET4:
3067 : case TET10:
3068 : case TET14:
3069 0 : return TET4;
3070 0 : case HEX8:
3071 : case HEX27:
3072 : case HEX20:
3073 0 : return HEX8;
3074 480 : case PRISM6:
3075 : case PRISM15:
3076 : case PRISM18:
3077 : case PRISM20:
3078 : case PRISM21:
3079 480 : return PRISM6;
3080 0 : case PYRAMID5:
3081 : case PYRAMID13:
3082 : case PYRAMID14:
3083 : case PYRAMID18:
3084 0 : return PYRAMID5;
3085 :
3086 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
3087 :
3088 16 : case INFEDGE2:
3089 16 : return INFEDGE2;
3090 72 : case INFQUAD4:
3091 : case INFQUAD6:
3092 72 : return INFQUAD4;
3093 0 : case INFHEX8:
3094 : case INFHEX16:
3095 : case INFHEX18:
3096 0 : return INFHEX8;
3097 0 : case INFPRISM6:
3098 : case INFPRISM12:
3099 0 : return INFPRISM6;
3100 :
3101 : #endif
3102 :
3103 0 : default:
3104 : // unknown element
3105 0 : return INVALID_ELEM;
3106 : }
3107 : }
3108 :
3109 :
3110 :
3111 2826009 : ElemType Elem::second_order_equivalent_type (const ElemType et,
3112 : const bool full_ordered)
3113 : {
3114 2826009 : switch (et)
3115 : {
3116 2 : case NODEELEM:
3117 2 : return NODEELEM;
3118 3806 : case EDGE2:
3119 : case EDGE3:
3120 : {
3121 : // full_ordered not relevant
3122 3806 : return EDGE3;
3123 : }
3124 :
3125 0 : case EDGE4:
3126 : {
3127 : // full_ordered not relevant
3128 0 : return EDGE4;
3129 : }
3130 :
3131 6359 : case TRI3:
3132 : case TRI6:
3133 : {
3134 : // full_ordered not relevant
3135 6359 : return TRI6;
3136 : }
3137 :
3138 0 : case TRI7:
3139 0 : return TRI7;
3140 :
3141 : // Currently there is no TRISHELL6, so similarly to other types
3142 : // where this is the case, we just return the input.
3143 0 : case TRISHELL3:
3144 0 : return TRISHELL3;
3145 :
3146 18948 : case QUAD4:
3147 : case QUAD8:
3148 : {
3149 18948 : if (full_ordered)
3150 1624 : return QUAD9;
3151 : else
3152 466 : return QUAD8;
3153 : }
3154 :
3155 8802 : case QUADSHELL4:
3156 : case QUADSHELL8:
3157 : {
3158 8802 : if (full_ordered)
3159 896 : return QUADSHELL9;
3160 : else
3161 448 : return QUADSHELL8;
3162 : }
3163 :
3164 250 : case QUAD9:
3165 : {
3166 : // full_ordered not relevant
3167 250 : return QUAD9;
3168 : }
3169 :
3170 2 : case QUADSHELL9:
3171 : {
3172 : // full_ordered not relevant
3173 2 : return QUADSHELL9;
3174 : }
3175 :
3176 1362240 : case TET4:
3177 : case TET10:
3178 : {
3179 : // full_ordered not relevant
3180 1362240 : return TET10;
3181 : }
3182 :
3183 0 : case TET14:
3184 0 : return TET14;
3185 :
3186 599508 : case HEX8:
3187 : case HEX20:
3188 : {
3189 : // see below how this correlates with INFHEX8
3190 599508 : if (full_ordered)
3191 39324 : return HEX27;
3192 : else
3193 17090 : return HEX20;
3194 : }
3195 :
3196 460 : case HEX27:
3197 : {
3198 : // full_ordered not relevant
3199 460 : return HEX27;
3200 : }
3201 :
3202 456152 : case PRISM6:
3203 : case PRISM15:
3204 : {
3205 456152 : if (full_ordered)
3206 29598 : return PRISM18;
3207 : else
3208 14786 : return PRISM15;
3209 : }
3210 :
3211 : // full_ordered not relevant, already fully second order
3212 12 : case PRISM18:
3213 12 : return PRISM18;
3214 0 : case PRISM20:
3215 0 : return PRISM20;
3216 0 : case PRISM21:
3217 0 : return PRISM21;
3218 0 : case PYRAMID18:
3219 0 : return PYRAMID18;
3220 :
3221 360964 : case PYRAMID5:
3222 : case PYRAMID13:
3223 : {
3224 360964 : if (full_ordered)
3225 5084 : return PYRAMID14;
3226 : else
3227 10168 : return PYRAMID13;
3228 : }
3229 :
3230 0 : case PYRAMID14:
3231 : {
3232 : // full_ordered not relevant
3233 0 : return PYRAMID14;
3234 : }
3235 :
3236 :
3237 :
3238 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
3239 :
3240 : // infinite elements
3241 0 : case INFEDGE2:
3242 : {
3243 0 : return INFEDGE2;
3244 : }
3245 :
3246 5 : case INFQUAD4:
3247 : case INFQUAD6:
3248 : {
3249 : // full_ordered not relevant
3250 5 : return INFQUAD6;
3251 : }
3252 :
3253 1032 : case INFHEX8:
3254 : case INFHEX16:
3255 : {
3256 : /*
3257 : * Note that this matches with \p Hex8:
3258 : * For full-ordered, \p InfHex18 and \p Hex27
3259 : * belong together, and for not full-ordered,
3260 : * \p InfHex16 and \p Hex20 belong together.
3261 : */
3262 1032 : if (full_ordered)
3263 452 : return INFHEX18;
3264 : else
3265 226 : return INFHEX16;
3266 : }
3267 :
3268 4 : case INFHEX18:
3269 : {
3270 : // full_ordered not relevant
3271 4 : return INFHEX18;
3272 : }
3273 :
3274 5 : case INFPRISM6:
3275 : case INFPRISM12:
3276 : {
3277 : // full_ordered not relevant
3278 5 : return INFPRISM12;
3279 : }
3280 :
3281 : #endif
3282 :
3283 :
3284 0 : default:
3285 : {
3286 : // what did we miss?
3287 0 : libmesh_error_msg("No second order equivalent element type for et = "
3288 : << Utility::enum_to_string(et));
3289 : }
3290 : }
3291 : }
3292 :
3293 :
3294 :
3295 4958683 : ElemType Elem::complete_order_equivalent_type (const ElemType et)
3296 : {
3297 4958683 : switch (et)
3298 : {
3299 0 : case NODEELEM:
3300 0 : return NODEELEM;
3301 1278 : case EDGE2:
3302 : case EDGE3:
3303 1278 : return EDGE3;
3304 :
3305 0 : case EDGE4:
3306 0 : return EDGE4;
3307 :
3308 5524 : case TRI3:
3309 : case TRI6:
3310 : case TRI7:
3311 5524 : return TRI7;
3312 :
3313 : // Currently there is no TRISHELL6, so similarly to other types
3314 : // where this is the case, we just return the input.
3315 0 : case TRISHELL3:
3316 0 : return TRISHELL3;
3317 :
3318 900 : case QUAD4:
3319 : case QUAD8:
3320 : case QUAD9:
3321 900 : return QUAD9;
3322 :
3323 0 : case QUADSHELL4:
3324 : case QUADSHELL8:
3325 : case QUADSHELL9:
3326 0 : return QUADSHELL9;
3327 :
3328 4776099 : case TET4:
3329 : case TET10:
3330 : case TET14:
3331 4776099 : return TET14;
3332 :
3333 0 : case HEX8:
3334 : case HEX20:
3335 : case HEX27:
3336 0 : return HEX27;
3337 :
3338 : // We don't strictly need the 21st node for DoFs, but some code
3339 : // depends on it for e.g. refinement patterns
3340 32953 : case PRISM6:
3341 : case PRISM15:
3342 : case PRISM18:
3343 : case PRISM20:
3344 : case PRISM21:
3345 32953 : return PRISM21;
3346 :
3347 141929 : case PYRAMID5:
3348 : case PYRAMID13:
3349 : case PYRAMID14:
3350 : case PYRAMID18:
3351 141929 : return PYRAMID18;
3352 :
3353 :
3354 :
3355 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
3356 0 : case INFEDGE2:
3357 0 : return INFEDGE2;
3358 :
3359 0 : case INFQUAD4:
3360 : case INFQUAD6:
3361 0 : return INFQUAD6;
3362 :
3363 0 : case INFHEX8:
3364 : case INFHEX16:
3365 : case INFHEX18:
3366 0 : return INFHEX18;
3367 :
3368 : // Probably ought to implement INFPRISM13 and/or 14; until then
3369 : // we'll just return the not-complete-order ElemType, in
3370 : // accordance which what seems like bad precedent...
3371 0 : case INFPRISM6:
3372 : case INFPRISM12:
3373 0 : return INFPRISM12;
3374 : #endif // LIBMESH_ENABLE_INFINITE_ELEMENTS
3375 :
3376 0 : default:
3377 : {
3378 : // what did we miss?
3379 0 : libmesh_error_msg("No complete order equivalent element type for et = "
3380 : << Utility::enum_to_string(et));
3381 : }
3382 : }
3383 : }
3384 :
3385 :
3386 :
3387 0 : Elem::side_iterator Elem::boundary_sides_begin()
3388 : {
3389 0 : Predicates::BoundarySide<SideIter> bsp;
3390 0 : return side_iterator(this->_first_side(), this->_last_side(), bsp);
3391 : }
3392 :
3393 :
3394 :
3395 :
3396 0 : Elem::side_iterator Elem::boundary_sides_end()
3397 : {
3398 0 : Predicates::BoundarySide<SideIter> bsp;
3399 0 : return side_iterator(this->_last_side(), this->_last_side(), bsp);
3400 : }
3401 :
3402 :
3403 :
3404 :
3405 8697 : Real Elem::volume () const
3406 : {
3407 : // The default implementation builds a finite element of the correct
3408 : // order and sums up the JxW contributions. This can be expensive,
3409 : // so the various element types can overload this method and compute
3410 : // the volume more efficiently.
3411 8697 : const FEFamily mapping_family = FEMap::map_fe_type(*this);
3412 8697 : const FEType fe_type(this->default_order(), mapping_family);
3413 :
3414 8697 : std::unique_ptr<FEBase> fe (FEBase::build(this->dim(),
3415 9017 : fe_type));
3416 :
3417 : // Use a map with a negative Jacobian tolerance, in case we're asked
3418 : // for the net volume of a tangled element
3419 320 : fe->get_fe_map().set_jacobian_tolerance(std::numeric_limits<Real>::lowest());
3420 :
3421 806 : const std::vector<Real> & JxW = fe->get_JxW();
3422 :
3423 : // The default quadrature rule should integrate the mass matrix,
3424 : // thus it should be plenty to compute the area
3425 8697 : QGauss qrule (this->dim(), fe_type.default_quadrature_order());
3426 :
3427 8697 : fe->attach_quadrature_rule(&qrule);
3428 :
3429 8697 : fe->reinit(this);
3430 :
3431 320 : Real vol=0.;
3432 164982 : for (auto jxw : JxW)
3433 156285 : vol += jxw;
3434 :
3435 9017 : return vol;
3436 :
3437 8067 : }
3438 :
3439 :
3440 :
3441 43166808 : BoundingBox Elem::loose_bounding_box () const
3442 : {
3443 43166808 : Point pmin = this->point(0);
3444 43166808 : Point pmax = pmin;
3445 :
3446 43166808 : unsigned int n_points = this->n_nodes();
3447 225387718 : for (unsigned int p=0; p != n_points; ++p)
3448 728883640 : for (unsigned d=0; d<LIBMESH_DIM; ++d)
3449 : {
3450 83243682 : const Point & pt = this->point(p);
3451 546662730 : if (pmin(d) > pt(d))
3452 40622849 : pmin(d) = pt(d);
3453 :
3454 546662730 : if (pmax(d) < pt(d))
3455 89207256 : pmax(d) = pt(d);
3456 : }
3457 :
3458 46602699 : return BoundingBox(pmin, pmax);
3459 : }
3460 :
3461 : Point
3462 710 : Elem::side_vertex_average_normal(const unsigned int s) const
3463 : {
3464 710 : unsigned int dim = this->dim();
3465 710 : const std::unique_ptr<const Elem> face = this->build_side_ptr(s);
3466 : std::unique_ptr<libMesh::FEBase> fe(
3467 730 : libMesh::FEBase::build(dim, libMesh::FEType(this->default_order())));
3468 50 : const std::vector<Point> & normals = fe->get_normals();
3469 710 : std::vector<Point> ref_side_vertex_average_v = {face->reference_elem()->vertex_average()};
3470 710 : fe->reinit(this, s, TOLERANCE, &ref_side_vertex_average_v);
3471 1420 : return normals[0];
3472 670 : }
3473 :
3474 3432141 : bool Elem::is_vertex_on_parent(unsigned int c,
3475 : unsigned int n) const
3476 : {
3477 : #ifdef LIBMESH_ENABLE_AMR
3478 :
3479 3432141 : unsigned int my_n_vertices = this->n_vertices();
3480 17139651 : for (unsigned int n_parent = 0; n_parent != my_n_vertices;
3481 : ++n_parent)
3482 16896454 : if (this->node_ptr(n_parent) == this->child_ptr(c)->node_ptr(n))
3483 55241 : return true;
3484 256923 : return false;
3485 :
3486 : #else
3487 :
3488 : // No AMR?
3489 : libmesh_ignore(c,n);
3490 : libmesh_error_msg("ERROR: AMR disabled, how did we get here?");
3491 : return true;
3492 :
3493 : #endif
3494 : }
3495 :
3496 :
3497 :
3498 0 : unsigned int Elem::opposite_side(const unsigned int /*s*/) const
3499 : {
3500 : // If the subclass didn't rederive this, using it is an error
3501 0 : libmesh_not_implemented();
3502 : }
3503 :
3504 :
3505 :
3506 0 : unsigned int Elem::opposite_node(const unsigned int /*n*/,
3507 : const unsigned int /*s*/) const
3508 : {
3509 : // If the subclass didn't rederive this, using it is an error
3510 0 : libmesh_not_implemented();
3511 : }
3512 :
3513 :
3514 37659 : unsigned int Elem::center_node_on_side(const unsigned short libmesh_dbg_var(side)) const
3515 : {
3516 3160 : libmesh_assert_less (side, this->n_sides());
3517 37659 : return invalid_uint;
3518 : }
3519 :
3520 :
3521 185667 : void Elem::swap2boundarysides(unsigned short s1,
3522 : unsigned short s2,
3523 : BoundaryInfo * boundary_info) const
3524 : {
3525 18780 : std::vector<boundary_id_type> ids1, ids2;
3526 185667 : boundary_info->boundary_ids(this, s1, ids1);
3527 185667 : boundary_info->boundary_ids(this, s2, ids2);
3528 185667 : boundary_info->remove_side(this, s1);
3529 185667 : boundary_info->remove_side(this, s2);
3530 185667 : if (!ids1.empty())
3531 7713 : boundary_info->add_side(this, s2, ids1);
3532 185667 : if (!ids2.empty())
3533 4392 : boundary_info->add_side(this, s1, ids2);
3534 185667 : }
3535 :
3536 :
3537 367644 : void Elem::swap2boundaryedges(unsigned short e1,
3538 : unsigned short e2,
3539 : BoundaryInfo * boundary_info) const
3540 : {
3541 37788 : std::vector<boundary_id_type> ids1, ids2;
3542 367644 : boundary_info->edge_boundary_ids(this, e1, ids1);
3543 367644 : boundary_info->edge_boundary_ids(this, e2, ids2);
3544 367644 : boundary_info->remove_edge(this, e1);
3545 367644 : boundary_info->remove_edge(this, e2);
3546 367644 : if (!ids1.empty())
3547 0 : boundary_info->add_edge(this, e2, ids1);
3548 367644 : if (!ids2.empty())
3549 0 : boundary_info->add_edge(this, e1, ids2);
3550 367644 : }
3551 :
3552 : bool
3553 172235 : Elem::is_internal(const unsigned int i) const
3554 : {
3555 172235 : switch (this->dim())
3556 : {
3557 1 : case 0:
3558 1 : return false;
3559 :
3560 216 : case 1:
3561 216 : return !this->is_vertex(i);
3562 :
3563 57654 : case 2:
3564 57654 : return !this->is_vertex(i) && !this->is_edge(i);
3565 :
3566 114353 : case 3:
3567 114353 : return !this->is_vertex(i) && !this->is_edge(i) && !this->is_face(i);
3568 :
3569 0 : default:
3570 0 : libmesh_error_msg("impossible element dimension " << std::to_string(this->dim()));
3571 : return 0;
3572 : }
3573 : }
3574 :
3575 : bool
3576 689829030 : Elem::positive_edge_orientation(const unsigned int i) const
3577 : {
3578 56442160 : libmesh_assert_less (i, this->n_edges());
3579 :
3580 689829030 : return this->point(this->local_edge_node(i, 0)) >
3581 746271190 : this->point(this->local_edge_node(i, 1));
3582 : }
3583 :
3584 : bool
3585 902687688 : Elem::positive_face_orientation(const unsigned int i) const
3586 : {
3587 75970818 : libmesh_assert_less (i, this->n_faces());
3588 :
3589 : // Get the number of vertices N of face i. Note that for 3d elements, i.e.
3590 : // elements for which this->n_faces() > 0, the number of vertices on any of
3591 : // its sides (or faces) is just the number of that face's sides (or edges).
3592 902687688 : auto side_i = this->side_ptr(i);
3593 902687688 : const unsigned int N = side_i->n_sides();
3594 :
3595 978658506 : const std::vector<unsigned int> nodes = this->nodes_on_side(i);
3596 :
3597 3221577036 : auto cmp = [&](const unsigned int & m, const unsigned int & n) -> bool
3598 4784738490 : { return this->point(m) < this->point(n); };
3599 :
3600 902687688 : const unsigned int V = std::distance(nodes.begin(),
3601 902687688 : std::min_element(nodes.begin(), nodes.begin() + N, cmp));
3602 :
3603 1957317012 : return cmp(nodes[(V + N - 1) % N], nodes[(V + 1) % N]);
3604 750746052 : }
3605 :
3606 : } // namespace libMesh
|
