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/remote_elem.h"
35 : #include "libmesh/tensor_value.h"
36 : #include "libmesh/threads.h"
37 : #include "libmesh/enum_to_string.h"
38 :
39 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
40 : #include "libmesh/inf_fe.h"
41 : #include "libmesh/fe_interface_macros.h"
42 : #endif
43 :
44 : namespace libMesh
45 : {
46 :
47 12861431 : FEAbstract::FEAbstract(const unsigned int d,
48 12861431 : const FEType & fet) :
49 11418504 : _fe_map( FEMap::build(fet) ),
50 11418504 : dim(d),
51 11418504 : calculations_started(false),
52 11418504 : calculate_dual(false),
53 11418504 : calculate_default_dual_coeff(true),
54 11418504 : calculate_nothing(false),
55 11418504 : calculate_map(false),
56 11418504 : calculate_phi(false),
57 11418504 : calculate_dphi(false),
58 11418504 : calculate_d2phi(false),
59 11418504 : calculate_curl_phi(false),
60 11418504 : calculate_div_phi(false),
61 11418504 : calculate_dphiref(false),
62 11418504 : fe_type(fet),
63 11418504 : _elem_type(INVALID_ELEM),
64 11418504 : _elem(nullptr),
65 11418504 : _elem_p_level(0),
66 11418504 : _p_level(0),
67 11418504 : qrule(nullptr),
68 11418504 : shapes_on_quadrature(false),
69 11418504 : _n_total_qp(0),
70 12861431 : _add_p_level_in_reinit(fet.p_refinement)
71 : {
72 12861431 : }
73 :
74 :
75 11418504 : FEAbstract::~FEAbstract() = default;
76 :
77 :
78 6731810 : std::unique_ptr<FEAbstract> FEAbstract::build(const unsigned int dim,
79 : const FEType & fet)
80 : {
81 6731810 : switch (dim)
82 : {
83 : // 0D
84 62220 : case 0:
85 : {
86 62220 : switch (fet.family)
87 : {
88 0 : case CLOUGH:
89 0 : return std::make_unique<FE<0,CLOUGH>>(fet);
90 :
91 0 : case HERMITE:
92 0 : return std::make_unique<FE<0,HERMITE>>(fet);
93 :
94 37340 : case LAGRANGE:
95 37340 : return std::make_unique<FE<0,LAGRANGE>>(fet);
96 :
97 0 : case LAGRANGE_VEC:
98 0 : return std::make_unique<FE<0,LAGRANGE_VEC>>(fet);
99 :
100 840 : case L2_LAGRANGE:
101 840 : return std::make_unique<FE<0,L2_LAGRANGE>>(fet);
102 :
103 0 : case L2_LAGRANGE_GLL:
104 0 : return std::make_unique<FE<0,L2_LAGRANGE_GLL>>(fet);
105 :
106 840 : case L2_LAGRANGE_VEC:
107 840 : return std::make_unique<FE<0,L2_LAGRANGE_VEC>>(fet);
108 :
109 0 : case HIERARCHIC_VEC:
110 0 : return std::make_unique<FE<0,HIERARCHIC_VEC>>(fet);
111 :
112 0 : case HIERARCHIC:
113 0 : return std::make_unique<FE<0,HIERARCHIC>>(fet);
114 :
115 1136 : case L2_HIERARCHIC:
116 1136 : return std::make_unique<FE<0,L2_HIERARCHIC>>(fet);
117 :
118 0 : case L2_HIERARCHIC_VEC:
119 0 : return std::make_unique<FE<0,L2_HIERARCHIC_VEC>>(fet);
120 :
121 840 : case SIDE_HIERARCHIC:
122 840 : return std::make_unique<FE<0,SIDE_HIERARCHIC>>(fet);
123 :
124 3408 : case MONOMIAL:
125 3408 : return std::make_unique<FE<0,MONOMIAL>>(fet);
126 :
127 0 : case MONOMIAL_VEC:
128 0 : return std::make_unique<FE<0,MONOMIAL_VEC>>(fet);
129 :
130 : #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
131 0 : case SZABAB:
132 0 : return std::make_unique<FE<0,SZABAB>>(fet);
133 :
134 0 : case BERNSTEIN:
135 0 : return std::make_unique<FE<0,BERNSTEIN>>(fet);
136 :
137 11012 : case RATIONAL_BERNSTEIN:
138 11012 : return std::make_unique<FE<0,RATIONAL_BERNSTEIN>>(fet);
139 : #endif
140 :
141 0 : case XYZ:
142 0 : return std::make_unique<FEXYZ<0>>(fet);
143 :
144 3964 : case SCALAR:
145 3964 : return std::make_unique<FEScalar<0>>(fet);
146 :
147 0 : case NEDELEC_ONE:
148 0 : return std::make_unique<FENedelecOne<0>>(fet);
149 :
150 2840 : case RAVIART_THOMAS:
151 2840 : return std::make_unique<FERaviartThomas<0>>(fet);
152 :
153 0 : case L2_RAVIART_THOMAS:
154 0 : return std::make_unique<FEL2RaviartThomas<0>>(fet);
155 :
156 0 : case SUBDIVISION:
157 0 : return std::make_unique<FESubdivision>(fet);
158 :
159 0 : default:
160 0 : libmesh_error_msg("ERROR: Bad FEType.family= " << Utility::enum_to_string(fet.family));
161 : }
162 : }
163 : // 1D
164 738346 : case 1:
165 : {
166 738346 : switch (fet.family)
167 : {
168 0 : case CLOUGH:
169 0 : return std::make_unique<FE<1,CLOUGH>>(fet);
170 :
171 141252 : case HERMITE:
172 141252 : return std::make_unique<FE<1,HERMITE>>(fet);
173 :
174 68654 : case LAGRANGE:
175 68654 : return std::make_unique<FE<1,LAGRANGE>>(fet);
176 :
177 0 : case LAGRANGE_VEC:
178 0 : return std::make_unique<FE<1,LAGRANGE_VEC>>(fet);
179 :
180 28626 : case L2_LAGRANGE:
181 28626 : return std::make_unique<FE<1,L2_LAGRANGE>>(fet);
182 :
183 114504 : case L2_LAGRANGE_GLL:
184 114504 : return std::make_unique<FE<1,L2_LAGRANGE_GLL>>(fet);
185 :
186 0 : case L2_LAGRANGE_VEC:
187 0 : return std::make_unique<FE<1,L2_LAGRANGE_VEC>>(fet);
188 :
189 0 : case HIERARCHIC_VEC:
190 0 : return std::make_unique<FE<1,HIERARCHIC_VEC>>(fet);
191 :
192 116246 : case HIERARCHIC:
193 116246 : return std::make_unique<FE<1,HIERARCHIC>>(fet);
194 :
195 47710 : case L2_HIERARCHIC:
196 47710 : return std::make_unique<FE<1,L2_HIERARCHIC>>(fet);
197 :
198 0 : case L2_HIERARCHIC_VEC:
199 0 : return std::make_unique<FE<1,L2_HIERARCHIC_VEC>>(fet);
200 :
201 24222 : case SIDE_HIERARCHIC:
202 24222 : return std::make_unique<FE<1,SIDE_HIERARCHIC>>(fet);
203 :
204 63544 : case MONOMIAL:
205 63544 : return std::make_unique<FE<1,MONOMIAL>>(fet);
206 :
207 0 : case MONOMIAL_VEC:
208 0 : return std::make_unique<FE<1,MONOMIAL_VEC>>(fet);
209 :
210 : #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
211 38168 : case SZABAB:
212 38168 : return std::make_unique<FE<1,SZABAB>>(fet);
213 :
214 38168 : case BERNSTEIN:
215 38168 : return std::make_unique<FE<1,BERNSTEIN>>(fet);
216 :
217 19084 : case RATIONAL_BERNSTEIN:
218 19084 : return std::make_unique<FE<1,RATIONAL_BERNSTEIN>>(fet);
219 : #endif
220 :
221 38168 : case XYZ:
222 38168 : return std::make_unique<FEXYZ<1>>(fet);
223 :
224 0 : case SCALAR:
225 0 : return std::make_unique<FEScalar<1>>(fet);
226 :
227 0 : case NEDELEC_ONE:
228 0 : return std::make_unique<FENedelecOne<1>>(fet);
229 :
230 0 : case RAVIART_THOMAS:
231 0 : return std::make_unique<FERaviartThomas<1>>(fet);
232 :
233 0 : case L2_RAVIART_THOMAS:
234 0 : return std::make_unique<FEL2RaviartThomas<1>>(fet);
235 :
236 0 : case SUBDIVISION:
237 0 : return std::make_unique<FESubdivision>(fet);
238 :
239 0 : default:
240 0 : libmesh_error_msg("ERROR: Bad FEType.family= " << Utility::enum_to_string(fet.family));
241 : }
242 : }
243 :
244 :
245 : // 2D
246 3257939 : case 2:
247 : {
248 3257939 : switch (fet.family)
249 : {
250 67506 : case CLOUGH:
251 67506 : return std::make_unique<FE<2,CLOUGH>>(fet);
252 :
253 60548 : case HERMITE:
254 60548 : return std::make_unique<FE<2,HERMITE>>(fet);
255 :
256 1488785 : case LAGRANGE:
257 1488785 : return std::make_unique<FE<2,LAGRANGE>>(fet);
258 :
259 61344 : case LAGRANGE_VEC:
260 61344 : return std::make_unique<FE<2,LAGRANGE_VEC>>(fet);
261 :
262 100508 : case L2_LAGRANGE:
263 100508 : return std::make_unique<FE<2,L2_LAGRANGE>>(fet);
264 :
265 114504 : case L2_LAGRANGE_GLL:
266 114504 : return std::make_unique<FE<2,L2_LAGRANGE_GLL>>(fet);
267 :
268 4236 : case L2_LAGRANGE_VEC:
269 4236 : return std::make_unique<FE<2,L2_LAGRANGE_VEC>>(fet);
270 :
271 0 : case HIERARCHIC_VEC:
272 0 : return std::make_unique<FE<2,HIERARCHIC_VEC>>(fet);
273 :
274 473148 : case HIERARCHIC:
275 473148 : return std::make_unique<FE<2,HIERARCHIC>>(fet);
276 :
277 133134 : case L2_HIERARCHIC:
278 133134 : return std::make_unique<FE<2,L2_HIERARCHIC>>(fet);
279 :
280 3408 : case L2_HIERARCHIC_VEC:
281 3408 : return std::make_unique<FE<2,L2_HIERARCHIC_VEC>>(fet);
282 :
283 86824 : case SIDE_HIERARCHIC:
284 86824 : return std::make_unique<FE<2,SIDE_HIERARCHIC>>(fet);
285 :
286 241010 : case MONOMIAL:
287 241010 : return std::make_unique<FE<2,MONOMIAL>>(fet);
288 :
289 1704 : case MONOMIAL_VEC:
290 1704 : return std::make_unique<FE<2,MONOMIAL_VEC>>(fet);
291 :
292 : #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
293 106430 : case SZABAB:
294 106430 : return std::make_unique<FE<2,SZABAB>>(fet);
295 :
296 104962 : case BERNSTEIN:
297 104962 : return std::make_unique<FE<2,BERNSTEIN>>(fet);
298 :
299 38126 : case RATIONAL_BERNSTEIN:
300 38126 : return std::make_unique<FE<2,RATIONAL_BERNSTEIN>>(fet);
301 : #endif
302 :
303 157312 : case XYZ:
304 157312 : return std::make_unique<FEXYZ<2>>(fet);
305 :
306 3396 : case SCALAR:
307 3396 : return std::make_unique<FEScalar<2>>(fet);
308 :
309 3714 : case NEDELEC_ONE:
310 3714 : return std::make_unique<FENedelecOne<2>>(fet);
311 :
312 5640 : case RAVIART_THOMAS:
313 5640 : return std::make_unique<FERaviartThomas<2>>(fet);
314 :
315 1700 : case L2_RAVIART_THOMAS:
316 1700 : return std::make_unique<FEL2RaviartThomas<2>>(fet);
317 :
318 0 : case SUBDIVISION:
319 0 : return std::make_unique<FESubdivision>(fet);
320 :
321 0 : default:
322 0 : libmesh_error_msg("ERROR: Bad FEType.family= " << Utility::enum_to_string(fet.family));
323 : }
324 : }
325 :
326 :
327 : // 3D
328 2673305 : case 3:
329 : {
330 2673305 : switch (fet.family)
331 : {
332 0 : case CLOUGH:
333 0 : libmesh_error_msg("ERROR: Clough-Tocher elements currently only support 1D and 2D");
334 :
335 14313 : case HERMITE:
336 14313 : return std::make_unique<FE<3,HERMITE>>(fet);
337 :
338 554953 : case LAGRANGE:
339 554953 : return std::make_unique<FE<3,LAGRANGE>>(fet);
340 :
341 92868 : case LAGRANGE_VEC:
342 92868 : return std::make_unique<FE<3,LAGRANGE_VEC>>(fet);
343 :
344 206914 : case L2_LAGRANGE:
345 206914 : return std::make_unique<FE<3,L2_LAGRANGE>>(fet);
346 :
347 171756 : case L2_LAGRANGE_GLL:
348 171756 : return std::make_unique<FE<3,L2_LAGRANGE_GLL>>(fet);
349 :
350 2556 : case L2_LAGRANGE_VEC:
351 2556 : return std::make_unique<FE<3,L2_LAGRANGE_VEC>>(fet);
352 :
353 0 : case HIERARCHIC_VEC:
354 0 : return std::make_unique<FE<3,HIERARCHIC_VEC>>(fet);
355 :
356 267465 : case HIERARCHIC:
357 267465 : return std::make_unique<FE<3,HIERARCHIC>>(fet);
358 :
359 281106 : case L2_HIERARCHIC:
360 281106 : return std::make_unique<FE<3,L2_HIERARCHIC>>(fet);
361 :
362 1278 : case L2_HIERARCHIC_VEC:
363 1278 : return std::make_unique<FE<3,L2_HIERARCHIC_VEC>>(fet);
364 :
365 121110 : case SIDE_HIERARCHIC:
366 121110 : return std::make_unique<FE<3,SIDE_HIERARCHIC>>(fet);
367 :
368 445674 : case MONOMIAL:
369 445674 : return std::make_unique<FE<3,MONOMIAL>>(fet);
370 :
371 0 : case MONOMIAL_VEC:
372 0 : return std::make_unique<FE<3,MONOMIAL_VEC>>(fet);
373 :
374 : #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
375 0 : case SZABAB:
376 0 : return std::make_unique<FE<3,SZABAB>>(fet);
377 :
378 100191 : case BERNSTEIN:
379 100191 : return std::make_unique<FE<3,BERNSTEIN>>(fet);
380 :
381 45207 : case RATIONAL_BERNSTEIN:
382 45207 : return std::make_unique<FE<3,RATIONAL_BERNSTEIN>>(fet);
383 : #endif
384 :
385 360381 : case XYZ:
386 360381 : return std::make_unique<FEXYZ<3>>(fet);
387 :
388 852 : case SCALAR:
389 852 : return std::make_unique<FEScalar<3>>(fet);
390 :
391 3306 : case NEDELEC_ONE:
392 3306 : return std::make_unique<FENedelecOne<3>>(fet);
393 :
394 2100 : case RAVIART_THOMAS:
395 2100 : return std::make_unique<FERaviartThomas<3>>(fet);
396 :
397 1275 : case L2_RAVIART_THOMAS:
398 1275 : return std::make_unique<FEL2RaviartThomas<3>>(fet);
399 :
400 0 : default:
401 0 : libmesh_error_msg("ERROR: Bad FEType.family= " << Utility::enum_to_string(fet.family));
402 : }
403 : }
404 :
405 0 : default:
406 0 : libmesh_error_msg("Invalid dimension dim = " << dim);
407 : }
408 : }
409 :
410 :
411 :
412 34153238 : void FEAbstract::get_refspace_nodes(const ElemType itemType, std::vector<Point> & nodes)
413 : {
414 34153238 : const unsigned int n_nodes = Elem::type_to_n_nodes_map[itemType];
415 34153238 : if (n_nodes == invalid_uint)
416 0 : libmesh_error_msg("Number of nodes is not well-defined for " <<
417 : Utility::enum_to_string(itemType));
418 :
419 34153238 : nodes.resize(n_nodes);
420 34153238 : switch(itemType)
421 : {
422 28512 : case NODEELEM:
423 : {
424 114048 : nodes[0] = Point (0.,0.,0.);
425 114048 : return;
426 : }
427 6673 : case EDGE3:
428 : {
429 80374 : nodes[2] = Point (0.,0.,0.);
430 : libmesh_fallthrough();
431 : }
432 6996760 : case EDGE2:
433 : {
434 6996760 : nodes[0] = Point (-1.,0.,0.);
435 6996760 : nodes[1] = Point (1.,0.,0.);
436 6996760 : return;
437 : }
438 4 : case EDGE4: // not nested with EDGE3
439 : {
440 48 : nodes[0] = Point (-1.,0.,0.);
441 48 : nodes[1] = Point (1.,0.,0.);
442 48 : nodes[2] = Point (-1./3.,0.,0.);
443 48 : nodes[3] = Point (1./3.,0.,0.);
444 48 : return;
445 : }
446 124229 : case TRI7:
447 : {
448 1505023 : nodes[6] = Point (1./3.,1./3.,0.);
449 : libmesh_fallthrough();
450 : }
451 5930802 : case TRI6:
452 : {
453 5930802 : nodes[3] = Point (.5,0.,0.);
454 5930802 : nodes[4] = Point (.5,.5,0.);
455 5930802 : nodes[5] = Point (0.,.5,0.);
456 : libmesh_fallthrough();
457 : }
458 8267834 : case TRI3:
459 : case TRISHELL3:
460 : {
461 8267834 : nodes[0] = Point (0.,0.,0.);
462 8267834 : nodes[1] = Point (1.,0.,0.);
463 8267834 : nodes[2] = Point (0.,1.,0.);
464 8267834 : return;
465 : }
466 194738 : case QUAD9:
467 : case QUADSHELL9:
468 : {
469 2085466 : nodes[8] = Point (0.,0.,0.);
470 : libmesh_fallthrough();
471 : }
472 2210244 : case QUAD8:
473 : case QUADSHELL8:
474 : {
475 2210244 : nodes[4] = Point (0.,-1.,0.);
476 2210244 : nodes[5] = Point (1.,0.,0.);
477 2210244 : nodes[6] = Point (0.,1.,0.);
478 2210244 : nodes[7] = Point (-1.,0.,0.);
479 : libmesh_fallthrough();
480 : }
481 4887580 : case QUAD4:
482 : case QUADSHELL4:
483 : {
484 4887580 : nodes[0] = Point (-1.,-1.,0.);
485 4887580 : nodes[1] = Point (1.,-1.,0.);
486 4887580 : nodes[2] = Point (1.,1.,0.);
487 4887580 : nodes[3] = Point (-1.,1.,0.);
488 4887580 : return;
489 : }
490 997830 : case TET14:
491 : {
492 12069846 : nodes[10] = Point (1/Real(3),1/Real(3),0.);
493 12069846 : nodes[11] = Point (1/Real(3),0.,1/Real(3));
494 12069846 : nodes[12] = Point (1/Real(3),1/Real(3),1/Real(3));
495 12069846 : nodes[13] = Point (0.,1/Real(3),1/Real(3));
496 : libmesh_fallthrough();
497 : }
498 12131450 : case TET10:
499 : {
500 12131450 : nodes[4] = Point (.5,0.,0.);
501 12131450 : nodes[5] = Point (.5,.5,0.);
502 12131450 : nodes[6] = Point (0.,.5,0.);
503 12131450 : nodes[7] = Point (0.,0.,.5);
504 12131450 : nodes[8] = Point (.5,0.,.5);
505 12131450 : nodes[9] = Point (0.,.5,.5);
506 : libmesh_fallthrough();
507 : }
508 12133694 : case TET4:
509 : {
510 12133694 : nodes[0] = Point (0.,0.,0.);
511 12133694 : nodes[1] = Point (1.,0.,0.);
512 12133694 : nodes[2] = Point (0.,1.,0.);
513 12133694 : nodes[3] = Point (0.,0.,1.);
514 12133694 : return;
515 : }
516 35028 : case HEX27:
517 : {
518 481763 : nodes[20] = Point (0.,0.,-1.);
519 481763 : nodes[21] = Point (0.,-1.,0.);
520 481763 : nodes[22] = Point (1.,0.,0.);
521 481763 : nodes[23] = Point (0.,1.,0.);
522 481763 : nodes[24] = Point (-1.,0.,0.);
523 481763 : nodes[25] = Point (0.,0.,1.);
524 481763 : nodes[26] = Point (0.,0.,0.);
525 : libmesh_fallthrough();
526 : }
527 487517 : case HEX20:
528 : {
529 487517 : nodes[8] = Point (0.,-1.,-1.);
530 487517 : nodes[9] = Point (1.,0.,-1.);
531 487517 : nodes[10] = Point (0.,1.,-1.);
532 487517 : nodes[11] = Point (-1.,0.,-1.);
533 487517 : nodes[12] = Point (-1.,-1.,0.);
534 487517 : nodes[13] = Point (1.,-1.,0.);
535 487517 : nodes[14] = Point (1.,1.,0.);
536 487517 : nodes[15] = Point (-1.,1.,0.);
537 487517 : nodes[16] = Point (0.,-1.,1.);
538 487517 : nodes[17] = Point (1.,0.,1.);
539 487517 : nodes[18] = Point (0.,1.,1.);
540 487517 : nodes[19] = Point (-1.,0.,1.);
541 : libmesh_fallthrough();
542 : }
543 1527283 : case HEX8:
544 : {
545 1527283 : nodes[0] = Point (-1.,-1.,-1.);
546 1527283 : nodes[1] = Point (1.,-1.,-1.);
547 1527283 : nodes[2] = Point (1.,1.,-1.);
548 1527283 : nodes[3] = Point (-1.,1.,-1.);
549 1527283 : nodes[4] = Point (-1.,-1.,1.);
550 1527283 : nodes[5] = Point (1.,-1.,1.);
551 1527283 : nodes[6] = Point (1.,1.,1.);
552 1527283 : nodes[7] = Point (-1.,1.,1.);
553 1527283 : return;
554 : }
555 7976 : case PRISM21:
556 : {
557 97556 : nodes[20] = Point (1/Real(3),1/Real(3),0);
558 : libmesh_fallthrough();
559 : }
560 154895 : case PRISM20:
561 : {
562 154895 : nodes[18] = Point (1/Real(3),1/Real(3),-1);
563 154895 : nodes[19] = Point (1/Real(3),1/Real(3),1);
564 : libmesh_fallthrough();
565 : }
566 198903 : case PRISM18:
567 : {
568 198903 : nodes[15] = Point (.5,0.,0.);
569 198903 : nodes[16] = Point (.5,.5,0.);
570 198903 : nodes[17] = Point (0.,.5,0.);
571 : libmesh_fallthrough();
572 : }
573 199478 : case PRISM15:
574 : {
575 199478 : nodes[6] = Point (.5,0.,-1.);
576 199478 : nodes[7] = Point (.5,.5,-1.);
577 199478 : nodes[8] = Point (0.,.5,-1.);
578 199478 : nodes[9] = Point (0.,0.,0.);
579 199478 : nodes[10] = Point (1.,0.,0.);
580 199478 : nodes[11] = Point (0.,1.,0.);
581 199478 : nodes[12] = Point (.5,0.,1.);
582 199478 : nodes[13] = Point (.5,.5,1.);
583 199478 : nodes[14] = Point (0.,.5,1.);
584 : libmesh_fallthrough();
585 : }
586 224696 : case PRISM6:
587 : {
588 224696 : nodes[0] = Point (0.,0.,-1.);
589 224696 : nodes[1] = Point (1.,0.,-1.);
590 224696 : nodes[2] = Point (0.,1.,-1.);
591 224696 : nodes[3] = Point (0.,0.,1.);
592 224696 : nodes[4] = Point (1.,0.,1.);
593 224696 : nodes[5] = Point (0.,1.,1.);
594 224696 : return;
595 : }
596 14 : case PYRAMID18:
597 : {
598 : // triangle centers
599 447 : nodes[14] = Point (0., -2/Real(3),1/Real(3));
600 447 : nodes[15] = Point (2/Real(3),0.,1/Real(3));
601 447 : nodes[16] = Point (0.,2/Real(3),1/Real(3));
602 447 : nodes[17] = Point (-2/Real(3),0.,1/Real(3));
603 :
604 : libmesh_fallthrough();
605 : }
606 848 : case PYRAMID14:
607 : {
608 : // base center
609 848 : nodes[13] = Point (0.,0.,0.);
610 :
611 : libmesh_fallthrough();
612 : }
613 1249 : case PYRAMID13:
614 : {
615 : // base midedge
616 1249 : nodes[5] = Point (0.,-1.,0.);
617 1249 : nodes[6] = Point (1.,0.,0.);
618 1249 : nodes[7] = Point (0.,1.,0.);
619 1249 : nodes[8] = Point (-1,0.,0.);
620 :
621 : // lateral midedge
622 1249 : nodes[9] = Point (-.5,-.5,.5);
623 1249 : nodes[10] = Point (.5,-.5,.5);
624 1249 : nodes[11] = Point (.5,.5,.5);
625 1249 : nodes[12] = Point (-.5,.5,.5);
626 :
627 : libmesh_fallthrough();
628 : }
629 1295 : case PYRAMID5:
630 : {
631 : // base corners
632 1295 : nodes[0] = Point (-1.,-1.,0.);
633 1295 : nodes[1] = Point (1.,-1.,0.);
634 1295 : nodes[2] = Point (1.,1.,0.);
635 1295 : nodes[3] = Point (-1.,1.,0.);
636 : // apex
637 1295 : nodes[4] = Point (0.,0.,1.);
638 1295 : return;
639 : }
640 :
641 0 : default:
642 0 : libmesh_error_msg("ERROR: Unknown element type " << Utility::enum_to_string(itemType));
643 : }
644 : }
645 :
646 :
647 :
648 : #ifdef LIBMESH_ENABLE_DEPRECATED
649 0 : bool FEAbstract::on_reference_element(const Point & p, const ElemType t, const Real eps)
650 : {
651 : // Use Elem::on_reference_element() instead
652 : libmesh_deprecated();
653 :
654 0 : libmesh_assert_greater_equal (eps, 0.);
655 :
656 0 : const Real xi = p(0);
657 : #if LIBMESH_DIM > 1
658 0 : const Real eta = p(1);
659 : #else
660 : const Real eta = 0.;
661 : #endif
662 : #if LIBMESH_DIM > 2
663 0 : const Real zeta = p(2);
664 : #else
665 : const Real zeta = 0.;
666 : #endif
667 :
668 0 : switch (t)
669 : {
670 0 : case NODEELEM:
671 : {
672 0 : return (!xi && !eta && !zeta);
673 : }
674 0 : case EDGE2:
675 : case EDGE3:
676 : case EDGE4:
677 : {
678 : // The reference 1D element is [-1,1].
679 0 : if ((xi >= -1.-eps) &&
680 0 : (xi <= 1.+eps))
681 0 : return true;
682 :
683 0 : return false;
684 : }
685 :
686 :
687 0 : case TRI3:
688 : case TRISHELL3:
689 : case TRI6:
690 : case TRI7:
691 : {
692 : // The reference triangle is isosceles
693 : // and is bound by xi=0, eta=0, and xi+eta=1.
694 0 : if ((xi >= 0.-eps) &&
695 0 : (eta >= 0.-eps) &&
696 0 : ((xi + eta) <= 1.+eps))
697 0 : return true;
698 :
699 0 : return false;
700 : }
701 :
702 :
703 0 : case QUAD4:
704 : case QUADSHELL4:
705 : case QUAD8:
706 : case QUADSHELL8:
707 : case QUAD9:
708 : case QUADSHELL9:
709 : {
710 : // The reference quadrilateral element is [-1,1]^2.
711 0 : if ((xi >= -1.-eps) &&
712 0 : (xi <= 1.+eps) &&
713 0 : (eta >= -1.-eps) &&
714 0 : (eta <= 1.+eps))
715 0 : return true;
716 :
717 0 : return false;
718 : }
719 :
720 :
721 0 : case TET4:
722 : case TET10:
723 : case TET14:
724 : {
725 : // The reference tetrahedral is isosceles
726 : // and is bound by xi=0, eta=0, zeta=0,
727 : // and xi+eta+zeta=1.
728 0 : if ((xi >= 0.-eps) &&
729 0 : (eta >= 0.-eps) &&
730 0 : (zeta >= 0.-eps) &&
731 0 : ((xi + eta + zeta) <= 1.+eps))
732 0 : return true;
733 :
734 0 : return false;
735 : }
736 :
737 :
738 0 : case HEX8:
739 : case HEX20:
740 : case HEX27:
741 : {
742 : /*
743 : if ((xi >= -1.) &&
744 : (xi <= 1.) &&
745 : (eta >= -1.) &&
746 : (eta <= 1.) &&
747 : (zeta >= -1.) &&
748 : (zeta <= 1.))
749 : return true;
750 : */
751 :
752 : // The reference hexahedral element is [-1,1]^3.
753 0 : if ((xi >= -1.-eps) &&
754 0 : (xi <= 1.+eps) &&
755 0 : (eta >= -1.-eps) &&
756 0 : (eta <= 1.+eps) &&
757 0 : (zeta >= -1.-eps) &&
758 0 : (zeta <= 1.+eps))
759 : {
760 : // libMesh::out << "Strange Point:\n";
761 : // p.print();
762 0 : return true;
763 : }
764 :
765 0 : return false;
766 : }
767 :
768 0 : case PRISM6:
769 : case PRISM15:
770 : case PRISM18:
771 : case PRISM20:
772 : case PRISM21:
773 : {
774 : // Figure this one out...
775 : // inside the reference triangle with zeta in [-1,1]
776 0 : if ((xi >= 0.-eps) &&
777 0 : (eta >= 0.-eps) &&
778 0 : (zeta >= -1.-eps) &&
779 0 : (zeta <= 1.+eps) &&
780 0 : ((xi + eta) <= 1.+eps))
781 0 : return true;
782 :
783 0 : return false;
784 : }
785 :
786 :
787 0 : case PYRAMID5:
788 : case PYRAMID13:
789 : case PYRAMID14:
790 : case PYRAMID18:
791 : {
792 : // Check that the point is on the same side of all the faces
793 : // by testing whether:
794 : //
795 : // n_i.(x - x_i) <= 0
796 : //
797 : // for each i, where:
798 : // n_i is the outward normal of face i,
799 : // x_i is a point on face i.
800 0 : if ((-eta - 1. + zeta <= 0.+eps) &&
801 0 : ( xi - 1. + zeta <= 0.+eps) &&
802 0 : ( eta - 1. + zeta <= 0.+eps) &&
803 0 : ( -xi - 1. + zeta <= 0.+eps) &&
804 0 : ( zeta >= 0.-eps))
805 0 : return true;
806 :
807 0 : return false;
808 : }
809 :
810 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
811 0 : case INFHEX8:
812 : case INFHEX16:
813 : case INFHEX18:
814 : {
815 : // The reference infhex8 is a [-1,1]^3.
816 0 : if ((xi >= -1.-eps) &&
817 0 : (xi <= 1.+eps) &&
818 0 : (eta >= -1.-eps) &&
819 0 : (eta <= 1.+eps) &&
820 0 : (zeta >= -1.-eps) &&
821 0 : (zeta <= 1.+eps))
822 : {
823 0 : return true;
824 : }
825 0 : return false;
826 : }
827 :
828 0 : case INFPRISM6:
829 : case INFPRISM12:
830 : {
831 : // inside the reference triangle with zeta in [-1,1]
832 0 : if ((xi >= 0.-eps) &&
833 0 : (eta >= 0.-eps) &&
834 0 : (zeta >= -1.-eps) &&
835 0 : (zeta <= 1.+eps) &&
836 0 : ((xi + eta) <= 1.+eps))
837 : {
838 0 : return true;
839 : }
840 :
841 0 : return false;
842 : }
843 : #endif
844 :
845 0 : default:
846 0 : libmesh_error_msg("ERROR: Unknown element type " << Utility::enum_to_string(t));
847 : }
848 :
849 : // If we get here then the point is _not_ in the
850 : // reference element. Better return false.
851 :
852 : return false;
853 : }
854 : #endif // LIBMESH_ENABLE_DEPRECATED
855 :
856 :
857 :
858 0 : void FEAbstract::print_JxW(std::ostream & os) const
859 : {
860 0 : this->_fe_map->print_JxW(os);
861 0 : }
862 :
863 :
864 :
865 0 : void FEAbstract::print_xyz(std::ostream & os) const
866 : {
867 0 : this->_fe_map->print_xyz(os);
868 0 : }
869 :
870 :
871 0 : void FEAbstract::print_info(std::ostream & os) const
872 : {
873 0 : os << "phi[i][j]: Shape function i at quadrature pt. j" << std::endl;
874 0 : this->print_phi(os);
875 :
876 0 : os << "dphi[i][j]: Shape function i's gradient at quadrature pt. j" << std::endl;
877 0 : this->print_dphi(os);
878 :
879 0 : os << "XYZ locations of the quadrature pts." << std::endl;
880 0 : this->print_xyz(os);
881 :
882 0 : os << "Values of JxW at the quadrature pts." << std::endl;
883 0 : this->print_JxW(os);
884 0 : }
885 :
886 :
887 0 : std::ostream & operator << (std::ostream & os, const FEAbstract & fe)
888 : {
889 0 : fe.print_info(os);
890 0 : return os;
891 : }
892 :
893 :
894 :
895 : #ifdef LIBMESH_ENABLE_AMR
896 :
897 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
898 1038978 : void FEAbstract::compute_node_constraints (NodeConstraints & constraints,
899 : const Elem * elem)
900 : {
901 519489 : libmesh_assert(elem);
902 :
903 1038978 : const unsigned int Dim = elem->dim();
904 :
905 : // Only constrain elements in 2,3D.
906 1038978 : if (Dim == 1)
907 96308 : return;
908 :
909 : // Only constrain active and ancestor elements
910 1024646 : if (elem->subactive())
911 39932 : return;
912 :
913 :
914 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
915 944782 : if (elem->infinite())
916 : {
917 2112 : const FEType fe_t(elem->default_order(), FEMap::map_fe_type(*elem));
918 :
919 : // expand the infinite_compute_constraint in its template-arguments.
920 2112 : switch(Dim)
921 : {
922 0 : case 2:
923 : {
924 0 : inf_fe_family_mapping_switch(2, inf_compute_node_constraints (constraints, elem) , ,; break;);
925 0 : break;
926 : }
927 1056 : case 3:
928 : {
929 2112 : inf_fe_family_mapping_switch(3, inf_compute_node_constraints (constraints, elem) , ,; break;);
930 1056 : break;
931 : }
932 0 : default:
933 0 : libmesh_error_msg("Invalid dim = " << Dim);
934 : }
935 1056 : return;
936 : }
937 :
938 : #endif
939 942670 : const FEFamily mapping_family = FEMap::map_fe_type(*elem);
940 942670 : const FEType fe_type(elem->default_side_order(), mapping_family);
941 :
942 : // Pull objects out of the loop to reduce heap operations
943 942670 : std::vector<const Node *> my_nodes, parent_nodes;
944 942670 : std::unique_ptr<const Elem> my_side, parent_side;
945 :
946 : // Look at the element faces. Check to see if we need to
947 : // build constraints.
948 4689836 : for (auto s : elem->side_index_range())
949 7369643 : if (elem->neighbor_ptr(s) != nullptr &&
950 3497788 : elem->neighbor_ptr(s) != remote_elem)
951 3497788 : if (elem->neighbor_ptr(s)->level() < elem->level()) // constrain dofs shared between
952 : { // this element and ones coarser
953 : // than this element.
954 : // Get pointers to the elements of interest and its parent.
955 117000 : const Elem * parent = elem->parent();
956 :
957 : // This can't happen... Only level-0 elements have nullptr
958 : // parents, and no level-0 elements can be at a higher
959 : // level than their neighbors!
960 58500 : libmesh_assert(parent);
961 :
962 117000 : elem->build_side_ptr(my_side, s);
963 117000 : parent->build_side_ptr(parent_side, s);
964 :
965 117000 : const unsigned int n_side_nodes = my_side->n_nodes();
966 :
967 58500 : my_nodes.clear();
968 117000 : my_nodes.reserve (n_side_nodes);
969 58500 : parent_nodes.clear();
970 117000 : parent_nodes.reserve (n_side_nodes);
971 :
972 437380 : for (unsigned int n=0; n != n_side_nodes; ++n)
973 480570 : my_nodes.push_back(my_side->node_ptr(n));
974 :
975 437380 : for (unsigned int n=0; n != n_side_nodes; ++n)
976 480570 : parent_nodes.push_back(parent_side->node_ptr(n));
977 :
978 378880 : for (unsigned int my_side_n=0;
979 437380 : my_side_n < n_side_nodes;
980 : my_side_n++)
981 : {
982 : // We can have an FE type that supports an order
983 : // partially, such that sides do not support the same
984 : // order. E.g. we say that a LAGRANGE PRISM21 supports
985 : // "third" order to distinguish its shape functions from
986 : // a PRISM18, but the QUAD9 sides will still only
987 : // support second order.
988 320380 : FEType side_fe_type = fe_type;
989 : const int side_max_order =
990 320380 : FEInterface::max_order(fe_type, my_side->type());
991 :
992 320380 : if ((int)fe_type.order > side_max_order)
993 2736 : side_fe_type.order = side_max_order;
994 :
995 : // Do not use the p_level(), if any, that is inherited by the side.
996 160190 : libmesh_assert_less
997 : (my_side_n,
998 : FEInterface::n_dofs(side_fe_type, /*extra_order=*/0,
999 : my_side.get()));
1000 :
1001 320380 : const Node * my_node = my_nodes[my_side_n];
1002 :
1003 : // The support point of the DOF
1004 320380 : const Point & support_point = *my_node;
1005 :
1006 : // Figure out where my node lies on their reference element.
1007 : const Point mapped_point = FEMap::inverse_map(Dim-1,
1008 : parent_side.get(),
1009 320380 : support_point);
1010 :
1011 : // Compute the parent's side shape function values.
1012 1152266 : for (unsigned int their_side_n=0;
1013 1312456 : their_side_n < n_side_nodes;
1014 : their_side_n++)
1015 : {
1016 : // Do not use the p_level(), if any, that is inherited by the side.
1017 496038 : libmesh_assert_less
1018 : (their_side_n,
1019 : FEInterface::n_dofs(side_fe_type,
1020 : /*extra_order=*/0,
1021 : parent_side.get()));
1022 :
1023 1488114 : const Node * their_node = parent_nodes[their_side_n];
1024 496038 : libmesh_assert(their_node);
1025 :
1026 : // Do not use the p_level(), if any, that is inherited by the side.
1027 992076 : const Real their_value = FEInterface::shape(side_fe_type,
1028 : /*extra_order=*/0,
1029 : parent_side.get(),
1030 : their_side_n,
1031 992076 : mapped_point);
1032 :
1033 496038 : const Real their_mag = std::abs(their_value);
1034 : #ifdef DEBUG
1035 : // Protect for the case u_i ~= u_j,
1036 : // in which case i better equal j.
1037 496038 : if (their_mag > 0.999)
1038 : {
1039 81714 : libmesh_assert_equal_to (my_node, their_node);
1040 81714 : libmesh_assert_less (std::abs(their_value - 1.), 0.001);
1041 : }
1042 : else
1043 : #endif
1044 : // To make nodal constraints useful for constructing
1045 : // sparsity patterns faster, we need to get EVERY
1046 : // POSSIBLE constraint coupling identified, even if
1047 : // there is no coupling in the isoparametric
1048 : // Lagrange case.
1049 910362 : if (their_mag < 1.e-5)
1050 : {
1051 : // since we may be running this method concurrently
1052 : // on multiple threads we need to acquire a lock
1053 : // before modifying the shared constraint_row object.
1054 210986 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1055 :
1056 : // A reference to the constraint row.
1057 421972 : NodeConstraintRow & constraint_row = constraints[my_node].first;
1058 :
1059 421972 : constraint_row.emplace(their_node, 0.);
1060 : }
1061 : // To get nodal coordinate constraints right, only
1062 : // add non-zero and non-identity values for Lagrange
1063 : // basis functions.
1064 : else // (1.e-5 <= their_mag <= .999)
1065 : {
1066 : // since we may be running this method concurrently
1067 : // on multiple threads we need to acquire a lock
1068 : // before modifying the shared constraint_row object.
1069 406676 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1070 :
1071 : // A reference to the constraint row.
1072 488390 : NodeConstraintRow & constraint_row = constraints[my_node].first;
1073 :
1074 203338 : constraint_row.emplace(their_node, their_value);
1075 : }
1076 : }
1077 : }
1078 : }
1079 : }
1080 :
1081 : #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
1082 :
1083 : #endif // #ifdef LIBMESH_ENABLE_AMR
1084 :
1085 :
1086 :
1087 : #ifdef LIBMESH_ENABLE_PERIODIC
1088 :
1089 : #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
1090 129724 : void FEAbstract::compute_periodic_node_constraints (NodeConstraints & constraints,
1091 : const PeriodicBoundaries & boundaries,
1092 : const MeshBase & mesh,
1093 : const PointLocatorBase * point_locator,
1094 : const Elem * elem)
1095 : {
1096 : // Only bother if we truly have periodic boundaries
1097 129724 : if (boundaries.empty())
1098 38982 : return;
1099 :
1100 64862 : libmesh_assert(elem);
1101 :
1102 : // Only constrain active elements with this method
1103 64862 : if (!elem->active())
1104 19491 : return;
1105 :
1106 90742 : const unsigned int Dim = elem->dim();
1107 :
1108 90742 : const FEFamily mapping_family = FEMap::map_fe_type(*elem);
1109 90742 : const FEType fe_type(elem->default_side_order(), mapping_family);
1110 :
1111 : // Pull objects out of the loop to reduce heap operations
1112 90742 : std::vector<const Node *> my_nodes, neigh_nodes;
1113 90742 : std::unique_ptr<const Elem> my_side, neigh_side;
1114 :
1115 : // Look at the element faces. Check to see if we need to
1116 : // build constraints.
1117 90742 : std::vector<boundary_id_type> bc_ids;
1118 441726 : for (auto s : elem->side_index_range())
1119 : {
1120 526476 : if (elem->neighbor_ptr(s))
1121 172162 : continue;
1122 :
1123 6660 : mesh.get_boundary_info().boundary_ids (elem, s, bc_ids);
1124 12272 : for (const auto & boundary_id : bc_ids)
1125 : {
1126 5612 : const PeriodicBoundaryBase * periodic = boundaries.boundary(boundary_id);
1127 5612 : if (periodic)
1128 : {
1129 2644 : libmesh_assert(point_locator);
1130 :
1131 : // Get pointers to the element's neighbor.
1132 : unsigned int s_neigh;
1133 5288 : const Elem * neigh = boundaries.neighbor(boundary_id, *point_locator, elem, s, &s_neigh);
1134 :
1135 5288 : libmesh_error_msg_if
1136 : (!neigh, "PeriodicBoundaries can't find a periodic neighbor for element " <<
1137 : elem->id() << " side " << s);
1138 :
1139 : // h refinement constraints:
1140 : // constrain dofs shared between
1141 : // this element and ones as coarse
1142 : // as or coarser than this element.
1143 5288 : if (neigh->level() <= elem->level())
1144 : {
1145 : #ifdef LIBMESH_ENABLE_AMR
1146 2196 : libmesh_assert(neigh->active());
1147 : #endif // #ifdef LIBMESH_ENABLE_AMR
1148 :
1149 4392 : elem->build_side_ptr(my_side, s);
1150 4392 : neigh->build_side_ptr(neigh_side, s_neigh);
1151 :
1152 4392 : const unsigned int n_side_nodes = my_side->n_nodes();
1153 :
1154 2196 : my_nodes.clear();
1155 4392 : my_nodes.reserve (n_side_nodes);
1156 2196 : neigh_nodes.clear();
1157 4392 : neigh_nodes.reserve (n_side_nodes);
1158 :
1159 14000 : for (unsigned int n=0; n != n_side_nodes; ++n)
1160 14412 : my_nodes.push_back(my_side->node_ptr(n));
1161 :
1162 14000 : for (unsigned int n=0; n != n_side_nodes; ++n)
1163 14412 : neigh_nodes.push_back(neigh_side->node_ptr(n));
1164 :
1165 : // Make sure we're not adding recursive constraints
1166 : // due to the redundancy in the way we add periodic
1167 : // boundary constraints, or adding constraints to
1168 : // nodes that already have AMR constraints
1169 6588 : std::vector<bool> skip_constraint(n_side_nodes, false);
1170 :
1171 11804 : for (unsigned int my_side_n=0;
1172 14000 : my_side_n < n_side_nodes;
1173 : my_side_n++)
1174 : {
1175 : // Do not use the p_level(), if any, that is inherited by the side.
1176 4804 : libmesh_assert_less (my_side_n, FEInterface::n_dofs(fe_type, /*extra_order=*/0, my_side.get()));
1177 :
1178 14412 : const Node * my_node = my_nodes[my_side_n];
1179 :
1180 : // If we've already got a constraint on this
1181 : // node, then the periodic constraint is
1182 : // redundant
1183 : {
1184 4804 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1185 :
1186 4804 : if (constraints.count(my_node))
1187 : {
1188 3896 : skip_constraint[my_side_n] = true;
1189 1948 : continue;
1190 : }
1191 : }
1192 :
1193 : // Compute the neighbors's side shape function values.
1194 16496 : for (unsigned int their_side_n=0;
1195 19352 : their_side_n < n_side_nodes;
1196 : their_side_n++)
1197 : {
1198 : // Do not use the p_level(), if any, that is inherited by the side.
1199 6820 : libmesh_assert_less (their_side_n, FEInterface::n_dofs(fe_type, /*extra_order=*/0, neigh_side.get()));
1200 :
1201 20460 : const Node * their_node = neigh_nodes[their_side_n];
1202 :
1203 : // If there's a constraint on an opposing node,
1204 : // we need to see if it's constrained by
1205 : // *our side* making any periodic constraint
1206 : // on us recursive
1207 : {
1208 6820 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1209 :
1210 6820 : if (!constraints.count(their_node))
1211 4152 : continue;
1212 :
1213 : const NodeConstraintRow & their_constraint_row =
1214 5336 : constraints[their_node].first;
1215 :
1216 19684 : for (unsigned int orig_side_n=0;
1217 22352 : orig_side_n < n_side_nodes;
1218 : orig_side_n++)
1219 : {
1220 : // Do not use the p_level(), if any, that is inherited by the side.
1221 8508 : libmesh_assert_less (orig_side_n, FEInterface::n_dofs(fe_type, /*extra_order=*/0, my_side.get()));
1222 :
1223 25524 : const Node * orig_node = my_nodes[orig_side_n];
1224 :
1225 8508 : if (their_constraint_row.count(orig_node))
1226 10732 : skip_constraint[orig_side_n] = true;
1227 : }
1228 : }
1229 : }
1230 : }
1231 11804 : for (unsigned int my_side_n=0;
1232 14000 : my_side_n < n_side_nodes;
1233 : my_side_n++)
1234 : {
1235 : // Do not use the p_level(), if any, that is inherited by the side.
1236 4804 : libmesh_assert_less (my_side_n, FEInterface::n_dofs(fe_type, /*extra_order=*/0, my_side.get()));
1237 :
1238 14412 : if (skip_constraint[my_side_n])
1239 5596 : continue;
1240 :
1241 4012 : const Node * my_node = my_nodes[my_side_n];
1242 :
1243 : // Figure out where my node lies on their reference element.
1244 4012 : const Point neigh_point = periodic->get_corresponding_pos(*my_node);
1245 :
1246 : // Figure out where my node lies on their reference element.
1247 : const Point mapped_point =
1248 : FEMap::inverse_map(Dim-1, neigh_side.get(),
1249 4012 : neigh_point);
1250 :
1251 10590 : for (unsigned int their_side_n=0;
1252 12596 : their_side_n < n_side_nodes;
1253 : their_side_n++)
1254 : {
1255 : // Do not use the p_level(), if any, that is inherited by the side.
1256 4292 : libmesh_assert_less (their_side_n, FEInterface::n_dofs(fe_type, /*extra_order=*/0, neigh_side.get()));
1257 :
1258 12876 : const Node * their_node = neigh_nodes[their_side_n];
1259 4292 : libmesh_assert(their_node);
1260 :
1261 : // Do not use the p_level(), if any, that is inherited by the side.
1262 8584 : const Real their_value = FEInterface::shape(fe_type,
1263 : /*extra_order=*/0,
1264 : neigh_side.get(),
1265 : their_side_n,
1266 8584 : mapped_point);
1267 :
1268 : // since we may be running this method concurrently
1269 : // on multiple threads we need to acquire a lock
1270 : // before modifying the shared constraint_row object.
1271 : {
1272 8584 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1273 :
1274 : NodeConstraintRow & constraint_row =
1275 8584 : constraints[my_node].first;
1276 :
1277 4292 : constraint_row.emplace(their_node, their_value);
1278 : }
1279 : }
1280 : }
1281 : }
1282 : }
1283 : }
1284 : }
1285 : }
1286 : #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
1287 :
1288 : #endif // LIBMESH_ENABLE_PERIODIC
1289 :
1290 :
1291 2156421 : unsigned int FEAbstract::n_quadrature_points () const
1292 : {
1293 240304 : if (this->shapes_on_quadrature)
1294 : {
1295 240288 : libmesh_assert(this->qrule);
1296 240288 : libmesh_assert_equal_to(this->qrule->n_points(),
1297 : this->_n_total_qp);
1298 : }
1299 2156421 : return this->_n_total_qp;
1300 : }
1301 :
1302 : } // namespace libMesh
|
