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