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/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 42872832 : 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 42872832 : const FEFamily mapping_family = FEMap::map_fe_type(*elem);
40 42872832 : 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 42872832 : FEInterface::n_shape_functions (map_fe_type, /*extra_order=*/0, elem);
46 :
47 42872832 : static const Point dofpt[2] = {Point(-1,-1,-1), Point(1,1,1)};
48 :
49 : FEInterface::shape_deriv_ptr shape_deriv_ptr =
50 42872832 : FEInterface::shape_deriv_function(map_fe_type, elem);
51 :
52 128618496 : for (int p = 0; p != 2; ++p)
53 : {
54 85745664 : dxdxi[0][p] = 0;
55 85745664 : dxdxi[1][p] = 0;
56 92891136 : dxdxi[2][p] = 0;
57 : #ifdef DEBUG
58 7145472 : dydxi[p] = 0;
59 7145472 : dzdeta[p] = 0;
60 7145472 : dxdzeta[p] = 0;
61 7145472 : dzdxi[p] = 0;
62 7145472 : dxdeta[p] = 0;
63 7145472 : dydzeta[p] = 0;
64 : #endif
65 2400878592 : for (int i = 0; i != n_mapping_shape_functions; ++i)
66 : {
67 : const Real ddxi = shape_deriv_ptr
68 2315132928 : (map_fe_type, elem, i, 0, dofpt[p], /*add_p_level=*/false);
69 : const Real ddeta = shape_deriv_ptr
70 2315132928 : (map_fe_type, elem, i, 1, dofpt[p], /*add_p_level=*/false);
71 : const Real ddzeta = shape_deriv_ptr
72 2315132928 : (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 385855488 : const Point & point_i = elem->point(i);
77 2315132928 : dxdxi[0][p] += point_i(0) * ddxi;
78 2315132928 : dxdxi[1][p] += point_i(1) * ddeta;
79 2508060672 : dxdxi[2][p] += point_i(2) * ddzeta;
80 : #ifdef DEBUG
81 192927744 : dydxi[p] += point_i(1) * ddxi;
82 192927744 : dzdeta[p] += point_i(2) * ddeta;
83 192927744 : dxdzeta[p] += point_i(0) * ddzeta;
84 192927744 : dzdxi[p] += point_i(2) * ddxi;
85 192927744 : dxdeta[p] += point_i(0) * ddeta;
86 192927744 : dydzeta[p] += point_i(1) * ddzeta;
87 : #endif
88 : }
89 :
90 : // No singular elements!
91 7145472 : libmesh_assert(dxdxi[0][p]);
92 7145472 : libmesh_assert(dxdxi[1][p]);
93 7145472 : libmesh_assert(dxdxi[2][p]);
94 : // No non-rectilinear or non-axis-aligned elements!
95 : #ifdef DEBUG
96 7145472 : libmesh_assert_less (std::abs(dydxi[p]), TOLERANCE);
97 7145472 : libmesh_assert_less (std::abs(dzdeta[p]), TOLERANCE);
98 7145472 : libmesh_assert_less (std::abs(dxdzeta[p]), TOLERANCE);
99 7145472 : libmesh_assert_less (std::abs(dzdxi[p]), TOLERANCE);
100 7145472 : libmesh_assert_less (std::abs(dxdeta[p]), TOLERANCE);
101 7145472 : libmesh_assert_less (std::abs(dydzeta[p]), TOLERANCE);
102 : #endif
103 : }
104 42872832 : }
105 :
106 :
107 :
108 42872832 : 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 3572736 : bases1D.clear();
114 42872832 : bases1D.resize(3,0);
115 3572736 : Real coef = 1.0;
116 :
117 42872832 : unsigned int e = o-2;
118 :
119 : // Nodes
120 42872832 : if (i < 64)
121 : {
122 42872832 : switch (i / 8)
123 : {
124 446592 : case 0:
125 446592 : break;
126 893184 : case 1:
127 5359104 : bases1D[0] = 1;
128 5359104 : break;
129 893184 : case 2:
130 5359104 : bases1D[0] = 1;
131 5359104 : bases1D[1] = 1;
132 5359104 : break;
133 893184 : case 3:
134 5359104 : bases1D[1] = 1;
135 5359104 : break;
136 893184 : case 4:
137 5359104 : bases1D[2] = 1;
138 5359104 : break;
139 893184 : case 5:
140 5359104 : bases1D[0] = 1;
141 5359104 : bases1D[2] = 1;
142 5359104 : break;
143 893184 : case 6:
144 5359104 : bases1D[0] = 1;
145 5359104 : bases1D[1] = 1;
146 5359104 : bases1D[2] = 1;
147 5359104 : break;
148 893184 : case 7:
149 5359104 : bases1D[1] = 1;
150 5359104 : bases1D[2] = 1;
151 5359104 : break;
152 0 : default:
153 0 : libmesh_error_msg("Invalid basis node = " << i/8);
154 : }
155 :
156 42872832 : unsigned int basisnum = i%8;
157 42872832 : switch (basisnum) // DoF type
158 : {
159 446592 : case 0: // DoF = value at node
160 446592 : coef = 1.0;
161 446592 : break;
162 893184 : case 1: // DoF = x derivative at node
163 5359104 : coef = dxdxi[0][bases1D[0]];
164 5359104 : bases1D[0] += 2; break;
165 893184 : case 2: // DoF = y derivative at node
166 5359104 : coef = dxdxi[1][bases1D[1]];
167 5359104 : bases1D[1] += 2; break;
168 893184 : case 3: // DoF = xy derivative at node
169 5805696 : coef = dxdxi[0][bases1D[0]] * dxdxi[1][bases1D[1]];
170 5359104 : bases1D[0] += 2; bases1D[1] += 2; break;
171 893184 : case 4: // DoF = z derivative at node
172 5359104 : coef = dxdxi[2][bases1D[2]];
173 5359104 : bases1D[2] += 2; break;
174 893184 : case 5: // DoF = xz derivative at node
175 5805696 : coef = dxdxi[0][bases1D[0]] * dxdxi[2][bases1D[2]];
176 5359104 : bases1D[0] += 2; bases1D[2] += 2; break;
177 893184 : case 6: // DoF = yz derivative at node
178 5805696 : coef = dxdxi[1][bases1D[1]] * dxdxi[2][bases1D[2]];
179 5359104 : bases1D[1] += 2; bases1D[2] += 2; break;
180 893184 : case 7: // DoF = xyz derivative at node
181 6252288 : coef = dxdxi[0][bases1D[0]] * dxdxi[1][bases1D[1]] * dxdxi[2][bases1D[2]];
182 5359104 : 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 3572736 : libmesh_assert(coef);
368 42872832 : return coef;
369 : }
370 :
371 :
372 : } // end anonymous namespace
373 :
374 :
375 : namespace libMesh
376 : {
377 :
378 :
379 16678932 : LIBMESH_DEFAULT_VECTORIZED_FE(3,HERMITE)
380 :
381 :
382 : template <>
383 4331520 : 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 360960 : libmesh_assert(elem);
390 :
391 5414400 : std::vector<std::vector<Real>> dxdxi(3, std::vector<Real>(2, 0));
392 :
393 : #ifdef DEBUG
394 1443840 : std::vector<Real> dydxi(2), dzdeta(2), dxdzeta(2);
395 1443840 : std::vector<Real> dzdxi(2), dxdeta(2), dydzeta(2);
396 : #endif //DEBUG
397 :
398 4331520 : hermite_compute_coefs(elem, dxdxi
399 : #ifdef DEBUG
400 : , dydxi, dzdeta, dxdzeta, dzdxi, dxdeta, dydzeta
401 : #endif
402 : );
403 :
404 4331520 : const ElemType type = elem->type();
405 :
406 : const Order totalorder =
407 4692480 : order + add_p_level*elem->p_level();
408 :
409 4331520 : switch (totalorder)
410 : {
411 : // 3rd-order tricubic Hermite functions
412 4331520 : case THIRD:
413 : {
414 4331520 : switch (type)
415 : {
416 360960 : case HEX8:
417 : case HEX20:
418 : case HEX27:
419 : {
420 360960 : libmesh_assert_less (i, 64);
421 :
422 360960 : std::vector<unsigned int> bases1D;
423 :
424 4331520 : Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i);
425 :
426 4331520 : return coef *
427 4331520 : FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
428 4331520 : FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
429 8663040 : 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 3609600 : }
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 12847104 : 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 1070592 : libmesh_assert(elem);
473 1070592 : libmesh_assert (j == 0 || j == 1 || j == 2);
474 :
475 16058880 : std::vector<std::vector<Real>> dxdxi(3, std::vector<Real>(2, 0));
476 :
477 : #ifdef DEBUG
478 4282368 : std::vector<Real> dydxi(2), dzdeta(2), dxdzeta(2);
479 4282368 : std::vector<Real> dzdxi(2), dxdeta(2), dydzeta(2);
480 : #endif //DEBUG
481 :
482 12847104 : hermite_compute_coefs(elem, dxdxi
483 : #ifdef DEBUG
484 : , dydxi, dzdeta, dxdzeta, dzdxi, dxdeta, dydzeta
485 : #endif
486 : );
487 :
488 12847104 : const ElemType type = elem->type();
489 :
490 : const Order totalorder =
491 13917696 : order + add_p_level*elem->p_level();
492 :
493 12847104 : switch (totalorder)
494 : {
495 : // 3rd-order tricubic Hermite functions
496 12847104 : case THIRD:
497 : {
498 12847104 : switch (type)
499 : {
500 1070592 : case HEX8:
501 : case HEX20:
502 : case HEX27:
503 : {
504 1070592 : libmesh_assert_less (i, 64);
505 :
506 2141184 : std::vector<unsigned int> bases1D;
507 :
508 12847104 : Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i);
509 :
510 12847104 : switch (j) // Derivative type
511 : {
512 4282368 : case 0:
513 4282368 : return coef *
514 4282368 : FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) *
515 4282368 : FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
516 4282368 : FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
517 : break;
518 4282368 : case 1:
519 4282368 : return coef *
520 4282368 : FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
521 4282368 : FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1)) *
522 4282368 : FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
523 : break;
524 4282368 : case 2:
525 4282368 : return coef *
526 4282368 : FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
527 4282368 : FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
528 4282368 : 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 10705920 : }
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 25694208 : 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 2141184 : libmesh_assert(elem);
583 :
584 32117760 : std::vector<std::vector<Real>> dxdxi(3, std::vector<Real>(2, 0));
585 :
586 : #ifdef DEBUG
587 8564736 : std::vector<Real> dydxi(2), dzdeta(2), dxdzeta(2);
588 8564736 : std::vector<Real> dzdxi(2), dxdeta(2), dydzeta(2);
589 : #endif //DEBUG
590 :
591 25694208 : hermite_compute_coefs(elem, dxdxi
592 : #ifdef DEBUG
593 : , dydxi, dzdeta, dxdzeta, dzdxi, dxdeta, dydzeta
594 : #endif
595 : );
596 :
597 25694208 : const ElemType type = elem->type();
598 :
599 : const Order totalorder =
600 27835392 : order + add_p_level*elem->p_level();
601 :
602 25694208 : switch (totalorder)
603 : {
604 : // 3rd-order tricubic Hermite functions
605 25694208 : case THIRD:
606 : {
607 25694208 : switch (type)
608 : {
609 2141184 : case HEX8:
610 : case HEX20:
611 : case HEX27:
612 : {
613 2141184 : libmesh_assert_less (i, 64);
614 :
615 4282368 : std::vector<unsigned int> bases1D;
616 :
617 25694208 : Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i);
618 :
619 25694208 : switch (j) // Derivative type
620 : {
621 4282368 : case 0:
622 4282368 : return coef *
623 4282368 : FEHermite<1>::hermite_raw_shape_second_deriv(bases1D[0],p(0)) *
624 4282368 : FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
625 4282368 : FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
626 : break;
627 4282368 : case 1:
628 4282368 : return coef *
629 4282368 : FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) *
630 4282368 : FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1)) *
631 4282368 : FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
632 : break;
633 4282368 : case 2:
634 4282368 : return coef *
635 4282368 : FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
636 4282368 : FEHermite<1>::hermite_raw_shape_second_deriv(bases1D[1],p(1)) *
637 4282368 : FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
638 : break;
639 4282368 : case 3:
640 4282368 : return coef *
641 4282368 : FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) *
642 4282368 : FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
643 4282368 : FEHermite<1>::hermite_raw_shape_deriv(bases1D[2],p(2));
644 : break;
645 4282368 : case 4:
646 4282368 : return coef *
647 4282368 : FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
648 4282368 : FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1)) *
649 4282368 : FEHermite<1>::hermite_raw_shape_deriv(bases1D[2],p(2));
650 : break;
651 4282368 : case 5:
652 4282368 : return coef *
653 4282368 : FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
654 4282368 : FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
655 4282368 : 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 21411840 : }
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
|
