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 : // Local includes
20 : #include "libmesh/fe.h"
21 : #include "libmesh/elem.h"
22 : #include "libmesh/fe_interface.h"
23 : #include "libmesh/number_lookups.h"
24 : #include "libmesh/enum_to_string.h"
25 :
26 : namespace
27 : {
28 : using namespace libMesh;
29 :
30 : // Compute the static coefficients for an element
31 39923712 : void hermite_compute_coefs(const Elem * elem, std::vector<std::vector<Real>> & dxdxi
32 :
33 : #ifdef DEBUG
34 : , std::vector<Real> & dydxi, std::vector<Real> & dzdeta, std::vector<Real> & dxdzeta,
35 : std::vector<Real> & dzdxi, std::vector<Real> & dxdeta, std::vector<Real> & dydzeta
36 : #endif
37 : )
38 : {
39 39923712 : const FEFamily mapping_family = FEMap::map_fe_type(*elem);
40 39923712 : const FEType map_fe_type(elem->default_order(), mapping_family);
41 :
42 : // Note: we explicitly don't consider the elem->p_level() when
43 : // computing the number of mapping shape functions.
44 : const int n_mapping_shape_functions =
45 39923712 : FEInterface::n_shape_functions (map_fe_type, /*extra_order=*/0, elem);
46 :
47 39923712 : static const Point dofpt[2] = {Point(-1,-1,-1), Point(1,1,1)};
48 :
49 : FEInterface::shape_deriv_ptr shape_deriv_ptr =
50 39923712 : FEInterface::shape_deriv_function(map_fe_type, elem);
51 :
52 119771136 : for (int p = 0; p != 2; ++p)
53 : {
54 79847424 : dxdxi[0][p] = 0;
55 79847424 : dxdxi[1][p] = 0;
56 86501376 : dxdxi[2][p] = 0;
57 : #ifdef DEBUG
58 6653952 : dydxi[p] = 0;
59 6653952 : dzdeta[p] = 0;
60 6653952 : dxdzeta[p] = 0;
61 6653952 : dzdxi[p] = 0;
62 6653952 : dxdeta[p] = 0;
63 6653952 : dydzeta[p] = 0;
64 : #endif
65 2235727872 : for (int i = 0; i != n_mapping_shape_functions; ++i)
66 : {
67 : const Real ddxi = shape_deriv_ptr
68 2155880448 : (map_fe_type, elem, i, 0, dofpt[p], /*add_p_level=*/false);
69 : const Real ddeta = shape_deriv_ptr
70 2155880448 : (map_fe_type, elem, i, 1, dofpt[p], /*add_p_level=*/false);
71 : const Real ddzeta = shape_deriv_ptr
72 2155880448 : (map_fe_type, elem, i, 2, dofpt[p], /*add_p_level=*/false);
73 :
74 : // dxdeta, dxdzeta, dydxi, dydzeta, dzdxi, dzdeta should all
75 : // be 0!
76 359313408 : const Point & point_i = elem->point(i);
77 2155880448 : dxdxi[0][p] += point_i(0) * ddxi;
78 2155880448 : dxdxi[1][p] += point_i(1) * ddeta;
79 2335537152 : dxdxi[2][p] += point_i(2) * ddzeta;
80 : #ifdef DEBUG
81 179656704 : dydxi[p] += point_i(1) * ddxi;
82 179656704 : dzdeta[p] += point_i(2) * ddeta;
83 179656704 : dxdzeta[p] += point_i(0) * ddzeta;
84 179656704 : dzdxi[p] += point_i(2) * ddxi;
85 179656704 : dxdeta[p] += point_i(0) * ddeta;
86 179656704 : dydzeta[p] += point_i(1) * ddzeta;
87 : #endif
88 : }
89 :
90 : // No singular elements!
91 6653952 : libmesh_assert(dxdxi[0][p]);
92 6653952 : libmesh_assert(dxdxi[1][p]);
93 6653952 : libmesh_assert(dxdxi[2][p]);
94 : // No non-rectilinear or non-axis-aligned elements!
95 : #ifdef DEBUG
96 6653952 : libmesh_assert_less (std::abs(dydxi[p]), TOLERANCE);
97 6653952 : libmesh_assert_less (std::abs(dzdeta[p]), TOLERANCE);
98 6653952 : libmesh_assert_less (std::abs(dxdzeta[p]), TOLERANCE);
99 6653952 : libmesh_assert_less (std::abs(dzdxi[p]), TOLERANCE);
100 6653952 : libmesh_assert_less (std::abs(dxdeta[p]), TOLERANCE);
101 6653952 : libmesh_assert_less (std::abs(dydzeta[p]), TOLERANCE);
102 : #endif
103 : }
104 39923712 : }
105 :
106 :
107 :
108 39923712 : Real hermite_bases_3D (std::vector<unsigned int> & bases1D,
109 : const std::vector<std::vector<Real>> & dxdxi,
110 : const Order & o,
111 : unsigned int i)
112 : {
113 3326976 : bases1D.clear();
114 39923712 : bases1D.resize(3,0);
115 3326976 : Real coef = 1.0;
116 :
117 39923712 : unsigned int e = o-2;
118 :
119 : // Nodes
120 39923712 : if (i < 64)
121 : {
122 39923712 : switch (i / 8)
123 : {
124 415872 : case 0:
125 415872 : break;
126 831744 : case 1:
127 4990464 : bases1D[0] = 1;
128 4990464 : break;
129 831744 : case 2:
130 4990464 : bases1D[0] = 1;
131 4990464 : bases1D[1] = 1;
132 4990464 : break;
133 831744 : case 3:
134 4990464 : bases1D[1] = 1;
135 4990464 : break;
136 831744 : case 4:
137 4990464 : bases1D[2] = 1;
138 4990464 : break;
139 831744 : case 5:
140 4990464 : bases1D[0] = 1;
141 4990464 : bases1D[2] = 1;
142 4990464 : break;
143 831744 : case 6:
144 4990464 : bases1D[0] = 1;
145 4990464 : bases1D[1] = 1;
146 4990464 : bases1D[2] = 1;
147 4990464 : break;
148 831744 : case 7:
149 4990464 : bases1D[1] = 1;
150 4990464 : bases1D[2] = 1;
151 4990464 : break;
152 0 : default:
153 0 : libmesh_error_msg("Invalid basis node = " << i/8);
154 : }
155 :
156 39923712 : unsigned int basisnum = i%8;
157 39923712 : switch (basisnum) // DoF type
158 : {
159 415872 : case 0: // DoF = value at node
160 415872 : coef = 1.0;
161 415872 : break;
162 831744 : case 1: // DoF = x derivative at node
163 4990464 : coef = dxdxi[0][bases1D[0]];
164 4990464 : bases1D[0] += 2; break;
165 831744 : case 2: // DoF = y derivative at node
166 4990464 : coef = dxdxi[1][bases1D[1]];
167 4990464 : bases1D[1] += 2; break;
168 831744 : case 3: // DoF = xy derivative at node
169 5406336 : coef = dxdxi[0][bases1D[0]] * dxdxi[1][bases1D[1]];
170 4990464 : bases1D[0] += 2; bases1D[1] += 2; break;
171 831744 : case 4: // DoF = z derivative at node
172 4990464 : coef = dxdxi[2][bases1D[2]];
173 4990464 : bases1D[2] += 2; break;
174 831744 : case 5: // DoF = xz derivative at node
175 5406336 : coef = dxdxi[0][bases1D[0]] * dxdxi[2][bases1D[2]];
176 4990464 : bases1D[0] += 2; bases1D[2] += 2; break;
177 831744 : case 6: // DoF = yz derivative at node
178 5406336 : coef = dxdxi[1][bases1D[1]] * dxdxi[2][bases1D[2]];
179 4990464 : bases1D[1] += 2; bases1D[2] += 2; break;
180 831744 : case 7: // DoF = xyz derivative at node
181 5822208 : coef = dxdxi[0][bases1D[0]] * dxdxi[1][bases1D[1]] * dxdxi[2][bases1D[2]];
182 4990464 : bases1D[0] += 2; bases1D[1] += 2; bases1D[2] += 2; break;
183 0 : default:
184 0 : libmesh_error_msg("Invalid basisnum = " << basisnum);
185 : }
186 : }
187 : // Edges
188 0 : else if (i < 64 + 12*4*e)
189 : {
190 0 : unsigned int basisnum = (i - 64) % (4*e);
191 0 : switch ((i - 64) / (4*e))
192 : {
193 0 : case 0:
194 0 : bases1D[0] = basisnum / 4 + 4;
195 0 : bases1D[1] = basisnum % 4 / 2 * 2;
196 0 : bases1D[2] = basisnum % 2 * 2;
197 0 : if (basisnum % 4 / 2)
198 0 : coef *= dxdxi[1][0];
199 0 : if (basisnum % 2)
200 0 : coef *= dxdxi[2][0];
201 0 : break;
202 0 : case 1:
203 0 : bases1D[0] = basisnum % 4 / 2 * 2 + 1;
204 0 : bases1D[1] = basisnum / 4 + 4;
205 0 : bases1D[2] = basisnum % 2 * 2;
206 0 : if (basisnum % 4 / 2)
207 0 : coef *= dxdxi[0][1];
208 0 : if (basisnum % 2)
209 0 : coef *= dxdxi[2][0];
210 0 : break;
211 0 : case 2:
212 0 : bases1D[0] = basisnum / 4 + 4;
213 0 : bases1D[1] = basisnum % 4 / 2 * 2 + 1;
214 0 : bases1D[2] = basisnum % 2 * 2;
215 0 : if (basisnum % 4 / 2)
216 0 : coef *= dxdxi[1][1];
217 0 : if (basisnum % 2)
218 0 : coef *= dxdxi[2][0];
219 0 : break;
220 0 : case 3:
221 0 : bases1D[0] = basisnum % 4 / 2 * 2;
222 0 : bases1D[1] = basisnum / 4 + 4;
223 0 : bases1D[2] = basisnum % 2 * 2;
224 0 : if (basisnum % 4 / 2)
225 0 : coef *= dxdxi[0][0];
226 0 : if (basisnum % 2)
227 0 : coef *= dxdxi[2][0];
228 0 : break;
229 0 : case 4:
230 0 : bases1D[0] = basisnum % 4 / 2 * 2;
231 0 : bases1D[1] = basisnum % 2 * 2;
232 0 : bases1D[2] = basisnum / 4 + 4;
233 0 : if (basisnum % 4 / 2)
234 0 : coef *= dxdxi[0][0];
235 0 : if (basisnum % 2)
236 0 : coef *= dxdxi[1][0];
237 0 : break;
238 0 : case 5:
239 0 : bases1D[0] = basisnum % 4 / 2 * 2 + 1;
240 0 : bases1D[1] = basisnum % 2 * 2;
241 0 : bases1D[2] = basisnum / 4 + 4;
242 0 : if (basisnum % 4 / 2)
243 0 : coef *= dxdxi[0][1];
244 0 : if (basisnum % 2)
245 0 : coef *= dxdxi[1][0];
246 0 : break;
247 0 : case 6:
248 0 : bases1D[0] = basisnum % 4 / 2 * 2 + 1;
249 0 : bases1D[1] = basisnum % 2 * 2 + 1;
250 0 : bases1D[2] = basisnum / 4 + 4;
251 0 : if (basisnum % 4 / 2)
252 0 : coef *= dxdxi[0][1];
253 0 : if (basisnum % 2)
254 0 : coef *= dxdxi[1][1];
255 0 : break;
256 0 : case 7:
257 0 : bases1D[0] = basisnum % 4 / 2 * 2;
258 0 : bases1D[1] = basisnum % 2 * 2 + 1;
259 0 : bases1D[2] = basisnum / 4 + 4;
260 0 : if (basisnum % 4 / 2)
261 0 : coef *= dxdxi[0][0];
262 0 : if (basisnum % 2)
263 0 : coef *= dxdxi[1][1];
264 0 : break;
265 0 : case 8:
266 0 : bases1D[0] = basisnum / 4 + 4;
267 0 : bases1D[1] = basisnum % 4 / 2 * 2;
268 0 : bases1D[2] = basisnum % 2 * 2 + 1;
269 0 : if (basisnum % 4 / 2)
270 0 : coef *= dxdxi[1][0];
271 0 : if (basisnum % 2)
272 0 : coef *= dxdxi[2][1];
273 0 : break;
274 0 : case 9:
275 0 : bases1D[0] = basisnum % 4 / 2 * 2 + 1;
276 0 : bases1D[1] = basisnum / 4 + 4;
277 0 : bases1D[2] = basisnum % 2 * 2;
278 0 : if (basisnum % 4 / 2)
279 0 : coef *= dxdxi[0][1];
280 0 : if (basisnum % 2)
281 0 : coef *= dxdxi[2][1];
282 0 : break;
283 0 : case 10:
284 0 : bases1D[0] = basisnum / 4 + 4;
285 0 : bases1D[1] = basisnum % 4 / 2 * 2 + 1;
286 0 : bases1D[2] = basisnum % 2 * 2 + 1;
287 0 : if (basisnum % 4 / 2)
288 0 : coef *= dxdxi[1][1];
289 0 : if (basisnum % 2)
290 0 : coef *= dxdxi[2][1];
291 0 : break;
292 0 : case 11:
293 0 : bases1D[0] = basisnum % 4 / 2 * 2;
294 0 : bases1D[1] = basisnum / 4 + 4;
295 0 : bases1D[2] = basisnum % 2 * 2 + 1;
296 0 : if (basisnum % 4 / 2)
297 0 : coef *= dxdxi[0][0];
298 0 : if (basisnum % 2)
299 0 : coef *= dxdxi[2][1];
300 0 : break;
301 0 : default:
302 0 : libmesh_error_msg("Invalid basis node = " << (i - 64) / (4*e));
303 : }
304 : }
305 : // Faces
306 0 : else if (i < 64 + 12*4*e + 6*2*e*e)
307 : {
308 0 : unsigned int basisnum = (i - 64 - 12*4*e) % (2*e*e);
309 0 : switch ((i - 64 - 12*4*e) / (2*e*e))
310 : {
311 0 : case 0:
312 0 : bases1D[0] = square_number_column[basisnum / 2];
313 0 : bases1D[1] = square_number_row[basisnum / 2];
314 0 : bases1D[2] = basisnum % 2 * 2;
315 0 : if (basisnum % 2)
316 0 : coef *= dxdxi[2][0];
317 0 : break;
318 0 : case 1:
319 0 : bases1D[0] = square_number_column[basisnum / 2];
320 0 : bases1D[1] = basisnum % 2 * 2;
321 0 : bases1D[2] = square_number_row[basisnum / 2];
322 0 : if (basisnum % 2)
323 0 : coef *= dxdxi[1][0];
324 0 : break;
325 0 : case 2:
326 0 : bases1D[0] = basisnum % 2 * 2 + 1;
327 0 : bases1D[1] = square_number_column[basisnum / 2];
328 0 : bases1D[2] = square_number_row[basisnum / 2];
329 0 : if (basisnum % 2)
330 0 : coef *= dxdxi[0][1];
331 0 : break;
332 0 : case 3:
333 0 : bases1D[0] = square_number_column[basisnum / 2];
334 0 : bases1D[1] = basisnum % 2 * 2 + 1;
335 0 : bases1D[2] = square_number_row[basisnum / 2];
336 0 : if (basisnum % 2)
337 0 : coef *= dxdxi[1][1];
338 0 : break;
339 0 : case 4:
340 0 : bases1D[0] = basisnum % 2 * 2;
341 0 : bases1D[1] = square_number_column[basisnum / 2];
342 0 : bases1D[2] = square_number_row[basisnum / 2];
343 0 : if (basisnum % 2)
344 0 : coef *= dxdxi[0][0];
345 0 : break;
346 0 : case 5:
347 0 : bases1D[0] = square_number_column[basisnum / 2];
348 0 : bases1D[1] = square_number_row[basisnum / 2];
349 0 : bases1D[2] = basisnum % 2 * 2 + 1;
350 0 : if (basisnum % 2)
351 0 : coef *= dxdxi[2][1];
352 0 : break;
353 0 : default:
354 0 : libmesh_error_msg("Invalid basis node = " << (i - 64 - 12*4*e) / (2*e*e));
355 : }
356 : }
357 : // Interior
358 : else
359 : {
360 0 : unsigned int basisnum = i - 64 - 12*4*e;
361 0 : bases1D[0] = cube_number_column[basisnum] + 4;
362 0 : bases1D[1] = cube_number_row[basisnum] + 4;
363 0 : bases1D[2] = cube_number_page[basisnum] + 4;
364 : }
365 :
366 : // No singular elements
367 3326976 : libmesh_assert(coef);
368 39923712 : return coef;
369 : }
370 :
371 :
372 : } // end anonymous namespace
373 :
374 :
375 : namespace libMesh
376 : {
377 :
378 :
379 15490068 : LIBMESH_DEFAULT_VECTORIZED_FE(3,HERMITE)
380 :
381 :
382 : template <>
383 4036608 : Real FE<3,HERMITE>::shape(const Elem * elem,
384 : const Order order,
385 : const unsigned int i,
386 : const Point & p,
387 : const bool add_p_level)
388 : {
389 336384 : libmesh_assert(elem);
390 :
391 5045760 : std::vector<std::vector<Real>> dxdxi(3, std::vector<Real>(2, 0));
392 :
393 : #ifdef DEBUG
394 1345536 : std::vector<Real> dydxi(2), dzdeta(2), dxdzeta(2);
395 1345536 : std::vector<Real> dzdxi(2), dxdeta(2), dydzeta(2);
396 : #endif //DEBUG
397 :
398 4036608 : hermite_compute_coefs(elem, dxdxi
399 : #ifdef DEBUG
400 : , dydxi, dzdeta, dxdzeta, dzdxi, dxdeta, dydzeta
401 : #endif
402 : );
403 :
404 4036608 : const ElemType type = elem->type();
405 :
406 : const Order totalorder =
407 4372992 : order + add_p_level*elem->p_level();
408 :
409 4036608 : switch (totalorder)
410 : {
411 : // 3rd-order tricubic Hermite functions
412 4036608 : case THIRD:
413 : {
414 4036608 : switch (type)
415 : {
416 336384 : case HEX8:
417 : case HEX20:
418 : case HEX27:
419 : {
420 336384 : libmesh_assert_less (i, 64);
421 :
422 336384 : std::vector<unsigned int> bases1D;
423 :
424 4036608 : Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i);
425 :
426 4036608 : return coef *
427 4036608 : FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
428 4036608 : FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
429 8073216 : FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
430 : }
431 0 : default:
432 0 : libmesh_error_msg("ERROR: Unsupported element type " << Utility::enum_to_string(type));
433 : }
434 : }
435 : // by default throw an error
436 0 : default:
437 0 : libmesh_error_msg("ERROR: Unsupported polynomial order " << totalorder);
438 : }
439 3363840 : }
440 :
441 :
442 : template <>
443 0 : Real FE<3,HERMITE>::shape(const ElemType,
444 : const Order,
445 : const unsigned int,
446 : const Point &)
447 : {
448 0 : libmesh_error_msg("Hermite elements require the real element \nto construct gradient-based degrees of freedom.");
449 : return 0.;
450 : }
451 :
452 : template <>
453 0 : Real FE<3,HERMITE>::shape(const FEType fet,
454 : const Elem * elem,
455 : const unsigned int i,
456 : const Point & p,
457 : const bool add_p_level)
458 : {
459 0 : return FE<3,HERMITE>::shape(elem, fet.order, i, p, add_p_level);
460 : }
461 :
462 :
463 :
464 : template <>
465 11962368 : Real FE<3,HERMITE>::shape_deriv(const Elem * elem,
466 : const Order order,
467 : const unsigned int i,
468 : const unsigned int j,
469 : const Point & p,
470 : const bool add_p_level)
471 : {
472 996864 : libmesh_assert(elem);
473 996864 : libmesh_assert (j == 0 || j == 1 || j == 2);
474 :
475 14952960 : std::vector<std::vector<Real>> dxdxi(3, std::vector<Real>(2, 0));
476 :
477 : #ifdef DEBUG
478 3987456 : std::vector<Real> dydxi(2), dzdeta(2), dxdzeta(2);
479 3987456 : std::vector<Real> dzdxi(2), dxdeta(2), dydzeta(2);
480 : #endif //DEBUG
481 :
482 11962368 : hermite_compute_coefs(elem, dxdxi
483 : #ifdef DEBUG
484 : , dydxi, dzdeta, dxdzeta, dzdxi, dxdeta, dydzeta
485 : #endif
486 : );
487 :
488 11962368 : const ElemType type = elem->type();
489 :
490 : const Order totalorder =
491 12959232 : order + add_p_level*elem->p_level();
492 :
493 11962368 : switch (totalorder)
494 : {
495 : // 3rd-order tricubic Hermite functions
496 11962368 : case THIRD:
497 : {
498 11962368 : switch (type)
499 : {
500 996864 : case HEX8:
501 : case HEX20:
502 : case HEX27:
503 : {
504 996864 : libmesh_assert_less (i, 64);
505 :
506 1993728 : std::vector<unsigned int> bases1D;
507 :
508 11962368 : Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i);
509 :
510 11962368 : switch (j) // Derivative type
511 : {
512 3987456 : case 0:
513 3987456 : return coef *
514 3987456 : FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) *
515 3987456 : FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
516 3987456 : FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
517 : break;
518 3987456 : case 1:
519 3987456 : return coef *
520 3987456 : FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
521 3987456 : FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1)) *
522 3987456 : FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
523 : break;
524 3987456 : case 2:
525 3987456 : return coef *
526 3987456 : FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
527 3987456 : FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
528 3987456 : FEHermite<1>::hermite_raw_shape_deriv(bases1D[2],p(2));
529 : break;
530 0 : default:
531 0 : libmesh_error_msg("Invalid shape function derivative j = " << j);
532 : }
533 :
534 : }
535 0 : default:
536 0 : libmesh_error_msg("ERROR: Unsupported element type " << Utility::enum_to_string(type));
537 : }
538 : }
539 : // by default throw an error
540 0 : default:
541 0 : libmesh_error_msg("ERROR: Unsupported polynomial order " << totalorder);
542 : }
543 9968640 : }
544 :
545 :
546 : template <>
547 0 : Real FE<3,HERMITE>::shape_deriv(const ElemType,
548 : const Order,
549 : const unsigned int,
550 : const unsigned int,
551 : const Point &)
552 : {
553 0 : libmesh_error_msg("Hermite elements require the real element \nto construct gradient-based degrees of freedom.");
554 : return 0.;
555 : }
556 :
557 :
558 : template <>
559 0 : Real FE<3,HERMITE>::shape_deriv(const FEType fet,
560 : const Elem * elem,
561 : const unsigned int i,
562 : const unsigned int j,
563 : const Point & p,
564 : const bool add_p_level)
565 : {
566 0 : return FE<3,HERMITE>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
567 : }
568 :
569 :
570 :
571 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
572 :
573 :
574 : template <>
575 23924736 : Real FE<3,HERMITE>::shape_second_deriv(const Elem * elem,
576 : const Order order,
577 : const unsigned int i,
578 : const unsigned int j,
579 : const Point & p,
580 : const bool add_p_level)
581 : {
582 1993728 : libmesh_assert(elem);
583 :
584 29905920 : std::vector<std::vector<Real>> dxdxi(3, std::vector<Real>(2, 0));
585 :
586 : #ifdef DEBUG
587 7974912 : std::vector<Real> dydxi(2), dzdeta(2), dxdzeta(2);
588 7974912 : std::vector<Real> dzdxi(2), dxdeta(2), dydzeta(2);
589 : #endif //DEBUG
590 :
591 23924736 : hermite_compute_coefs(elem, dxdxi
592 : #ifdef DEBUG
593 : , dydxi, dzdeta, dxdzeta, dzdxi, dxdeta, dydzeta
594 : #endif
595 : );
596 :
597 23924736 : const ElemType type = elem->type();
598 :
599 : const Order totalorder =
600 25918464 : order + add_p_level*elem->p_level();
601 :
602 23924736 : switch (totalorder)
603 : {
604 : // 3rd-order tricubic Hermite functions
605 23924736 : case THIRD:
606 : {
607 23924736 : switch (type)
608 : {
609 1993728 : case HEX8:
610 : case HEX20:
611 : case HEX27:
612 : {
613 1993728 : libmesh_assert_less (i, 64);
614 :
615 3987456 : std::vector<unsigned int> bases1D;
616 :
617 23924736 : Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i);
618 :
619 23924736 : switch (j) // Derivative type
620 : {
621 3987456 : case 0:
622 3987456 : return coef *
623 3987456 : FEHermite<1>::hermite_raw_shape_second_deriv(bases1D[0],p(0)) *
624 3987456 : FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
625 3987456 : FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
626 : break;
627 3987456 : case 1:
628 3987456 : return coef *
629 3987456 : FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) *
630 3987456 : FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1)) *
631 3987456 : FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
632 : break;
633 3987456 : case 2:
634 3987456 : return coef *
635 3987456 : FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
636 3987456 : FEHermite<1>::hermite_raw_shape_second_deriv(bases1D[1],p(1)) *
637 3987456 : FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
638 : break;
639 3987456 : case 3:
640 3987456 : return coef *
641 3987456 : FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) *
642 3987456 : FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
643 3987456 : FEHermite<1>::hermite_raw_shape_deriv(bases1D[2],p(2));
644 : break;
645 3987456 : case 4:
646 3987456 : return coef *
647 3987456 : FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
648 3987456 : FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1)) *
649 3987456 : FEHermite<1>::hermite_raw_shape_deriv(bases1D[2],p(2));
650 : break;
651 3987456 : case 5:
652 3987456 : return coef *
653 3987456 : FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
654 3987456 : FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
655 3987456 : FEHermite<1>::hermite_raw_shape_second_deriv(bases1D[2],p(2));
656 : break;
657 0 : default:
658 0 : libmesh_error_msg("Invalid shape function derivative j = " << j);
659 : }
660 :
661 : }
662 0 : default:
663 0 : libmesh_error_msg("ERROR: Unsupported element type " << Utility::enum_to_string(type));
664 : }
665 : }
666 : // by default throw an error
667 0 : default:
668 0 : libmesh_error_msg("ERROR: Unsupported polynomial order " << totalorder);
669 : }
670 19937280 : }
671 :
672 :
673 : template <>
674 0 : Real FE<3,HERMITE>::shape_second_deriv(const ElemType,
675 : const Order,
676 : const unsigned int,
677 : const unsigned int,
678 : const Point &)
679 : {
680 0 : libmesh_error_msg("Hermite elements require the real element \nto construct gradient-based degrees of freedom.");
681 : return 0.;
682 : }
683 :
684 : template <>
685 0 : Real FE<3,HERMITE>::shape_second_deriv(const FEType fet,
686 : const Elem * elem,
687 : const unsigned int i,
688 : const unsigned int j,
689 : const Point & p,
690 : const bool add_p_level)
691 : {
692 0 : return FE<3,HERMITE>::shape_second_deriv(elem, fet.order, i, j, p, add_p_level);
693 : }
694 :
695 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
696 :
697 : } // namespace libMesh
|
