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