Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 : // libmesh includes
19 : #include "libmesh/fe.h"
20 : #include "libmesh/libmesh_logging.h"
21 : #include "libmesh/enum_elem_type.h"
22 : #include "libmesh/boundary_info.h"
23 : #include "libmesh/mesh_base.h"
24 : #include "libmesh/dense_matrix.h"
25 : #include "libmesh/dense_vector.h"
26 : #include "libmesh/dof_map.h"
27 : #include "libmesh/elem.h"
28 : #include "libmesh/fe_interface.h"
29 : #include "libmesh/numeric_vector.h"
30 : #include "libmesh/periodic_boundaries.h"
31 : #include "libmesh/periodic_boundary.h"
32 : #include "libmesh/quadrature.h"
33 : #include "libmesh/quadrature_gauss.h"
34 : #include "libmesh/reference_elem.h"
35 : #include "libmesh/remote_elem.h"
36 : #include "libmesh/tensor_value.h"
37 : #include "libmesh/threads.h"
38 : #include "libmesh/enum_to_string.h"
39 :
40 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
41 : #include "libmesh/inf_fe.h"
42 : #include "libmesh/fe_interface_macros.h"
43 : #endif
44 :
45 : namespace libMesh
46 : {
47 :
48 16954969 : FEAbstract::FEAbstract(const unsigned int d,
49 16954969 : const FEType & fet) :
50 15266874 : _fe_map( FEMap::build(fet) ),
51 15266874 : dim(d),
52 15266874 : calculations_started(false),
53 15266874 : calculate_dual(false),
54 15266874 : calculate_default_dual_coeff(true),
55 15266874 : calculate_nothing(false),
56 15266874 : calculate_map(false),
57 15266874 : calculate_phi(false),
58 15266874 : calculate_dphi(false),
59 15266874 : calculate_d2phi(false),
60 15266874 : calculate_curl_phi(false),
61 15266874 : calculate_div_phi(false),
62 15266874 : calculate_dphiref(false),
63 15266874 : fe_type(fet),
64 15266874 : _elem_type(INVALID_ELEM),
65 15266874 : _elem(nullptr),
66 15266874 : _elem_p_level(0),
67 15266874 : _p_level(0),
68 15266874 : qrule(nullptr),
69 15266874 : shapes_on_quadrature(false),
70 15266874 : _n_total_qp(0),
71 16954969 : _add_p_level_in_reinit(true)
72 : {
73 16954969 : }
74 :
75 :
76 15266874 : FEAbstract::~FEAbstract() = default;
77 :
78 :
79 11091653 : std::unique_ptr<FEAbstract> FEAbstract::build(const unsigned int dim,
80 : const FEType & fet)
81 : {
82 11091653 : switch (dim)
83 : {
84 : // 0D
85 124440 : case 0:
86 : {
87 124440 : switch (fet.family)
88 : {
89 0 : case CLOUGH:
90 0 : return std::make_unique<FE<0,CLOUGH>>(fet);
91 :
92 0 : case HERMITE:
93 0 : return std::make_unique<FE<0,HERMITE>>(fet);
94 :
95 74680 : case LAGRANGE:
96 74680 : return std::make_unique<FE<0,LAGRANGE>>(fet);
97 :
98 0 : case LAGRANGE_VEC:
99 0 : return std::make_unique<FE<0,LAGRANGE_VEC>>(fet);
100 :
101 1680 : case L2_LAGRANGE:
102 1680 : return std::make_unique<FE<0,L2_LAGRANGE>>(fet);
103 :
104 1680 : case L2_LAGRANGE_VEC:
105 1680 : return std::make_unique<FE<0,L2_LAGRANGE_VEC>>(fet);
106 :
107 0 : case HIERARCHIC_VEC:
108 0 : return std::make_unique<FE<0,HIERARCHIC_VEC>>(fet);
109 :
110 0 : case HIERARCHIC:
111 0 : return std::make_unique<FE<0,HIERARCHIC>>(fet);
112 :
113 2272 : case L2_HIERARCHIC:
114 2272 : return std::make_unique<FE<0,L2_HIERARCHIC>>(fet);
115 :
116 0 : case L2_HIERARCHIC_VEC:
117 0 : return std::make_unique<FE<0,L2_HIERARCHIC_VEC>>(fet);
118 :
119 1680 : case SIDE_HIERARCHIC:
120 1680 : return std::make_unique<FE<0,SIDE_HIERARCHIC>>(fet);
121 :
122 6816 : case MONOMIAL:
123 6816 : return std::make_unique<FE<0,MONOMIAL>>(fet);
124 :
125 0 : case MONOMIAL_VEC:
126 0 : return std::make_unique<FE<0,MONOMIAL_VEC>>(fet);
127 :
128 : #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
129 0 : case SZABAB:
130 0 : return std::make_unique<FE<0,SZABAB>>(fet);
131 :
132 0 : case BERNSTEIN:
133 0 : return std::make_unique<FE<0,BERNSTEIN>>(fet);
134 :
135 22024 : case RATIONAL_BERNSTEIN:
136 22024 : return std::make_unique<FE<0,RATIONAL_BERNSTEIN>>(fet);
137 : #endif
138 :
139 0 : case XYZ:
140 0 : return std::make_unique<FEXYZ<0>>(fet);
141 :
142 7928 : case SCALAR:
143 7928 : return std::make_unique<FEScalar<0>>(fet);
144 :
145 0 : case NEDELEC_ONE:
146 0 : return std::make_unique<FENedelecOne<0>>(fet);
147 :
148 5680 : case RAVIART_THOMAS:
149 5680 : return std::make_unique<FERaviartThomas<0>>(fet);
150 :
151 0 : case L2_RAVIART_THOMAS:
152 0 : return std::make_unique<FEL2RaviartThomas<0>>(fet);
153 :
154 0 : case SUBDIVISION:
155 0 : return std::make_unique<FESubdivision>(fet);
156 :
157 0 : default:
158 0 : libmesh_error_msg("ERROR: Bad FEType.family= " << Utility::enum_to_string(fet.family));
159 : }
160 : }
161 : // 1D
162 1093544 : case 1:
163 : {
164 1093544 : switch (fet.family)
165 : {
166 0 : case CLOUGH:
167 0 : return std::make_unique<FE<1,CLOUGH>>(fet);
168 :
169 273696 : case HERMITE:
170 273696 : return std::make_unique<FE<1,HERMITE>>(fet);
171 :
172 128500 : case LAGRANGE:
173 128500 : return std::make_unique<FE<1,LAGRANGE>>(fet);
174 :
175 0 : case LAGRANGE_VEC:
176 0 : return std::make_unique<FE<1,LAGRANGE_VEC>>(fet);
177 :
178 44040 : case L2_LAGRANGE:
179 44040 : return std::make_unique<FE<1,L2_LAGRANGE>>(fet);
180 :
181 0 : case L2_LAGRANGE_VEC:
182 0 : return std::make_unique<FE<1,L2_LAGRANGE_VEC>>(fet);
183 :
184 0 : case HIERARCHIC_VEC:
185 0 : return std::make_unique<FE<1,HIERARCHIC_VEC>>(fet);
186 :
187 214876 : case HIERARCHIC:
188 214876 : return std::make_unique<FE<1,HIERARCHIC>>(fet);
189 :
190 73400 : case L2_HIERARCHIC:
191 73400 : return std::make_unique<FE<1,L2_HIERARCHIC>>(fet);
192 :
193 0 : case L2_HIERARCHIC_VEC:
194 0 : return std::make_unique<FE<1,L2_HIERARCHIC_VEC>>(fet);
195 :
196 48444 : case SIDE_HIERARCHIC:
197 48444 : return std::make_unique<FE<1,SIDE_HIERARCHIC>>(fet);
198 :
199 105068 : case MONOMIAL:
200 203936 : return std::make_unique<FE<1,MONOMIAL>>(fet);
201 :
202 0 : case MONOMIAL_VEC:
203 0 : return std::make_unique<FE<1,MONOMIAL_VEC>>(fet);
204 :
205 : #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
206 58720 : case SZABAB:
207 58720 : return std::make_unique<FE<1,SZABAB>>(fet);
208 :
209 58720 : case BERNSTEIN:
210 58720 : return std::make_unique<FE<1,BERNSTEIN>>(fet);
211 :
212 29360 : case RATIONAL_BERNSTEIN:
213 29360 : return std::make_unique<FE<1,RATIONAL_BERNSTEIN>>(fet);
214 : #endif
215 :
216 58720 : case XYZ:
217 58720 : return std::make_unique<FEXYZ<1>>(fet);
218 :
219 0 : case SCALAR:
220 0 : return std::make_unique<FEScalar<1>>(fet);
221 :
222 0 : case NEDELEC_ONE:
223 0 : return std::make_unique<FENedelecOne<1>>(fet);
224 :
225 0 : case RAVIART_THOMAS:
226 0 : return std::make_unique<FERaviartThomas<1>>(fet);
227 :
228 0 : case L2_RAVIART_THOMAS:
229 0 : return std::make_unique<FEL2RaviartThomas<1>>(fet);
230 :
231 0 : case SUBDIVISION:
232 0 : return std::make_unique<FESubdivision>(fet);
233 :
234 0 : default:
235 0 : libmesh_error_msg("ERROR: Bad FEType.family= " << Utility::enum_to_string(fet.family));
236 : }
237 : }
238 :
239 :
240 : // 2D
241 5785963 : case 2:
242 : {
243 5785963 : switch (fet.family)
244 : {
245 126204 : case CLOUGH:
246 126204 : return std::make_unique<FE<2,CLOUGH>>(fet);
247 :
248 112288 : case HERMITE:
249 112288 : return std::make_unique<FE<2,HERMITE>>(fet);
250 :
251 2910007 : case LAGRANGE:
252 2910007 : return std::make_unique<FE<2,LAGRANGE>>(fet);
253 :
254 122688 : case LAGRANGE_VEC:
255 122688 : return std::make_unique<FE<2,LAGRANGE_VEC>>(fet);
256 :
257 156976 : case L2_LAGRANGE:
258 156976 : return std::make_unique<FE<2,L2_LAGRANGE>>(fet);
259 :
260 8472 : case L2_LAGRANGE_VEC:
261 8472 : return std::make_unique<FE<2,L2_LAGRANGE_VEC>>(fet);
262 :
263 0 : case HIERARCHIC_VEC:
264 0 : return std::make_unique<FE<2,HIERARCHIC_VEC>>(fet);
265 :
266 897852 : case HIERARCHIC:
267 897852 : return std::make_unique<FE<2,HIERARCHIC>>(fet);
268 :
269 209016 : case L2_HIERARCHIC:
270 209016 : return std::make_unique<FE<2,L2_HIERARCHIC>>(fet);
271 :
272 6816 : case L2_HIERARCHIC_VEC:
273 6816 : return std::make_unique<FE<2,L2_HIERARCHIC_VEC>>(fet);
274 :
275 173648 : case SIDE_HIERARCHIC:
276 173648 : return std::make_unique<FE<2,SIDE_HIERARCHIC>>(fet);
277 :
278 396592 : case MONOMIAL:
279 396592 : return std::make_unique<FE<2,MONOMIAL>>(fet);
280 :
281 3408 : case MONOMIAL_VEC:
282 3408 : return std::make_unique<FE<2,MONOMIAL_VEC>>(fet);
283 :
284 : #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
285 164416 : case SZABAB:
286 164416 : return std::make_unique<FE<2,SZABAB>>(fet);
287 :
288 161480 : case BERNSTEIN:
289 161480 : return std::make_unique<FE<2,BERNSTEIN>>(fet);
290 :
291 63040 : case RATIONAL_BERNSTEIN:
292 63040 : return std::make_unique<FE<2,RATIONAL_BERNSTEIN>>(fet);
293 : #endif
294 :
295 244160 : case XYZ:
296 244160 : return std::make_unique<FEXYZ<2>>(fet);
297 :
298 6792 : case SCALAR:
299 6792 : return std::make_unique<FEScalar<2>>(fet);
300 :
301 7428 : case NEDELEC_ONE:
302 7428 : return std::make_unique<FENedelecOne<2>>(fet);
303 :
304 11280 : case RAVIART_THOMAS:
305 11280 : return std::make_unique<FERaviartThomas<2>>(fet);
306 :
307 3400 : case L2_RAVIART_THOMAS:
308 3400 : return std::make_unique<FEL2RaviartThomas<2>>(fet);
309 :
310 0 : case SUBDIVISION:
311 0 : return std::make_unique<FESubdivision>(fet);
312 :
313 0 : default:
314 0 : libmesh_error_msg("ERROR: Bad FEType.family= " << Utility::enum_to_string(fet.family));
315 : }
316 : }
317 :
318 :
319 : // 3D
320 4087706 : case 3:
321 : {
322 4087706 : switch (fet.family)
323 : {
324 0 : case CLOUGH:
325 0 : libmesh_error_msg("ERROR: Clough-Tocher elements currently only support 1D and 2D");
326 :
327 22020 : case HERMITE:
328 22020 : return std::make_unique<FE<3,HERMITE>>(fet);
329 :
330 954346 : case LAGRANGE:
331 954346 : return std::make_unique<FE<3,LAGRANGE>>(fet);
332 :
333 185736 : case LAGRANGE_VEC:
334 185736 : return std::make_unique<FE<3,LAGRANGE_VEC>>(fet);
335 :
336 321202 : case L2_LAGRANGE:
337 321202 : return std::make_unique<FE<3,L2_LAGRANGE>>(fet);
338 :
339 5112 : case L2_LAGRANGE_VEC:
340 5112 : return std::make_unique<FE<3,L2_LAGRANGE_VEC>>(fet);
341 :
342 0 : case HIERARCHIC_VEC:
343 0 : return std::make_unique<FE<3,HIERARCHIC_VEC>>(fet);
344 :
345 422628 : case HIERARCHIC:
346 422628 : return std::make_unique<FE<3,HIERARCHIC>>(fet);
347 :
348 436698 : case L2_HIERARCHIC:
349 436698 : return std::make_unique<FE<3,L2_HIERARCHIC>>(fet);
350 :
351 2556 : case L2_HIERARCHIC_VEC:
352 2556 : return std::make_unique<FE<3,L2_HIERARCHIC_VEC>>(fet);
353 :
354 242220 : case SIDE_HIERARCHIC:
355 242220 : return std::make_unique<FE<3,SIDE_HIERARCHIC>>(fet);
356 :
357 699774 : case MONOMIAL:
358 699774 : return std::make_unique<FE<3,MONOMIAL>>(fet);
359 :
360 0 : case MONOMIAL_VEC:
361 0 : return std::make_unique<FE<3,MONOMIAL_VEC>>(fet);
362 :
363 : #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
364 0 : case SZABAB:
365 0 : return std::make_unique<FE<3,SZABAB>>(fet);
366 :
367 154140 : case BERNSTEIN:
368 154140 : return std::make_unique<FE<3,BERNSTEIN>>(fet);
369 :
370 70596 : case RATIONAL_BERNSTEIN:
371 70596 : return std::make_unique<FE<3,RATIONAL_BERNSTEIN>>(fet);
372 : #endif
373 :
374 555612 : case XYZ:
375 555612 : return std::make_unique<FEXYZ<3>>(fet);
376 :
377 1704 : case SCALAR:
378 1704 : return std::make_unique<FEScalar<3>>(fet);
379 :
380 6612 : case NEDELEC_ONE:
381 6612 : return std::make_unique<FENedelecOne<3>>(fet);
382 :
383 4200 : case RAVIART_THOMAS:
384 4200 : return std::make_unique<FERaviartThomas<3>>(fet);
385 :
386 2550 : case L2_RAVIART_THOMAS:
387 2550 : return std::make_unique<FEL2RaviartThomas<3>>(fet);
388 :
389 0 : default:
390 0 : libmesh_error_msg("ERROR: Bad FEType.family= " << Utility::enum_to_string(fet.family));
391 : }
392 : }
393 :
394 0 : default:
395 0 : libmesh_error_msg("Invalid dimension dim = " << dim);
396 : }
397 : }
398 :
399 :
400 :
401 34116366 : void FEAbstract::get_refspace_nodes(const ElemType itemType, std::vector<Point> & nodes)
402 : {
403 34116366 : const unsigned int n_nodes = Elem::type_to_n_nodes_map[itemType];
404 34116366 : if (n_nodes == invalid_uint)
405 0 : libmesh_error_msg("Number of nodes is not well-defined for " <<
406 : Utility::enum_to_string(itemType));
407 :
408 34116366 : nodes.resize(n_nodes);
409 :
410 : // For any element with a well-defined n_nodes we should also have a
411 : // well-defined reference element.
412 34116366 : const Elem & ref_elem = ReferenceElem::get(itemType);
413 :
414 317539688 : for (auto i : make_range(n_nodes))
415 332262538 : nodes[i] = ref_elem.point(i);
416 34116366 : }
417 :
418 :
419 :
420 : #ifdef LIBMESH_ENABLE_DEPRECATED
421 0 : bool FEAbstract::on_reference_element(const Point & p, const ElemType t, const Real eps)
422 : {
423 : // Use Elem::on_reference_element() instead
424 : libmesh_deprecated();
425 :
426 0 : libmesh_assert_greater_equal (eps, 0.);
427 :
428 0 : const Real xi = p(0);
429 : #if LIBMESH_DIM > 1
430 0 : const Real eta = p(1);
431 : #else
432 : const Real eta = 0.;
433 : #endif
434 : #if LIBMESH_DIM > 2
435 0 : const Real zeta = p(2);
436 : #else
437 : const Real zeta = 0.;
438 : #endif
439 :
440 0 : switch (t)
441 : {
442 0 : case NODEELEM:
443 : {
444 0 : return (!xi && !eta && !zeta);
445 : }
446 0 : case EDGE2:
447 : case EDGE3:
448 : case EDGE4:
449 : {
450 : // The reference 1D element is [-1,1].
451 0 : if ((xi >= -1.-eps) &&
452 0 : (xi <= 1.+eps))
453 0 : return true;
454 :
455 0 : return false;
456 : }
457 :
458 :
459 0 : case TRI3:
460 : case TRISHELL3:
461 : case TRI6:
462 : case TRI7:
463 : {
464 : // The reference triangle is isosceles
465 : // and is bound by xi=0, eta=0, and xi+eta=1.
466 0 : if ((xi >= 0.-eps) &&
467 0 : (eta >= 0.-eps) &&
468 0 : ((xi + eta) <= 1.+eps))
469 0 : return true;
470 :
471 0 : return false;
472 : }
473 :
474 :
475 0 : case QUAD4:
476 : case QUADSHELL4:
477 : case QUAD8:
478 : case QUADSHELL8:
479 : case QUAD9:
480 : case QUADSHELL9:
481 : {
482 : // The reference quadrilateral element is [-1,1]^2.
483 0 : if ((xi >= -1.-eps) &&
484 0 : (xi <= 1.+eps) &&
485 0 : (eta >= -1.-eps) &&
486 0 : (eta <= 1.+eps))
487 0 : return true;
488 :
489 0 : return false;
490 : }
491 :
492 :
493 0 : case TET4:
494 : case TET10:
495 : case TET14:
496 : {
497 : // The reference tetrahedral is isosceles
498 : // and is bound by xi=0, eta=0, zeta=0,
499 : // and xi+eta+zeta=1.
500 0 : if ((xi >= 0.-eps) &&
501 0 : (eta >= 0.-eps) &&
502 0 : (zeta >= 0.-eps) &&
503 0 : ((xi + eta + zeta) <= 1.+eps))
504 0 : return true;
505 :
506 0 : return false;
507 : }
508 :
509 :
510 0 : case HEX8:
511 : case HEX20:
512 : case HEX27:
513 : {
514 : /*
515 : if ((xi >= -1.) &&
516 : (xi <= 1.) &&
517 : (eta >= -1.) &&
518 : (eta <= 1.) &&
519 : (zeta >= -1.) &&
520 : (zeta <= 1.))
521 : return true;
522 : */
523 :
524 : // The reference hexahedral element is [-1,1]^3.
525 0 : if ((xi >= -1.-eps) &&
526 0 : (xi <= 1.+eps) &&
527 0 : (eta >= -1.-eps) &&
528 0 : (eta <= 1.+eps) &&
529 0 : (zeta >= -1.-eps) &&
530 0 : (zeta <= 1.+eps))
531 : {
532 : // libMesh::out << "Strange Point:\n";
533 : // p.print();
534 0 : return true;
535 : }
536 :
537 0 : return false;
538 : }
539 :
540 0 : case PRISM6:
541 : case PRISM15:
542 : case PRISM18:
543 : case PRISM20:
544 : case PRISM21:
545 : {
546 : // Figure this one out...
547 : // inside the reference triangle with zeta in [-1,1]
548 0 : if ((xi >= 0.-eps) &&
549 0 : (eta >= 0.-eps) &&
550 0 : (zeta >= -1.-eps) &&
551 0 : (zeta <= 1.+eps) &&
552 0 : ((xi + eta) <= 1.+eps))
553 0 : return true;
554 :
555 0 : return false;
556 : }
557 :
558 :
559 0 : case PYRAMID5:
560 : case PYRAMID13:
561 : case PYRAMID14:
562 : case PYRAMID18:
563 : {
564 : // Check that the point is on the same side of all the faces
565 : // by testing whether:
566 : //
567 : // n_i.(x - x_i) <= 0
568 : //
569 : // for each i, where:
570 : // n_i is the outward normal of face i,
571 : // x_i is a point on face i.
572 0 : if ((-eta - 1. + zeta <= 0.+eps) &&
573 0 : ( xi - 1. + zeta <= 0.+eps) &&
574 0 : ( eta - 1. + zeta <= 0.+eps) &&
575 0 : ( -xi - 1. + zeta <= 0.+eps) &&
576 0 : ( zeta >= 0.-eps))
577 0 : return true;
578 :
579 0 : return false;
580 : }
581 :
582 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
583 0 : case INFHEX8:
584 : case INFHEX16:
585 : case INFHEX18:
586 : {
587 : // The reference infhex8 is a [-1,1]^3.
588 0 : if ((xi >= -1.-eps) &&
589 0 : (xi <= 1.+eps) &&
590 0 : (eta >= -1.-eps) &&
591 0 : (eta <= 1.+eps) &&
592 0 : (zeta >= -1.-eps) &&
593 0 : (zeta <= 1.+eps))
594 : {
595 0 : return true;
596 : }
597 0 : return false;
598 : }
599 :
600 0 : case INFPRISM6:
601 : case INFPRISM12:
602 : {
603 : // inside the reference triangle with zeta in [-1,1]
604 0 : if ((xi >= 0.-eps) &&
605 0 : (eta >= 0.-eps) &&
606 0 : (zeta >= -1.-eps) &&
607 0 : (zeta <= 1.+eps) &&
608 0 : ((xi + eta) <= 1.+eps))
609 : {
610 0 : return true;
611 : }
612 :
613 0 : return false;
614 : }
615 : #endif
616 :
617 0 : default:
618 0 : libmesh_error_msg("ERROR: Unknown element type " << Utility::enum_to_string(t));
619 : }
620 :
621 : // If we get here then the point is _not_ in the
622 : // reference element. Better return false.
623 :
624 : return false;
625 : }
626 : #endif // LIBMESH_ENABLE_DEPRECATED
627 :
628 :
629 :
630 0 : void FEAbstract::print_JxW(std::ostream & os) const
631 : {
632 0 : this->_fe_map->print_JxW(os);
633 0 : }
634 :
635 :
636 :
637 0 : void FEAbstract::print_xyz(std::ostream & os) const
638 : {
639 0 : this->_fe_map->print_xyz(os);
640 0 : }
641 :
642 :
643 0 : void FEAbstract::print_info(std::ostream & os) const
644 : {
645 0 : os << "phi[i][j]: Shape function i at quadrature pt. j" << std::endl;
646 0 : this->print_phi(os);
647 :
648 0 : os << "dphi[i][j]: Shape function i's gradient at quadrature pt. j" << std::endl;
649 0 : this->print_dphi(os);
650 :
651 0 : os << "XYZ locations of the quadrature pts." << std::endl;
652 0 : this->print_xyz(os);
653 :
654 0 : os << "Values of JxW at the quadrature pts." << std::endl;
655 0 : this->print_JxW(os);
656 0 : }
657 :
658 :
659 0 : std::ostream & operator << (std::ostream & os, const FEAbstract & fe)
660 : {
661 0 : fe.print_info(os);
662 0 : return os;
663 : }
664 :
665 :
666 :
667 : #ifdef LIBMESH_ENABLE_AMR
668 :
669 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
670 1029068 : void FEAbstract::compute_node_constraints (NodeConstraints & constraints,
671 : const Elem * elem)
672 : {
673 514534 : libmesh_assert(elem);
674 :
675 1029068 : const unsigned int Dim = elem->dim();
676 :
677 : // Only constrain elements in 2,3D.
678 1029068 : if (Dim == 1)
679 96308 : return;
680 :
681 : // Only constrain active and ancestor elements
682 1014736 : if (elem->subactive())
683 39932 : return;
684 :
685 :
686 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
687 934872 : if (elem->infinite())
688 : {
689 2112 : const FEType fe_t(elem->default_order(), FEMap::map_fe_type(*elem));
690 :
691 : // expand the infinite_compute_constraint in its template-arguments.
692 2112 : switch(Dim)
693 : {
694 0 : case 2:
695 : {
696 0 : inf_fe_family_mapping_switch(2, inf_compute_node_constraints (constraints, elem) , ,; break;);
697 0 : break;
698 : }
699 1056 : case 3:
700 : {
701 2112 : inf_fe_family_mapping_switch(3, inf_compute_node_constraints (constraints, elem) , ,; break;);
702 1056 : break;
703 : }
704 0 : default:
705 0 : libmesh_error_msg("Invalid dim = " << Dim);
706 : }
707 1056 : return;
708 : }
709 :
710 : #endif
711 932760 : const FEFamily mapping_family = FEMap::map_fe_type(*elem);
712 932760 : const FEType fe_type(elem->default_side_order(), mapping_family);
713 :
714 : // Pull objects out of the loop to reduce heap operations
715 932760 : std::vector<const Node *> my_nodes, parent_nodes;
716 932760 : std::unique_ptr<const Elem> my_side, parent_side;
717 :
718 : // Look at the element faces. Check to see if we need to
719 : // build constraints.
720 4637580 : for (auto s : elem->side_index_range())
721 7293068 : if (elem->neighbor_ptr(s) != nullptr &&
722 3471676 : elem->neighbor_ptr(s) != remote_elem)
723 3471676 : if (elem->neighbor_ptr(s)->level() < elem->level()) // constrain dofs shared between
724 : { // this element and ones coarser
725 : // than this element.
726 : // Get pointers to the elements of interest and its parent.
727 117000 : const Elem * parent = elem->parent();
728 :
729 : // This can't happen... Only level-0 elements have nullptr
730 : // parents, and no level-0 elements can be at a higher
731 : // level than their neighbors!
732 58500 : libmesh_assert(parent);
733 :
734 117000 : elem->build_side_ptr(my_side, s);
735 117000 : parent->build_side_ptr(parent_side, s);
736 :
737 117000 : const unsigned int n_side_nodes = my_side->n_nodes();
738 :
739 58500 : my_nodes.clear();
740 117000 : my_nodes.reserve (n_side_nodes);
741 58500 : parent_nodes.clear();
742 117000 : parent_nodes.reserve (n_side_nodes);
743 :
744 437380 : for (unsigned int n=0; n != n_side_nodes; ++n)
745 480570 : my_nodes.push_back(my_side->node_ptr(n));
746 :
747 437380 : for (unsigned int n=0; n != n_side_nodes; ++n)
748 480570 : parent_nodes.push_back(parent_side->node_ptr(n));
749 :
750 378880 : for (unsigned int my_side_n=0;
751 437380 : my_side_n < n_side_nodes;
752 : my_side_n++)
753 : {
754 : // We can have an FE type that supports an order
755 : // partially, such that sides do not support the same
756 : // order. E.g. we say that a LAGRANGE PRISM21 supports
757 : // "third" order to distinguish its shape functions from
758 : // a PRISM18, but the QUAD9 sides will still only
759 : // support second order.
760 320380 : FEType side_fe_type = fe_type;
761 : const int side_max_order =
762 320380 : FEInterface::max_order(fe_type, my_side->type());
763 :
764 320380 : if ((int)fe_type.order > side_max_order)
765 2736 : side_fe_type.order = side_max_order;
766 :
767 : // Do not use the p_level(), if any, that is inherited by the side.
768 160190 : libmesh_assert_less
769 : (my_side_n,
770 : FEInterface::n_dofs(side_fe_type, /*extra_order=*/0,
771 : my_side.get()));
772 :
773 320380 : const Node * my_node = my_nodes[my_side_n];
774 :
775 : // The support point of the DOF
776 320380 : const Point & support_point = *my_node;
777 :
778 : // Figure out where my node lies on their reference element.
779 : const Point mapped_point = FEMap::inverse_map(Dim-1,
780 : parent_side.get(),
781 320380 : support_point);
782 :
783 : // Compute the parent's side shape function values.
784 1152266 : for (unsigned int their_side_n=0;
785 1312456 : their_side_n < n_side_nodes;
786 : their_side_n++)
787 : {
788 : // Do not use the p_level(), if any, that is inherited by the side.
789 496038 : libmesh_assert_less
790 : (their_side_n,
791 : FEInterface::n_dofs(side_fe_type,
792 : /*extra_order=*/0,
793 : parent_side.get()));
794 :
795 1488114 : const Node * their_node = parent_nodes[their_side_n];
796 496038 : libmesh_assert(their_node);
797 :
798 : // Do not use the p_level(), if any, that is inherited by the side.
799 992076 : const Real their_value = FEInterface::shape(side_fe_type,
800 : /*extra_order=*/0,
801 : parent_side.get(),
802 : their_side_n,
803 992076 : mapped_point);
804 :
805 496038 : const Real their_mag = std::abs(their_value);
806 : #ifdef DEBUG
807 : // Protect for the case u_i ~= u_j,
808 : // in which case i better equal j.
809 496038 : if (their_mag > 0.999)
810 : {
811 81714 : libmesh_assert_equal_to (my_node, their_node);
812 81714 : libmesh_assert_less (std::abs(their_value - 1.), 0.001);
813 : }
814 : else
815 : #endif
816 : // To make nodal constraints useful for constructing
817 : // sparsity patterns faster, we need to get EVERY
818 : // POSSIBLE constraint coupling identified, even if
819 : // there is no coupling in the isoparametric
820 : // Lagrange case.
821 910362 : if (their_mag < 1.e-5)
822 : {
823 : // since we may be running this method concurrently
824 : // on multiple threads we need to acquire a lock
825 : // before modifying the shared constraint_row object.
826 210986 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
827 :
828 : // A reference to the constraint row.
829 421972 : NodeConstraintRow & constraint_row = constraints[my_node].first;
830 :
831 421972 : constraint_row.emplace(their_node, 0.);
832 : }
833 : // To get nodal coordinate constraints right, only
834 : // add non-zero and non-identity values for Lagrange
835 : // basis functions.
836 : else // (1.e-5 <= their_mag <= .999)
837 : {
838 : // since we may be running this method concurrently
839 : // on multiple threads we need to acquire a lock
840 : // before modifying the shared constraint_row object.
841 406676 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
842 :
843 : // A reference to the constraint row.
844 488390 : NodeConstraintRow & constraint_row = constraints[my_node].first;
845 :
846 203338 : constraint_row.emplace(their_node, their_value);
847 : }
848 : }
849 : }
850 : }
851 : }
852 :
853 : #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
854 :
855 : #endif // #ifdef LIBMESH_ENABLE_AMR
856 :
857 :
858 :
859 : #ifdef LIBMESH_ENABLE_PERIODIC
860 :
861 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
862 129724 : void FEAbstract::compute_periodic_node_constraints (NodeConstraints & constraints,
863 : const PeriodicBoundaries & boundaries,
864 : const MeshBase & mesh,
865 : const PointLocatorBase * point_locator,
866 : const Elem * elem)
867 : {
868 : // Only bother if we truly have periodic boundaries
869 129724 : if (boundaries.empty())
870 38982 : return;
871 :
872 64862 : libmesh_assert(elem);
873 :
874 : // Only constrain active elements with this method
875 64862 : if (!elem->active())
876 19491 : return;
877 :
878 90742 : const unsigned int Dim = elem->dim();
879 :
880 90742 : const FEFamily mapping_family = FEMap::map_fe_type(*elem);
881 90742 : const FEType fe_type(elem->default_side_order(), mapping_family);
882 :
883 : // Pull objects out of the loop to reduce heap operations
884 90742 : std::vector<const Node *> my_nodes, neigh_nodes;
885 90742 : std::unique_ptr<const Elem> my_side, neigh_side;
886 :
887 : // Look at the element faces. Check to see if we need to
888 : // build constraints.
889 90742 : std::vector<boundary_id_type> bc_ids;
890 441726 : for (auto s : elem->side_index_range())
891 : {
892 526476 : if (elem->neighbor_ptr(s))
893 172162 : continue;
894 :
895 6660 : mesh.get_boundary_info().boundary_ids (elem, s, bc_ids);
896 12272 : for (const auto & boundary_id : bc_ids)
897 : {
898 5612 : const PeriodicBoundaryBase * periodic = boundaries.boundary(boundary_id);
899 5612 : if (periodic)
900 : {
901 2644 : libmesh_assert(point_locator);
902 :
903 : // Get pointers to the element's neighbor.
904 : unsigned int s_neigh;
905 5288 : const Elem * neigh = boundaries.neighbor(boundary_id, *point_locator, elem, s, &s_neigh);
906 :
907 5288 : libmesh_error_msg_if
908 : (!neigh, "PeriodicBoundaries can't find a periodic neighbor for element " <<
909 : elem->id() << " side " << s);
910 :
911 : // h refinement constraints:
912 : // constrain dofs shared between
913 : // this element and ones as coarse
914 : // as or coarser than this element.
915 5288 : if (neigh->level() <= elem->level())
916 : {
917 : #ifdef LIBMESH_ENABLE_AMR
918 2196 : libmesh_assert(neigh->active());
919 : #endif // #ifdef LIBMESH_ENABLE_AMR
920 :
921 4392 : elem->build_side_ptr(my_side, s);
922 4392 : neigh->build_side_ptr(neigh_side, s_neigh);
923 :
924 4392 : const unsigned int n_side_nodes = my_side->n_nodes();
925 :
926 2196 : my_nodes.clear();
927 4392 : my_nodes.reserve (n_side_nodes);
928 2196 : neigh_nodes.clear();
929 4392 : neigh_nodes.reserve (n_side_nodes);
930 :
931 14000 : for (unsigned int n=0; n != n_side_nodes; ++n)
932 14412 : my_nodes.push_back(my_side->node_ptr(n));
933 :
934 14000 : for (unsigned int n=0; n != n_side_nodes; ++n)
935 14412 : neigh_nodes.push_back(neigh_side->node_ptr(n));
936 :
937 : // Make sure we're not adding recursive constraints
938 : // due to the redundancy in the way we add periodic
939 : // boundary constraints, or adding constraints to
940 : // nodes that already have AMR constraints
941 6588 : std::vector<bool> skip_constraint(n_side_nodes, false);
942 :
943 11804 : for (unsigned int my_side_n=0;
944 14000 : my_side_n < n_side_nodes;
945 : my_side_n++)
946 : {
947 : // Do not use the p_level(), if any, that is inherited by the side.
948 4804 : libmesh_assert_less (my_side_n, FEInterface::n_dofs(fe_type, /*extra_order=*/0, my_side.get()));
949 :
950 14412 : const Node * my_node = my_nodes[my_side_n];
951 :
952 : // If we've already got a constraint on this
953 : // node, then the periodic constraint is
954 : // redundant
955 : {
956 4804 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
957 :
958 4804 : if (constraints.count(my_node))
959 : {
960 3896 : skip_constraint[my_side_n] = true;
961 1948 : continue;
962 : }
963 : }
964 :
965 : // Compute the neighbors's side shape function values.
966 16496 : for (unsigned int their_side_n=0;
967 19352 : their_side_n < n_side_nodes;
968 : their_side_n++)
969 : {
970 : // Do not use the p_level(), if any, that is inherited by the side.
971 6820 : libmesh_assert_less (their_side_n, FEInterface::n_dofs(fe_type, /*extra_order=*/0, neigh_side.get()));
972 :
973 20460 : const Node * their_node = neigh_nodes[their_side_n];
974 :
975 : // If there's a constraint on an opposing node,
976 : // we need to see if it's constrained by
977 : // *our side* making any periodic constraint
978 : // on us recursive
979 : {
980 6820 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
981 :
982 6820 : if (!constraints.count(their_node))
983 4152 : continue;
984 :
985 : const NodeConstraintRow & their_constraint_row =
986 5336 : constraints[their_node].first;
987 :
988 19684 : for (unsigned int orig_side_n=0;
989 22352 : orig_side_n < n_side_nodes;
990 : orig_side_n++)
991 : {
992 : // Do not use the p_level(), if any, that is inherited by the side.
993 8508 : libmesh_assert_less (orig_side_n, FEInterface::n_dofs(fe_type, /*extra_order=*/0, my_side.get()));
994 :
995 25524 : const Node * orig_node = my_nodes[orig_side_n];
996 :
997 8508 : if (their_constraint_row.count(orig_node))
998 10732 : skip_constraint[orig_side_n] = true;
999 : }
1000 : }
1001 : }
1002 : }
1003 11804 : for (unsigned int my_side_n=0;
1004 14000 : my_side_n < n_side_nodes;
1005 : my_side_n++)
1006 : {
1007 : // Do not use the p_level(), if any, that is inherited by the side.
1008 4804 : libmesh_assert_less (my_side_n, FEInterface::n_dofs(fe_type, /*extra_order=*/0, my_side.get()));
1009 :
1010 14412 : if (skip_constraint[my_side_n])
1011 5596 : continue;
1012 :
1013 4012 : const Node * my_node = my_nodes[my_side_n];
1014 :
1015 : // Figure out where my node lies on their reference element.
1016 4012 : const Point neigh_point = periodic->get_corresponding_pos(*my_node);
1017 :
1018 : // Figure out where my node lies on their reference element.
1019 : const Point mapped_point =
1020 : FEMap::inverse_map(Dim-1, neigh_side.get(),
1021 4012 : neigh_point);
1022 :
1023 10590 : for (unsigned int their_side_n=0;
1024 12596 : their_side_n < n_side_nodes;
1025 : their_side_n++)
1026 : {
1027 : // Do not use the p_level(), if any, that is inherited by the side.
1028 4292 : libmesh_assert_less (their_side_n, FEInterface::n_dofs(fe_type, /*extra_order=*/0, neigh_side.get()));
1029 :
1030 12876 : const Node * their_node = neigh_nodes[their_side_n];
1031 4292 : libmesh_assert(their_node);
1032 :
1033 : // Do not use the p_level(), if any, that is inherited by the side.
1034 8584 : const Real their_value = FEInterface::shape(fe_type,
1035 : /*extra_order=*/0,
1036 : neigh_side.get(),
1037 : their_side_n,
1038 8584 : mapped_point);
1039 :
1040 : // since we may be running this method concurrently
1041 : // on multiple threads we need to acquire a lock
1042 : // before modifying the shared constraint_row object.
1043 : {
1044 8584 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1045 :
1046 : NodeConstraintRow & constraint_row =
1047 8584 : constraints[my_node].first;
1048 :
1049 4292 : constraint_row.emplace(their_node, their_value);
1050 : }
1051 : }
1052 : }
1053 : }
1054 : }
1055 : }
1056 : }
1057 : }
1058 : #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
1059 :
1060 : #endif // LIBMESH_ENABLE_PERIODIC
1061 :
1062 :
1063 2156439 : unsigned int FEAbstract::n_quadrature_points () const
1064 : {
1065 240304 : if (this->shapes_on_quadrature)
1066 : {
1067 240288 : libmesh_assert(this->qrule);
1068 240288 : libmesh_assert_equal_to(this->qrule->n_points(),
1069 : this->_n_total_qp);
1070 : }
1071 2156439 : return this->_n_total_qp;
1072 : }
1073 :
1074 : } // namespace libMesh
|
