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