Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // Local includes
21 : #include "libmesh/dof_map.h"
22 : #include "libmesh/elem.h"
23 : #include "libmesh/enum_to_string.h"
24 : #include "libmesh/fe.h"
25 : #include "libmesh/fe_interface.h"
26 : #include "libmesh/fe_macro.h"
27 : #include "libmesh/remote_elem.h"
28 : #include "libmesh/threads.h"
29 :
30 :
31 : #ifdef LIBMESH_ENABLE_AMR
32 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
33 : // to have the macro 'inf_fe_switch' available.
34 : #include "libmesh/fe_interface_macros.h"
35 : #include "libmesh/inf_fe.h"
36 : #endif
37 : #endif
38 :
39 : namespace libMesh
40 : {
41 :
42 : // global helper function
43 8596910 : void lagrange_nodal_soln(const Elem * elem,
44 : const Order order,
45 : const std::vector<Number> & elem_soln,
46 : std::vector<Number> & nodal_soln,
47 : const bool add_p_level)
48 : {
49 8596910 : const unsigned int n_nodes = elem->n_nodes();
50 8596910 : const ElemType type = elem->type();
51 :
52 9362368 : const Order totalorder = order + add_p_level*elem->p_level();
53 :
54 8596910 : nodal_soln.resize(n_nodes);
55 :
56 :
57 :
58 8596910 : switch (totalorder)
59 : {
60 : // linear Lagrange shape functions
61 7713227 : case FIRST:
62 : {
63 7713227 : switch (type)
64 : {
65 328 : case EDGE3:
66 : {
67 164 : libmesh_assert_equal_to (elem_soln.size(), 2);
68 164 : libmesh_assert_equal_to (nodal_soln.size(), 3);
69 :
70 1968 : nodal_soln[0] = elem_soln[0];
71 1968 : nodal_soln[1] = elem_soln[1];
72 1968 : nodal_soln[2] = .5*(elem_soln[0] + elem_soln[1]);
73 :
74 1968 : return;
75 : }
76 :
77 0 : case EDGE4:
78 : {
79 0 : libmesh_assert_equal_to (elem_soln.size(), 2);
80 0 : libmesh_assert_equal_to (nodal_soln.size(), 4);
81 :
82 0 : nodal_soln[0] = elem_soln[0];
83 0 : nodal_soln[1] = elem_soln[1];
84 0 : nodal_soln[2] = (2.*elem_soln[0] + elem_soln[1])/3.;
85 0 : nodal_soln[3] = (elem_soln[0] + 2.*elem_soln[1])/3.;
86 :
87 0 : return;
88 : }
89 :
90 :
91 1350 : case TRI7:
92 450 : libmesh_assert_equal_to (nodal_soln.size(), 7);
93 5400 : nodal_soln[6] = (elem_soln[0] + elem_soln[1] + elem_soln[2])/3.;
94 : libmesh_fallthrough();
95 15120 : case TRI6:
96 : {
97 1292 : libmesh_assert (type == TRI7 || nodal_soln.size() == 6);
98 1292 : libmesh_assert_equal_to (elem_soln.size(), 3);
99 :
100 15120 : nodal_soln[0] = elem_soln[0];
101 15120 : nodal_soln[1] = elem_soln[1];
102 15120 : nodal_soln[2] = elem_soln[2];
103 15120 : nodal_soln[3] = .5*(elem_soln[0] + elem_soln[1]);
104 15120 : nodal_soln[4] = .5*(elem_soln[1] + elem_soln[2]);
105 15120 : nodal_soln[5] = .5*(elem_soln[2] + elem_soln[0]);
106 :
107 15120 : return;
108 : }
109 :
110 :
111 22647 : case QUAD8:
112 : case QUAD9:
113 : {
114 22647 : libmesh_assert_equal_to (elem_soln.size(), 4);
115 :
116 22647 : if (type == QUAD8)
117 0 : libmesh_assert_equal_to (nodal_soln.size(), 8);
118 : else
119 22647 : libmesh_assert_equal_to (nodal_soln.size(), 9);
120 :
121 :
122 269869 : nodal_soln[0] = elem_soln[0];
123 269869 : nodal_soln[1] = elem_soln[1];
124 269869 : nodal_soln[2] = elem_soln[2];
125 269869 : nodal_soln[3] = elem_soln[3];
126 269869 : nodal_soln[4] = .5*(elem_soln[0] + elem_soln[1]);
127 269869 : nodal_soln[5] = .5*(elem_soln[1] + elem_soln[2]);
128 269869 : nodal_soln[6] = .5*(elem_soln[2] + elem_soln[3]);
129 269869 : nodal_soln[7] = .5*(elem_soln[3] + elem_soln[0]);
130 :
131 269869 : if (type == QUAD9)
132 269869 : nodal_soln[8] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3]);
133 :
134 22647 : return;
135 : }
136 :
137 :
138 4608 : case TET14:
139 1536 : libmesh_assert_equal_to (nodal_soln.size(), 14);
140 18432 : nodal_soln[10] = (elem_soln[0] + elem_soln[1] + elem_soln[2])/3.;
141 18432 : nodal_soln[11] = (elem_soln[0] + elem_soln[1] + elem_soln[3])/3.;
142 18432 : nodal_soln[12] = (elem_soln[1] + elem_soln[2] + elem_soln[3])/3.;
143 18432 : nodal_soln[13] = (elem_soln[0] + elem_soln[2] + elem_soln[3])/3.;
144 : libmesh_fallthrough();
145 18432 : case TET10:
146 : {
147 1536 : libmesh_assert_equal_to (elem_soln.size(), 4);
148 1536 : libmesh_assert (type == TET14 || nodal_soln.size() == 10);
149 :
150 18432 : nodal_soln[0] = elem_soln[0];
151 18432 : nodal_soln[1] = elem_soln[1];
152 18432 : nodal_soln[2] = elem_soln[2];
153 18432 : nodal_soln[3] = elem_soln[3];
154 18432 : nodal_soln[4] = .5*(elem_soln[0] + elem_soln[1]);
155 18432 : nodal_soln[5] = .5*(elem_soln[1] + elem_soln[2]);
156 18432 : nodal_soln[6] = .5*(elem_soln[2] + elem_soln[0]);
157 18432 : nodal_soln[7] = .5*(elem_soln[3] + elem_soln[0]);
158 18432 : nodal_soln[8] = .5*(elem_soln[3] + elem_soln[1]);
159 18432 : nodal_soln[9] = .5*(elem_soln[3] + elem_soln[2]);
160 :
161 18432 : return;
162 : }
163 :
164 :
165 1936 : case HEX20:
166 : case HEX27:
167 : {
168 1936 : libmesh_assert_equal_to (elem_soln.size(), 8);
169 :
170 1936 : if (type == HEX20)
171 0 : libmesh_assert_equal_to (nodal_soln.size(), 20);
172 : else
173 1936 : libmesh_assert_equal_to (nodal_soln.size(), 27);
174 :
175 23112 : nodal_soln[0] = elem_soln[0];
176 23112 : nodal_soln[1] = elem_soln[1];
177 23112 : nodal_soln[2] = elem_soln[2];
178 23112 : nodal_soln[3] = elem_soln[3];
179 23112 : nodal_soln[4] = elem_soln[4];
180 23112 : nodal_soln[5] = elem_soln[5];
181 23112 : nodal_soln[6] = elem_soln[6];
182 23112 : nodal_soln[7] = elem_soln[7];
183 23112 : nodal_soln[8] = .5*(elem_soln[0] + elem_soln[1]);
184 23112 : nodal_soln[9] = .5*(elem_soln[1] + elem_soln[2]);
185 23112 : nodal_soln[10] = .5*(elem_soln[2] + elem_soln[3]);
186 23112 : nodal_soln[11] = .5*(elem_soln[3] + elem_soln[0]);
187 23112 : nodal_soln[12] = .5*(elem_soln[0] + elem_soln[4]);
188 23112 : nodal_soln[13] = .5*(elem_soln[1] + elem_soln[5]);
189 23112 : nodal_soln[14] = .5*(elem_soln[2] + elem_soln[6]);
190 23112 : nodal_soln[15] = .5*(elem_soln[3] + elem_soln[7]);
191 23112 : nodal_soln[16] = .5*(elem_soln[4] + elem_soln[5]);
192 23112 : nodal_soln[17] = .5*(elem_soln[5] + elem_soln[6]);
193 23112 : nodal_soln[18] = .5*(elem_soln[6] + elem_soln[7]);
194 23112 : nodal_soln[19] = .5*(elem_soln[4] + elem_soln[7]);
195 :
196 23112 : if (type == HEX27)
197 : {
198 23112 : nodal_soln[20] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3]);
199 23112 : nodal_soln[21] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[4] + elem_soln[5]);
200 23112 : nodal_soln[22] = .25*(elem_soln[1] + elem_soln[2] + elem_soln[5] + elem_soln[6]);
201 23112 : nodal_soln[23] = .25*(elem_soln[2] + elem_soln[3] + elem_soln[6] + elem_soln[7]);
202 23112 : nodal_soln[24] = .25*(elem_soln[3] + elem_soln[0] + elem_soln[7] + elem_soln[4]);
203 23112 : nodal_soln[25] = .25*(elem_soln[4] + elem_soln[5] + elem_soln[6] + elem_soln[7]);
204 :
205 26984 : nodal_soln[26] = .125*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3] +
206 21296 : elem_soln[4] + elem_soln[5] + elem_soln[6] + elem_soln[7]);
207 : }
208 :
209 1936 : return;
210 : }
211 :
212 :
213 0 : case PRISM21:
214 0 : nodal_soln[20] = (elem_soln[9] + elem_soln[10] + elem_soln[11])/Real(3);
215 : libmesh_fallthrough();
216 0 : case PRISM20:
217 0 : if (type == PRISM20)
218 0 : libmesh_assert_equal_to (nodal_soln.size(), 20);
219 0 : nodal_soln[18] = (elem_soln[0] + elem_soln[1] + elem_soln[2])/Real(3);
220 0 : nodal_soln[19] = (elem_soln[3] + elem_soln[4] + elem_soln[5])/Real(3);
221 : libmesh_fallthrough();
222 0 : case PRISM18:
223 0 : if (type == PRISM18)
224 0 : libmesh_assert_equal_to (nodal_soln.size(), 18);
225 0 : nodal_soln[15] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[4] + elem_soln[3]);
226 0 : nodal_soln[16] = .25*(elem_soln[1] + elem_soln[2] + elem_soln[5] + elem_soln[4]);
227 0 : nodal_soln[17] = .25*(elem_soln[2] + elem_soln[0] + elem_soln[3] + elem_soln[5]);
228 : libmesh_fallthrough();
229 0 : case PRISM15:
230 : {
231 0 : libmesh_assert_equal_to (elem_soln.size(), 6);
232 :
233 0 : if (type == PRISM15)
234 0 : libmesh_assert_equal_to (nodal_soln.size(), 15);
235 :
236 0 : nodal_soln[0] = elem_soln[0];
237 0 : nodal_soln[1] = elem_soln[1];
238 0 : nodal_soln[2] = elem_soln[2];
239 0 : nodal_soln[3] = elem_soln[3];
240 0 : nodal_soln[4] = elem_soln[4];
241 0 : nodal_soln[5] = elem_soln[5];
242 0 : nodal_soln[6] = .5*(elem_soln[0] + elem_soln[1]);
243 0 : nodal_soln[7] = .5*(elem_soln[1] + elem_soln[2]);
244 0 : nodal_soln[8] = .5*(elem_soln[0] + elem_soln[2]);
245 0 : nodal_soln[9] = .5*(elem_soln[0] + elem_soln[3]);
246 0 : nodal_soln[10] = .5*(elem_soln[1] + elem_soln[4]);
247 0 : nodal_soln[11] = .5*(elem_soln[2] + elem_soln[5]);
248 0 : nodal_soln[12] = .5*(elem_soln[3] + elem_soln[4]);
249 0 : nodal_soln[13] = .5*(elem_soln[4] + elem_soln[5]);
250 0 : nodal_soln[14] = .5*(elem_soln[3] + elem_soln[5]);
251 :
252 0 : return;
253 : }
254 :
255 0 : case PYRAMID18:
256 : {
257 0 : libmesh_assert_equal_to (nodal_soln.size(), 18);
258 :
259 0 : nodal_soln[14] = (elem_soln[0] + elem_soln[1] + elem_soln[4])/Real(3);
260 0 : nodal_soln[15] = (elem_soln[1] + elem_soln[2] + elem_soln[4])/Real(3);
261 0 : nodal_soln[16] = (elem_soln[2] + elem_soln[3] + elem_soln[4])/Real(3);
262 0 : nodal_soln[17] = (elem_soln[0] + elem_soln[3] + elem_soln[4])/Real(3);
263 :
264 : libmesh_fallthrough();
265 : }
266 :
267 0 : case PYRAMID14:
268 : {
269 0 : if (type == PYRAMID14)
270 0 : libmesh_assert_equal_to (nodal_soln.size(), 14);
271 :
272 0 : nodal_soln[13] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3]);
273 :
274 : libmesh_fallthrough();
275 : }
276 :
277 0 : case PYRAMID13:
278 : {
279 0 : libmesh_assert_equal_to (elem_soln.size(), 5);
280 :
281 0 : if (type == PYRAMID13)
282 0 : libmesh_assert_equal_to (nodal_soln.size(), 13);
283 :
284 0 : nodal_soln[0] = elem_soln[0];
285 0 : nodal_soln[1] = elem_soln[1];
286 0 : nodal_soln[2] = elem_soln[2];
287 0 : nodal_soln[3] = elem_soln[3];
288 0 : nodal_soln[4] = elem_soln[4];
289 0 : nodal_soln[5] = .5*(elem_soln[0] + elem_soln[1]);
290 0 : nodal_soln[6] = .5*(elem_soln[1] + elem_soln[2]);
291 0 : nodal_soln[7] = .5*(elem_soln[2] + elem_soln[3]);
292 0 : nodal_soln[8] = .5*(elem_soln[3] + elem_soln[0]);
293 0 : nodal_soln[9] = .5*(elem_soln[0] + elem_soln[4]);
294 0 : nodal_soln[10] = .5*(elem_soln[1] + elem_soln[4]);
295 0 : nodal_soln[11] = .5*(elem_soln[2] + elem_soln[4]);
296 0 : nodal_soln[12] = .5*(elem_soln[3] + elem_soln[4]);
297 :
298 0 : return;
299 : }
300 7384726 : default:
301 : {
302 : // By default the element solution _is_ nodal,
303 : // so just copy it.
304 7384726 : nodal_soln = elem_soln;
305 :
306 7384726 : return;
307 : }
308 : }
309 : }
310 :
311 881163 : case SECOND:
312 : {
313 881163 : switch (type)
314 : {
315 0 : case EDGE4:
316 : {
317 0 : libmesh_assert_equal_to (elem_soln.size(), 3);
318 0 : libmesh_assert_equal_to (nodal_soln.size(), 4);
319 :
320 : // Project quadratic solution onto cubic element nodes
321 0 : nodal_soln[0] = elem_soln[0];
322 0 : nodal_soln[1] = elem_soln[1];
323 0 : nodal_soln[2] = (2.*elem_soln[0] - elem_soln[1] +
324 0 : 8.*elem_soln[2])/9.;
325 0 : nodal_soln[3] = (-elem_soln[0] + 2.*elem_soln[1] +
326 0 : 8.*elem_soln[2])/9.;
327 0 : return;
328 : }
329 :
330 450 : case TRI7:
331 : {
332 450 : libmesh_assert_equal_to (elem_soln.size(), 6);
333 450 : libmesh_assert_equal_to (nodal_soln.size(), 7);
334 :
335 37800 : for (int i=0; i != 6; ++i)
336 37800 : nodal_soln[i] = elem_soln[i];
337 :
338 5850 : nodal_soln[6] = -1./9. * (elem_soln[0] + elem_soln[1] + elem_soln[2])
339 4950 : +4./9. * (elem_soln[3] + elem_soln[4] + elem_soln[5]);
340 :
341 5400 : return;
342 : }
343 :
344 1536 : case TET14:
345 : {
346 1536 : libmesh_assert_equal_to (elem_soln.size(), 10);
347 1536 : libmesh_assert_equal_to (nodal_soln.size(), 14);
348 :
349 202752 : for (int i=0; i != 10; ++i)
350 215040 : nodal_soln[i] = elem_soln[i];
351 :
352 19968 : nodal_soln[10] = -1./9. * (elem_soln[0] + elem_soln[1] + elem_soln[2])
353 16896 : +4./9. * (elem_soln[4] + elem_soln[5] + elem_soln[6]);
354 19968 : nodal_soln[11] = -1./9. * (elem_soln[0] + elem_soln[1] + elem_soln[3])
355 16896 : +4./9. * (elem_soln[4] + elem_soln[7] + elem_soln[8]);
356 19968 : nodal_soln[12] = -1./9. * (elem_soln[1] + elem_soln[2] + elem_soln[3])
357 16896 : +4./9. * (elem_soln[5] + elem_soln[8] + elem_soln[9]);
358 19968 : nodal_soln[13] = -1./9. * (elem_soln[0] + elem_soln[2] + elem_soln[3])
359 16896 : +4./9. * (elem_soln[6] + elem_soln[7] + elem_soln[9]);
360 :
361 18432 : return;
362 : }
363 :
364 606 : case PRISM21:
365 : {
366 2424 : nodal_soln[20] = (elem_soln[9] + elem_soln[10] + elem_soln[11])/Real(3);
367 : libmesh_fallthrough();
368 : }
369 202 : case PRISM20:
370 : {
371 202 : if (type == PRISM20)
372 0 : libmesh_assert_equal_to (nodal_soln.size(), 20);
373 :
374 46056 : for (int i=0; i != 18; ++i)
375 50904 : nodal_soln[i] = elem_soln[i];
376 :
377 2424 : nodal_soln[18] = (elem_soln[0] + elem_soln[1] + elem_soln[2])/Real(3);
378 2424 : nodal_soln[19] = (elem_soln[3] + elem_soln[4] + elem_soln[5])/Real(3);
379 2424 : return;
380 : }
381 :
382 0 : case PYRAMID18:
383 : {
384 0 : libmesh_assert_equal_to (nodal_soln.size(), 18);
385 :
386 0 : for (int i=0; i != 14; ++i)
387 0 : nodal_soln[i] = elem_soln[i];
388 :
389 0 : nodal_soln[14] = (elem_soln[0] + elem_soln[1] + elem_soln[4])/Real(3);
390 0 : nodal_soln[15] = (elem_soln[1] + elem_soln[2] + elem_soln[4])/Real(3);
391 0 : nodal_soln[16] = (elem_soln[2] + elem_soln[3] + elem_soln[4])/Real(3);
392 0 : nodal_soln[17] = (elem_soln[0] + elem_soln[3] + elem_soln[4])/Real(3);
393 :
394 : libmesh_fallthrough();
395 : }
396 :
397 854907 : default:
398 : {
399 : // By default the element solution _is_ nodal, so just
400 : // copy the portion relevant to the nodal solution.
401 : // (this is the whole nodal solution for true Lagrange
402 : // elements, but a smaller part for "L2 Lagrange"
403 75879 : libmesh_assert_less_equal(nodal_soln.size(), elem_soln.size());
404 8558616 : for (const auto i : index_range(nodal_soln))
405 9070433 : nodal_soln[i] = elem_soln[i];
406 :
407 151758 : return;
408 : }
409 : }
410 : }
411 :
412 :
413 :
414 :
415 210 : default:
416 : {
417 : // By default the element solution _is_ nodal, so just copy
418 : // the portion relevant to the nodal solution. (this is the
419 : // whole nodal solution for true Lagrange elements, but a
420 : // smaller part for "L2 Lagrange"
421 210 : libmesh_assert_less_equal(nodal_soln.size(), elem_soln.size());
422 49392 : for (const auto i : index_range(nodal_soln))
423 54684 : nodal_soln[i] = elem_soln[i];
424 :
425 420 : return;
426 : }
427 : }
428 : }
429 :
430 :
431 : // Anonymous namespace for local helper functions
432 : namespace {
433 1238611585 : unsigned int lagrange_n_dofs(const ElemType t, const Elem * e, const Order o)
434 : {
435 108235564 : libmesh_assert(!e || e->type() == t);
436 :
437 1238611585 : switch (o)
438 : {
439 : // lagrange can only be constant on a single node
440 146021 : case CONSTANT:
441 : {
442 146021 : switch (t)
443 : {
444 12157 : case NODEELEM:
445 12157 : return 1;
446 :
447 0 : default:
448 0 : libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
449 : }
450 : }
451 :
452 : // linear Lagrange shape functions
453 289996873 : case FIRST:
454 : {
455 289996873 : switch (t)
456 : {
457 1352 : case NODEELEM:
458 1352 : return 1;
459 :
460 14590287 : case EDGE2:
461 : case EDGE3:
462 : case EDGE4:
463 14590287 : return 2;
464 :
465 6633714 : case TRI3:
466 : case TRISHELL3:
467 : case TRI3SUBDIVISION:
468 : case TRI6:
469 : case TRI7:
470 6633714 : return 3;
471 :
472 23461617 : case QUAD4:
473 : case QUADSHELL4:
474 : case QUAD8:
475 : case QUADSHELL8:
476 : case QUAD9:
477 : case QUADSHELL9:
478 23461617 : return 4;
479 :
480 1450476 : case TET4:
481 : case TET10:
482 : case TET14:
483 1450476 : return 4;
484 :
485 16721010 : case HEX8:
486 : case HEX20:
487 : case HEX27:
488 16721010 : return 8;
489 :
490 1754306 : case PRISM6:
491 : case PRISM15:
492 : case PRISM18:
493 : case PRISM20:
494 : case PRISM21:
495 1754306 : return 6;
496 :
497 1338847 : case PYRAMID5:
498 : case PYRAMID13:
499 : case PYRAMID14:
500 : case PYRAMID18:
501 1338847 : return 5;
502 :
503 0 : case INVALID_ELEM:
504 0 : return 0;
505 :
506 3714952 : case C0POLYGON:
507 : case C0POLYHEDRON:
508 : // Polygons and polyhedra require using newer FE APIs
509 3714952 : if (!e)
510 0 : libmesh_error_msg("Code (see stack trace) used an outdated FE function overload.\n"
511 631722 : "n_dofs() on a polygon or polyhedron is not defined by ElemType alone.");
512 3714952 : return e->n_nodes();
513 :
514 0 : default:
515 0 : libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
516 : }
517 : }
518 :
519 :
520 : // quadratic Lagrange shape functions
521 723048389 : case SECOND:
522 : {
523 723048389 : switch (t)
524 : {
525 195 : case NODEELEM:
526 195 : return 1;
527 :
528 5327004 : case EDGE3:
529 5327004 : return 3;
530 :
531 47285778 : case TRI6:
532 : case TRI7:
533 47285778 : return 6;
534 :
535 1321957 : case QUAD8:
536 : case QUADSHELL8:
537 1321957 : return 8;
538 :
539 21937851 : case QUAD9:
540 : case QUADSHELL9:
541 21937851 : return 9;
542 :
543 70566205 : case TET10:
544 : case TET14:
545 70566205 : return 10;
546 :
547 1380349 : case HEX20:
548 1380349 : return 20;
549 :
550 117730463 : case HEX27:
551 117730463 : return 27;
552 :
553 12649134 : case PRISM15:
554 12649134 : return 15;
555 :
556 3076907 : case PRISM18:
557 : case PRISM20:
558 : case PRISM21:
559 3076907 : return 18;
560 :
561 5209427 : case PYRAMID13:
562 5209427 : return 13;
563 :
564 353845 : case PYRAMID14:
565 : case PYRAMID18:
566 353845 : return 14;
567 :
568 0 : case INVALID_ELEM:
569 0 : return 0;
570 :
571 0 : default:
572 0 : libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
573 : }
574 : }
575 :
576 225420302 : case THIRD:
577 : {
578 225420302 : switch (t)
579 : {
580 0 : case NODEELEM:
581 0 : return 1;
582 :
583 616 : case EDGE4:
584 616 : return 4;
585 :
586 2274767 : case PRISM20:
587 2274767 : return 20;
588 :
589 3098398 : case PRISM21:
590 3098398 : return 21;
591 :
592 496071 : case PYRAMID18:
593 496071 : return 18;
594 :
595 82562892 : case TRI7:
596 82562892 : return 7;
597 :
598 6903944 : case TET14:
599 6903944 : return 14;
600 :
601 0 : case INVALID_ELEM:
602 0 : return 0;
603 :
604 0 : default:
605 0 : libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
606 : }
607 : }
608 :
609 0 : default:
610 0 : libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for LAGRANGE FE family!");
611 : }
612 : }
613 :
614 :
615 :
616 1142628466 : unsigned int lagrange_n_dofs_at_node(const ElemType t,
617 : const Order o,
618 : const unsigned int n)
619 : {
620 1142628466 : switch (o)
621 : {
622 : // lagrange can only be constant on a single node
623 0 : case CONSTANT:
624 : {
625 0 : switch (t)
626 : {
627 0 : case NODEELEM:
628 0 : return 1;
629 :
630 0 : default:
631 0 : libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
632 : }
633 : }
634 :
635 : // linear Lagrange shape functions
636 831007147 : case FIRST:
637 : {
638 831007147 : switch (t)
639 : {
640 1648 : case NODEELEM:
641 1648 : return 1;
642 :
643 26848902 : case EDGE2:
644 : case EDGE3:
645 : case EDGE4:
646 : {
647 26848902 : switch (n)
648 : {
649 1314237 : case 0:
650 : case 1:
651 1314237 : return 1;
652 :
653 1208462 : default:
654 1208462 : return 0;
655 : }
656 : }
657 :
658 74388071 : case TRI3:
659 : case TRISHELL3:
660 : case TRI3SUBDIVISION:
661 : case TRI6:
662 : case TRI7:
663 : {
664 74388071 : switch (n)
665 : {
666 4834230 : case 0:
667 : case 1:
668 : case 2:
669 4834230 : return 1;
670 :
671 20237 : default:
672 20237 : return 0;
673 : }
674 : }
675 :
676 491925633 : case QUAD4:
677 : case QUADSHELL4:
678 : case QUAD8:
679 : case QUADSHELL8:
680 : case QUAD9:
681 : case QUADSHELL9:
682 : {
683 491925633 : switch (n)
684 : {
685 42859705 : case 0:
686 : case 1:
687 : case 2:
688 : case 3:
689 42859705 : return 1;
690 :
691 6348191 : default:
692 6348191 : return 0;
693 : }
694 : }
695 :
696 2322 : case C0POLYGON:
697 : case C0POLYHEDRON:
698 2322 : return 1;
699 :
700 :
701 71102891 : case TET4:
702 : case TET10:
703 : case TET14:
704 : {
705 71102891 : switch (n)
706 : {
707 4126622 : case 0:
708 : case 1:
709 : case 2:
710 : case 3:
711 4126622 : return 1;
712 :
713 0 : default:
714 0 : return 0;
715 : }
716 : }
717 :
718 165289507 : case HEX8:
719 : case HEX20:
720 : case HEX27:
721 : {
722 165289507 : switch (n)
723 : {
724 11925167 : case 0:
725 : case 1:
726 : case 2:
727 : case 3:
728 : case 4:
729 : case 5:
730 : case 6:
731 : case 7:
732 11925167 : return 1;
733 :
734 2603770 : default:
735 2603770 : return 0;
736 : }
737 : }
738 :
739 1319257 : case PRISM6:
740 : case PRISM15:
741 : case PRISM18:
742 : case PRISM20:
743 : case PRISM21:
744 : {
745 1319257 : switch (n)
746 : {
747 81391 : case 0:
748 : case 1:
749 : case 2:
750 : case 3:
751 : case 4:
752 : case 5:
753 81391 : return 1;
754 :
755 0 : default:
756 0 : return 0;
757 : }
758 : }
759 :
760 70092 : case PYRAMID5:
761 : case PYRAMID13:
762 : case PYRAMID14:
763 : case PYRAMID18:
764 : {
765 70092 : switch (n)
766 : {
767 3582 : case 0:
768 : case 1:
769 : case 2:
770 : case 3:
771 : case 4:
772 3582 : return 1;
773 :
774 0 : default:
775 0 : return 0;
776 : }
777 : }
778 :
779 0 : case INVALID_ELEM:
780 0 : return 0;
781 :
782 0 : default:
783 0 : libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
784 : }
785 : }
786 :
787 : // quadratic Lagrange shape functions
788 303776555 : case SECOND:
789 : {
790 303776555 : switch (t)
791 : {
792 : // quadratic lagrange has one dof at each node
793 26012262 : case NODEELEM:
794 : case EDGE3:
795 : case TRI6:
796 : case QUAD8:
797 : case QUADSHELL8:
798 : case QUAD9:
799 : case QUADSHELL9:
800 : case TET10:
801 : case HEX20:
802 : case HEX27:
803 : case PRISM15:
804 : case PRISM18:
805 : case PYRAMID13:
806 : case PYRAMID14:
807 26012262 : return 1;
808 :
809 3645483 : case PRISM20:
810 : case PRISM21:
811 3645483 : return (n < 18);
812 :
813 0 : case PYRAMID18:
814 0 : return (n < 14);
815 :
816 29406 : case TRI7:
817 29406 : return (n < 6);
818 :
819 782361 : case TET14:
820 782361 : return (n < 10);
821 :
822 0 : case INVALID_ELEM:
823 0 : return 0;
824 :
825 0 : default:
826 0 : libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
827 : }
828 : }
829 :
830 7844764 : case THIRD:
831 : {
832 7293470 : switch (t)
833 : {
834 551056 : case NODEELEM:
835 : case EDGE4:
836 : case PRISM20:
837 : case PRISM21:
838 : case PYRAMID18:
839 : case TRI7:
840 : case TET14:
841 551056 : return 1;
842 :
843 0 : case INVALID_ELEM:
844 0 : return 0;
845 :
846 0 : default:
847 0 : libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
848 : }
849 : }
850 :
851 0 : default:
852 0 : libmesh_error_msg("Unsupported order: " << o );
853 : }
854 : }
855 :
856 :
857 :
858 87646174 : unsigned int lagrange_n_dofs_at_node(const Elem & e,
859 : const Order o,
860 : const unsigned int n)
861 : {
862 94307864 : return lagrange_n_dofs_at_node(e.type(), o, n);
863 : }
864 :
865 :
866 : #ifdef LIBMESH_ENABLE_AMR
867 4343507 : void lagrange_compute_constraints (DofConstraints & constraints,
868 : DofMap & dof_map,
869 : const unsigned int variable_number,
870 : const Elem * elem,
871 : const unsigned Dim)
872 : {
873 : // Only constrain elements in 2,3D.
874 4343507 : if (Dim == 1)
875 380840 : return;
876 :
877 414697 : libmesh_assert(elem);
878 :
879 : // Only constrain active and ancestor elements
880 4343507 : if (elem->subactive())
881 36396 : return;
882 :
883 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
884 1051705 : if (elem->infinite())
885 : {
886 3212 : const FEType fe_t = dof_map.variable_type(variable_number);
887 :
888 : // expand the infinite_compute_constraint in its template-arguments.
889 3212 : switch(Dim)
890 : {
891 0 : case 2:
892 : {
893 0 : inf_fe_family_mapping_switch(2, inf_compute_constraints (constraints, dof_map, variable_number, elem) , ,; break;);
894 0 : break;
895 : }
896 3212 : case 3:
897 : {
898 3212 : inf_fe_family_mapping_switch(3, inf_compute_constraints (constraints, dof_map, variable_number, elem) , ,; break;);
899 1056 : break;
900 : }
901 0 : default:
902 0 : libmesh_error_msg("Invalid dim = " << Dim);
903 : }
904 1056 : return;
905 : }
906 : #endif
907 :
908 3962667 : const Variable & var = dof_map.variable(variable_number);
909 3962667 : const FEType fe_type = var.type();
910 3962667 : const bool add_p_level = dof_map.should_p_refine_var(variable_number);
911 :
912 : // Pull objects out of the loop to reduce heap operations
913 754490 : std::vector<dof_id_type> my_dof_indices, parent_dof_indices;
914 4256670 : std::unique_ptr<const Elem> my_side, parent_side;
915 :
916 : // Look at the element faces. Check to see if we need to
917 : // build constraints.
918 20446443 : for (auto s : elem->side_index_range())
919 18045061 : if (const Elem * const neigh = elem->neighbor_ptr(s);
920 15478633 : neigh && neigh != remote_elem)
921 : {
922 : // constrain dofs shared between this element and ones coarser
923 15509827 : if (neigh->level() < elem->level() &&
924 318852 : var.active_on_subdomain(neigh->subdomain_id()))
925 : {
926 : // than this element.
927 : // Get pointers to the elements of interest and its parent.
928 72316 : const Elem * parent = elem->parent();
929 :
930 : // This can't happen... Only level-0 elements have nullptr
931 : // parents, and no level-0 elements can be at a higher
932 : // level than their neighbors!
933 36158 : libmesh_assert(parent);
934 :
935 318510 : elem->build_side_ptr(my_side, s);
936 318510 : parent->build_side_ptr(parent_side, s);
937 :
938 : // We have some elements with interior bubble function bases
939 : // and lower-order sides. We need to only ask for the lower
940 : // order in those cases.
941 318510 : FEType side_fe_type = fe_type;
942 36158 : int extra_order = 0;
943 318510 : if (my_side->default_order() < fe_type.order)
944 : {
945 1104 : side_fe_type.order = my_side->default_order();
946 1196 : extra_order = (int)(side_fe_type.order) -
947 92 : (int)(fe_type.order);
948 : }
949 :
950 : // This function gets called element-by-element, so there
951 : // will be a lot of memory allocation going on. We can
952 : // at least minimize this for the case of the dof indices
953 : // by efficiently preallocating the requisite storage.
954 318510 : my_dof_indices.reserve (my_side->n_nodes());
955 318510 : parent_dof_indices.reserve (parent_side->n_nodes());
956 :
957 318510 : dof_map.dof_indices (my_side.get(), my_dof_indices,
958 72316 : variable_number, extra_order);
959 318510 : dof_map.dof_indices (parent_side.get(), parent_dof_indices,
960 72316 : variable_number, extra_order);
961 :
962 : const unsigned int n_side_dofs =
963 318510 : FEInterface::n_dofs(side_fe_type, my_side.get());
964 : const unsigned int n_parent_side_dofs =
965 318510 : FEInterface::n_dofs(side_fe_type, parent_side.get());
966 1167050 : for (unsigned int my_dof=0; my_dof != n_side_dofs; my_dof++)
967 : {
968 93762 : libmesh_assert_less (my_dof, my_side->n_nodes());
969 :
970 : // My global dof index.
971 942302 : const dof_id_type my_dof_g = my_dof_indices[my_dof];
972 :
973 : // Hunt for "constraining against myself" cases before
974 : // we bother creating a constraint row
975 93762 : bool self_constraint = false;
976 2743303 : for (unsigned int their_dof=0;
977 2837065 : their_dof != n_parent_side_dofs; their_dof++)
978 : {
979 255853 : libmesh_assert_less (their_dof, parent_side->n_nodes());
980 :
981 : // Their global dof index.
982 : const dof_id_type their_dof_g =
983 2344019 : parent_dof_indices[their_dof];
984 :
985 2344019 : if (their_dof_g == my_dof_g)
986 : {
987 40858 : self_constraint = true;
988 40858 : break;
989 : }
990 : }
991 :
992 848540 : if (self_constraint)
993 635005 : continue;
994 :
995 : DofConstraintRow * constraint_row;
996 :
997 : // we may be running constraint methods concurrently
998 : // on multiple threads, so we need a lock to
999 : // ensure that this constraint is "ours"
1000 : {
1001 52904 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1002 :
1003 493046 : if (dof_map.is_constrained_dof(my_dof_g))
1004 29643 : continue;
1005 :
1006 213535 : constraint_row = &(constraints[my_dof_g]);
1007 23261 : libmesh_assert(constraint_row->empty());
1008 : }
1009 :
1010 : // The support point of the DOF
1011 46522 : const Point & support_point = my_side->point(my_dof);
1012 :
1013 : // Figure out where my node lies on their reference element.
1014 167013 : const Point mapped_point = FEMap::inverse_map(Dim-1,
1015 : parent_side.get(),
1016 46522 : support_point);
1017 :
1018 : // Compute the parent's side shape function values.
1019 702354 : for (unsigned int their_dof=0;
1020 892628 : their_dof != n_parent_side_dofs; their_dof++)
1021 : {
1022 73275 : libmesh_assert_less (their_dof, parent_side->n_nodes());
1023 :
1024 : // Their global dof index.
1025 : const dof_id_type their_dof_g =
1026 752368 : parent_dof_indices[their_dof];
1027 :
1028 679093 : const Real their_dof_value = FEInterface::shape(side_fe_type,
1029 : parent_side.get(),
1030 : their_dof,
1031 679093 : mapped_point);
1032 :
1033 : // Only add non-zero and non-identity values
1034 : // for Lagrange basis functions.
1035 737191 : if ((std::abs(their_dof_value) > 1.e-5) &&
1036 58098 : (std::abs(their_dof_value) < .999))
1037 : {
1038 58098 : constraint_row->emplace(their_dof_g, their_dof_value);
1039 : }
1040 : #ifdef DEBUG
1041 : // Protect for the case u_i = 0.999 u_j,
1042 : // in which case i better equal j.
1043 15177 : else if (their_dof_value >= .999)
1044 0 : libmesh_assert_equal_to (my_dof_g, their_dof_g);
1045 : #endif
1046 : }
1047 : }
1048 : }
1049 :
1050 13601486 : if (elem->active() && add_p_level)
1051 : {
1052 : // p refinement constraints:
1053 : // constrain dofs shared between
1054 : // active elements and neighbors with
1055 : // lower polynomial degrees
1056 : const unsigned int min_p_level =
1057 14640664 : neigh->min_p_level_by_neighbor(elem, elem->p_level());
1058 14640664 : if (min_p_level < elem->p_level())
1059 0 : libmesh_error_msg(
1060 : "Mismatched p-refinement levels for LAGRANGE is not currently supported");
1061 : }
1062 : }
1063 2914174 : } // lagrange_compute_constraints()
1064 : #endif // #ifdef LIBMESH_ENABLE_AMR
1065 :
1066 : } // anonymous namespace
1067 :
1068 :
1069 : // Instantiate (side_) nodal_soln() function for every dimension
1070 8483027 : LIBMESH_FE_NODAL_SOLN(LAGRANGE, lagrange_nodal_soln)
1071 1584 : LIBMESH_FE_SIDE_NODAL_SOLN(LAGRANGE)
1072 :
1073 :
1074 : // Do full-specialization for every dimension, instead
1075 : // of explicit instantiation at the end of this function.
1076 : // This could be macro-ified.
1077 0 : template <> unsigned int FE<0,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, nullptr, o); }
1078 0 : template <> unsigned int FE<1,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, nullptr, o); }
1079 1590 : template <> unsigned int FE<2,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, nullptr, o); }
1080 0 : template <> unsigned int FE<3,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, nullptr, o); }
1081 :
1082 150384 : template <> unsigned int FE<0,LAGRANGE>::n_dofs(const Elem * e, const Order o) { return lagrange_n_dofs(e->type(), e, o); }
1083 113029705 : template <> unsigned int FE<1,LAGRANGE>::n_dofs(const Elem * e, const Order o) { return lagrange_n_dofs(e->type(), e, o); }
1084 630760768 : template <> unsigned int FE<2,LAGRANGE>::n_dofs(const Elem * e, const Order o) { return lagrange_n_dofs(e->type(), e, o); }
1085 494669138 : template <> unsigned int FE<3,LAGRANGE>::n_dofs(const Elem * e, const Order o) { return lagrange_n_dofs(e->type(), e, o); }
1086 :
1087 :
1088 : // Do full-specialization for every dimension, instead
1089 : // of explicit instantiation at the end of this function.
1090 21208 : template <> unsigned int FE<0,LAGRANGE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(t, o, n); }
1091 27016723 : template <> unsigned int FE<1,LAGRANGE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(t, o, n); }
1092 775037175 : template <> unsigned int FE<2,LAGRANGE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(t, o, n); }
1093 246245496 : template <> unsigned int FE<3,LAGRANGE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(t, o, n); }
1094 :
1095 23152 : template <> unsigned int FE<0,LAGRANGE>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(e, o, n); }
1096 98853 : template <> unsigned int FE<1,LAGRANGE>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(e, o, n); }
1097 66730735 : template <> unsigned int FE<2,LAGRANGE>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(e, o, n); }
1098 27455124 : template <> unsigned int FE<3,LAGRANGE>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(e, o, n); }
1099 :
1100 :
1101 : // Lagrange elements have no dofs per element
1102 : // (just at the nodes)
1103 0 : template <> unsigned int FE<0,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1104 0 : template <> unsigned int FE<1,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1105 0 : template <> unsigned int FE<2,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1106 0 : template <> unsigned int FE<3,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1107 :
1108 32784 : template <> unsigned int FE<0,LAGRANGE>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1109 9172944 : template <> unsigned int FE<1,LAGRANGE>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1110 166466859 : template <> unsigned int FE<2,LAGRANGE>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1111 36284706 : template <> unsigned int FE<3,LAGRANGE>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1112 :
1113 : // Lagrange FEMs are always C^0 continuous
1114 0 : template <> FEContinuity FE<0,LAGRANGE>::get_continuity() const { return C_ZERO; }
1115 14895 : template <> FEContinuity FE<1,LAGRANGE>::get_continuity() const { return C_ZERO; }
1116 333351 : template <> FEContinuity FE<2,LAGRANGE>::get_continuity() const { return C_ZERO; }
1117 133253 : template <> FEContinuity FE<3,LAGRANGE>::get_continuity() const { return C_ZERO; }
1118 :
1119 : // Lagrange FEMs are not hierarchic
1120 0 : template <> bool FE<0,LAGRANGE>::is_hierarchic() const { return false; }
1121 24 : template <> bool FE<1,LAGRANGE>::is_hierarchic() const { return false; }
1122 180 : template <> bool FE<2,LAGRANGE>::is_hierarchic() const { return false; }
1123 656 : template <> bool FE<3,LAGRANGE>::is_hierarchic() const { return false; }
1124 :
1125 : // Lagrange FEM shapes do not need reinit (is this always true?)
1126 3084 : template <> bool FE<0,LAGRANGE>::shapes_need_reinit() const { return false; }
1127 30254 : template <> bool FE<1,LAGRANGE>::shapes_need_reinit() const { return false; }
1128 137062338 : template <> bool FE<2,LAGRANGE>::shapes_need_reinit() const { return false; }
1129 5867132 : template <> bool FE<3,LAGRANGE>::shapes_need_reinit() const { return false; }
1130 :
1131 : // Methods for computing Lagrange constraints. Note: we pass the
1132 : // dimension as the last argument to the anonymous helper function.
1133 : // Also note: we only need instantiations of this function for
1134 : // Dim==2 and 3.
1135 : #ifdef LIBMESH_ENABLE_AMR
1136 : template <>
1137 3159585 : void FE<2,LAGRANGE>::compute_constraints (DofConstraints & constraints,
1138 : DofMap & dof_map,
1139 : const unsigned int variable_number,
1140 : const Elem * elem)
1141 3159585 : { lagrange_compute_constraints(constraints, dof_map, variable_number, elem, /*Dim=*/2); }
1142 :
1143 : template <>
1144 1183922 : void FE<3,LAGRANGE>::compute_constraints (DofConstraints & constraints,
1145 : DofMap & dof_map,
1146 : const unsigned int variable_number,
1147 : const Elem * elem)
1148 1183922 : { lagrange_compute_constraints(constraints, dof_map, variable_number, elem, /*Dim=*/3); }
1149 : #endif // LIBMESH_ENABLE_AMR
1150 :
1151 : } // namespace libMesh
|
