Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 : #include "libmesh/libmesh_config.h"
20 :
21 : #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
22 :
23 : // libmesh includes
24 : #include "libmesh/fe.h"
25 : #include "libmesh/elem.h"
26 : #include "libmesh/number_lookups.h"
27 : #include "libmesh/utility.h"
28 : #include "libmesh/enum_to_string.h"
29 :
30 :
31 : namespace {
32 :
33 : using namespace libMesh;
34 :
35 : /*
36 : * Helper function for finding indices, when computing quad shape
37 : * functions and derivatives via a tensor-product of 1D functions and
38 : * derivatives.
39 : */
40 850351704 : std::pair<unsigned int, unsigned int> quad_i0_i1 (const unsigned int i,
41 : const Order totalorder,
42 : const Elem & elem)
43 : {
44 78807455 : libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
45 :
46 : unsigned int i0, i1;
47 :
48 : // Vertex DoFs
49 156484636 : if (i == 0)
50 8761999 : { i0 = 0; i1 = 0; }
51 139086224 : else if (i == 1)
52 8761999 : { i0 = 1; i1 = 0; }
53 121687812 : else if (i == 2)
54 8761999 : { i0 = 1; i1 = 1; }
55 104289400 : else if (i == 3)
56 8761999 : { i0 = 0; i1 = 1; }
57 :
58 :
59 : // Edge DoFs
60 472148088 : else if (i < totalorder + 3u)
61 94366872 : { i0 = i - 2; i1 = 0; }
62 377781216 : else if (i < 2u*totalorder + 2)
63 94366872 : { i0 = 1; i1 = i - totalorder - 1; }
64 283414344 : else if (i < 3u*totalorder + 1)
65 94366872 : { i0 = i - 2u*totalorder; i1 = 1; }
66 189047472 : else if (i < 4u*totalorder)
67 94366872 : { i0 = 0; i1 = i - 3u*totalorder + 1; }
68 : // Interior DoFs
69 : else
70 : {
71 94680600 : unsigned int basisnum = i - 4*totalorder;
72 94680600 : i0 = square_number_column[basisnum] + 2;
73 94680600 : i1 = square_number_row[basisnum] + 2;
74 : }
75 :
76 :
77 : // Flip odd degree of freedom values if necessary
78 : // to keep continuity on sides
79 850351704 : if ((i>= 4 && i<= 4+ totalorder-2u) && elem.point(0) > elem.point(1)) i0=totalorder+2-i0;
80 758671443 : else if ((i>= 4+ totalorder-1u && i<= 4+2*totalorder-3u) && elem.point(1) > elem.point(2)) i1=totalorder+2-i1;
81 666664644 : else if ((i>= 4+2*totalorder-2u && i<= 4+3*totalorder-4u) && elem.point(3) > elem.point(2)) i0=totalorder+2-i0;
82 574991673 : else if ((i>= 4+3*totalorder-3u && i<= 4+4*totalorder-5u) && elem.point(0) > elem.point(3)) i1=totalorder+2-i1;
83 :
84 929159159 : return std::make_pair(i0, i1);
85 : }
86 :
87 : }
88 :
89 : namespace libMesh
90 : {
91 :
92 :
93 15336884 : LIBMESH_DEFAULT_VECTORIZED_FE(2,BERNSTEIN)
94 :
95 :
96 : template <>
97 508028966 : Real FE<2,BERNSTEIN>::shape(const Elem * elem,
98 : const Order order,
99 : const unsigned int i,
100 : const Point & p,
101 : const bool add_p_level)
102 : {
103 47887451 : libmesh_assert(elem);
104 :
105 508028966 : const ElemType type = elem->type();
106 :
107 : const Order totalorder =
108 554786143 : order + add_p_level*elem->p_level();
109 :
110 : // Declare that we are using our own special power function
111 : // from the Utility namespace. This saves typing later.
112 : using Utility::pow;
113 :
114 461271789 : switch (type)
115 : {
116 : // Hierarchic shape functions on the quadrilateral.
117 498373206 : case QUAD4:
118 : case QUADSHELL4:
119 : case QUAD9:
120 : case QUADSHELL9:
121 : {
122 : // Compute quad shape functions as a tensor-product
123 498373206 : auto [i0, i1] = quad_i0_i1(i, totalorder, *elem);
124 :
125 544272650 : return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0, p(0))*
126 544272650 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1, p(1)));
127 : }
128 : // handle serendipity QUAD8 element separately
129 1761792 : case QUAD8:
130 : case QUADSHELL8:
131 : {
132 146816 : libmesh_assert_less (totalorder, 3);
133 :
134 1761792 : const Real xi = p(0);
135 1761792 : const Real eta = p(1);
136 :
137 146816 : libmesh_assert_less (i, 8);
138 :
139 : // 0 1 2 3 4 5 6 7 8
140 : static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
141 : static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
142 : static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5};
143 :
144 : //B_t,i0(i)|xi * B_s,i1(i)|eta + scal(i) * B_t,2|xi * B_t,2|eta
145 1761792 : return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)*
146 1761792 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta)
147 1908608 : +scal[i]*
148 1761792 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[8], xi)*
149 1761792 : FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[8], eta));
150 :
151 : }
152 :
153 7334638 : case TRI3:
154 : case TRISHELL3:
155 151587 : libmesh_assert_less (totalorder, 2);
156 : libmesh_fallthrough();
157 : case TRI6:
158 : case TRI7:
159 7893968 : switch (totalorder)
160 : {
161 1669752 : case FIRST:
162 : {
163 1669752 : const Real x=p(0);
164 1669752 : const Real y=p(1);
165 1669752 : const Real r=1.-x-y;
166 :
167 151587 : libmesh_assert_less (i, 3);
168 :
169 1669752 : switch(i)
170 : {
171 50529 : case 0: return r; //f0,0,1
172 556584 : case 1: return x; //f0,1,1
173 556584 : case 2: return y; //f1,0,1
174 :
175 0 : default:
176 0 : libmesh_error_msg("Invalid shape function index i = " << i);
177 : }
178 : }
179 1059216 : case SECOND:
180 : {
181 1059216 : const Real x=p(0);
182 1059216 : const Real y=p(1);
183 1059216 : const Real r=1.-x-y;
184 :
185 95340 : libmesh_assert_less (i, 6);
186 :
187 1059216 : switch(i)
188 : {
189 176536 : case 0: return r*r;
190 176536 : case 1: return x*x;
191 176536 : case 2: return y*y;
192 :
193 176536 : case 3: return 2.*x*r;
194 176536 : case 4: return 2.*x*y;
195 176536 : case 5: return 2.*r*y;
196 :
197 0 : default:
198 0 : libmesh_error_msg("Invalid shape function index i = " << i);
199 : }
200 : }
201 1944200 : case THIRD:
202 : {
203 1944200 : const Real x=p(0);
204 1944200 : const Real y=p(1);
205 1944200 : const Real r=1.-x-y;
206 174700 : libmesh_assert_less (i, 10);
207 :
208 174700 : unsigned int shape=i;
209 :
210 :
211 1944200 : if ((i==3||i==4) && elem->positive_edge_orientation(0)) shape=7-i;
212 1944200 : if ((i==5||i==6) && elem->positive_edge_orientation(1)) shape=11-i;
213 1944200 : if ((i==7||i==8) && !elem->positive_edge_orientation(2)) shape=15-i;
214 :
215 1944200 : switch(shape)
216 : {
217 194420 : case 0: return r*r*r;
218 194420 : case 1: return x*x*x;
219 194420 : case 2: return y*y*y;
220 :
221 194420 : case 3: return 3.*x*r*r;
222 194420 : case 4: return 3.*x*x*r;
223 :
224 194420 : case 5: return 3.*y*x*x;
225 194420 : case 6: return 3.*y*y*x;
226 :
227 194420 : case 7: return 3.*y*r*r;
228 194420 : case 8: return 3.*y*y*r;
229 :
230 194420 : case 9: return 6.*x*y*r;
231 :
232 0 : default:
233 0 : libmesh_error_msg("Invalid shape function index shape = " << shape);
234 : }
235 : }
236 3220800 : case FOURTH:
237 : {
238 3220800 : const Real x=p(0);
239 3220800 : const Real y=p(1);
240 3220800 : const Real r=1-x-y;
241 289290 : unsigned int shape=i;
242 :
243 289290 : libmesh_assert_less (i, 15);
244 :
245 3220800 : if ((i==3||i== 5) && elem->positive_edge_orientation(0)) shape=8-i;
246 3220800 : if ((i==6||i== 8) && elem->positive_edge_orientation(1)) shape=14-i;
247 3220800 : if ((i==9||i==11) && !elem->positive_edge_orientation(2)) shape=20-i;
248 :
249 :
250 3220800 : switch(shape)
251 : {
252 : // point functions
253 214720 : case 0: return r*r*r*r;
254 214720 : case 1: return x*x*x*x;
255 214720 : case 2: return y*y*y*y;
256 :
257 : // edge functions
258 214720 : case 3: return 4.*x*r*r*r;
259 214720 : case 4: return 6.*x*x*r*r;
260 214720 : case 5: return 4.*x*x*x*r;
261 :
262 214720 : case 6: return 4.*y*x*x*x;
263 214720 : case 7: return 6.*y*y*x*x;
264 214720 : case 8: return 4.*y*y*y*x;
265 :
266 214720 : case 9: return 4.*y*r*r*r;
267 214720 : case 10: return 6.*y*y*r*r;
268 214720 : case 11: return 4.*y*y*y*r;
269 :
270 : // inner functions
271 214720 : case 12: return 12.*x*y*r*r;
272 214720 : case 13: return 12.*x*x*y*r;
273 214720 : case 14: return 12.*x*y*y*r;
274 :
275 0 : default:
276 0 : libmesh_error_msg("Invalid shape function index shape = " << shape);
277 : }
278 : }
279 0 : case FIFTH:
280 : {
281 0 : const Real x=p(0);
282 0 : const Real y=p(1);
283 0 : const Real r=1-x-y;
284 0 : unsigned int shape=i;
285 :
286 0 : libmesh_assert_less (i, 21);
287 :
288 0 : if ((i>= 3&&i<= 6) && elem->positive_edge_orientation(0)) shape=9-i;
289 0 : if ((i>= 7&&i<=10) && elem->positive_edge_orientation(1)) shape=17-i;
290 0 : if ((i>=11&&i<=14) && !elem->positive_edge_orientation(2)) shape=25-i;
291 :
292 0 : switch(shape)
293 : {
294 : //point functions
295 0 : case 0: return pow<5>(r);
296 0 : case 1: return pow<5>(x);
297 0 : case 2: return pow<5>(y);
298 :
299 : //edge functions
300 0 : case 3: return 5.*x *pow<4>(r);
301 0 : case 4: return 10.*pow<2>(x)*pow<3>(r);
302 0 : case 5: return 10.*pow<3>(x)*pow<2>(r);
303 0 : case 6: return 5.*pow<4>(x)*r;
304 :
305 0 : case 7: return 5.*y *pow<4>(x);
306 0 : case 8: return 10.*pow<2>(y)*pow<3>(x);
307 0 : case 9: return 10.*pow<3>(y)*pow<2>(x);
308 0 : case 10: return 5.*pow<4>(y)*x;
309 :
310 0 : case 11: return 5.*y *pow<4>(r);
311 0 : case 12: return 10.*pow<2>(y)*pow<3>(r);
312 0 : case 13: return 10.*pow<3>(y)*pow<2>(r);
313 0 : case 14: return 5.*pow<4>(y)*r;
314 :
315 : //inner functions
316 0 : case 15: return 20.*x*y*pow<3>(r);
317 0 : case 16: return 30.*x*pow<2>(y)*pow<2>(r);
318 0 : case 17: return 30.*pow<2>(x)*y*pow<2>(r);
319 0 : case 18: return 20.*x*pow<3>(y)*r;
320 0 : case 19: return 20.*pow<3>(x)*y*r;
321 0 : case 20: return 30.*pow<2>(x)*pow<2>(y)*r;
322 :
323 0 : default:
324 0 : libmesh_error_msg("Invalid shape function index shape = " << shape);
325 : }
326 : }
327 0 : case SIXTH:
328 : {
329 0 : const Real x=p(0);
330 0 : const Real y=p(1);
331 0 : const Real r=1-x-y;
332 0 : unsigned int shape=i;
333 :
334 0 : libmesh_assert_less (i, 28);
335 :
336 0 : if ((i>= 3&&i<= 7) && elem->positive_edge_orientation(0)) shape=10-i;
337 0 : if ((i>= 8&&i<=12) && elem->positive_edge_orientation(1)) shape=20-i;
338 0 : if ((i>=13&&i<=17) && !elem->positive_edge_orientation(2)) shape=30-i;
339 :
340 0 : switch(shape)
341 : {
342 : //point functions
343 0 : case 0: return pow<6>(r);
344 0 : case 1: return pow<6>(x);
345 0 : case 2: return pow<6>(y);
346 :
347 : //edge functions
348 0 : case 3: return 6.*x *pow<5>(r);
349 0 : case 4: return 15.*pow<2>(x)*pow<4>(r);
350 0 : case 5: return 20.*pow<3>(x)*pow<3>(r);
351 0 : case 6: return 15.*pow<4>(x)*pow<2>(r);
352 0 : case 7: return 6.*pow<5>(x)*r;
353 :
354 0 : case 8: return 6.*y *pow<5>(x);
355 0 : case 9: return 15.*pow<2>(y)*pow<4>(x);
356 0 : case 10: return 20.*pow<3>(y)*pow<3>(x);
357 0 : case 11: return 15.*pow<4>(y)*pow<2>(x);
358 0 : case 12: return 6.*pow<5>(y)*x;
359 :
360 0 : case 13: return 6.*y *pow<5>(r);
361 0 : case 14: return 15.*pow<2>(y)*pow<4>(r);
362 0 : case 15: return 20.*pow<3>(y)*pow<3>(r);
363 0 : case 16: return 15.*pow<4>(y)*pow<2>(r);
364 0 : case 17: return 6.*pow<5>(y)*r;
365 :
366 : //inner functions
367 0 : case 18: return 30.*x*y*pow<4>(r);
368 0 : case 19: return 60.*x*pow<2>(y)*pow<3>(r);
369 0 : case 20: return 60.* pow<2>(x)*y*pow<3>(r);
370 0 : case 21: return 60.*x*pow<3>(y)*pow<2>(r);
371 0 : case 22: return 60.*pow<3>(x)*y*pow<2>(r);
372 0 : case 23: return 90.*pow<2>(x)*pow<2>(y)*pow<2>(r);
373 0 : case 24: return 30.*x*pow<4>(y)*r;
374 0 : case 25: return 60.*pow<2>(x)*pow<3>(y)*r;
375 0 : case 26: return 60.*pow<3>(x)*pow<2>(y)*r;
376 0 : case 27: return 30.*pow<4>(x)*y*r;
377 :
378 0 : default:
379 0 : libmesh_error_msg("Invalid shape function index shape = " << shape);
380 : } // switch shape
381 : }
382 0 : default:
383 0 : libmesh_error_msg("Invalid totalorder = " << totalorder);
384 : } // switch order
385 :
386 0 : default:
387 0 : libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
388 : } // switch type
389 : }
390 :
391 :
392 : template <>
393 0 : Real FE<2,BERNSTEIN>::shape(const ElemType,
394 : const Order,
395 : const unsigned int,
396 : const Point &)
397 : {
398 0 : libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge orientation is needed.");
399 : return 0.;
400 : }
401 :
402 :
403 : template <>
404 0 : Real FE<2,BERNSTEIN>::shape(const FEType fet,
405 : const Elem * elem,
406 : const unsigned int i,
407 : const Point & p,
408 : const bool add_p_level)
409 : {
410 0 : return FE<2,BERNSTEIN>::shape(elem, fet.order, i, p, add_p_level);
411 : }
412 :
413 :
414 :
415 : template <>
416 353827214 : Real FE<2,BERNSTEIN>::shape_deriv(const Elem * elem,
417 : const Order order,
418 : const unsigned int i,
419 : const unsigned int j,
420 : const Point & p,
421 : const bool add_p_level)
422 : {
423 31956046 : libmesh_assert(elem);
424 :
425 353827214 : const ElemType type = elem->type();
426 :
427 : const Order totalorder =
428 385783260 : order + add_p_level*elem->p_level();
429 :
430 321871168 : switch (type)
431 : {
432 : // Hierarchic shape functions on the quadrilateral.
433 349244190 : case QUAD4:
434 : case QUAD9:
435 : case QUADSHELL9:
436 : {
437 : // Compute quad shape functions as a tensor-product
438 349244190 : auto [i0, i1] = quad_i0_i1(i, totalorder, *elem);
439 :
440 349244190 : switch (j)
441 : {
442 : // d()/dxi
443 174622095 : case 0:
444 190397034 : return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0, 0, p(0))*
445 190397034 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1, p(1)));
446 :
447 : // d()/deta
448 174622095 : case 1:
449 190397034 : return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0, p(0))*
450 190397034 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1, 0, p(1)));
451 :
452 0 : default:
453 0 : libmesh_error_msg("Invalid shape function derivative j = " << j);
454 : }
455 : }
456 :
457 : // Bernstein shape functions on the 8-noded quadrilateral
458 : // is handled separately.
459 1101312 : case QUAD8:
460 : case QUADSHELL8:
461 : {
462 91776 : libmesh_assert_less (totalorder, 3);
463 :
464 1101312 : const Real xi = p(0);
465 1101312 : const Real eta = p(1);
466 :
467 91776 : libmesh_assert_less (i, 8);
468 :
469 : // 0 1 2 3 4 5 6 7 8
470 : static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
471 : static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
472 : static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5};
473 1101312 : switch (j)
474 : {
475 : // d()/dxi
476 45888 : case 0:
477 550656 : return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
478 550656 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)
479 596544 : +scal[i]*
480 550656 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[8], 0, xi)*
481 550656 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[8], eta));
482 :
483 : // d()/deta
484 45888 : case 1:
485 550656 : return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
486 550656 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)
487 596544 : +scal[i]*
488 550656 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[8], xi)*
489 550656 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[8], 0, eta));
490 :
491 0 : default:
492 0 : libmesh_error_msg("Invalid shape function derivative j = " << j);
493 : }
494 : }
495 :
496 3232960 : case TRI3:
497 : case TRISHELL3:
498 65640 : libmesh_assert_less (totalorder, 2);
499 : libmesh_fallthrough();
500 : case TRI6:
501 : case TRI7:
502 : {
503 3481712 : return fe_fdm_deriv(elem, order, i, j, p, add_p_level, FE<2,BERNSTEIN>::shape);
504 : }
505 :
506 0 : default:
507 0 : libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
508 : }
509 : }
510 :
511 :
512 :
513 :
514 : template <>
515 0 : Real FE<2,BERNSTEIN>::shape_deriv(const ElemType,
516 : const Order,
517 : const unsigned int,
518 : const unsigned int,
519 : const Point &)
520 : {
521 0 : libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge orientation is needed.");
522 : return 0.;
523 : }
524 :
525 : template <>
526 0 : Real FE<2,BERNSTEIN>::shape_deriv(const FEType fet,
527 : const Elem * elem,
528 : const unsigned int i,
529 : const unsigned int j,
530 : const Point & p,
531 : const bool add_p_level)
532 : {
533 0 : return FE<2,BERNSTEIN>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
534 : }
535 :
536 :
537 :
538 :
539 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
540 :
541 :
542 :
543 : template <>
544 3983124 : Real FE<2,BERNSTEIN>::shape_second_deriv(const Elem * elem,
545 : const Order order,
546 : const unsigned int i,
547 : const unsigned int j,
548 : const Point & p,
549 : const bool add_p_level)
550 : {
551 340590 : libmesh_assert(elem);
552 :
553 3983124 : const ElemType type = elem->type();
554 :
555 : const Order totalorder =
556 4323714 : order + add_p_level*elem->p_level();
557 :
558 3642534 : switch (type)
559 : {
560 : // Hierarchic shape functions on the quadrilateral.
561 2734308 : case QUAD4:
562 : case QUAD9:
563 : case QUADSHELL9:
564 : {
565 : // Compute quad shape functions as a tensor-product
566 2734308 : auto [i0, i1] = quad_i0_i1(i, totalorder, *elem);
567 :
568 2734308 : switch (j)
569 : {
570 : // d^2() / dxi^2
571 911436 : case 0:
572 987389 : return (FE<1,BERNSTEIN>::shape_second_deriv(EDGE3, totalorder, i0, 0, p(0))*
573 987389 : FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1, p(1)));
574 :
575 : // d^2() / dxi deta
576 911436 : case 1:
577 987389 : return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0, 0, p(0))*
578 987389 : FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1, 0, p(1)));
579 :
580 : // d^2() / deta^2
581 911436 : case 2:
582 987389 : return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0, p(0))*
583 987389 : FE<1,BERNSTEIN>::shape_second_deriv(EDGE3, totalorder, i1, 0, p(1)));
584 :
585 0 : default:
586 0 : libmesh_error_msg("Invalid shape function derivative j = " << j);
587 : }
588 : }
589 :
590 : // Going to be lazy again about the hard cases.
591 1155534 : case TRI3:
592 : case TRISHELL3:
593 19449 : libmesh_assert_less (totalorder, 2);
594 : libmesh_fallthrough();
595 : case QUAD8:
596 : case QUADSHELL8:
597 : case TRI6:
598 : case TRI7:
599 : {
600 1248816 : return fe_fdm_second_deriv(elem, order, i, j, p, add_p_level,
601 1248816 : FE<2,BERNSTEIN>::shape_deriv);
602 : }
603 :
604 0 : default:
605 0 : libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
606 : }
607 : }
608 :
609 :
610 : template <>
611 0 : Real FE<2,BERNSTEIN>::shape_second_deriv(const ElemType,
612 : const Order,
613 : const unsigned int,
614 : const unsigned int,
615 : const Point &)
616 : {
617 0 : libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge orientation is needed.");
618 : return 0.;
619 : }
620 :
621 :
622 : template <>
623 0 : Real FE<2,BERNSTEIN>::shape_second_deriv(const FEType fet,
624 : const Elem * elem,
625 : const unsigned int i,
626 : const unsigned int j,
627 : const Point & p,
628 : const bool add_p_level)
629 : {
630 0 : libmesh_assert(elem);
631 0 : return FE<2,BERNSTEIN>::shape_second_deriv(elem, fet.order, i, j, p, add_p_level);
632 : }
633 :
634 :
635 :
636 : #endif
637 :
638 : } // namespace libMesh
639 :
640 :
641 : #endif// LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
|
