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 : // Local includes
20 : #include "libmesh/fe.h"
21 : #include "libmesh/elem.h"
22 : #include "libmesh/number_lookups.h"
23 : #include "libmesh/enum_to_string.h"
24 : #include "libmesh/cell_tet4.h" // We need edge_nodes_map
25 : #include "libmesh/face_tri3.h" // Faster to construct these on the stack
26 : #include "libmesh/face_quad4.h"
27 :
28 : // Anonymous namespace for functions shared by HIERARCHIC and
29 : // L2_HIERARCHIC implementations. Implementations appear at the bottom
30 : // of this file.
31 : namespace
32 : {
33 : using namespace libMesh;
34 :
35 : unsigned int cube_side(const Point & p);
36 :
37 : Point cube_side_point(unsigned int sidenum, const Point & interior_point);
38 :
39 : template <FEFamily T>
40 : Real fe_hierarchic_3D_shape(const Elem * elem,
41 : const Order order,
42 : const unsigned int i,
43 : const Point & p,
44 : const bool add_p_level);
45 :
46 : template <FEFamily T>
47 : Real fe_hierarchic_3D_shape_deriv(const Elem * elem,
48 : const Order order,
49 : const unsigned int i,
50 : const unsigned int j,
51 : const Point & p,
52 : const bool add_p_level);
53 :
54 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
55 :
56 : template <FEFamily T>
57 : Real fe_hierarchic_3D_shape_second_deriv(const Elem * elem,
58 : const Order order,
59 : const unsigned int i,
60 : const unsigned int j,
61 : const Point & p,
62 : const bool add_p_level);
63 :
64 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
65 :
66 : #if LIBMESH_DIM > 2
67 : // Remap non-face-nodes based on point ordering
68 : template <unsigned int N_nodes>
69 60725016 : unsigned int remap_node(unsigned int n,
70 : const Elem & elem,
71 : unsigned int nodebegin)
72 : {
73 : std::array<const Point *, N_nodes> points;
74 :
75 303625080 : for (auto i : IntRange<unsigned int>(0, N_nodes))
76 263088224 : points[i] = &elem.point(nodebegin+i);
77 :
78 5047040 : std::sort(points.begin(), points.end(),
79 28448840 : [](const Point * a, const Point * b)
80 313884136 : { return *a < *b; });
81 :
82 60725016 : const Point * pn = points[n-nodebegin];
83 :
84 151812540 : for (auto i : IntRange<unsigned int>(nodebegin, nodebegin+N_nodes))
85 164430140 : if (pn == &elem.point(i))
86 5047040 : return i;
87 :
88 0 : libmesh_assert(false);
89 0 : return libMesh::invalid_uint;
90 : }
91 :
92 :
93 88992834 : void cube_remap(unsigned int & side_i,
94 : const Elem & elem,
95 : unsigned int sidenum,
96 : const Elem & side,
97 : unsigned int totalorder,
98 : Point & sidep)
99 : {
100 : // "vertex" nodes are now decoupled from vertices, so we have
101 : // to order them consistently otherwise
102 88992834 : if (side_i < 4)
103 17967096 : side_i = remap_node<4>(side_i, side, 0);
104 :
105 : // And "edge" nodes are decoupled from edges, so we have to
106 : // reorder them too!
107 71025738 : else if (side_i < 4u*totalorder)
108 : {
109 42757920 : unsigned int side_node = (side_i - 4)/(totalorder-1)+4;
110 42757920 : side_node = remap_node<4>(side_node, side, 4);
111 49867040 : side_i = ((side_i - 4) % (totalorder - 1)) // old local edge_i
112 42757920 : + 4 + (side_node-4)*(totalorder-1);
113 : }
114 :
115 : // Interior dofs in 2D don't care about where xi/eta point in
116 : // physical space, but here we need them to match from both
117 : // sides of a face!
118 : else
119 : {
120 : // The side numbers its nodes as the side node map of elem does.
121 : // Rotate the least node to the origin, then flip about the
122 : // diagonal through it if its lesser neighbor is its predecessor
123 28267818 : const unsigned int orientation = elem.face_orientation(sidenum);
124 28267818 : const unsigned int min_side_node = orientation / 2;
125 28267818 : const bool flip = orientation % 2;
126 :
127 28267818 : switch (min_side_node) {
128 3752036 : case 0:
129 3752036 : if (flip)
130 157319 : std::swap(sidep(0), sidep(1));
131 312044 : break;
132 6531006 : case 1:
133 6531006 : sidep(0) = -sidep(0);
134 6531006 : if (!flip)
135 250609 : std::swap(sidep(0), sidep(1));
136 543947 : break;
137 7038408 : case 2:
138 7038408 : sidep(0) = -sidep(0);
139 7038408 : sidep(1) = -sidep(1);
140 7038408 : if (flip)
141 197284 : std::swap(sidep(0), sidep(1));
142 585520 : break;
143 10946368 : case 3:
144 10946368 : sidep(1) = -sidep(1);
145 10946368 : if (!flip)
146 387266 : std::swap(sidep(0), sidep(1));
147 909033 : break;
148 0 : default:
149 0 : libmesh_error();
150 : }
151 : }
152 88992834 : }
153 :
154 :
155 974545994 : void cube_indices(const Elem * elem,
156 : const unsigned int totalorder,
157 : const unsigned int i,
158 : Real & xi, Real & eta, Real & zeta,
159 : unsigned int & i0,
160 : unsigned int & i1,
161 : unsigned int & i2)
162 : {
163 : // The only way to make any sense of this
164 : // is to look at the mgflo/mg2/mgf documentation
165 : // and make the cut-out cube!
166 : // Example i0 and i1 values for totalorder = 3:
167 : // FIXME - these examples are incorrect now that we've got truly
168 : // hierarchic basis functions
169 : // Nodes 0 1 2 3 4 5 6 7 8 8 9 9 10 10 11 11 12 12 13 13 14 14 15 15 16 16 17 17 18 18 19 19 20 20 20 20 21 21 21 21 22 22 22 22 23 23 23 23 24 24 24 24 25 25 25 25 26 26 26 26 26 26 26 26
170 : // DOFS 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 18 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 60 62 63
171 : // static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 1, 1, 1, 1, 2, 3, 2, 3, 0, 0, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3};
172 : // static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3, 0, 0, 0, 0, 2, 3, 2, 3, 1, 1, 1, 1, 2, 3, 2, 3, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3};
173 : // static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3};
174 :
175 : // the number of DoFs per edge appears everywhere:
176 974545994 : const unsigned int e = totalorder - 1u;
177 :
178 78858204 : libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u)*(totalorder+1u));
179 :
180 974545994 : Real xi_saved = xi, eta_saved = eta, zeta_saved = zeta;
181 :
182 : // Vertices:
183 974545994 : if (i == 0)
184 : {
185 18358646 : i0 = 0;
186 18358646 : i1 = 0;
187 18358646 : i2 = 0;
188 : }
189 154755322 : else if (i == 1)
190 : {
191 18357244 : i0 = 1;
192 18357244 : i1 = 0;
193 18357244 : i2 = 0;
194 : }
195 151793708 : else if (i == 2)
196 : {
197 18357738 : i0 = 1;
198 18357738 : i1 = 1;
199 18357738 : i2 = 0;
200 : }
201 148831866 : else if (i == 3)
202 : {
203 18359140 : i0 = 0;
204 18359140 : i1 = 1;
205 18359140 : i2 = 0;
206 : }
207 145870552 : else if (i == 4)
208 : {
209 18360094 : i0 = 0;
210 18360094 : i1 = 0;
211 18360094 : i2 = 1;
212 : }
213 142908462 : else if (i == 5)
214 : {
215 18358692 : i0 = 1;
216 18358692 : i1 = 0;
217 18358692 : i2 = 1;
218 : }
219 139945844 : else if (i == 6)
220 : {
221 18359186 : i0 = 1;
222 18359186 : i1 = 1;
223 18359186 : i2 = 1;
224 : }
225 136982998 : else if (i == 7)
226 : {
227 18360588 : i0 = 0;
228 18360588 : i1 = 1;
229 18360588 : i2 = 1;
230 : }
231 : // Edge 0
232 827674666 : else if (i < 8 + e)
233 : {
234 27552752 : i0 = i - 6;
235 27552752 : i1 = 0;
236 27552752 : i2 = 0;
237 27552752 : if (elem->positive_edge_orientation(0))
238 400914 : xi = -xi_saved;
239 : }
240 : // Edge 1
241 800121914 : else if (i < 8 + 2*e)
242 : {
243 27551036 : i0 = 1;
244 27551036 : i1 = i - e - 6;
245 27551036 : i2 = 0;
246 27551036 : if (elem->positive_edge_orientation(1))
247 19533788 : eta = -eta_saved;
248 : }
249 : // Edge 2
250 772570878 : else if (i < 8 + 3*e)
251 : {
252 27554082 : i0 = i - 2*e - 6;
253 27554082 : i1 = 1;
254 27554082 : i2 = 0;
255 27554082 : if (!elem->positive_edge_orientation(2))
256 403174 : xi = -xi_saved;
257 : }
258 : // Edge 3
259 745016796 : else if (i < 8 + 4*e)
260 : {
261 27555658 : i0 = 0;
262 27555658 : i1 = i - 3*e - 6;
263 27555658 : i2 = 0;
264 27555658 : if (elem->positive_edge_orientation(3))
265 19528120 : eta = -eta_saved;
266 : }
267 : // Edge 4
268 717461138 : else if (i < 8 + 5*e)
269 : {
270 27556818 : i0 = 0;
271 27556818 : i1 = 0;
272 27556818 : i2 = i - 4*e - 6;
273 27556818 : if (elem->positive_edge_orientation(4))
274 19562650 : zeta = -zeta_saved;
275 : }
276 : // Edge 5
277 689904320 : else if (i < 8 + 6*e)
278 : {
279 27552196 : i0 = 1;
280 27552196 : i1 = 0;
281 27552196 : i2 = i - 5*e - 6;
282 27552196 : if (elem->positive_edge_orientation(5))
283 19571400 : zeta = -zeta_saved;
284 : }
285 : // Edge 6
286 662352124 : else if (i < 8 + 7*e)
287 : {
288 27553526 : i0 = 1;
289 27553526 : i1 = 1;
290 27553526 : i2 = i - 6*e - 6;
291 27553526 : if (elem->positive_edge_orientation(6))
292 19579660 : zeta = -zeta_saved;
293 : }
294 : // Edge 7
295 634798598 : else if (i < 8 + 8*e)
296 : {
297 27558148 : i0 = 0;
298 27558148 : i1 = 1;
299 27558148 : i2 = i - 7*e - 6;
300 27558148 : if (elem->positive_edge_orientation(7))
301 19570910 : zeta = -zeta_saved;
302 : }
303 : // Edge 8
304 607240450 : else if (i < 8 + 9*e)
305 : {
306 27557592 : i0 = i - 8*e - 6;
307 27557592 : i1 = 0;
308 27557592 : i2 = 1;
309 27557592 : if (elem->positive_edge_orientation(8))
310 409730 : xi = -xi_saved;
311 : }
312 : // Edge 9
313 579682858 : else if (i < 8 + 10*e)
314 : {
315 27555876 : i0 = 1;
316 27555876 : i1 = i - 9*e - 6;
317 27555876 : i2 = 1;
318 27555876 : if (elem->positive_edge_orientation(9))
319 19538912 : eta = -eta_saved;
320 : }
321 : // Edge 10
322 552126982 : else if (i < 8 + 11*e)
323 : {
324 27558922 : i0 = i - 10*e - 6;
325 27558922 : i1 = 1;
326 27558922 : i2 = 1;
327 27558922 : if (!elem->positive_edge_orientation(10))
328 411990 : xi = -xi_saved;
329 : }
330 : // Edge 11
331 524568060 : else if (i < 8 + 12*e)
332 : {
333 27560498 : i0 = 0;
334 27560498 : i1 = i - 11*e - 6;
335 27560498 : i2 = 1;
336 27560498 : if (elem->positive_edge_orientation(11))
337 19533244 : eta = -eta_saved;
338 : }
339 : // Faces
340 497007562 : else if (i < 8 + 12*e + 6*e*e)
341 : {
342 355220604 : const unsigned int face = (i - 8 - 12*e) / (e*e);
343 355220604 : const unsigned int basisnum = (i - 8 - 12*e) % (e*e);
344 :
345 : // The face's own frame has its origin at the least vertex and its
346 : // first axis toward that vertex's lesser neighbor, so every element
347 : // sharing the face agrees on it. The origin sits at -1 along both
348 : // axes, so the in-face coordinates are projections onto them.
349 : const std::vector<unsigned int> v =
350 384011484 : elem->oriented_face_vertices(face);
351 355220604 : const Point origin = elem->master_point(v[0]);
352 355220604 : const Point first = (elem->master_point(v[1]) - origin) / 2;
353 355220604 : const Point second = (elem->master_point(v[3]) - origin) / 2;
354 28790880 : const Point p_saved(xi_saved, eta_saved, zeta_saved);
355 :
356 355220604 : Real * const coord[3] = {&xi, &eta, &zeta};
357 355220604 : unsigned int * const index[3] = {&i0, &i1, &i2};
358 :
359 1420882416 : for (const auto d : make_range(3u))
360 1065661812 : if (first(d) != 0)
361 : {
362 355220604 : *index[d] = square_number_row[basisnum] + 2;
363 384011484 : *coord[d] = p_saved * first;
364 : }
365 710441208 : else if (second(d) != 0)
366 : {
367 355220604 : *index[d] = square_number_column[basisnum] + 2;
368 384011484 : *coord[d] = p_saved * second;
369 : }
370 : else
371 : // The vertex mode that is 1 on this face
372 355220604 : *index[d] = (origin(d) > 0);
373 : }
374 :
375 : // Internal DoFs
376 : else
377 : {
378 141786958 : unsigned int basisnum = i - 8 - 12*e - 6*e*e;
379 141786958 : i0 = cube_number_column[basisnum] + 2;
380 141786958 : i1 = cube_number_row[basisnum] + 2;
381 141786958 : i2 = cube_number_page[basisnum] + 2;
382 : }
383 974545994 : }
384 :
385 :
386 : // Order the barycentric coordinates of triangular prism face \p face_num
387 : // by the face's vertices. The triangle interior basis is not symmetric
388 : // in them, so both elements sharing the face must order them alike.
389 57351816 : void orient_triangle_coords(const Elem & elem,
390 : const unsigned int face_num,
391 : const Point & xi_eta_saved,
392 : Point & xi_eta)
393 : {
394 : const std::vector<unsigned int> face_vertex =
395 57351816 : elem.oriented_face_vertices(face_num);
396 :
397 57351816 : const Real barycentric[3] = {1 - xi_eta_saved(0) - xi_eta_saved(1),
398 10242216 : xi_eta_saved(0),
399 57351816 : xi_eta_saved(1)};
400 :
401 : // The face's least vertex takes the place of the triangle's vertex 0, so xi and eta are the
402 : // barycentric coordinates of the next two vertices in the face's order. Face 4's vertices
403 : // 3-5 lie above vertices 0-2, so % 3 maps the vertices of either face to the triangle's.
404 57351816 : xi_eta(0) = barycentric[face_vertex[1] % 3];
405 57351816 : xi_eta(1) = barycentric[face_vertex[2] % 3];
406 57351816 : }
407 :
408 :
409 955370748 : void prism_indices(const Elem * elem,
410 : const unsigned int totalorder,
411 : const unsigned int i,
412 : Point & xi_eta, Real & zeta,
413 : unsigned int & i01,
414 : unsigned int & i2)
415 : {
416 : // the number of DoFs per edge appears everywhere:
417 955370748 : const unsigned int e = totalorder - 1u;
418 :
419 85295454 : libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u)*(totalorder+2u)/2u);
420 :
421 955370748 : Point xi_eta_saved = xi_eta;
422 955370748 : Real zeta_saved = zeta;
423 :
424 : // Vertices:
425 955370748 : if (i == 0)
426 : {
427 25320466 : i01 = 0;
428 25320466 : i2 = 0;
429 : }
430 166070340 : else if (i == 1)
431 : {
432 25318256 : i01 = 1;
433 25318256 : i2 = 0;
434 : }
435 161549876 : else if (i == 2)
436 : {
437 25316670 : i01 = 2;
438 25316670 : i2 = 0;
439 : }
440 157029828 : else if (i == 3)
441 : {
442 25320706 : i01 = 0;
443 25320706 : i2 = 1;
444 : }
445 152509260 : else if (i == 4)
446 : {
447 25318496 : i01 = 1;
448 25318496 : i2 = 1;
449 : }
450 147988796 : else if (i == 5)
451 : {
452 25316910 : i01 = 2;
453 25316910 : i2 = 1;
454 : }
455 : // Edges 0,1,2 (vertices 6,7,8)
456 803459244 : else if (i < 6 + 3*e)
457 : {
458 : // The TRI code will handle any flips here
459 123395400 : i01 = i - 3;
460 123395400 : i2 = 0;
461 : }
462 : // Edge 3,4,5 (vertices 9,10,11)
463 680063844 : else if (i < 6 + 6*e)
464 : {
465 123505776 : i01 = (i - 6 - 3*e)/e; // which tri DoF are we?
466 123505776 : i2 = (i - 6 - 3*e)%e+2; // edge DoF? +2 to skip endpoints
467 : // EDGE evaluations don't flip, so handle that here
468 123505776 : if (elem->positive_edge_orientation(i01+3))
469 123175096 : zeta = -zeta;
470 : }
471 : // Edge 6,7,8 (vertices 12,13,14)
472 556558068 : else if (i < 6 + 9*e)
473 : {
474 : // The TRI code will handle any flips here
475 123397272 : i01 = i - 3 - 6*e;
476 123397272 : i2 = 1;
477 : }
478 : // Faces 1, 2, 3 (nodes 15, 16, 17), the quadrilaterals
479 433160796 : else if (i < 6 + 9*e + 3*e*e)
480 : {
481 295542216 : const unsigned int face = (i - 6 - 9*e) / (e*e) + 1;
482 295542216 : const unsigned int basisnum = (i - 6 - 9*e) % (e*e);
483 :
484 : // indexes in edge numbering
485 295542216 : const unsigned int s0 = square_number_row[basisnum] + 2;
486 295542216 : const unsigned int s1 = square_number_column[basisnum] + 2;
487 :
488 : // The face's side node map runs along triangle edge face-1 from
489 : // n0 to n1, then up to the top triangle and back. Its own frame has
490 : // its origin at the least vertex and its first axis toward that
491 : // vertex's lesser neighbor, so every element sharing the face
492 : // agrees on it.
493 295542216 : const unsigned int n0 = elem->local_side_node(face, 0),
494 295542216 : n1 = elem->local_side_node(face, 1);
495 : const std::vector<unsigned int> v =
496 321929412 : elem->oriented_face_vertices(face);
497 :
498 : // The first axis runs along zeta: the edge takes the column mode
499 295542216 : const bool swap = (v[0] < 3) != (v[1] < 3);
500 295542216 : i01 = (swap ? s1 : s0) + 1 + (face-1)*e; // edge to triangle side numbering
501 295542216 : i2 = swap ? s0 : s1;
502 :
503 : // The origin is on the top triangle
504 295542216 : if (v[0] >= 3)
505 294895420 : zeta = -zeta_saved;
506 :
507 : // The origin is at the far end of the edge: reflect the triangle by
508 : // exchanging the barycentric coordinates of the edge's vertices
509 295542216 : if (v[0] % 3 == n1 % 3)
510 : {
511 160502448 : Real barycentric[3] = {1 - xi_eta_saved(0) - xi_eta_saved(1),
512 27924800 : xi_eta_saved(0),
513 160502448 : xi_eta_saved(1)};
514 160502448 : std::swap(barycentric[n0 % 3], barycentric[n1 % 3]);
515 160502448 : xi_eta(0) = barycentric[1];
516 160502448 : xi_eta(1) = barycentric[2];
517 : }
518 : }
519 : // Face 0, node 18 - node order due to hierarchic numbering
520 137618580 : else if (i < 6 + 9*e + 3*e*e + e*(e-1)/2)
521 : {
522 28675644 : i01 = i - 3 - 6*e - 3*e*e;
523 28675644 : i2 = 0;
524 28675644 : orient_triangle_coords(*elem, 0, xi_eta_saved, xi_eta);
525 : }
526 : // Face 4
527 108942936 : else if (i < 6 + 9*e + 3*e*e + e*(e-1))
528 : {
529 28676172 : i01 = i - 3 - 6*e - 3*e*e - e*(e-1)/2;
530 28676172 : i2 = 1;
531 28676172 : orient_triangle_coords(*elem, 4, xi_eta_saved, xi_eta);
532 : }
533 : // Internal DoFs
534 : else
535 : {
536 : // We won't bother with any internal DoF reordering / flipping;
537 : // that's fine unless we ever get to 4D.
538 80266764 : unsigned int basisnum = i - 6 - 9*e - 3*e*e - e*(e-1);
539 80266764 : i01 = prism_number_triangle[basisnum] + 3 + 3*e;
540 80266764 : i2 = prism_number_page[basisnum] + 2;
541 : }
542 955370748 : }
543 :
544 : #endif // LIBMESH_DIM > 2
545 :
546 : } // end anonymous namespace
547 :
548 :
549 :
550 : namespace libMesh
551 : {
552 :
553 :
554 124758451 : LIBMESH_DEFAULT_VECTORIZED_FE(3,HIERARCHIC)
555 134368094 : LIBMESH_DEFAULT_VECTORIZED_FE(3,L2_HIERARCHIC)
556 97033343 : LIBMESH_DEFAULT_VECTORIZED_FE(3,SIDE_HIERARCHIC)
557 :
558 :
559 : template <>
560 0 : Real FE<3,HIERARCHIC>::shape(const ElemType,
561 : const Order,
562 : const unsigned int,
563 : const Point &)
564 : {
565 0 : libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
566 : return 0.;
567 : }
568 :
569 :
570 :
571 : template <>
572 0 : Real FE<3,L2_HIERARCHIC>::shape(const ElemType,
573 : const Order,
574 : const unsigned int,
575 : const Point &)
576 : {
577 0 : libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
578 : return 0.;
579 : }
580 :
581 :
582 :
583 : template <>
584 0 : Real FE<3,SIDE_HIERARCHIC>::shape(const ElemType,
585 : const Order,
586 : const unsigned int,
587 : const Point &)
588 : {
589 0 : libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
590 : return 0.;
591 : }
592 :
593 :
594 :
595 : template <>
596 1140188348 : Real FE<3,HIERARCHIC>::shape(const Elem * elem,
597 : const Order order,
598 : const unsigned int i,
599 : const Point & p,
600 : const bool add_p_level)
601 : {
602 1140188348 : return fe_hierarchic_3D_shape<HIERARCHIC>(elem, order, i, p, add_p_level);
603 : }
604 :
605 :
606 : template <>
607 0 : Real FE<3,HIERARCHIC>::shape(const FEType fet,
608 : const Elem * elem,
609 : const unsigned int i,
610 : const Point & p,
611 : const bool add_p_level)
612 : {
613 0 : return fe_hierarchic_3D_shape<HIERARCHIC>(elem, fet.order, i, p, add_p_level);
614 : }
615 :
616 :
617 :
618 :
619 : template <>
620 1520145288 : Real FE<3,L2_HIERARCHIC>::shape(const Elem * elem,
621 : const Order order,
622 : const unsigned int i,
623 : const Point & p,
624 : const bool add_p_level)
625 : {
626 1520145288 : return fe_hierarchic_3D_shape<L2_HIERARCHIC>(elem, order, i, p, add_p_level);
627 : }
628 :
629 :
630 : template <>
631 0 : Real FE<3,L2_HIERARCHIC>::shape(const FEType fet,
632 : const Elem * elem,
633 : const unsigned int i,
634 : const Point & p,
635 : const bool add_p_level)
636 : {
637 0 : return fe_hierarchic_3D_shape<L2_HIERARCHIC>(elem, fet.order, i, p, add_p_level);
638 : }
639 :
640 :
641 :
642 : template <>
643 4529357029 : Real FE<3,SIDE_HIERARCHIC>::shape(const Elem * elem,
644 : const Order order,
645 : const unsigned int i,
646 : const Point & p,
647 : const bool add_p_level)
648 : {
649 : #if LIBMESH_DIM == 3
650 376761768 : libmesh_assert(elem);
651 4529357029 : const ElemType type = elem->type();
652 :
653 4906118797 : const Order totalorder = order + add_p_level*elem->p_level();
654 :
655 4529357029 : switch (type)
656 : {
657 41575248 : case HEX27:
658 : {
659 500333274 : const unsigned int dofs_per_side = (totalorder+1u)*(totalorder+1u);
660 41575248 : libmesh_assert_less(i, 6*dofs_per_side);
661 :
662 500333274 : const unsigned int sidenum = cube_side(p);
663 500333274 : if (sidenum > 5)
664 25920 : return std::numeric_limits<Real>::quiet_NaN();
665 :
666 499903326 : const unsigned int dof_offset = sidenum * dofs_per_side;
667 :
668 499903326 : if (i < dof_offset) // i is on a previous side
669 17119013 : return 0;
670 :
671 293517049 : if (i >= dof_offset + dofs_per_side) // i is on a later side
672 17227987 : return 0;
673 :
674 86646501 : if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
675 14264 : return 1;
676 :
677 86478594 : unsigned int side_i = i - dof_offset;
678 :
679 86478594 : std::unique_ptr<const Elem> side = elem->build_side_ptr(sidenum);
680 :
681 86478594 : Point sidep = cube_side_point(sidenum, p);
682 :
683 86478594 : cube_remap(side_i, *elem, sidenum, *side, totalorder, sidep);
684 :
685 86478594 : return FE<2,HIERARCHIC>::shape(side.get(), order, side_i, sidep, add_p_level);
686 72102466 : }
687 :
688 240081920 : case TET14:
689 : {
690 2884736640 : const unsigned int dofs_per_side = (totalorder+1u)*(totalorder+2u)/2u;
691 240081920 : libmesh_assert_less(i, 4*dofs_per_side);
692 :
693 2884736640 : const Real zeta[4] = { Real(1.) - p(0) - p(1) - p(2), p(0), p(1), p(2) };
694 :
695 240081920 : unsigned int face_num = 0;
696 2884736640 : if (zeta[0] > zeta[3] &&
697 992688136 : zeta[1] > zeta[3] &&
698 79217916 : zeta[2] > zeta[3])
699 : {
700 60007536 : face_num = 0;
701 : }
702 2163466560 : else if (zeta[0] > zeta[2] &&
703 750599560 : zeta[1] > zeta[2] &&
704 60010416 : zeta[3] > zeta[2])
705 : {
706 60010416 : face_num = 1;
707 : }
708 1442858576 : else if (zeta[1] > zeta[0] &&
709 721207504 : zeta[2] > zeta[0] &&
710 59987424 : zeta[3] > zeta[0])
711 : {
712 59987424 : face_num = 2;
713 : }
714 : else
715 : {
716 : // We'd better not be right between two faces
717 60076544 : libmesh_assert (zeta[0] > zeta[1] &&
718 : zeta[2] > zeta[1] &&
719 : zeta[3] > zeta[1]);
720 60076544 : face_num = 3;
721 : }
722 :
723 2884736640 : if (i < face_num * dofs_per_side ||
724 1802742588 : i >= (face_num+1) * dofs_per_side)
725 180061440 : return 0;
726 :
727 721184160 : if (totalorder == 0)
728 178112 : return 1;
729 :
730 : const std::vector<unsigned int> face_vertex =
731 778887144 : elem->oriented_face_vertices(face_num);
732 :
733 : // We only need a Tri3 to evaluate L2_HIERARCHIC on the affine
734 : // master element
735 119684736 : Tri3 side;
736 :
737 : // We pinky swear not to modify these nodes
738 59842368 : Elem & e = const_cast<Elem &>(*elem);
739 719044776 : side.set_node(0, e.node_ptr(face_vertex[0]));
740 719044776 : side.set_node(1, e.node_ptr(face_vertex[1]));
741 719044776 : side.set_node(2, e.node_ptr(face_vertex[2]));
742 :
743 719044776 : const unsigned int basisnum = i - face_num*dofs_per_side;
744 :
745 719044776 : Point sidep {zeta[face_vertex[1]], zeta[face_vertex[2]]};
746 :
747 719044776 : return FE<2,L2_HIERARCHIC>::shape(&side, totalorder,
748 59842368 : basisnum, sidep, false);
749 : }
750 :
751 95104600 : case PRISM20:
752 : case PRISM21:
753 : {
754 1144287115 : const unsigned int dofs_per_quad = (totalorder+1u)*(totalorder+1u);
755 1144287115 : const unsigned int dofs_per_tri = (totalorder+1u)*(totalorder+2u)/2u;
756 95104600 : libmesh_assert_less(i, 3*dofs_per_quad + 2*dofs_per_tri);
757 :
758 : // We only need a Tri3 or Quad4 to evaluate L2_HIERARCHIC on
759 : // the affine master element
760 190209200 : Tri3 tri;
761 190209200 : Quad4 quad;
762 95104600 : Elem * side = &quad;
763 95104600 : unsigned int dofs_on_side = dofs_per_quad;
764 :
765 : // We pinky swear not to modify the nodes we'll point to
766 95104600 : Elem & e = const_cast<Elem &>(*elem);
767 :
768 95104600 : Point sidep;
769 :
770 : // Face number calculation is tricky - the ordering of side
771 : // nodes on Prisms does *not* match the ordering of sides!
772 : // (the mid-triangle side nodes were added "later")
773 : // Here face_num will be the numbering that matches the side
774 : // number, but i_offset will have to consider the nodal
775 : // ordering.
776 95104600 : unsigned int face_num = 0;
777 95104600 : unsigned int i_offset = 0;
778 :
779 : // Triangular coordinates
780 1144287115 : const Real zeta[3] = { Real(1.) - p(0) - p(1), p(0), p(1) };
781 :
782 : // Closeness to midplane
783 1144287115 : const Real zmid = 1 - std::abs(p(2));
784 :
785 1144287115 : if (zeta[1] > zeta[2] && zeta[0] > zeta[2] &&
786 377613345 : zmid > 3*zeta[2]) // face 1, quad
787 : {
788 21012280 : face_num = 1;
789 21012280 : i_offset = 0;
790 : }
791 891497900 : else if (zeta[1] > zeta[0] && zeta[2] > zeta[0] &&
792 378143710 : zmid > 3*zeta[0]) // face 2, quad
793 : {
794 21012280 : face_num = 2;
795 21012280 : i_offset = dofs_per_quad;
796 : }
797 638853620 : else if (zeta[0] > zeta[1] && zeta[2] > zeta[1] &&
798 376661480 : zmid > 3*zeta[1]) // face 3, quad
799 : {
800 21012280 : face_num = 3;
801 252934150 : i_offset = 2*dofs_per_quad;
802 : }
803 586705960 : else if (p(2) + 1 < 3*zeta[0] &&
804 402016570 : p(2) + 1 < 3*zeta[1] &&
805 193260030 : p(2) + 1 < 3*zeta[2]) // face 0, tri
806 : {
807 16097100 : face_num = 0;
808 193260030 : i_offset = 3*dofs_per_quad;
809 16097100 : dofs_on_side = dofs_per_tri;
810 16097100 : side = &tri;
811 : }
812 369348220 : else if (1 - p(2) < 3*zeta[0] &&
813 208630100 : 1 - p(2) < 3*zeta[1] &&
814 192659440 : 1 - p(2) < 3*zeta[2]) // face 4, tri
815 : {
816 15970660 : face_num = 4;
817 192659440 : i_offset = dofs_per_tri + 3*dofs_per_quad;
818 15970660 : dofs_on_side = dofs_per_tri;
819 15970660 : side = &tri;
820 : }
821 : else
822 : {
823 0 : libmesh_error_msg("Evaluating SIDE_HIERARCHIC right between two Prism faces?");
824 : }
825 :
826 1144287115 : if (i < i_offset ||
827 663248625 : i >= i_offset + dofs_on_side)
828 75535040 : return 0;
829 :
830 235450955 : if (totalorder == 0)
831 1400 : return 1;
832 :
833 : const std::vector<unsigned int> face_vertex =
834 255001500 : elem->oriented_face_vertices(face_num);
835 :
836 255001500 : side->set_node(0, e.node_ptr(face_vertex[0]));
837 255001500 : side->set_node(1, e.node_ptr(face_vertex[1]));
838 255001500 : side->set_node(2, e.node_ptr(face_vertex[2]));
839 255001500 : if (face_vertex.size() > 3)
840 194216130 : side->set_node(3, e.node_ptr(face_vertex[3]));
841 :
842 235433340 : if (face_num == 0 || face_num == 4)
843 56121930 : sidep = {zeta[face_vertex[1]%3], zeta[face_vertex[2]%3]};
844 : else
845 : {
846 : // Transform a coordinate from the master prism to the
847 : // master quad, based on two vertex indices defining the
848 : // coordinate's direction
849 358622820 : auto coord_val = [p](int v1, int v2){
850 358622820 : if (v2-v1 == 3)
851 78513840 : return p(2);
852 280108980 : else if (v2-v1 == -3)
853 100797570 : return -p(2);
854 179311410 : else if (v1%3 == 0 && v2%3 == 1)
855 31965930 : return 2*p(0)-1;
856 147345480 : else if (v2%3 == 0 && v1%3 == 1)
857 27804540 : return 1-2*p(0);
858 119540940 : else if (v1%3 == 1 && v2%3 == 2)
859 31432920 : return p(1)-p(0);
860 88108020 : else if (v2%3 == 1 && v1%3 == 2)
861 28303320 : return p(0)-p(1);
862 59804700 : else if (v1%3 == 2 && v2%3 == 0)
863 20748270 : return 1-2*p(1);
864 39056430 : else if (v2%3 == 2 && v1%3 == 0)
865 39056430 : return 2*p(1)-1;
866 : else
867 0 : libmesh_error();
868 179311410 : };
869 :
870 209120850 : sidep = {coord_val(face_vertex[0], face_vertex[1]),
871 29809440 : coord_val(face_vertex[0], face_vertex[3])};
872 : }
873 :
874 235433340 : const unsigned int basisnum = i - i_offset;
875 :
876 235433340 : return FE<2,L2_HIERARCHIC>::shape(side, totalorder,
877 19568160 : basisnum, sidep, false);
878 : }
879 :
880 :
881 0 : default:
882 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
883 : }
884 :
885 : #else // LIBMESH_DIM != 3
886 : libmesh_ignore(elem, order, i, p, add_p_level);
887 : libmesh_not_implemented();
888 : #endif
889 : }
890 :
891 :
892 : template <>
893 0 : Real FE<3,SIDE_HIERARCHIC>::shape(const FEType fet,
894 : const Elem * elem,
895 : const unsigned int i,
896 : const Point & p,
897 : const bool add_p_level)
898 : {
899 0 : return FE<3,SIDE_HIERARCHIC>::shape(elem,fet.order, i, p, add_p_level);
900 : }
901 :
902 :
903 : template <>
904 0 : Real FE<3,HIERARCHIC>::shape_deriv(const ElemType,
905 : const Order,
906 : const unsigned int,
907 : const unsigned int,
908 : const Point & )
909 : {
910 0 : libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
911 : return 0.;
912 : }
913 :
914 :
915 :
916 : template <>
917 0 : Real FE<3,L2_HIERARCHIC>::shape_deriv(const ElemType,
918 : const Order,
919 : const unsigned int,
920 : const unsigned int,
921 : const Point & )
922 : {
923 0 : libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
924 : return 0.;
925 : }
926 :
927 :
928 :
929 : template <>
930 0 : Real FE<3,SIDE_HIERARCHIC>::shape_deriv(const ElemType,
931 : const Order,
932 : const unsigned int,
933 : const unsigned int,
934 : const Point & )
935 : {
936 0 : libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
937 : return 0.;
938 : }
939 :
940 :
941 :
942 : template <>
943 421855626 : Real FE<3,HIERARCHIC>::shape_deriv(const Elem * elem,
944 : const Order order,
945 : const unsigned int i,
946 : const unsigned int j,
947 : const Point & p,
948 : const bool add_p_level)
949 : {
950 808746099 : return fe_hierarchic_3D_shape_deriv<HIERARCHIC>(elem, order, i, j, p, add_p_level);
951 : }
952 :
953 :
954 : template <>
955 0 : Real FE<3,HIERARCHIC>::shape_deriv(const FEType fet,
956 : const Elem * elem,
957 : const unsigned int i,
958 : const unsigned int j,
959 : const Point & p,
960 : const bool add_p_level)
961 : {
962 0 : return fe_hierarchic_3D_shape_deriv<HIERARCHIC>(elem, fet.order, i, j, p, add_p_level);
963 : }
964 :
965 :
966 :
967 : template <>
968 695091552 : Real FE<3,L2_HIERARCHIC>::shape_deriv(const Elem * elem,
969 : const Order order,
970 : const unsigned int i,
971 : const unsigned int j,
972 : const Point & p,
973 : const bool add_p_level)
974 : {
975 1333177182 : return fe_hierarchic_3D_shape_deriv<L2_HIERARCHIC>(elem, order, i, j, p, add_p_level);
976 : }
977 :
978 :
979 : template <>
980 0 : Real FE<3,L2_HIERARCHIC>::shape_deriv(const FEType fet,
981 : const Elem * elem,
982 : const unsigned int i,
983 : const unsigned int j,
984 : const Point & p,
985 : const bool add_p_level)
986 : {
987 0 : return fe_hierarchic_3D_shape_deriv<L2_HIERARCHIC>(elem, fet.order, i, j, p, add_p_level);
988 : }
989 :
990 :
991 :
992 : template <>
993 2042311680 : Real FE<3,SIDE_HIERARCHIC>::shape_deriv(const Elem * elem,
994 : const Order order,
995 : const unsigned int i,
996 : const unsigned int j,
997 : const Point & p,
998 : const bool add_p_level)
999 : {
1000 : #if LIBMESH_DIM == 3
1001 170192640 : libmesh_assert(elem);
1002 2042311680 : const ElemType type = elem->type();
1003 :
1004 2212504320 : const Order totalorder = order + add_p_level*elem->p_level();
1005 :
1006 2042311680 : if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
1007 62400 : return 0; // constants have zero derivative
1008 :
1009 2041562880 : switch (type)
1010 : {
1011 19448640 : case HEX27:
1012 : {
1013 : // I need to debug the p>2 case here...
1014 233383680 : if (totalorder > 2)
1015 228355200 : return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<3,SIDE_HIERARCHIC>::shape);
1016 :
1017 5028480 : const unsigned int dofs_per_side = (totalorder+1u)*(totalorder+1u);
1018 419040 : libmesh_assert_less(i, 6*dofs_per_side);
1019 :
1020 5028480 : const unsigned int sidenum = cube_side(p);
1021 5028480 : if (sidenum > 5)
1022 0 : return std::numeric_limits<Real>::quiet_NaN();
1023 :
1024 5028480 : const unsigned int dof_offset = sidenum * dofs_per_side;
1025 :
1026 5028480 : if (i < dof_offset) // i is on a previous side
1027 174600 : return 0;
1028 :
1029 2933280 : if (i >= dof_offset + dofs_per_side) // i is on a later side
1030 174600 : return 0;
1031 :
1032 838080 : unsigned int side_i = i - dof_offset;
1033 :
1034 838080 : std::unique_ptr<const Elem> side = elem->build_side_ptr(sidenum);
1035 :
1036 838080 : Point sidep = cube_side_point(sidenum, p);
1037 :
1038 838080 : cube_remap(side_i, *elem, sidenum, *side, totalorder, sidep);
1039 :
1040 : // What direction on the side corresponds to the derivative
1041 : // direction we want?
1042 69840 : unsigned int sidej = 100;
1043 :
1044 : // Do we need a -1 here to flip that direction?
1045 69840 : Real f = 1.;
1046 :
1047 : switch (j)
1048 : {
1049 279360 : case 0: // d()/dxi
1050 : {
1051 : switch (sidenum)
1052 : {
1053 3880 : case 0:
1054 3880 : sidej = 1;
1055 3880 : break;
1056 7760 : case 1:
1057 3880 : sidej = 0;
1058 7760 : break;
1059 3880 : case 2:
1060 3880 : return 0;
1061 46560 : case 3:
1062 3880 : sidej = 0;
1063 3880 : f = -1;
1064 46560 : break;
1065 3880 : case 4:
1066 3880 : return 0;
1067 7760 : case 5:
1068 3880 : sidej = 0;
1069 7760 : break;
1070 0 : default:
1071 0 : libmesh_error();
1072 : }
1073 15520 : break;
1074 : }
1075 279360 : case 1: // d()/deta
1076 : {
1077 : switch (sidenum)
1078 : {
1079 3880 : case 0:
1080 3880 : sidej = 0;
1081 3880 : break;
1082 3880 : case 1:
1083 3880 : return 0;
1084 3880 : case 2:
1085 3880 : sidej = 0;
1086 3880 : break;
1087 3880 : case 3:
1088 3880 : return 0;
1089 46560 : case 4:
1090 3880 : sidej = 0;
1091 3880 : f = -1;
1092 46560 : break;
1093 7760 : case 5:
1094 3880 : sidej = 1;
1095 7760 : break;
1096 0 : default:
1097 0 : libmesh_error();
1098 : }
1099 15520 : break;
1100 : }
1101 279360 : case 2: // d()/dzeta
1102 : {
1103 : switch (sidenum)
1104 : {
1105 3880 : case 0:
1106 3880 : return 0;
1107 15520 : case 1:
1108 : case 2:
1109 : case 3:
1110 : case 4:
1111 15520 : sidej = 1;
1112 15520 : break;
1113 3880 : case 5:
1114 3880 : return 0;
1115 0 : default:
1116 0 : libmesh_error();
1117 : }
1118 15520 : break;
1119 : }
1120 :
1121 0 : default:
1122 0 : libmesh_error_msg("Invalid derivative index j = " << j);
1123 : }
1124 :
1125 558720 : return f * FE<2,HIERARCHIC>::shape_deriv(side.get(), order,
1126 : side_i, sidej, sidep,
1127 558720 : add_p_level);
1128 698400 : }
1129 :
1130 1808179200 : case TET14:
1131 : case PRISM20:
1132 : case PRISM21:
1133 : {
1134 1808179200 : return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<3,SIDE_HIERARCHIC>::shape);
1135 : }
1136 :
1137 0 : default:
1138 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
1139 : }
1140 :
1141 : #else // LIBMESH_DIM != 3
1142 : libmesh_ignore(elem, order, i, j, p, add_p_level);
1143 : libmesh_not_implemented();
1144 : #endif
1145 : }
1146 :
1147 :
1148 : template <>
1149 0 : Real FE<3,SIDE_HIERARCHIC>::shape_deriv(const FEType fet,
1150 : const Elem * elem,
1151 : const unsigned int i,
1152 : const unsigned int j,
1153 : const Point & p,
1154 : const bool add_p_level)
1155 : {
1156 0 : return FE<3,SIDE_HIERARCHIC>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
1157 : }
1158 :
1159 :
1160 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1161 :
1162 : template <>
1163 0 : Real FE<3,HIERARCHIC>::shape_second_deriv(const ElemType,
1164 : const Order,
1165 : const unsigned int,
1166 : const unsigned int,
1167 : const Point & )
1168 : {
1169 0 : libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
1170 : return 0.;
1171 : }
1172 :
1173 :
1174 :
1175 : template <>
1176 0 : Real FE<3,L2_HIERARCHIC>::shape_second_deriv(const ElemType,
1177 : const Order,
1178 : const unsigned int,
1179 : const unsigned int,
1180 : const Point & )
1181 : {
1182 0 : libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
1183 : return 0.;
1184 : }
1185 :
1186 :
1187 :
1188 : template <>
1189 0 : Real FE<3,SIDE_HIERARCHIC>::shape_second_deriv(const ElemType,
1190 : const Order,
1191 : const unsigned int,
1192 : const unsigned int,
1193 : const Point & )
1194 : {
1195 0 : libmesh_error_msg("Hierarchic shape functions require an Elem for edge/face orientation.");
1196 : return 0.;
1197 : }
1198 :
1199 :
1200 :
1201 : template <>
1202 156002052 : Real FE<3,HIERARCHIC>::shape_second_deriv(const Elem * elem,
1203 : const Order order,
1204 : const unsigned int i,
1205 : const unsigned int j,
1206 : const Point & p,
1207 : const bool add_p_level)
1208 : {
1209 299079726 : return fe_hierarchic_3D_shape_second_deriv<HIERARCHIC>(elem, order, i, j, p, add_p_level);
1210 : }
1211 :
1212 :
1213 :
1214 : template <>
1215 0 : Real FE<3,HIERARCHIC>::shape_second_deriv(const FEType fet,
1216 : const Elem * elem,
1217 : const unsigned int i,
1218 : const unsigned int j,
1219 : const Point & p,
1220 : const bool add_p_level)
1221 : {
1222 0 : return fe_hierarchic_3D_shape_second_deriv<HIERARCHIC>(elem, fet.order, i, j, p, add_p_level);
1223 : }
1224 :
1225 :
1226 :
1227 : template <>
1228 252148476 : Real FE<3,L2_HIERARCHIC>::shape_second_deriv(const Elem * elem,
1229 : const Order order,
1230 : const unsigned int i,
1231 : const unsigned int j,
1232 : const Point & p,
1233 : const bool add_p_level)
1234 : {
1235 483628464 : return fe_hierarchic_3D_shape_second_deriv<L2_HIERARCHIC>(elem, order, i, j, p, add_p_level);
1236 : }
1237 :
1238 :
1239 : template <>
1240 0 : Real FE<3,L2_HIERARCHIC>::shape_second_deriv(const FEType fet,
1241 : const Elem * elem,
1242 : const unsigned int i,
1243 : const unsigned int j,
1244 : const Point & p,
1245 : const bool add_p_level)
1246 : {
1247 0 : return fe_hierarchic_3D_shape_second_deriv<L2_HIERARCHIC>(elem, fet.order, i, j, p, add_p_level);
1248 : }
1249 :
1250 :
1251 : template <>
1252 826168320 : Real FE<3,SIDE_HIERARCHIC>::shape_second_deriv(const Elem * elem,
1253 : const Order order,
1254 : const unsigned int i,
1255 : const unsigned int j,
1256 : const Point & p,
1257 : const bool add_p_level)
1258 : {
1259 : #if LIBMESH_DIM == 3
1260 68847360 : libmesh_assert(elem);
1261 826168320 : const ElemType type = elem->type();
1262 :
1263 895015680 : const Order totalorder = order + add_p_level*elem->p_level();
1264 :
1265 826168320 : if (totalorder == 0) // special case since raw HIERARCHIC lacks CONSTANTs
1266 124800 : return 0; // constants have zero derivative
1267 :
1268 824670720 : switch (type)
1269 : {
1270 8449920 : case HEX27:
1271 : {
1272 : // I need to debug the p>2 case here...
1273 101399040 : if (totalorder > 2)
1274 91342080 : return fe_fdm_second_deriv(elem, order, i, j, p, add_p_level,
1275 7611840 : FE<3,SIDE_HIERARCHIC>::shape_deriv);
1276 :
1277 10056960 : const unsigned int dofs_per_side = (totalorder+1u)*(totalorder+1u);
1278 838080 : libmesh_assert_less(i, 6*dofs_per_side);
1279 :
1280 10056960 : const unsigned int sidenum = cube_side(p);
1281 10056960 : if (sidenum > 5)
1282 0 : return std::numeric_limits<Real>::quiet_NaN();
1283 :
1284 10056960 : const unsigned int dof_offset = sidenum * dofs_per_side;
1285 :
1286 10056960 : if (i < dof_offset) // i is on a previous side
1287 349200 : return 0;
1288 :
1289 5866560 : if (i >= dof_offset + dofs_per_side) // i is on a later side
1290 349200 : return 0;
1291 :
1292 1676160 : unsigned int side_i = i - dof_offset;
1293 :
1294 1676160 : std::unique_ptr<const Elem> side = elem->build_side_ptr(sidenum);
1295 :
1296 1676160 : Point sidep = cube_side_point(sidenum, p);
1297 :
1298 1676160 : cube_remap(side_i, *elem, sidenum, *side, totalorder, sidep);
1299 :
1300 : // What second derivative or mixed derivative on the side
1301 : // corresponds to the xi/eta/zeta mix we were asked for?
1302 139680 : unsigned int sidej = 100;
1303 :
1304 : // Do we need a -1 here to flip the final derivative value?
1305 139680 : Real f = 1.;
1306 :
1307 : switch (j)
1308 : {
1309 279360 : case 0: // d^2()/dxi^2
1310 : {
1311 : switch (sidenum)
1312 : {
1313 3880 : case 0:
1314 3880 : sidej = 2;
1315 3880 : break;
1316 7760 : case 1:
1317 3880 : sidej = 0;
1318 7760 : break;
1319 3880 : case 2:
1320 3880 : return 0;
1321 7760 : case 3:
1322 3880 : sidej = 0;
1323 7760 : break;
1324 3880 : case 4:
1325 3880 : return 0;
1326 7760 : case 5:
1327 3880 : sidej = 0;
1328 7760 : break;
1329 0 : default:
1330 0 : libmesh_error();
1331 : }
1332 15520 : break;
1333 : }
1334 279360 : case 1: // d^2()/dxideta
1335 : {
1336 : switch (sidenum)
1337 : {
1338 3880 : case 0:
1339 3880 : sidej = 1;
1340 3880 : break;
1341 15520 : case 1:
1342 : case 2:
1343 : case 3:
1344 : case 4:
1345 15520 : return 0;
1346 3880 : case 5:
1347 3880 : sidej = 1;
1348 3880 : break;
1349 0 : default:
1350 0 : libmesh_error();
1351 : }
1352 7760 : break;
1353 : }
1354 279360 : case 2: // d^2()/deta^2
1355 : {
1356 : switch (sidenum)
1357 : {
1358 3880 : case 0:
1359 3880 : sidej = 0;
1360 3880 : break;
1361 3880 : case 1:
1362 3880 : return 0;
1363 3880 : case 2:
1364 3880 : sidej = 0;
1365 3880 : break;
1366 3880 : case 3:
1367 3880 : return 0;
1368 3880 : case 4:
1369 3880 : sidej = 0;
1370 3880 : break;
1371 3880 : case 5:
1372 3880 : sidej = 2;
1373 3880 : break;
1374 0 : default:
1375 0 : libmesh_error();
1376 : }
1377 15520 : break;
1378 : }
1379 279360 : case 3: // d^2()/dxidzeta
1380 : {
1381 : switch (sidenum)
1382 : {
1383 3880 : case 0:
1384 3880 : return 0;
1385 3880 : case 1:
1386 3880 : sidej = 1;
1387 3880 : break;
1388 3880 : case 2:
1389 3880 : return 0;
1390 3880 : case 3:
1391 3880 : sidej = 1;
1392 3880 : f = -1;
1393 3880 : break;
1394 7760 : case 4:
1395 : case 5:
1396 7760 : return 0;
1397 0 : default:
1398 0 : libmesh_error();
1399 : }
1400 7760 : break;
1401 : }
1402 279360 : case 4: // d^2()/detadzeta
1403 : {
1404 : switch (sidenum)
1405 : {
1406 7760 : case 0:
1407 : case 1:
1408 7760 : return 0;
1409 3880 : case 2:
1410 3880 : sidej = 1;
1411 3880 : break;
1412 3880 : case 3:
1413 3880 : return 0;
1414 3880 : case 4:
1415 3880 : sidej = 1;
1416 3880 : f = -1;
1417 3880 : break;
1418 3880 : case 5:
1419 3880 : return 0;
1420 0 : default:
1421 0 : libmesh_error();
1422 : }
1423 7760 : break;
1424 : }
1425 279360 : case 5: // d^2()/dzeta^2
1426 : {
1427 : switch (sidenum)
1428 : {
1429 3880 : case 0:
1430 3880 : return 0;
1431 15520 : case 1:
1432 : case 2:
1433 : case 3:
1434 : case 4:
1435 15520 : sidej = 2;
1436 15520 : break;
1437 3880 : case 5:
1438 3880 : return 0;
1439 0 : default:
1440 0 : libmesh_error();
1441 : }
1442 15520 : break;
1443 : }
1444 :
1445 0 : default:
1446 0 : libmesh_error_msg("Invalid derivative index j = " << j);
1447 : }
1448 :
1449 838080 : return f * FE<2,HIERARCHIC>::shape_second_deriv(side.get(),
1450 : order, side_i,
1451 : sidej, sidep,
1452 838080 : add_p_level);
1453 1396800 : }
1454 :
1455 723271680 : case TET14:
1456 : case PRISM20:
1457 : case PRISM21:
1458 : {
1459 723271680 : return fe_fdm_second_deriv(elem, order, i, j, p, add_p_level,
1460 723271680 : FE<3,SIDE_HIERARCHIC>::shape_deriv);
1461 : }
1462 :
1463 0 : default:
1464 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
1465 : }
1466 :
1467 : #else // LIBMESH_DIM != 3
1468 : libmesh_ignore(elem, order, i, j, p, add_p_level);
1469 : libmesh_not_implemented();
1470 : #endif
1471 : }
1472 :
1473 :
1474 : template <>
1475 0 : Real FE<3,SIDE_HIERARCHIC>::shape_second_deriv(const FEType fet,
1476 : const Elem * elem,
1477 : const unsigned int i,
1478 : const unsigned int j,
1479 : const Point & p,
1480 : const bool add_p_level)
1481 : {
1482 0 : return FE<3,SIDE_HIERARCHIC>::shape_second_deriv(elem, fet.order, i, j, p, add_p_level);
1483 : }
1484 :
1485 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
1486 :
1487 : } // namespace libMesh
1488 :
1489 :
1490 :
1491 : namespace
1492 : {
1493 : using namespace libMesh;
1494 :
1495 :
1496 515418714 : unsigned int cube_side (const Point & p)
1497 : {
1498 515418714 : const Real xi = p(0), eta = p(1), zeta = p(2);
1499 515418714 : const Real absxi = std::abs(xi),
1500 515418714 : abseta = std::abs(eta),
1501 515418714 : abszeta = std::abs(zeta);
1502 515418714 : const Real maxabs_xi_eta = std::max(absxi, abseta),
1503 515418714 : maxabs_xi_zeta = std::max(absxi, abszeta),
1504 515418714 : maxabs_eta_zeta = std::max(abseta, abszeta);
1505 :
1506 515418714 : if (zeta < -maxabs_xi_eta)
1507 7151148 : return 0;
1508 429516196 : else if (eta < -maxabs_xi_zeta)
1509 7171126 : return 1;
1510 343477936 : else if (xi > maxabs_eta_zeta)
1511 7188558 : return 2;
1512 257348620 : else if (eta > maxabs_xi_zeta)
1513 7158032 : return 3;
1514 171597400 : else if (xi < -maxabs_eta_zeta)
1515 7006032 : return 4;
1516 86307298 : else if (zeta > maxabs_xi_eta)
1517 85877350 : return 5;
1518 :
1519 : // We need to be able to return invalid values for cases where
1520 : // mixed FE are being evaluated together on edges and vertices
1521 25920 : return 65535;
1522 : }
1523 :
1524 :
1525 :
1526 88992834 : Point cube_side_point(unsigned int sidenum, const Point & p)
1527 : {
1528 7397584 : Point sidep;
1529 :
1530 88992834 : switch (sidenum)
1531 : {
1532 14834316 : case 0:
1533 14834316 : sidep(0) = p(1);
1534 14834316 : sidep(1) = p(0);
1535 14834316 : break;
1536 14866568 : case 1:
1537 14866568 : sidep(0) = p(0);
1538 14866568 : sidep(1) = p(2);
1539 14866568 : break;
1540 14873064 : case 2:
1541 14873064 : sidep(0) = p(1);
1542 14873064 : sidep(1) = p(2);
1543 14873064 : break;
1544 14818788 : case 3:
1545 14818788 : sidep(0) = -p(0);
1546 14818788 : sidep(1) = p(2);
1547 14818788 : break;
1548 14750670 : case 4:
1549 14750670 : sidep(0) = -p(1);
1550 14750670 : sidep(1) = p(2);
1551 14750670 : break;
1552 14849428 : case 5:
1553 14849428 : sidep(0) = p(0);
1554 14849428 : sidep(1) = p(1);
1555 14849428 : break;
1556 0 : default:
1557 0 : libmesh_error();
1558 : }
1559 :
1560 88992834 : return sidep;
1561 : }
1562 :
1563 :
1564 :
1565 : template <FEFamily T>
1566 2660333636 : Real fe_hierarchic_3D_shape(const Elem * elem,
1567 : const Order order,
1568 : const unsigned int i,
1569 : const Point & p,
1570 : const bool add_p_level)
1571 : {
1572 : #if LIBMESH_DIM == 3
1573 :
1574 218716468 : libmesh_assert(elem);
1575 2660333636 : const ElemType type = elem->type();
1576 :
1577 2879050104 : const Order totalorder = order + add_p_level*elem->p_level();
1578 :
1579 2441617168 : switch (type)
1580 : {
1581 912816426 : case HEX8:
1582 : case HEX20:
1583 17128636 : libmesh_assert (T == L2_HIERARCHIC || totalorder < 2);
1584 : libmesh_fallthrough();
1585 : case HEX27:
1586 : {
1587 78858204 : libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u)*(totalorder+1u));
1588 :
1589 : // Compute hex shape functions as a tensor-product
1590 974545994 : Real xi = p(0);
1591 974545994 : Real eta = p(1);
1592 974545994 : Real zeta = p(2);
1593 :
1594 : unsigned int i0, i1, i2;
1595 :
1596 974545994 : cube_indices(elem, totalorder, i, xi, eta, zeta, i0, i1, i2);
1597 :
1598 1053404198 : return (FE<1,T>::shape(EDGE3, totalorder, i0, xi)*
1599 1053404198 : FE<1,T>::shape(EDGE3, totalorder, i1, eta)*
1600 1053404198 : FE<1,T>::shape(EDGE3, totalorder, i2, zeta));
1601 : }
1602 :
1603 886578588 : case PRISM6:
1604 : case PRISM15:
1605 16503294 : libmesh_assert (T == L2_HIERARCHIC || totalorder < 2);
1606 : libmesh_fallthrough();
1607 : case PRISM18:
1608 19362150 : libmesh_assert (T == L2_HIERARCHIC || totalorder < 3);
1609 : libmesh_fallthrough();
1610 : case PRISM20:
1611 : case PRISM21:
1612 : {
1613 85295454 : libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u)*(totalorder+2u)/2u);
1614 :
1615 : // Compute prism shape functions as a tensor-product.
1616 : // Non-const here, because prism_indices might need to do some
1617 : // flips before evaluating edge or face DoFs
1618 955370748 : Point xi_eta {p(0),p(1)};
1619 955370748 : Real zeta = p(2);
1620 :
1621 : unsigned int i01, i2;
1622 :
1623 955370748 : prism_indices(elem, totalorder, i, xi_eta, zeta, i01, i2);
1624 :
1625 : // We'll use the 2D Tri to handle any basis function flipping
1626 : // needed in xi/eta for triangle face+edge bases.
1627 85295454 : Tri3 tri;
1628 :
1629 : // We pinky swear not to modify these nodes
1630 85295454 : Elem & e = const_cast<Elem &>(*elem);
1631 955370748 : if (i2 == 0)
1632 : {
1633 40714860 : tri.set_node(0, e.node_ptr(0));
1634 20357430 : tri.set_node(1, e.node_ptr(1));
1635 20357430 : tri.set_node(2, e.node_ptr(2));
1636 : }
1637 727344312 : else if (i2 == 1)
1638 : {
1639 40714860 : tri.set_node(0, e.node_ptr(3));
1640 20357430 : tri.set_node(1, e.node_ptr(4));
1641 20357430 : tri.set_node(2, e.node_ptr(5));
1642 : }
1643 : else
1644 : {
1645 : // For interior DoFs, no flipping is necessary or done; we
1646 : // can just evaluate on any triangle ... but *not* the
1647 : // obvious 9,10,11 triangle, because that might not exist
1648 : // if we have L2_HIERARCHIC on Prism6.
1649 89161188 : tri.set_node(0, e.node_ptr(0));
1650 44580594 : tri.set_node(1, e.node_ptr(1));
1651 44580594 : tri.set_node(2, e.node_ptr(2));
1652 :
1653 : // For square face DoFs, prism_indices handles flipping,
1654 : // and we *can't* override that in the tri shape call.
1655 499314756 : if (i01 > 2 && i01 < 3u*totalorder)
1656 : {
1657 : // %(p-1) to find the edge number, %2 for even vs odd
1658 295542216 : const bool odd_basis = ((i01-3)%(totalorder-1))%2;
1659 295542216 : if (odd_basis)
1660 : {
1661 103152204 : const int tri_edge = (i01-3)/(totalorder-1);
1662 : // Flip nodes now to avoid triggering a shape
1663 : // function flip later
1664 121571448 : if (tri.point(tri_edge) > tri.point((tri_edge+1)%3))
1665 : {
1666 9745492 : Node * n = tri.node_ptr(tri_edge);
1667 4872746 : tri.set_node(tri_edge, tri.node_ptr((tri_edge+1)%3));
1668 4872746 : tri.set_node((tri_edge+1)%3, n);
1669 : }
1670 : }
1671 : }
1672 : }
1673 :
1674 955370748 : return (FE<2,L2_HIERARCHIC>::shape(&tri, totalorder, i01, xi_eta)*
1675 1040666202 : FE<1,L2_HIERARCHIC>::shape(EDGE2, totalorder, i2, zeta));
1676 : }
1677 :
1678 676997932 : case TET4:
1679 1143848 : libmesh_assert (T == L2_HIERARCHIC || totalorder < 2);
1680 : libmesh_fallthrough();
1681 : case TET10:
1682 3275051 : libmesh_assert (T == L2_HIERARCHIC || totalorder < 3);
1683 : libmesh_fallthrough();
1684 : case TET14:
1685 : {
1686 730416894 : const Real zeta[4] = { 1 - p(0) - p(1) - p(2), p(0), p(1), p(2) };
1687 :
1688 : // Nodal DoFs
1689 730416894 : if (i < 4)
1690 395715748 : return zeta[i];
1691 :
1692 : // Edge DoFs
1693 334701146 : else if (i < 6u*totalorder - 2u)
1694 : {
1695 283191414 : const unsigned int edge_num = (i - 4) / (totalorder - 1u);
1696 : // const int edge_node = edge_num + 4;
1697 283191414 : const unsigned int basisorder = i - 2 - ((totalorder - 1u) * edge_num);
1698 :
1699 283191414 : const unsigned int edgevertex0 = Tet4::edge_nodes_map[edge_num][0],
1700 283191414 : edgevertex1 = Tet4::edge_nodes_map[edge_num][1];
1701 :
1702 : // Get factors to account for edge-flipping
1703 20248452 : Real flip = 1;
1704 319450356 : if (basisorder%2 &&
1705 36258942 : elem->positive_edge_orientation(edge_num))
1706 937446 : flip = -1;
1707 :
1708 283191414 : const Real crossval = zeta[edgevertex0] + zeta[edgevertex1];
1709 283191414 : const Real edgenumerator = zeta[edgevertex1] - zeta[edgevertex0];
1710 :
1711 283191414 : if (crossval == 0.) // Yes, exact comparison; we seem numerically stable otherwise
1712 : {
1713 57913 : unsigned int basisfactorial = 1.;
1714 2574916 : for (unsigned int n=2; n <= basisorder; ++n)
1715 1522573 : basisfactorial *= n;
1716 :
1717 1052343 : return flip * std::pow(edgenumerator, basisorder) / basisfactorial;
1718 : }
1719 :
1720 282139071 : const Real edgeval = edgenumerator / crossval;
1721 20190539 : const Real crossfunc = std::pow(crossval, basisorder);
1722 :
1723 282139071 : return flip * crossfunc *
1724 282139071 : FE<1,HIERARCHIC>::shape(EDGE3, totalorder,
1725 282139071 : basisorder, edgeval);
1726 : }
1727 :
1728 : // Face DoFs
1729 51509732 : else if (i < 2u*totalorder*totalorder + 2u)
1730 : {
1731 48515412 : const int dofs_per_face = (totalorder - 1u) * (totalorder - 2u) / 2;
1732 48515412 : const int face_num = (i - (6u*totalorder - 2u)) / dofs_per_face;
1733 :
1734 48515412 : const std::vector<unsigned int> face_vertex =
1735 : elem->oriented_face_vertices(face_num);
1736 48515412 : const Real zeta0 = zeta[face_vertex[0]],
1737 48515412 : zeta1 = zeta[face_vertex[1]],
1738 48515412 : zeta2 = zeta[face_vertex[2]];
1739 :
1740 48515412 : const unsigned int basisnum =
1741 51011736 : i - 4 -
1742 51011736 : (totalorder - 1u) * /*n_edges*/6 -
1743 48515412 : (dofs_per_face * face_num);
1744 :
1745 48515412 : const unsigned int exp0 = triangular_number_column[basisnum] + 1;
1746 48515412 : const unsigned int exp1 = triangular_number_row[basisnum] + 1 -
1747 : triangular_number_column[basisnum];
1748 :
1749 2496324 : Real returnval = 1;
1750 110172776 : for (unsigned int n = 0; n != exp0; ++n)
1751 61657364 : returnval *= zeta0;
1752 110172776 : for (unsigned int n = 0; n != exp1; ++n)
1753 61657364 : returnval *= zeta1;
1754 48515412 : returnval *= zeta2;
1755 2496324 : return returnval;
1756 : }
1757 :
1758 : // Interior DoFs
1759 : else
1760 : {
1761 2994320 : const unsigned int basisnum = i - 2u*totalorder*totalorder - 2u;
1762 2994320 : const unsigned int exp0 = tetrahedral_number_column[basisnum] + 1;
1763 3290116 : const unsigned int exp1 = tetrahedral_number_row[basisnum] + 1 -
1764 2994320 : tetrahedral_number_column[basisnum] -
1765 2698524 : tetrahedral_number_page[basisnum];
1766 2994320 : const unsigned int exp2 = tetrahedral_number_page[basisnum] + 1;
1767 :
1768 147898 : Real returnval = 1;
1769 5988640 : for (unsigned int n = 0; n != exp0; ++n)
1770 2994320 : returnval *= zeta[0];
1771 5988640 : for (unsigned int n = 0; n != exp1; ++n)
1772 2994320 : returnval *= zeta[1];
1773 5988640 : for (unsigned int n = 0; n != exp2; ++n)
1774 2994320 : returnval *= zeta[2];
1775 2994320 : returnval *= zeta[3];
1776 2994320 : return returnval;
1777 : }
1778 : }
1779 :
1780 0 : default:
1781 0 : libmesh_error_msg("Invalid element type = " << Utility::enum_to_string(type));
1782 : }
1783 :
1784 : #else // LIBMESH_DIM != 3
1785 : libmesh_ignore(elem, order, i, p, add_p_level);
1786 : libmesh_not_implemented();
1787 : #endif
1788 : }
1789 :
1790 :
1791 :
1792 : template <FEFamily T>
1793 91971075 : Real fe_hierarchic_3D_shape_deriv(const Elem * elem,
1794 : const Order order,
1795 : const unsigned int i,
1796 : const unsigned int j,
1797 : const Point & p,
1798 : const bool add_p_level)
1799 : {
1800 1116947178 : return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<3,T>::shape);
1801 : }
1802 :
1803 :
1804 :
1805 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1806 :
1807 : template <FEFamily T>
1808 33592866 : Real fe_hierarchic_3D_shape_second_deriv(const Elem * elem,
1809 : const Order order,
1810 : const unsigned int i,
1811 : const unsigned int j,
1812 : const Point & p,
1813 : const bool add_p_level)
1814 : {
1815 408150528 : return fe_fdm_second_deriv(elem, order, i, j, p, add_p_level,
1816 33592866 : FE<3,T>::shape_deriv);
1817 : }
1818 :
1819 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
1820 :
1821 :
1822 : } // anonymous namespace
|
