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 2620388 : FEAbstract::FEAbstract(const unsigned int d,
48 2620388 : const FEType & fet) :
49 1175606 : _fe_map( FEMap::build(fet) ),
50 1175606 : dim(d),
51 1175606 : calculations_started(false),
52 1175606 : calculate_dual(false),
53 1175606 : calculate_default_dual_coeff(true),
54 1175606 : calculate_nothing(false),
55 1175606 : calculate_map(false),
56 1175606 : calculate_phi(false),
57 1175606 : calculate_dphi(false),
58 1175606 : calculate_d2phi(false),
59 1175606 : calculate_curl_phi(false),
60 1175606 : calculate_div_phi(false),
61 1175606 : calculate_dphiref(false),
62 1175606 : fe_type(fet),
63 1175606 : _elem_type(INVALID_ELEM),
64 1175606 : _elem(nullptr),
65 1175606 : _elem_p_level(0),
66 1175606 : _p_level(0),
67 1175606 : qrule(nullptr),
68 1175606 : shapes_on_quadrature(false),
69 1175606 : _n_total_qp(0),
70 2620388 : _add_p_level_in_reinit(fet.p_refinement)
71 : {
72 2620388 : }
73 :
74 :
75 1175606 : FEAbstract::~FEAbstract() = default;
76 :
77 :
78 725522 : std::unique_ptr<FEAbstract> FEAbstract::build(const unsigned int dim,
79 : const FEType & fet)
80 : {
81 725522 : switch (dim)
82 : {
83 : // 0D
84 6308 : case 0:
85 : {
86 6308 : 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 3676 : case LAGRANGE:
95 3676 : 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 72 : case L2_LAGRANGE:
101 72 : 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 72 : case L2_LAGRANGE_VEC:
107 72 : 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 112 : case L2_HIERARCHIC:
116 112 : 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 72 : case SIDE_HIERARCHIC:
122 72 : return std::make_unique<FE<0,SIDE_HIERARCHIC>>(fet);
123 :
124 336 : case MONOMIAL:
125 336 : 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 1308 : case RATIONAL_BERNSTEIN:
138 1308 : 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 380 : case SCALAR:
145 380 : return std::make_unique<FEScalar<0>>(fet);
146 :
147 0 : case NEDELEC_ONE:
148 0 : return std::make_unique<FENedelecOne<0>>(fet);
149 :
150 280 : case RAVIART_THOMAS:
151 280 : 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 77500 : case 1:
165 : {
166 77500 : switch (fet.family)
167 : {
168 0 : case CLOUGH:
169 0 : return std::make_unique<FE<1,CLOUGH>>(fet);
170 :
171 14052 : case HERMITE:
172 14052 : return std::make_unique<FE<1,HERMITE>>(fet);
173 :
174 6954 : case LAGRANGE:
175 6954 : 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 3042 : case L2_LAGRANGE:
181 3042 : return std::make_unique<FE<1,L2_LAGRANGE>>(fet);
182 :
183 12168 : case L2_LAGRANGE_GLL:
184 12168 : 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 11782 : case HIERARCHIC:
193 11782 : return std::make_unique<FE<1,HIERARCHIC>>(fet);
194 :
195 5070 : case L2_HIERARCHIC:
196 5070 : 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 2574 : case SIDE_HIERARCHIC:
202 2574 : return std::make_unique<FE<1,SIDE_HIERARCHIC>>(fet);
203 :
204 6648 : case MONOMIAL:
205 6648 : 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 4056 : case SZABAB:
212 4056 : return std::make_unique<FE<1,SZABAB>>(fet);
213 :
214 5070 : case BERNSTEIN:
215 5070 : return std::make_unique<FE<1,BERNSTEIN>>(fet);
216 :
217 2028 : case RATIONAL_BERNSTEIN:
218 2028 : return std::make_unique<FE<1,RATIONAL_BERNSTEIN>>(fet);
219 : #endif
220 :
221 4056 : case XYZ:
222 4056 : 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 365233 : case 2:
247 : {
248 365233 : switch (fet.family)
249 : {
250 6802 : case CLOUGH:
251 6802 : return std::make_unique<FE<2,CLOUGH>>(fet);
252 :
253 6116 : case HERMITE:
254 6116 : return std::make_unique<FE<2,HERMITE>>(fet);
255 :
256 185217 : case LAGRANGE:
257 185217 : return std::make_unique<FE<2,LAGRANGE>>(fet);
258 :
259 6048 : case LAGRANGE_VEC:
260 6048 : return std::make_unique<FE<2,LAGRANGE_VEC>>(fet);
261 :
262 10620 : case L2_LAGRANGE:
263 10620 : return std::make_unique<FE<2,L2_LAGRANGE>>(fet);
264 :
265 12168 : case L2_LAGRANGE_GLL:
266 12168 : return std::make_unique<FE<2,L2_LAGRANGE_GLL>>(fet);
267 :
268 396 : case L2_LAGRANGE_VEC:
269 396 : 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 44132 : case HIERARCHIC:
275 44132 : return std::make_unique<FE<2,HIERARCHIC>>(fet);
276 :
277 14078 : case L2_HIERARCHIC:
278 14078 : return std::make_unique<FE<2,L2_HIERARCHIC>>(fet);
279 :
280 336 : case L2_HIERARCHIC_VEC:
281 336 : return std::make_unique<FE<2,L2_HIERARCHIC_VEC>>(fet);
282 :
283 9192 : case SIDE_HIERARCHIC:
284 9192 : return std::make_unique<FE<2,SIDE_HIERARCHIC>>(fet);
285 :
286 25376 : case MONOMIAL:
287 25376 : return std::make_unique<FE<2,MONOMIAL>>(fet);
288 :
289 168 : case MONOMIAL_VEC:
290 168 : return std::make_unique<FE<2,MONOMIAL_VEC>>(fet);
291 :
292 : #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
293 11310 : case SZABAB:
294 11310 : return std::make_unique<FE<2,SZABAB>>(fet);
295 :
296 11154 : case BERNSTEIN:
297 11154 : return std::make_unique<FE<2,BERNSTEIN>>(fet);
298 :
299 4150 : case RATIONAL_BERNSTEIN:
300 4150 : return std::make_unique<FE<2,RATIONAL_BERNSTEIN>>(fet);
301 : #endif
302 :
303 16704 : case XYZ:
304 16704 : return std::make_unique<FEXYZ<2>>(fet);
305 :
306 324 : case SCALAR:
307 324 : return std::make_unique<FEScalar<2>>(fet);
308 :
309 258 : case NEDELEC_ONE:
310 258 : return std::make_unique<FENedelecOne<2>>(fet);
311 :
312 520 : case RAVIART_THOMAS:
313 520 : return std::make_unique<FERaviartThomas<2>>(fet);
314 :
315 164 : case L2_RAVIART_THOMAS:
316 164 : 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 276481 : case 3:
329 : {
330 276481 : 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 1521 : case HERMITE:
336 1521 : return std::make_unique<FE<3,HERMITE>>(fet);
337 :
338 52697 : case LAGRANGE:
339 52697 : return std::make_unique<FE<3,LAGRANGE>>(fet);
340 :
341 9156 : case LAGRANGE_VEC:
342 9156 : return std::make_unique<FE<3,LAGRANGE_VEC>>(fet);
343 :
344 21938 : case L2_LAGRANGE:
345 21938 : return std::make_unique<FE<3,L2_LAGRANGE>>(fet);
346 :
347 18252 : case L2_LAGRANGE_GLL:
348 18252 : return std::make_unique<FE<3,L2_LAGRANGE_GLL>>(fet);
349 :
350 252 : case L2_LAGRANGE_VEC:
351 252 : 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 28305 : case HIERARCHIC:
357 28305 : return std::make_unique<FE<3,HIERARCHIC>>(fet);
358 :
359 29802 : case L2_HIERARCHIC:
360 29802 : return std::make_unique<FE<3,L2_HIERARCHIC>>(fet);
361 :
362 126 : case L2_HIERARCHIC_VEC:
363 126 : return std::make_unique<FE<3,L2_HIERARCHIC_VEC>>(fet);
364 :
365 12870 : case SIDE_HIERARCHIC:
366 12870 : return std::make_unique<FE<3,SIDE_HIERARCHIC>>(fet);
367 :
368 47154 : case MONOMIAL:
369 47154 : 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 10647 : case BERNSTEIN:
379 10647 : return std::make_unique<FE<3,BERNSTEIN>>(fet);
380 :
381 4863 : case RATIONAL_BERNSTEIN:
382 4863 : return std::make_unique<FE<3,RATIONAL_BERNSTEIN>>(fet);
383 : #endif
384 :
385 38277 : case XYZ:
386 38277 : return std::make_unique<FEXYZ<3>>(fet);
387 :
388 84 : case SCALAR:
389 84 : return std::make_unique<FEScalar<3>>(fet);
390 :
391 234 : case NEDELEC_ONE:
392 234 : return std::make_unique<FENedelecOne<3>>(fet);
393 :
394 180 : case RAVIART_THOMAS:
395 180 : return std::make_unique<FERaviartThomas<3>>(fet);
396 :
397 123 : case L2_RAVIART_THOMAS:
398 123 : 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 11366086 : void FEAbstract::get_refspace_nodes(const ElemType itemType, std::vector<Point> & nodes)
413 : {
414 11366086 : const unsigned int n_nodes = Elem::type_to_n_nodes_map[itemType];
415 11366086 : 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 11366086 : nodes.resize(n_nodes);
420 11366086 : 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 26684 : nodes[2] = Point (0.,0.,0.);
430 : libmesh_fallthrough();
431 : }
432 1912982 : case EDGE2:
433 : {
434 1912982 : nodes[0] = Point (-1.,0.,0.);
435 1912982 : nodes[1] = Point (1.,0.,0.);
436 1912982 : return;
437 : }
438 5 : case EDGE4: // not nested with EDGE3
439 : {
440 20 : nodes[0] = Point (-1.,0.,0.);
441 20 : nodes[1] = Point (1.,0.,0.);
442 20 : nodes[2] = Point (-1./3.,0.,0.);
443 20 : nodes[3] = Point (1./3.,0.,0.);
444 20 : return;
445 : }
446 124229 : case TRI7:
447 : {
448 498168 : nodes[6] = Point (1./3.,1./3.,0.);
449 : libmesh_fallthrough();
450 : }
451 1913615 : case TRI6:
452 : {
453 1913615 : nodes[3] = Point (.5,0.,0.);
454 1913615 : nodes[4] = Point (.5,.5,0.);
455 1913615 : nodes[5] = Point (0.,.5,0.);
456 : libmesh_fallthrough();
457 : }
458 2692330 : case TRI3:
459 : case TRISHELL3:
460 : {
461 2692330 : nodes[0] = Point (0.,0.,0.);
462 2692330 : nodes[1] = Point (1.,0.,0.);
463 2692330 : nodes[2] = Point (0.,1.,0.);
464 2692330 : return;
465 : }
466 194754 : case QUAD9:
467 : case QUADSHELL9:
468 : {
469 689137 : nodes[8] = Point (0.,0.,0.);
470 : libmesh_fallthrough();
471 : }
472 727048 : case QUAD8:
473 : case QUADSHELL8:
474 : {
475 727048 : nodes[4] = Point (0.,-1.,0.);
476 727048 : nodes[5] = Point (1.,0.,0.);
477 727048 : nodes[6] = Point (0.,1.,0.);
478 727048 : nodes[7] = Point (-1.,0.,0.);
479 : libmesh_fallthrough();
480 : }
481 2062405 : case QUAD4:
482 : case QUADSHELL4:
483 : {
484 2062405 : nodes[0] = Point (-1.,-1.,0.);
485 2062405 : nodes[1] = Point (1.,-1.,0.);
486 2062405 : nodes[2] = Point (1.,1.,0.);
487 2062405 : nodes[3] = Point (-1.,1.,0.);
488 2062405 : return;
489 : }
490 997830 : case TET14:
491 : {
492 3998632 : nodes[10] = Point (1/Real(3),1/Real(3),0.);
493 3998632 : nodes[11] = Point (1/Real(3),0.,1/Real(3));
494 3998632 : nodes[12] = Point (1/Real(3),1/Real(3),1/Real(3));
495 3998632 : nodes[13] = Point (0.,1/Real(3),1/Real(3));
496 : libmesh_fallthrough();
497 : }
498 4016386 : case TET10:
499 : {
500 4016386 : nodes[4] = Point (.5,0.,0.);
501 4016386 : nodes[5] = Point (.5,.5,0.);
502 4016386 : nodes[6] = Point (0.,.5,0.);
503 4016386 : nodes[7] = Point (0.,0.,.5);
504 4016386 : nodes[8] = Point (.5,0.,.5);
505 4016386 : nodes[9] = Point (0.,.5,.5);
506 : libmesh_fallthrough();
507 : }
508 4016984 : case TET4:
509 : {
510 4016984 : nodes[0] = Point (0.,0.,0.);
511 4016984 : nodes[1] = Point (1.,0.,0.);
512 4016984 : nodes[2] = Point (0.,1.,0.);
513 4016984 : nodes[3] = Point (0.,0.,1.);
514 4016984 : return;
515 : }
516 35028 : case HEX27:
517 : {
518 143356 : nodes[20] = Point (0.,0.,-1.);
519 143356 : nodes[21] = Point (0.,-1.,0.);
520 143356 : nodes[22] = Point (1.,0.,0.);
521 143356 : nodes[23] = Point (0.,1.,0.);
522 143356 : nodes[24] = Point (-1.,0.,0.);
523 143356 : nodes[25] = Point (0.,0.,1.);
524 143356 : nodes[26] = Point (0.,0.,0.);
525 : libmesh_fallthrough();
526 : }
527 145046 : case HEX20:
528 : {
529 145046 : nodes[8] = Point (0.,-1.,-1.);
530 145046 : nodes[9] = Point (1.,0.,-1.);
531 145046 : nodes[10] = Point (0.,1.,-1.);
532 145046 : nodes[11] = Point (-1.,0.,-1.);
533 145046 : nodes[12] = Point (-1.,-1.,0.);
534 145046 : nodes[13] = Point (1.,-1.,0.);
535 145046 : nodes[14] = Point (1.,1.,0.);
536 145046 : nodes[15] = Point (-1.,1.,0.);
537 145046 : nodes[16] = Point (0.,-1.,1.);
538 145046 : nodes[17] = Point (1.,0.,1.);
539 145046 : nodes[18] = Point (0.,1.,1.);
540 145046 : nodes[19] = Point (-1.,0.,1.);
541 : libmesh_fallthrough();
542 : }
543 481920 : case HEX8:
544 : {
545 481920 : nodes[0] = Point (-1.,-1.,-1.);
546 481920 : nodes[1] = Point (1.,-1.,-1.);
547 481920 : nodes[2] = Point (1.,1.,-1.);
548 481920 : nodes[3] = Point (-1.,1.,-1.);
549 481920 : nodes[4] = Point (-1.,-1.,1.);
550 481920 : nodes[5] = Point (1.,-1.,1.);
551 481920 : nodes[6] = Point (1.,1.,1.);
552 481920 : nodes[7] = Point (-1.,1.,1.);
553 481920 : return;
554 : }
555 7976 : case PRISM21:
556 : {
557 31591 : nodes[20] = Point (1/Real(3),1/Real(3),0);
558 : libmesh_fallthrough();
559 : }
560 49769 : case PRISM20:
561 : {
562 49769 : nodes[18] = Point (1/Real(3),1/Real(3),-1);
563 49769 : nodes[19] = Point (1/Real(3),1/Real(3),1);
564 : libmesh_fallthrough();
565 : }
566 76749 : case PRISM18:
567 : {
568 76749 : nodes[15] = Point (.5,0.,0.);
569 76749 : nodes[16] = Point (.5,.5,0.);
570 76749 : nodes[17] = Point (0.,.5,0.);
571 : libmesh_fallthrough();
572 : }
573 76854 : case PRISM15:
574 : {
575 76854 : nodes[6] = Point (.5,0.,-1.);
576 76854 : nodes[7] = Point (.5,.5,-1.);
577 76854 : nodes[8] = Point (0.,.5,-1.);
578 76854 : nodes[9] = Point (0.,0.,0.);
579 76854 : nodes[10] = Point (1.,0.,0.);
580 76854 : nodes[11] = Point (0.,1.,0.);
581 76854 : nodes[12] = Point (.5,0.,1.);
582 76854 : nodes[13] = Point (.5,.5,1.);
583 76854 : nodes[14] = Point (0.,.5,1.);
584 : libmesh_fallthrough();
585 : }
586 85257 : case PRISM6:
587 : {
588 85257 : nodes[0] = Point (0.,0.,-1.);
589 85257 : nodes[1] = Point (1.,0.,-1.);
590 85257 : nodes[2] = Point (0.,1.,-1.);
591 85257 : nodes[3] = Point (0.,0.,1.);
592 85257 : nodes[4] = Point (1.,0.,1.);
593 85257 : nodes[5] = Point (0.,1.,1.);
594 85257 : return;
595 : }
596 14 : case PYRAMID18:
597 : {
598 : // triangle centers
599 49 : nodes[14] = Point (0., -2/Real(3),1/Real(3));
600 49 : nodes[15] = Point (2/Real(3),0.,1/Real(3));
601 49 : nodes[16] = Point (0.,2/Real(3),1/Real(3));
602 49 : nodes[17] = Point (-2/Real(3),0.,1/Real(3));
603 :
604 : libmesh_fallthrough();
605 : }
606 91 : case PYRAMID14:
607 : {
608 : // base center
609 91 : nodes[13] = Point (0.,0.,0.);
610 :
611 : libmesh_fallthrough();
612 : }
613 133 : case PYRAMID13:
614 : {
615 : // base midedge
616 133 : nodes[5] = Point (0.,-1.,0.);
617 133 : nodes[6] = Point (1.,0.,0.);
618 133 : nodes[7] = Point (0.,1.,0.);
619 133 : nodes[8] = Point (-1,0.,0.);
620 :
621 : // lateral midedge
622 133 : nodes[9] = Point (-.5,-.5,.5);
623 133 : nodes[10] = Point (.5,-.5,.5);
624 133 : nodes[11] = Point (.5,.5,.5);
625 133 : nodes[12] = Point (-.5,.5,.5);
626 :
627 : libmesh_fallthrough();
628 : }
629 140 : case PYRAMID5:
630 : {
631 : // base corners
632 140 : nodes[0] = Point (-1.,-1.,0.);
633 140 : nodes[1] = Point (1.,-1.,0.);
634 140 : nodes[2] = Point (1.,1.,0.);
635 140 : nodes[3] = Point (-1.,1.,0.);
636 : // apex
637 140 : nodes[4] = Point (0.,0.,1.);
638 140 : 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 1039738 : void FEAbstract::compute_node_constraints (NodeConstraints & constraints,
899 : const Elem * elem)
900 : {
901 519869 : libmesh_assert(elem);
902 :
903 1039738 : const unsigned int Dim = elem->dim();
904 :
905 : // Only constrain elements in 2,3D.
906 1039738 : if (Dim == 1)
907 96276 : return;
908 :
909 : // Only constrain active and ancestor elements
910 1025390 : if (elem->subactive())
911 39908 : return;
912 :
913 :
914 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
915 945574 : 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 943462 : const FEFamily mapping_family = FEMap::map_fe_type(*elem);
940 943462 : const FEType fe_type(elem->default_side_order(), mapping_family);
941 :
942 : // Pull objects out of the loop to reduce heap operations
943 943462 : std::vector<const Node *> my_nodes, parent_nodes;
944 943462 : 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 4693748 : for (auto s : elem->side_index_range())
949 7375533 : if (elem->neighbor_ptr(s) != nullptr &&
950 3500208 : elem->neighbor_ptr(s) != remote_elem)
951 3500208 : 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 117204 : 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 58602 : libmesh_assert(parent);
961 :
962 117204 : elem->build_side_ptr(my_side, s);
963 117204 : parent->build_side_ptr(parent_side, s);
964 :
965 117204 : const unsigned int n_side_nodes = my_side->n_nodes();
966 :
967 58602 : my_nodes.clear();
968 117204 : my_nodes.reserve (n_side_nodes);
969 58602 : parent_nodes.clear();
970 117204 : parent_nodes.reserve (n_side_nodes);
971 :
972 438196 : for (unsigned int n=0; n != n_side_nodes; ++n)
973 481488 : my_nodes.push_back(my_side->node_ptr(n));
974 :
975 438196 : for (unsigned int n=0; n != n_side_nodes; ++n)
976 481488 : parent_nodes.push_back(parent_side->node_ptr(n));
977 :
978 379594 : for (unsigned int my_side_n=0;
979 438196 : 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 320992 : FEType side_fe_type = fe_type;
989 : const int side_max_order =
990 320992 : FEInterface::max_order(fe_type, my_side->type());
991 :
992 320992 : 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 160496 : libmesh_assert_less
997 : (my_side_n,
998 : FEInterface::n_dofs(side_fe_type, /*extra_order=*/0,
999 : my_side.get()));
1000 :
1001 320992 : const Node * my_node = my_nodes[my_side_n];
1002 :
1003 : // The support point of the DOF
1004 320992 : 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 320992 : support_point);
1010 :
1011 : // Compute the parent's side shape function values.
1012 1154408 : for (unsigned int their_side_n=0;
1013 1314904 : 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 496956 : libmesh_assert_less
1018 : (their_side_n,
1019 : FEInterface::n_dofs(side_fe_type,
1020 : /*extra_order=*/0,
1021 : parent_side.get()));
1022 :
1023 1490868 : const Node * their_node = parent_nodes[their_side_n];
1024 496956 : libmesh_assert(their_node);
1025 :
1026 : // Do not use the p_level(), if any, that is inherited by the side.
1027 993912 : const Real their_value = FEInterface::shape(side_fe_type,
1028 : /*extra_order=*/0,
1029 : parent_side.get(),
1030 : their_side_n,
1031 993912 : mapped_point);
1032 :
1033 496956 : 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 496956 : if (their_mag > 0.999)
1038 : {
1039 81918 : libmesh_assert_equal_to (my_node, their_node);
1040 81918 : 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 911994 : 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 211394 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1055 :
1056 : // A reference to the constraint row.
1057 422788 : NodeConstraintRow & constraint_row = constraints[my_node].first;
1058 :
1059 422788 : 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 407288 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1070 :
1071 : // A reference to the constraint row.
1072 489206 : NodeConstraintRow & constraint_row = constraints[my_node].first;
1073 :
1074 203644 : 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 130068 : 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 130068 : if (boundaries.empty())
1098 38996 : return;
1099 :
1100 65034 : libmesh_assert(elem);
1101 :
1102 : // Only constrain active elements with this method
1103 65034 : if (!elem->active())
1104 19498 : return;
1105 :
1106 91072 : const unsigned int Dim = elem->dim();
1107 :
1108 91072 : const FEFamily mapping_family = FEMap::map_fe_type(*elem);
1109 91072 : const FEType fe_type(elem->default_side_order(), mapping_family);
1110 :
1111 : // Pull objects out of the loop to reduce heap operations
1112 91072 : std::vector<const Node *> my_nodes, neigh_nodes;
1113 91072 : 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 91072 : std::vector<boundary_id_type> bc_ids;
1118 443376 : for (auto s : elem->side_index_range())
1119 : {
1120 528456 : if (elem->neighbor_ptr(s))
1121 172686 : continue;
1122 :
1123 6932 : mesh.get_boundary_info().boundary_ids (elem, s, bc_ids);
1124 12816 : for (const auto & boundary_id : bc_ids)
1125 : {
1126 5884 : const PeriodicBoundaryBase * periodic = boundaries.boundary(boundary_id);
1127 5884 : if (periodic)
1128 : {
1129 2716 : libmesh_assert(point_locator);
1130 :
1131 : // Get pointers to the element's neighbor.
1132 : unsigned int s_neigh;
1133 5432 : const Elem * neigh = boundaries.neighbor(boundary_id, *point_locator, elem, s, &s_neigh);
1134 :
1135 5432 : 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 5432 : if (neigh->level() <= elem->level())
1144 : {
1145 : #ifdef LIBMESH_ENABLE_AMR
1146 2268 : libmesh_assert(neigh->active());
1147 : #endif // #ifdef LIBMESH_ENABLE_AMR
1148 :
1149 4536 : elem->build_side_ptr(my_side, s);
1150 4536 : neigh->build_side_ptr(neigh_side, s_neigh);
1151 :
1152 4536 : const unsigned int n_side_nodes = my_side->n_nodes();
1153 :
1154 2268 : my_nodes.clear();
1155 4536 : my_nodes.reserve (n_side_nodes);
1156 2268 : neigh_nodes.clear();
1157 4536 : neigh_nodes.reserve (n_side_nodes);
1158 :
1159 14560 : for (unsigned int n=0; n != n_side_nodes; ++n)
1160 15036 : my_nodes.push_back(my_side->node_ptr(n));
1161 :
1162 14560 : for (unsigned int n=0; n != n_side_nodes; ++n)
1163 15036 : 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 6804 : std::vector<bool> skip_constraint(n_side_nodes, false);
1170 :
1171 12292 : for (unsigned int my_side_n=0;
1172 14560 : 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 5012 : libmesh_assert_less (my_side_n, FEInterface::n_dofs(fe_type, /*extra_order=*/0, my_side.get()));
1177 :
1178 15036 : 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 5012 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1185 :
1186 5012 : if (constraints.count(my_node))
1187 : {
1188 3960 : skip_constraint[my_side_n] = true;
1189 1980 : continue;
1190 : }
1191 : }
1192 :
1193 : // Compute the neighbors's side shape function values.
1194 17702 : for (unsigned int their_side_n=0;
1195 20734 : 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 7335 : libmesh_assert_less (their_side_n, FEInterface::n_dofs(fe_type, /*extra_order=*/0, neigh_side.get()));
1200 :
1201 22005 : 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 7335 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1209 :
1210 7335 : if (!constraints.count(their_node))
1211 4505 : continue;
1212 :
1213 : const NodeConstraintRow & their_constraint_row =
1214 5660 : constraints[their_node].first;
1215 :
1216 20786 : for (unsigned int orig_side_n=0;
1217 23616 : 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 8978 : libmesh_assert_less (orig_side_n, FEInterface::n_dofs(fe_type, /*extra_order=*/0, my_side.get()));
1222 :
1223 26934 : const Node * orig_node = my_nodes[orig_side_n];
1224 :
1225 8978 : if (their_constraint_row.count(orig_node))
1226 11500 : skip_constraint[orig_side_n] = true;
1227 : }
1228 : }
1229 : }
1230 : }
1231 12292 : for (unsigned int my_side_n=0;
1232 14560 : 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 5012 : libmesh_assert_less (my_side_n, FEInterface::n_dofs(fe_type, /*extra_order=*/0, my_side.get()));
1237 :
1238 15036 : if (skip_constraint[my_side_n])
1239 5780 : continue;
1240 :
1241 4244 : const Node * my_node = my_nodes[my_side_n];
1242 :
1243 : // Figure out where my node lies on their reference element.
1244 4244 : 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 4244 : neigh_point);
1250 :
1251 11392 : for (unsigned int their_side_n=0;
1252 13514 : 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 4635 : libmesh_assert_less (their_side_n, FEInterface::n_dofs(fe_type, /*extra_order=*/0, neigh_side.get()));
1257 :
1258 13905 : const Node * their_node = neigh_nodes[their_side_n];
1259 4635 : libmesh_assert(their_node);
1260 :
1261 : // Do not use the p_level(), if any, that is inherited by the side.
1262 9270 : const Real their_value = FEInterface::shape(fe_type,
1263 : /*extra_order=*/0,
1264 : neigh_side.get(),
1265 : their_side_n,
1266 9270 : 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 9270 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1273 :
1274 : NodeConstraintRow & constraint_row =
1275 9270 : constraints[my_node].first;
1276 :
1277 4635 : 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 853261 : unsigned int FEAbstract::n_quadrature_points () const
1292 : {
1293 240269 : if (this->shapes_on_quadrature)
1294 : {
1295 240253 : libmesh_assert(this->qrule);
1296 240253 : libmesh_assert_equal_to(this->qrule->n_points(),
1297 : this->_n_total_qp);
1298 : }
1299 853261 : return this->_n_total_qp;
1300 : }
1301 :
1302 : } // namespace libMesh
|
