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/enum_to_string.h"
23 :
24 : namespace libMesh
25 : {
26 :
27 : template <>
28 10946952 : RealGradient FE<3,NEDELEC_ONE>::shape(const Elem * elem,
29 : const Order order,
30 : const unsigned int i,
31 : const Point & p,
32 : const bool add_p_level)
33 : {
34 : #if LIBMESH_DIM == 3
35 927312 : libmesh_assert(elem);
36 :
37 11874264 : const Order totalorder = order + add_p_level*elem->p_level();
38 927312 : libmesh_assert_less(i, n_dofs(elem->type(), totalorder));
39 :
40 10946952 : const char sign = elem->positive_edge_orientation(i) ? 1 : -1;
41 :
42 10946952 : const Real xi = p(0);
43 10946952 : const Real eta = p(1);
44 10946952 : const Real zeta = p(2);
45 :
46 10946952 : switch (totalorder)
47 : {
48 : // linear Nedelec (first kind) shape functions
49 10946952 : case FIRST:
50 : {
51 10946952 : switch (elem->type())
52 : {
53 1615752 : case HEX20:
54 : case HEX27:
55 : {
56 : // Even with a loose inverse_map tolerance we ought to
57 : // be nearly on the element interior in master
58 : // coordinates
59 125520 : libmesh_assert_less_equal ( std::fabs(xi), 1.0+10*TOLERANCE );
60 125520 : libmesh_assert_less_equal ( std::fabs(eta), 1.0+10*TOLERANCE );
61 125520 : libmesh_assert_less_equal ( std::fabs(zeta), 1.0+10*TOLERANCE );
62 :
63 1615752 : switch(i)
64 : {
65 134646 : case 0:
66 134646 : return sign * RealGradient( -0.125*(1.0-eta-zeta+eta*zeta), 0.0, 0.0 );
67 134646 : case 1:
68 134646 : return sign * RealGradient( 0.0, -0.125*(1.0+xi-zeta-xi*zeta), 0.0 );
69 134646 : case 2:
70 134646 : return sign * RealGradient( 0.125*(1.0+eta-zeta-eta*zeta), 0.0, 0.0 );
71 134646 : case 3:
72 134646 : return sign * RealGradient( 0.0, -0.125*(1.0-xi-zeta+xi*zeta), 0.0 );
73 134646 : case 4:
74 134646 : return sign * RealGradient( 0.0, 0.0, -0.125*(1.0-xi-eta+xi*eta) );
75 134646 : case 5:
76 134646 : return sign * RealGradient( 0.0, 0.0, -0.125*(1.0+xi-eta-xi*eta) );
77 134646 : case 6:
78 134646 : return sign * RealGradient( 0.0, 0.0, -0.125*(1.0+xi+eta+xi*eta) );
79 134646 : case 7:
80 134646 : return sign * RealGradient( 0.0, 0.0, -0.125*(1.0-xi+eta-xi*eta) );
81 134646 : case 8:
82 134646 : return sign * RealGradient( -0.125*(1.0-eta+zeta-eta*zeta), 0.0, 0.0 );
83 134646 : case 9:
84 134646 : return sign * RealGradient( 0.0, -0.125*(1.0+xi+zeta+xi*zeta), 0.0 );
85 134646 : case 10:
86 134646 : return sign * RealGradient( 0.125*(1.0+eta+zeta+eta*zeta), 0.0, 0.0 );
87 134646 : case 11:
88 134646 : return sign * RealGradient( 0.0, -0.125*(1.0-xi+zeta-xi*zeta), 0.0 );
89 :
90 0 : default:
91 0 : libmesh_error_msg("Invalid i = " << i);
92 : }
93 : }
94 :
95 9331200 : case TET10:
96 : case TET14:
97 : {
98 9331200 : switch(i)
99 : {
100 1555200 : case 0:
101 1555200 : return sign * RealGradient( -1.0+eta+zeta, -xi, -xi );
102 1555200 : case 1:
103 1555200 : return sign * RealGradient( eta, -xi, 0.0 );
104 1555200 : case 2:
105 1555200 : return sign * RealGradient( -eta, -1.0+xi+zeta, -eta );
106 1555200 : case 3:
107 1555200 : return sign * RealGradient( -zeta, -zeta, -1.0+xi+eta );
108 1555200 : case 4:
109 1555200 : return sign * RealGradient( zeta, 0.0, -xi );
110 1555200 : case 5:
111 1555200 : return sign * RealGradient( 0.0, zeta, -eta );
112 :
113 0 : default:
114 0 : libmesh_error_msg("Invalid i = " << i);
115 : }
116 : }
117 :
118 0 : default:
119 0 : libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(elem->type()));
120 : }
121 : }
122 :
123 : // unsupported order
124 0 : default:
125 0 : libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << totalorder);
126 : }
127 :
128 : #else // LIBMESH_DIM != 3
129 : libmesh_ignore(elem, order, i, p, add_p_level);
130 : libmesh_not_implemented();
131 : #endif
132 : }
133 :
134 :
135 : template <>
136 0 : RealGradient FE<3,NEDELEC_ONE>::shape(const ElemType,
137 : const Order,
138 : const unsigned int,
139 : const Point &)
140 : {
141 0 : libmesh_error_msg("Nedelec elements require the element type \nbecause edge orientation is needed.");
142 : return RealGradient();
143 : }
144 :
145 :
146 : template <>
147 0 : RealGradient FE<3,NEDELEC_ONE>::shape(const FEType fet,
148 : const Elem * elem,
149 : const unsigned int i,
150 : const Point & p,
151 : const bool add_p_level)
152 : {
153 0 : return FE<3,NEDELEC_ONE>::shape(elem, fet.order, i, p, add_p_level);
154 : }
155 :
156 :
157 : template <>
158 23136408 : RealGradient FE<3,NEDELEC_ONE>::shape_deriv(const Elem * elem,
159 : const Order order,
160 : const unsigned int i,
161 : const unsigned int j,
162 : const Point & p,
163 : const bool add_p_level)
164 : {
165 : #if LIBMESH_DIM == 3
166 1954800 : libmesh_assert(elem);
167 1954800 : libmesh_assert_less (j, 3);
168 :
169 25091208 : const Order totalorder = order + add_p_level*elem->p_level();
170 1954800 : libmesh_assert_less(i, n_dofs(elem->type(), totalorder));
171 :
172 23136408 : const char sign = elem->positive_edge_orientation(i) ? 1 : -1;
173 :
174 23136408 : const Real xi = p(0);
175 23136408 : const Real eta = p(1);
176 23136408 : const Real zeta = p(2);
177 :
178 23136408 : switch (totalorder)
179 : {
180 : // linear Nedelec (first kind) shape function first derivatives
181 23136408 : case FIRST:
182 : {
183 23136408 : switch (elem->type())
184 : {
185 686232 : case HEX20:
186 : case HEX27:
187 : {
188 : // Even with a loose inverse_map tolerance we ought to
189 : // be nearly on the element interior in master
190 : // coordinates
191 19440 : libmesh_assert_less_equal ( std::fabs(xi), 1.0+10*TOLERANCE );
192 19440 : libmesh_assert_less_equal ( std::fabs(eta), 1.0+10*TOLERANCE );
193 19440 : libmesh_assert_less_equal ( std::fabs(zeta), 1.0+10*TOLERANCE );
194 :
195 686232 : switch (j)
196 : {
197 : // d()/dxi
198 228744 : case 0:
199 : {
200 228744 : switch(i)
201 : {
202 2160 : case 0:
203 : case 2:
204 : case 8:
205 : case 10:
206 2160 : return RealGradient();
207 19062 : case 1:
208 19062 : return sign * RealGradient( 0.0, -0.125*(1.0-zeta) );
209 19062 : case 3:
210 19062 : return sign * RealGradient( 0.0, -0.125*(-1.0+zeta) );
211 19062 : case 4:
212 19062 : return sign * RealGradient( 0.0, 0.0, -0.125*(-1.0+eta) );
213 19062 : case 5:
214 19062 : return sign * RealGradient( 0.0, 0.0, -0.125*(1.0-eta) );
215 19062 : case 6:
216 19062 : return sign * RealGradient( 0.0, 0.0, -0.125*(1.0+eta) );
217 19062 : case 7:
218 19062 : return sign * RealGradient( 0.0, 0.0, -0.125*(-1.0-eta) );
219 19062 : case 9:
220 19062 : return sign * RealGradient( 0.0, -0.125*(1.0+zeta), 0.0 );
221 19062 : case 11:
222 19062 : return sign * RealGradient( 0.0, -0.125*(-1.0-zeta), 0.0 );
223 :
224 0 : default:
225 0 : libmesh_error_msg("Invalid i = " << i);
226 : } // switch(i)
227 :
228 : } // j = 0
229 :
230 : // d()/deta
231 228744 : case 1:
232 : {
233 228744 : switch(i)
234 : {
235 2160 : case 1:
236 : case 3:
237 : case 9:
238 : case 11:
239 2160 : return RealGradient();
240 19062 : case 0:
241 19062 : return sign * RealGradient( -0.125*(-1.0+zeta), 0.0, 0.0 );
242 19062 : case 2:
243 19062 : return sign * RealGradient( 0.125*(1.0-zeta), 0.0, 0.0 );
244 19062 : case 4:
245 19062 : return sign * RealGradient( 0.0, 0.0, -0.125*(-1.0+xi) );
246 19062 : case 5:
247 19062 : return sign * RealGradient( 0.0, 0.0, -0.125*(-1.0-xi) );
248 19062 : case 6:
249 19062 : return sign * RealGradient( 0.0, 0.0, -0.125*(1.0+xi) );
250 19062 : case 7:
251 19062 : return sign * RealGradient( 0.0, 0.0, -0.125*(1.0-xi) );
252 19062 : case 8:
253 19062 : return sign * RealGradient( -0.125*(-1.0-zeta), 0.0, 0.0 );
254 19062 : case 10:
255 19062 : return sign * RealGradient( 0.125*(1.0+zeta), 0.0, 0.0 );
256 :
257 0 : default:
258 0 : libmesh_error_msg("Invalid i = " << i);
259 : } // switch(i)
260 :
261 : } // j = 1
262 :
263 : // d()/dzeta
264 228744 : case 2:
265 : {
266 228744 : switch(i)
267 : {
268 2160 : case 4:
269 : case 5:
270 : case 6:
271 : case 7:
272 2160 : return RealGradient();
273 19062 : case 0:
274 19062 : return sign * RealGradient( -0.125*(-1.0+eta), 0.0, 0.0 );
275 19062 : case 1:
276 19062 : return sign * RealGradient( 0.0, -0.125*(-1.0-xi), 0.0 );
277 19062 : case 2:
278 19062 : return sign * RealGradient( 0.125*(-1.0-eta), 0.0, 0.0 );
279 19062 : case 3:
280 19062 : return sign * RealGradient( 0.0, -0.125*(-1.0+xi), 0.0 );
281 19062 : case 8:
282 19062 : return sign * RealGradient( -0.125*(1.0-eta), 0.0, 0.0 );
283 19062 : case 9:
284 19062 : return sign * RealGradient( 0.0, -0.125*(1.0+xi), 0.0 );
285 19062 : case 10:
286 19062 : return sign * RealGradient( 0.125*(1.0+eta), 0.0, 0.0 );
287 19062 : case 11:
288 19062 : return sign * RealGradient( 0.0, -0.125*(1.0-xi), 0.0 );
289 :
290 0 : default:
291 0 : libmesh_error_msg("Invalid i = " << i);
292 : } // switch(i)
293 :
294 : } // j = 2
295 :
296 0 : default:
297 0 : libmesh_error_msg("Invalid j = " << j);
298 : }
299 : }
300 :
301 22450176 : case TET10:
302 : case TET14:
303 : {
304 22450176 : switch (j)
305 : {
306 : // d()/dxi
307 7483392 : case 0:
308 : {
309 7483392 : switch(i)
310 : {
311 1247232 : case 0:
312 107520 : return sign * RealGradient( 0.0, -1.0, -1.0 );
313 1247232 : case 1:
314 107520 : return sign * RealGradient( 0.0, -1.0, 0.0 );
315 1247232 : case 2:
316 107520 : return sign * RealGradient( 0.0, 1.0, 0.0 );
317 1247232 : case 3:
318 107520 : return sign * RealGradient( 0.0, 0.0, 1.0 );
319 1247232 : case 4:
320 107520 : return sign * RealGradient( 0.0, 0.0, -1.0 );
321 107520 : case 5:
322 107520 : return RealGradient();
323 :
324 0 : default:
325 0 : libmesh_error_msg("Invalid i = " << i);
326 : } // switch(i)
327 :
328 : } // j = 0
329 :
330 : // d()/deta
331 7483392 : case 1:
332 : {
333 7483392 : switch(i)
334 : {
335 2494464 : case 0:
336 : case 1:
337 215040 : return sign * RealGradient( 1.0, 0.0, 0.0 );
338 1247232 : case 2:
339 107520 : return sign * RealGradient( -1.0, 0.0, -1.0 );
340 1247232 : case 3:
341 107520 : return sign * RealGradient( 0.0, 0.0, 1.0 );
342 107520 : case 4:
343 107520 : return RealGradient();
344 1247232 : case 5:
345 107520 : return sign * RealGradient( 0.0, 0.0, -1.0 );
346 :
347 0 : default:
348 0 : libmesh_error_msg("Invalid i = " << i);
349 : } // switch(i)
350 :
351 : } // j = 1
352 :
353 : // d()/dzeta
354 7483392 : case 2:
355 : {
356 7483392 : switch(i)
357 : {
358 1247232 : case 0:
359 107520 : return sign * RealGradient( 1.0, 0.0, 0.0 );
360 107520 : case 1:
361 107520 : return RealGradient();
362 2494464 : case 2:
363 : case 5:
364 215040 : return sign * RealGradient( 0.0, 1.0, 0.0 );
365 1247232 : case 3:
366 107520 : return sign * RealGradient( -1.0, -1.0, 0.0 );
367 1247232 : case 4:
368 107520 : return sign * RealGradient( 1.0, 0.0, 0.0 );
369 :
370 0 : default:
371 0 : libmesh_error_msg("Invalid i = " << i);
372 : } // switch(i)
373 :
374 : } // j = 2
375 :
376 0 : default:
377 0 : libmesh_error_msg("Invalid j = " << j);
378 : }
379 : }
380 :
381 0 : default:
382 0 : libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(elem->type()));
383 : }
384 : }
385 : // unsupported order
386 0 : default:
387 0 : libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << totalorder);
388 : }
389 :
390 : #else // LIBMESH_DIM != 3
391 : libmesh_ignore(elem, order, i, j, p, add_p_level);
392 : libmesh_not_implemented();
393 : #endif
394 : }
395 :
396 :
397 : template <>
398 0 : RealGradient FE<3,NEDELEC_ONE>::shape_deriv(const ElemType,
399 : const Order,
400 : const unsigned int,
401 : const unsigned int,
402 : const Point &)
403 : {
404 0 : libmesh_error_msg("Nedelec elements require the element type \nbecause edge orientation is needed.");
405 : return RealGradient();
406 : }
407 :
408 :
409 : template <>
410 0 : RealGradient FE<3,NEDELEC_ONE>::shape_deriv(const FEType fet,
411 : const Elem * elem,
412 : const unsigned int i,
413 : const unsigned int j,
414 : const Point & p,
415 : const bool add_p_level)
416 : {
417 0 : return FE<3,NEDELEC_ONE>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
418 : }
419 :
420 :
421 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
422 :
423 : template <>
424 0 : RealGradient FE<3,NEDELEC_ONE>::shape_second_deriv(const Elem * elem,
425 : const Order order,
426 : const unsigned int i,
427 : const unsigned int j,
428 : const Point &,
429 : const bool add_p_level)
430 : {
431 : #if LIBMESH_DIM == 3
432 :
433 0 : libmesh_assert(elem);
434 :
435 : // j = 0 ==> d^2 phi / dxi^2
436 : // j = 1 ==> d^2 phi / dxi deta
437 : // j = 2 ==> d^2 phi / deta^2
438 : // j = 3 ==> d^2 phi / dxi dzeta
439 : // j = 4 ==> d^2 phi / deta dzeta
440 : // j = 5 ==> d^2 phi / dzeta^2
441 0 : libmesh_assert_less (j, 6);
442 :
443 0 : const Order totalorder = order + add_p_level*elem->p_level();
444 0 : libmesh_assert_less(i, n_dofs(elem->type(), totalorder));
445 :
446 0 : const char sign = elem->positive_edge_orientation(i) ? 1 : -1;
447 :
448 0 : switch (totalorder)
449 : {
450 : // linear Nedelec (first kind) shape function second derivatives
451 0 : case FIRST:
452 : {
453 0 : switch (elem->type())
454 : {
455 0 : case HEX20:
456 : case HEX27:
457 : {
458 0 : switch (j)
459 : {
460 : // d^2()/dxi^2, d^2()/deta^2, d^2()/dzeta^2
461 0 : case 0:
462 : case 2:
463 : case 5:
464 0 : return RealGradient();
465 :
466 : // d^2()/dxideta
467 0 : case 1:
468 : {
469 0 : switch(i)
470 : {
471 0 : case 0:
472 : case 1:
473 : case 2:
474 : case 3:
475 : case 8:
476 : case 9:
477 : case 10:
478 : case 11:
479 0 : return RealGradient();
480 0 : case 4:
481 : case 6:
482 0 : return sign * RealGradient( 0.0, 0.0, -0.125 );
483 0 : case 5:
484 : case 7:
485 0 : return sign * RealGradient( 0.0, 0.0, 0.125 );
486 :
487 0 : default:
488 0 : libmesh_error_msg("Invalid i = " << i);
489 : } // switch(i)
490 :
491 : } // j = 1
492 :
493 : // d^2()/dxidzeta
494 0 : case 3:
495 : {
496 0 : switch(i)
497 : {
498 0 : case 0:
499 : case 2:
500 : case 4:
501 : case 5:
502 : case 6:
503 : case 7:
504 : case 8:
505 : case 10:
506 0 : return RealGradient();
507 0 : case 1:
508 : case 3:
509 : case 11:
510 0 : return sign * RealGradient( 0.0, 0.125 );
511 0 : case 9:
512 0 : return sign * RealGradient( 0.0, -0.125, 0.0 );
513 :
514 0 : default:
515 0 : libmesh_error_msg("Invalid i = " << i);
516 : } // switch(i)
517 :
518 : } // j = 3
519 :
520 : // d^2()/detadzeta
521 0 : case 4:
522 : {
523 0 : switch(i)
524 : {
525 0 : case 0:
526 0 : return sign * RealGradient( -0.125, 0.0, 0.0 );
527 0 : case 1:
528 : case 3:
529 : case 4:
530 : case 5:
531 : case 6:
532 : case 7:
533 : case 9:
534 : case 11:
535 0 : return RealGradient();
536 0 : case 2:
537 : case 8:
538 : case 10:
539 0 : return sign * RealGradient( 0.125, 0.0, 0.0 );
540 :
541 0 : default:
542 0 : libmesh_error_msg("Invalid i = " << i);
543 : } // switch(i)
544 :
545 : } // j = 4
546 :
547 0 : default:
548 0 : libmesh_error_msg("Invalid j = " << j);
549 : }
550 : }
551 :
552 : // All second derivatives for linear tets are zero.
553 0 : case TET10:
554 : case TET14:
555 0 : return RealGradient();
556 :
557 0 : default:
558 0 : libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(elem->type()));
559 :
560 : } //switch(type)
561 :
562 : }
563 :
564 : // unsupported order
565 0 : default:
566 0 : libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << totalorder);
567 : }
568 :
569 : #else // LIBMESH_DIM != 3
570 : libmesh_assert(true || p(0));
571 : libmesh_ignore(elem, order, i, j, add_p_level);
572 : libmesh_not_implemented();
573 : #endif
574 : }
575 :
576 :
577 : template <>
578 0 : RealGradient FE<3,NEDELEC_ONE>::shape_second_deriv(const ElemType,
579 : const Order,
580 : const unsigned int,
581 : const unsigned int,
582 : const Point &)
583 : {
584 0 : libmesh_error_msg("Nedelec elements require the element type \nbecause edge orientation is needed.");
585 : return RealGradient();
586 : }
587 :
588 :
589 : template <>
590 0 : RealGradient FE<3,NEDELEC_ONE>::shape_second_deriv(const FEType fet,
591 : const Elem * elem,
592 : const unsigned int i,
593 : const unsigned int j,
594 : const Point & p,
595 : const bool add_p_level)
596 : {
597 0 : return FE<3,NEDELEC_ONE>::shape_second_deriv(elem, fet.order, i, j, p, add_p_level);
598 : }
599 :
600 : #endif
601 :
602 : } // namespace libMesh
|
