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 :
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 8613254 : 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 8613254 : const unsigned int n_nodes = elem->n_nodes();
50 8613254 : const ElemType type = elem->type();
51 :
52 9379633 : const Order totalorder = order + add_p_level*elem->p_level();
53 :
54 8613254 : nodal_soln.resize(n_nodes);
55 :
56 :
57 :
58 8613254 : switch (totalorder)
59 : {
60 : // linear Lagrange shape functions
61 7714981 : case FIRST:
62 : {
63 7714981 : 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 1380 : case TRI7:
92 462 : libmesh_assert_equal_to (nodal_soln.size(), 7);
93 5532 : nodal_soln[6] = (elem_soln[0] + elem_soln[1] + elem_soln[2])/3.;
94 : libmesh_fallthrough();
95 15384 : case TRI6:
96 : {
97 1316 : libmesh_assert (type == TRI7 || nodal_soln.size() == 6);
98 1316 : libmesh_assert_equal_to (elem_soln.size(), 3);
99 :
100 15384 : nodal_soln[0] = elem_soln[0];
101 15384 : nodal_soln[1] = elem_soln[1];
102 15384 : nodal_soln[2] = elem_soln[2];
103 15384 : nodal_soln[3] = .5*(elem_soln[0] + elem_soln[1]);
104 15384 : nodal_soln[4] = .5*(elem_soln[1] + elem_soln[2]);
105 15384 : nodal_soln[5] = .5*(elem_soln[2] + elem_soln[0]);
106 :
107 15384 : 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 7386216 : default:
301 : {
302 : // By default the element solution _is_ nodal,
303 : // so just copy it.
304 7386216 : nodal_soln = elem_soln;
305 :
306 7386216 : return;
307 : }
308 : }
309 : }
310 :
311 890909 : case SECOND:
312 : {
313 890909 : 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 474 : case TRI7:
331 : {
332 474 : libmesh_assert_equal_to (elem_soln.size(), 6);
333 474 : libmesh_assert_equal_to (nodal_soln.size(), 7);
334 :
335 39648 : for (int i=0; i != 6; ++i)
336 39672 : nodal_soln[i] = elem_soln[i];
337 :
338 6138 : nodal_soln[6] = -1./9. * (elem_soln[0] + elem_soln[1] + elem_soln[2])
339 5202 : +4./9. * (elem_soln[3] + elem_soln[4] + elem_soln[5]);
340 :
341 5664 : return;
342 : }
343 :
344 1556 : case TET14:
345 : {
346 1556 : libmesh_assert_equal_to (elem_soln.size(), 10);
347 1556 : libmesh_assert_equal_to (nodal_soln.size(), 14);
348 :
349 210232 : for (int i=0; i != 10; ++i)
350 222240 : nodal_soln[i] = elem_soln[i];
351 :
352 20668 : nodal_soln[10] = -1./9. * (elem_soln[0] + elem_soln[1] + elem_soln[2])
353 17566 : +4./9. * (elem_soln[4] + elem_soln[5] + elem_soln[6]);
354 20668 : nodal_soln[11] = -1./9. * (elem_soln[0] + elem_soln[1] + elem_soln[3])
355 17566 : +4./9. * (elem_soln[4] + elem_soln[7] + elem_soln[8]);
356 20668 : nodal_soln[12] = -1./9. * (elem_soln[1] + elem_soln[2] + elem_soln[3])
357 17566 : +4./9. * (elem_soln[5] + elem_soln[8] + elem_soln[9]);
358 20668 : nodal_soln[13] = -1./9. * (elem_soln[0] + elem_soln[2] + elem_soln[3])
359 17566 : +4./9. * (elem_soln[6] + elem_soln[7] + elem_soln[9]);
360 :
361 19112 : return;
362 : }
363 :
364 548 : case PRISM21:
365 : {
366 3764 : nodal_soln[20] = -1./9. * (elem_soln[9] + elem_soln[10] + elem_soln[11])
367 2978 : +4./9. * (elem_soln[15] + elem_soln[16] + elem_soln[17]);
368 : libmesh_fallthrough();
369 : }
370 346 : case PRISM20:
371 : {
372 346 : libmesh_assert_equal_to (elem_soln.size(), 18);
373 346 : if (type == PRISM20)
374 72 : libmesh_assert_equal_to (nodal_soln.size(), 20);
375 :
376 76152 : for (int i=0; i != 18; ++i)
377 84600 : nodal_soln[i] = elem_soln[i];
378 :
379 4354 : nodal_soln[18] = -1./9. * (elem_soln[0] + elem_soln[1] + elem_soln[2])
380 3734 : +4./9. * (elem_soln[6] + elem_soln[7] + elem_soln[8]);
381 4354 : nodal_soln[19] = -1./9. * (elem_soln[3] + elem_soln[4] + elem_soln[5])
382 3734 : +4./9. * (elem_soln[12] + elem_soln[13] + elem_soln[14]);
383 4008 : return;
384 : }
385 :
386 28 : case PYRAMID18:
387 : {
388 28 : libmesh_assert_equal_to (elem_soln.size(), 14);
389 28 : libmesh_assert_equal_to (nodal_soln.size(), 18);
390 :
391 9660 : for (int i=0; i != 14; ++i)
392 9800 : nodal_soln[i] = elem_soln[i];
393 :
394 672 : nodal_soln[14] = -1./9. * (elem_soln[0] + elem_soln[1] + elem_soln[4])
395 630 : +4./9. * (elem_soln[5] + elem_soln[9] + elem_soln[10]);
396 672 : nodal_soln[15] = -1./9. * (elem_soln[1] + elem_soln[2] + elem_soln[4])
397 630 : +4./9. * (elem_soln[6] + elem_soln[10] + elem_soln[11]);
398 672 : nodal_soln[16] = -1./9. * (elem_soln[2] + elem_soln[3] + elem_soln[4])
399 630 : +4./9. * (elem_soln[7] + elem_soln[11] + elem_soln[12]);
400 672 : nodal_soln[17] = -1./9. * (elem_soln[0] + elem_soln[3] + elem_soln[4])
401 630 : +4./9. * (elem_soln[8] + elem_soln[9] + elem_soln[12]);
402 :
403 644 : return;
404 : }
405 :
406 76186 : default:
407 : {
408 : // By default the element solution _is_ nodal, so just
409 : // copy the portion relevant to the nodal solution.
410 : // (this is the whole nodal solution for true Lagrange
411 : // elements, but a smaller part for "L2 Lagrange"
412 76186 : libmesh_assert_less_equal(nodal_soln.size(), elem_soln.size());
413 8647024 : for (const auto i : index_range(nodal_soln))
414 9161423 : nodal_soln[i] = elem_soln[i];
415 :
416 152372 : return;
417 : }
418 : }
419 : }
420 :
421 :
422 :
423 :
424 494 : default:
425 : {
426 : // By default the element solution _is_ nodal, so just copy
427 : // the portion relevant to the nodal solution. (this is the
428 : // whole nodal solution for true Lagrange elements, but a
429 : // smaller part for "L2 Lagrange"
430 494 : libmesh_assert_less_equal(nodal_soln.size(), elem_soln.size());
431 134956 : for (const auto i : index_range(nodal_soln))
432 145388 : nodal_soln[i] = elem_soln[i];
433 :
434 988 : return;
435 : }
436 : }
437 : }
438 :
439 :
440 : // Anonymous namespace for local helper functions
441 : namespace {
442 1606823016 : unsigned int lagrange_n_dofs(const ElemType t, const Elem * e, const Order o)
443 : {
444 127635051 : libmesh_assert(!e || e->type() == t);
445 :
446 1606823016 : switch (o)
447 : {
448 : // lagrange can only be constant on a single node
449 147550 : case CONSTANT:
450 : {
451 147550 : switch (t)
452 : {
453 12196 : case NODEELEM:
454 12196 : return 1;
455 :
456 0 : default:
457 0 : libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
458 : }
459 : }
460 :
461 : // linear Lagrange shape functions
462 297940106 : case FIRST:
463 : {
464 297940106 : switch (t)
465 : {
466 1352 : case NODEELEM:
467 1352 : return 1;
468 :
469 14591527 : case EDGE2:
470 : case EDGE3:
471 : case EDGE4:
472 14591527 : return 2;
473 :
474 6634635 : case TRI3:
475 : case TRISHELL3:
476 : case TRI3SUBDIVISION:
477 : case TRI6:
478 : case TRI7:
479 6634635 : return 3;
480 :
481 23581175 : case QUAD4:
482 : case QUADSHELL4:
483 : case QUAD8:
484 : case QUADSHELL8:
485 : case QUAD9:
486 : case QUADSHELL9:
487 23581175 : return 4;
488 :
489 1676687 : case TET4:
490 : case TET10:
491 : case TET14:
492 1676687 : return 4;
493 :
494 17330426 : case HEX8:
495 : case HEX20:
496 : case HEX27:
497 17330426 : return 8;
498 :
499 1819462 : case PRISM6:
500 : case PRISM15:
501 : case PRISM18:
502 : case PRISM20:
503 : case PRISM21:
504 1819462 : return 6;
505 :
506 1889444 : case PYRAMID5:
507 : case PYRAMID13:
508 : case PYRAMID14:
509 : case PYRAMID18:
510 1889444 : return 5;
511 :
512 0 : case INVALID_ELEM:
513 0 : return 0;
514 :
515 6223964 : case C0POLYGON:
516 : case C0POLYHEDRON:
517 : // Polygons and polyhedra require using newer FE APIs
518 6223964 : if (!e)
519 0 : libmesh_error_msg("Code (see stack trace) used an outdated FE function overload.\n"
520 1058150 : "n_dofs() on a polygon or polyhedron is not defined by ElemType alone.");
521 6223964 : return e->n_nodes();
522 :
523 0 : default:
524 0 : libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
525 : }
526 : }
527 :
528 :
529 : // quadratic Lagrange shape functions
530 942027203 : case SECOND:
531 : {
532 942027203 : switch (t)
533 : {
534 195 : case NODEELEM:
535 195 : return 1;
536 :
537 6340729 : case EDGE3:
538 6340729 : return 3;
539 :
540 56011341 : case TRI6:
541 : case TRI7:
542 56011341 : return 6;
543 :
544 1376812 : case QUAD8:
545 : case QUADSHELL8:
546 1376812 : return 8;
547 :
548 31339842 : case QUAD9:
549 : case QUADSHELL9:
550 31339842 : return 9;
551 :
552 154526537 : case TET10:
553 : case TET14:
554 154526537 : return 10;
555 :
556 1977126 : case HEX20:
557 1977126 : return 20;
558 :
559 128449695 : case HEX27:
560 128449695 : return 27;
561 :
562 13615960 : case PRISM15:
563 13615960 : return 15;
564 :
565 3157577 : case PRISM18:
566 : case PRISM20:
567 : case PRISM21:
568 3157577 : return 18;
569 :
570 11529006 : case PYRAMID13:
571 11529006 : return 13;
572 :
573 681804 : case PYRAMID14:
574 : case PYRAMID18:
575 681804 : return 14;
576 :
577 0 : case INVALID_ELEM:
578 0 : return 0;
579 :
580 0 : default:
581 0 : libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
582 : }
583 : }
584 :
585 366708157 : case THIRD:
586 : {
587 366708157 : switch (t)
588 : {
589 0 : case NODEELEM:
590 0 : return 1;
591 :
592 5976 : case EDGE4:
593 5976 : return 4;
594 :
595 2351455 : case PRISM20:
596 2351455 : return 20;
597 :
598 3342625 : case PRISM21:
599 3342625 : return 21;
600 :
601 1267451 : case PYRAMID18:
602 1267451 : return 18;
603 :
604 84003293 : case TRI7:
605 84003293 : return 7;
606 :
607 14510843 : case TET14:
608 14510843 : return 14;
609 :
610 0 : case INVALID_ELEM:
611 0 : return 0;
612 :
613 0 : default:
614 0 : libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
615 : }
616 : }
617 :
618 0 : default:
619 0 : libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for LAGRANGE FE family!");
620 : }
621 : }
622 :
623 :
624 :
625 1167007868 : unsigned int lagrange_n_dofs_at_node(const ElemType t,
626 : const Order o,
627 : const unsigned int n)
628 : {
629 1167007868 : switch (o)
630 : {
631 : // lagrange can only be constant on a single node
632 0 : case CONSTANT:
633 : {
634 0 : switch (t)
635 : {
636 0 : case NODEELEM:
637 0 : return 1;
638 :
639 0 : default:
640 0 : libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
641 : }
642 : }
643 :
644 : // linear Lagrange shape functions
645 853483273 : case FIRST:
646 : {
647 853483273 : switch (t)
648 : {
649 1648 : case NODEELEM:
650 1648 : return 1;
651 :
652 26853902 : case EDGE2:
653 : case EDGE3:
654 : case EDGE4:
655 : {
656 26853902 : switch (n)
657 : {
658 1314585 : case 0:
659 : case 1:
660 1314585 : return 1;
661 :
662 1208438 : default:
663 1208438 : return 0;
664 : }
665 : }
666 :
667 74930166 : case TRI3:
668 : case TRISHELL3:
669 : case TRI3SUBDIVISION:
670 : case TRI6:
671 : case TRI7:
672 : {
673 74930166 : switch (n)
674 : {
675 4808583 : case 0:
676 : case 1:
677 : case 2:
678 4808583 : return 1;
679 :
680 20897 : default:
681 20897 : return 0;
682 : }
683 : }
684 :
685 492903837 : case QUAD4:
686 : case QUADSHELL4:
687 : case QUAD8:
688 : case QUADSHELL8:
689 : case QUAD9:
690 : case QUADSHELL9:
691 : {
692 492903837 : switch (n)
693 : {
694 42854626 : case 0:
695 : case 1:
696 : case 2:
697 : case 3:
698 42854626 : return 1;
699 :
700 6336891 : default:
701 6336891 : return 0;
702 : }
703 : }
704 :
705 3354 : case C0POLYGON:
706 : case C0POLYHEDRON:
707 3354 : return 1;
708 :
709 :
710 83978868 : case TET4:
711 : case TET10:
712 : case TET14:
713 : {
714 83978868 : switch (n)
715 : {
716 5021973 : case 0:
717 : case 1:
718 : case 2:
719 : case 3:
720 5021973 : return 1;
721 :
722 0 : default:
723 0 : return 0;
724 : }
725 : }
726 :
727 166699306 : case HEX8:
728 : case HEX20:
729 : case HEX27:
730 : {
731 166699306 : switch (n)
732 : {
733 11995643 : case 0:
734 : case 1:
735 : case 2:
736 : case 3:
737 : case 4:
738 : case 5:
739 : case 6:
740 : case 7:
741 11995643 : return 1;
742 :
743 2577436 : default:
744 2577436 : return 0;
745 : }
746 : }
747 :
748 4514051 : case PRISM6:
749 : case PRISM15:
750 : case PRISM18:
751 : case PRISM20:
752 : case PRISM21:
753 : {
754 4514051 : switch (n)
755 : {
756 325171 : case 0:
757 : case 1:
758 : case 2:
759 : case 3:
760 : case 4:
761 : case 5:
762 325171 : return 1;
763 :
764 0 : default:
765 0 : return 0;
766 : }
767 : }
768 :
769 3529224 : case PYRAMID5:
770 : case PYRAMID13:
771 : case PYRAMID14:
772 : case PYRAMID18:
773 : {
774 3529224 : switch (n)
775 : {
776 284414 : case 0:
777 : case 1:
778 : case 2:
779 : case 3:
780 : case 4:
781 284414 : return 1;
782 :
783 0 : default:
784 0 : return 0;
785 : }
786 : }
787 :
788 0 : case INVALID_ELEM:
789 0 : return 0;
790 :
791 0 : default:
792 0 : libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
793 : }
794 : }
795 :
796 : // quadratic Lagrange shape functions
797 302923934 : case SECOND:
798 : {
799 302923934 : switch (t)
800 : {
801 : // quadratic lagrange has one dof at each node
802 25773162 : case NODEELEM:
803 : case EDGE3:
804 : case TRI6:
805 : case QUAD8:
806 : case QUADSHELL8:
807 : case QUAD9:
808 : case QUADSHELL9:
809 : case TET10:
810 : case HEX20:
811 : case HEX27:
812 : case PRISM15:
813 : case PRISM18:
814 : case PYRAMID13:
815 : case PYRAMID14:
816 25773162 : return 1;
817 :
818 3815861 : case PRISM20:
819 : case PRISM21:
820 3815861 : return (n < 18);
821 :
822 752778 : case PYRAMID18:
823 752778 : return (n < 14);
824 :
825 41878 : case TRI7:
826 41878 : return (n < 6);
827 :
828 1130077 : case TET14:
829 1130077 : return (n < 10);
830 :
831 0 : case INVALID_ELEM:
832 0 : return 0;
833 :
834 0 : default:
835 0 : libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
836 : }
837 : }
838 :
839 10600661 : case THIRD:
840 : {
841 9859146 : switch (t)
842 : {
843 741025 : case NODEELEM:
844 : case EDGE4:
845 : case PRISM20:
846 : case PRISM21:
847 : case PYRAMID18:
848 : case TRI7:
849 : case TET14:
850 741025 : return 1;
851 :
852 0 : case INVALID_ELEM:
853 0 : return 0;
854 :
855 0 : default:
856 0 : libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
857 : }
858 : }
859 :
860 0 : default:
861 0 : libmesh_error_msg("Unsupported order: " << o );
862 : }
863 : }
864 :
865 :
866 :
867 90907129 : unsigned int lagrange_n_dofs_at_node(const Elem & e,
868 : const Order o,
869 : const unsigned int n)
870 : {
871 97563219 : return lagrange_n_dofs_at_node(e.type(), o, n);
872 : }
873 :
874 :
875 : #ifdef LIBMESH_ENABLE_AMR
876 4519284 : void lagrange_compute_constraints (DofConstraints & constraints,
877 : DofMap & dof_map,
878 : const unsigned int variable_number,
879 : const Elem * elem,
880 : const unsigned Dim)
881 : {
882 : // Only constrain elements in 2,3D.
883 4519284 : if (Dim == 1)
884 380840 : return;
885 :
886 422476 : libmesh_assert(elem);
887 :
888 : // Only constrain active and ancestor elements
889 4519284 : if (elem->subactive())
890 36396 : return;
891 :
892 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
893 1075098 : if (elem->infinite())
894 : {
895 3212 : const FEType fe_t = dof_map.variable_type(variable_number);
896 :
897 : // expand the infinite_compute_constraint in its template-arguments.
898 3212 : switch(Dim)
899 : {
900 0 : case 2:
901 : {
902 0 : inf_fe_family_mapping_switch(2, inf_compute_constraints (constraints, dof_map, variable_number, elem) , ,; break;);
903 0 : break;
904 : }
905 3212 : case 3:
906 : {
907 3212 : inf_fe_family_mapping_switch(3, inf_compute_constraints (constraints, dof_map, variable_number, elem) , ,; break;);
908 1056 : break;
909 : }
910 0 : default:
911 0 : libmesh_error_msg("Invalid dim = " << Dim);
912 : }
913 1056 : return;
914 : }
915 : #endif
916 :
917 4138444 : const Variable & var = dof_map.variable(variable_number);
918 4138444 : const FEType fe_type = var.type();
919 4138444 : const bool add_p_level = fe_type.p_refinement;
920 :
921 : // Pull objects out of the loop to reduce heap operations
922 770048 : std::vector<dof_id_type> my_dof_indices, parent_dof_indices;
923 4440282 : std::unique_ptr<const Elem> my_side, parent_side;
924 :
925 : // Look at the element faces. Check to see if we need to
926 : // build constraints.
927 21312890 : for (auto s : elem->side_index_range())
928 18767531 : if (const Elem * const neigh = elem->neighbor_ptr(s);
929 16110865 : neigh && neigh != remote_elem)
930 : {
931 : // constrain dofs shared between this element and ones coarser
932 16141919 : if (neigh->level() < elem->level() &&
933 318460 : var.active_on_subdomain(neigh->subdomain_id()))
934 : {
935 : // than this element.
936 : // Get pointers to the elements of interest and its parent.
937 72036 : const Elem * parent = elem->parent();
938 :
939 : // This can't happen... Only level-0 elements have nullptr
940 : // parents, and no level-0 elements can be at a higher
941 : // level than their neighbors!
942 36018 : libmesh_assert(parent);
943 :
944 318118 : elem->build_side_ptr(my_side, s);
945 318118 : parent->build_side_ptr(parent_side, s);
946 :
947 : // We have some elements with interior bubble function bases
948 : // and lower-order sides. We need to only ask for the lower
949 : // order in those cases.
950 318118 : FEType side_fe_type = fe_type;
951 36018 : int extra_order = 0;
952 318118 : if (my_side->default_order() < fe_type.order)
953 : {
954 1104 : side_fe_type.order = my_side->default_order();
955 1196 : extra_order = (int)(side_fe_type.order) -
956 92 : (int)(fe_type.order);
957 : }
958 :
959 : // This function gets called element-by-element, so there
960 : // will be a lot of memory allocation going on. We can
961 : // at least minimize this for the case of the dof indices
962 : // by efficiently preallocating the requisite storage.
963 318118 : my_dof_indices.reserve (my_side->n_nodes());
964 318118 : parent_dof_indices.reserve (parent_side->n_nodes());
965 :
966 318118 : dof_map.dof_indices (my_side.get(), my_dof_indices,
967 72036 : variable_number, extra_order);
968 318118 : dof_map.dof_indices (parent_side.get(), parent_dof_indices,
969 72036 : variable_number, extra_order);
970 :
971 : const unsigned int n_side_dofs =
972 318118 : FEInterface::n_dofs(side_fe_type, my_side.get());
973 : const unsigned int n_parent_side_dofs =
974 318118 : FEInterface::n_dofs(side_fe_type, parent_side.get());
975 1165482 : for (unsigned int my_dof=0; my_dof != n_side_dofs; my_dof++)
976 : {
977 93342 : libmesh_assert_less (my_dof, my_side->n_nodes());
978 :
979 : // My global dof index.
980 940706 : const dof_id_type my_dof_g = my_dof_indices[my_dof];
981 :
982 : // Hunt for "constraining against myself" cases before
983 : // we bother creating a constraint row
984 93342 : bool self_constraint = false;
985 2740391 : for (unsigned int their_dof=0;
986 2833733 : their_dof != n_parent_side_dofs; their_dof++)
987 : {
988 254803 : libmesh_assert_less (their_dof, parent_side->n_nodes());
989 :
990 : // Their global dof index.
991 : const dof_id_type their_dof_g =
992 2341079 : parent_dof_indices[their_dof];
993 :
994 2341079 : if (their_dof_g == my_dof_g)
995 : {
996 40578 : self_constraint = true;
997 40578 : break;
998 : }
999 : }
1000 :
1001 847364 : if (self_constraint)
1002 634219 : continue;
1003 :
1004 : DofConstraintRow * constraint_row;
1005 :
1006 : // we may be running constraint methods concurrently
1007 : // on multiple threads, so we need a lock to
1008 : // ensure that this constraint is "ours"
1009 : {
1010 52764 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1011 :
1012 492654 : if (dof_map.is_constrained_dof(my_dof_g))
1013 29643 : continue;
1014 :
1015 213145 : constraint_row = &(constraints[my_dof_g]);
1016 23121 : libmesh_assert(constraint_row->empty());
1017 : }
1018 :
1019 : // The support point of the DOF
1020 46242 : const Point & support_point = my_side->point(my_dof);
1021 :
1022 : // Figure out where my node lies on their reference element.
1023 166903 : const Point mapped_point = FEMap::inverse_map(Dim-1,
1024 : parent_side.get(),
1025 46242 : support_point);
1026 :
1027 : // Compute the parent's side shape function values.
1028 700991 : for (unsigned int their_dof=0;
1029 891015 : their_dof != n_parent_side_dofs; their_dof++)
1030 : {
1031 72838 : libmesh_assert_less (their_dof, parent_side->n_nodes());
1032 :
1033 : // Their global dof index.
1034 : const dof_id_type their_dof_g =
1035 750708 : parent_dof_indices[their_dof];
1036 :
1037 677870 : const Real their_dof_value = FEInterface::shape(side_fe_type,
1038 : parent_side.get(),
1039 : their_dof,
1040 677870 : mapped_point);
1041 :
1042 : // Only add non-zero and non-identity values
1043 : // for Lagrange basis functions.
1044 735548 : if ((std::abs(their_dof_value) > 1.e-5) &&
1045 57678 : (std::abs(their_dof_value) < .999))
1046 : {
1047 57678 : constraint_row->emplace(their_dof_g, their_dof_value);
1048 : }
1049 : #ifdef DEBUG
1050 : // Protect for the case u_i = 0.999 u_j,
1051 : // in which case i better equal j.
1052 15160 : else if (their_dof_value >= .999)
1053 0 : libmesh_assert_equal_to (my_dof_g, their_dof_g);
1054 : #endif
1055 : }
1056 : }
1057 : }
1058 :
1059 14233820 : if (elem->active() && add_p_level)
1060 : {
1061 : // p refinement constraints:
1062 : // constrain dofs shared between
1063 : // active elements and neighbors with
1064 : // lower polynomial degrees
1065 : const unsigned int min_p_level =
1066 15301446 : neigh->min_p_level_by_neighbor(elem, elem->p_level());
1067 15301446 : if (min_p_level < elem->p_level())
1068 0 : libmesh_error_msg(
1069 : "Mismatched p-refinement levels for LAGRANGE is not currently supported");
1070 : }
1071 : }
1072 3066558 : } // lagrange_compute_constraints()
1073 : #endif // #ifdef LIBMESH_ENABLE_AMR
1074 :
1075 : } // anonymous namespace
1076 :
1077 :
1078 : // Instantiate (side_) nodal_soln() function for every dimension
1079 8492345 : LIBMESH_FE_NODAL_SOLN(LAGRANGE, lagrange_nodal_soln)
1080 1584 : LIBMESH_FE_SIDE_NODAL_SOLN(LAGRANGE)
1081 :
1082 :
1083 : // Do full-specialization for every dimension, instead
1084 : // of explicit instantiation at the end of this function.
1085 : // This could be macro-ified.
1086 0 : template <> unsigned int FE<0,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, nullptr, o); }
1087 0 : template <> unsigned int FE<1,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, nullptr, o); }
1088 1587 : template <> unsigned int FE<2,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, nullptr, o); }
1089 0 : template <> unsigned int FE<3,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, nullptr, o); }
1090 :
1091 151913 : template <> unsigned int FE<0,LAGRANGE>::n_dofs(const Elem * e, const Order o) { return lagrange_n_dofs(e->type(), e, o); }
1092 116448380 : template <> unsigned int FE<1,LAGRANGE>::n_dofs(const Elem * e, const Order o) { return lagrange_n_dofs(e->type(), e, o); }
1093 734849773 : template <> unsigned int FE<2,LAGRANGE>::n_dofs(const Elem * e, const Order o) { return lagrange_n_dofs(e->type(), e, o); }
1094 755371363 : template <> unsigned int FE<3,LAGRANGE>::n_dofs(const Elem * e, const Order o) { return lagrange_n_dofs(e->type(), e, o); }
1095 :
1096 :
1097 : // Do full-specialization for every dimension, instead
1098 : // of explicit instantiation at the end of this function.
1099 21135 : 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); }
1100 27014673 : 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); }
1101 776266173 : 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); }
1102 266142668 : 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); }
1103 :
1104 23006 : 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); }
1105 97633 : 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); }
1106 67383137 : 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); }
1107 30059443 : 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); }
1108 :
1109 :
1110 : // Lagrange elements have no dofs per element
1111 : // (just at the nodes)
1112 0 : template <> unsigned int FE<0,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1113 0 : template <> unsigned int FE<1,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1114 0 : template <> unsigned int FE<2,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1115 0 : template <> unsigned int FE<3,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
1116 :
1117 32638 : template <> unsigned int FE<0,LAGRANGE>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1118 9172762 : template <> unsigned int FE<1,LAGRANGE>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1119 166726922 : template <> unsigned int FE<2,LAGRANGE>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1120 40926457 : template <> unsigned int FE<3,LAGRANGE>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
1121 :
1122 : // Lagrange FEMs are always C^0 continuous
1123 0 : template <> FEContinuity FE<0,LAGRANGE>::get_continuity() const { return C_ZERO; }
1124 17414 : template <> FEContinuity FE<1,LAGRANGE>::get_continuity() const { return C_ZERO; }
1125 338313 : template <> FEContinuity FE<2,LAGRANGE>::get_continuity() const { return C_ZERO; }
1126 161845 : template <> FEContinuity FE<3,LAGRANGE>::get_continuity() const { return C_ZERO; }
1127 :
1128 : // Lagrange FEMs are not hierarchic
1129 0 : template <> bool FE<0,LAGRANGE>::is_hierarchic() const { return false; }
1130 24 : template <> bool FE<1,LAGRANGE>::is_hierarchic() const { return false; }
1131 180 : template <> bool FE<2,LAGRANGE>::is_hierarchic() const { return false; }
1132 702 : template <> bool FE<3,LAGRANGE>::is_hierarchic() const { return false; }
1133 :
1134 : // Lagrange FEM shapes do not need reinit (is this always true?)
1135 3084 : template <> bool FE<0,LAGRANGE>::shapes_need_reinit() const { return false; }
1136 30479 : template <> bool FE<1,LAGRANGE>::shapes_need_reinit() const { return false; }
1137 137189902 : template <> bool FE<2,LAGRANGE>::shapes_need_reinit() const { return false; }
1138 7298177 : template <> bool FE<3,LAGRANGE>::shapes_need_reinit() const { return false; }
1139 :
1140 : // Methods for computing Lagrange constraints. Note: we pass the
1141 : // dimension as the last argument to the anonymous helper function.
1142 : // Also note: we only need instantiations of this function for
1143 : // Dim==2 and 3.
1144 : #ifdef LIBMESH_ENABLE_AMR
1145 : template <>
1146 3209263 : void FE<2,LAGRANGE>::compute_constraints (DofConstraints & constraints,
1147 : DofMap & dof_map,
1148 : const unsigned int variable_number,
1149 : const Elem * elem)
1150 3209263 : { lagrange_compute_constraints(constraints, dof_map, variable_number, elem, /*Dim=*/2); }
1151 :
1152 : template <>
1153 1310021 : void FE<3,LAGRANGE>::compute_constraints (DofConstraints & constraints,
1154 : DofMap & dof_map,
1155 : const unsigned int variable_number,
1156 : const Elem * elem)
1157 1310021 : { lagrange_compute_constraints(constraints, dof_map, variable_number, elem, /*Dim=*/3); }
1158 : #endif // LIBMESH_ENABLE_AMR
1159 :
1160 : } // namespace libMesh
|
