libMesh
Loading...
Searching...
No Matches
fe_nedelec_one_shape_3D.C
Go to the documentation of this file.
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/enum_to_string.h"
23
24namespace libMesh
25{
26
27template <>
28RealGradient 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 libmesh_assert(elem);
36
37 const Order totalorder = order + add_p_level*elem->p_level();
38 libmesh_assert_less(i, n_dofs(elem->type(), totalorder));
39
40 const char sign = elem->positive_edge_orientation(i) ? 1 : -1;
41
42 const Real xi = p(0);
43 const Real eta = p(1);
44 const Real zeta = p(2);
45
46 switch (totalorder)
47 {
48 // linear Nedelec (first kind) shape functions
49 case FIRST:
50 {
51 switch (elem->type())
52 {
53 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 libmesh_assert_less_equal ( std::fabs(xi), 1.0+10*TOLERANCE );
60 libmesh_assert_less_equal ( std::fabs(eta), 1.0+10*TOLERANCE );
61 libmesh_assert_less_equal ( std::fabs(zeta), 1.0+10*TOLERANCE );
62
63 switch(i)
64 {
65 case 0:
66 return sign * RealGradient( -0.125*(1.0-eta-zeta+eta*zeta), 0.0, 0.0 );
67 case 1:
68 return sign * RealGradient( 0.0, -0.125*(1.0+xi-zeta-xi*zeta), 0.0 );
69 case 2:
70 return sign * RealGradient( 0.125*(1.0+eta-zeta-eta*zeta), 0.0, 0.0 );
71 case 3:
72 return sign * RealGradient( 0.0, -0.125*(1.0-xi-zeta+xi*zeta), 0.0 );
73 case 4:
74 return sign * RealGradient( 0.0, 0.0, -0.125*(1.0-xi-eta+xi*eta) );
75 case 5:
76 return sign * RealGradient( 0.0, 0.0, -0.125*(1.0+xi-eta-xi*eta) );
77 case 6:
78 return sign * RealGradient( 0.0, 0.0, -0.125*(1.0+xi+eta+xi*eta) );
79 case 7:
80 return sign * RealGradient( 0.0, 0.0, -0.125*(1.0-xi+eta-xi*eta) );
81 case 8:
82 return sign * RealGradient( -0.125*(1.0-eta+zeta-eta*zeta), 0.0, 0.0 );
83 case 9:
84 return sign * RealGradient( 0.0, -0.125*(1.0+xi+zeta+xi*zeta), 0.0 );
85 case 10:
86 return sign * RealGradient( 0.125*(1.0+eta+zeta+eta*zeta), 0.0, 0.0 );
87 case 11:
88 return sign * RealGradient( 0.0, -0.125*(1.0-xi+zeta-xi*zeta), 0.0 );
89
90 default:
91 libmesh_error_msg("Invalid i = " << i);
92 }
93 }
94
95 case TET10:
96 case TET14:
97 {
98 switch(i)
99 {
100 case 0:
101 return sign * RealGradient( -1.0+eta+zeta, -xi, -xi );
102 case 1:
103 return sign * RealGradient( eta, -xi, 0.0 );
104 case 2:
105 return sign * RealGradient( -eta, -1.0+xi+zeta, -eta );
106 case 3:
107 return sign * RealGradient( -zeta, -zeta, -1.0+xi+eta );
108 case 4:
109 return sign * RealGradient( zeta, 0.0, -xi );
110 case 5:
111 return sign * RealGradient( 0.0, zeta, -eta );
112
113 default:
114 libmesh_error_msg("Invalid i = " << i);
115 }
116 }
117
118 default:
119 libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(elem->type()));
120 }
121 }
122
123 // unsupported order
124 default:
125 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
135template <>
136RealGradient FE<3,NEDELEC_ONE>::shape(const ElemType,
137 const Order,
138 const unsigned int,
139 const Point &)
140{
141 libmesh_error_msg("Nedelec elements require the element type \nbecause edge orientation is needed.");
142 return RealGradient();
143}
144
145
146template <>
147RealGradient 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 return FE<3,NEDELEC_ONE>::shape(elem, fet.order, i, p, add_p_level);
154}
155
156
157template <>
158RealGradient 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 libmesh_assert(elem);
167 libmesh_assert_less (j, 3);
168
169 const Order totalorder = order + add_p_level*elem->p_level();
170 libmesh_assert_less(i, n_dofs(elem->type(), totalorder));
171
172 const char sign = elem->positive_edge_orientation(i) ? 1 : -1;
173
174 const Real xi = p(0);
175 const Real eta = p(1);
176 const Real zeta = p(2);
177
178 switch (totalorder)
179 {
180 // linear Nedelec (first kind) shape function first derivatives
181 case FIRST:
182 {
183 switch (elem->type())
184 {
185 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 libmesh_assert_less_equal ( std::fabs(xi), 1.0+10*TOLERANCE );
192 libmesh_assert_less_equal ( std::fabs(eta), 1.0+10*TOLERANCE );
193 libmesh_assert_less_equal ( std::fabs(zeta), 1.0+10*TOLERANCE );
194
195 switch (j)
196 {
197 // d()/dxi
198 case 0:
199 {
200 switch(i)
201 {
202 case 0:
203 case 2:
204 case 8:
205 case 10:
206 return RealGradient();
207 case 1:
208 return sign * RealGradient( 0.0, -0.125*(1.0-zeta) );
209 case 3:
210 return sign * RealGradient( 0.0, -0.125*(-1.0+zeta) );
211 case 4:
212 return sign * RealGradient( 0.0, 0.0, -0.125*(-1.0+eta) );
213 case 5:
214 return sign * RealGradient( 0.0, 0.0, -0.125*(1.0-eta) );
215 case 6:
216 return sign * RealGradient( 0.0, 0.0, -0.125*(1.0+eta) );
217 case 7:
218 return sign * RealGradient( 0.0, 0.0, -0.125*(-1.0-eta) );
219 case 9:
220 return sign * RealGradient( 0.0, -0.125*(1.0+zeta), 0.0 );
221 case 11:
222 return sign * RealGradient( 0.0, -0.125*(-1.0-zeta), 0.0 );
223
224 default:
225 libmesh_error_msg("Invalid i = " << i);
226 } // switch(i)
227
228 } // j = 0
229
230 // d()/deta
231 case 1:
232 {
233 switch(i)
234 {
235 case 1:
236 case 3:
237 case 9:
238 case 11:
239 return RealGradient();
240 case 0:
241 return sign * RealGradient( -0.125*(-1.0+zeta), 0.0, 0.0 );
242 case 2:
243 return sign * RealGradient( 0.125*(1.0-zeta), 0.0, 0.0 );
244 case 4:
245 return sign * RealGradient( 0.0, 0.0, -0.125*(-1.0+xi) );
246 case 5:
247 return sign * RealGradient( 0.0, 0.0, -0.125*(-1.0-xi) );
248 case 6:
249 return sign * RealGradient( 0.0, 0.0, -0.125*(1.0+xi) );
250 case 7:
251 return sign * RealGradient( 0.0, 0.0, -0.125*(1.0-xi) );
252 case 8:
253 return sign * RealGradient( -0.125*(-1.0-zeta), 0.0, 0.0 );
254 case 10:
255 return sign * RealGradient( 0.125*(1.0+zeta), 0.0, 0.0 );
256
257 default:
258 libmesh_error_msg("Invalid i = " << i);
259 } // switch(i)
260
261 } // j = 1
262
263 // d()/dzeta
264 case 2:
265 {
266 switch(i)
267 {
268 case 4:
269 case 5:
270 case 6:
271 case 7:
272 return RealGradient();
273 case 0:
274 return sign * RealGradient( -0.125*(-1.0+eta), 0.0, 0.0 );
275 case 1:
276 return sign * RealGradient( 0.0, -0.125*(-1.0-xi), 0.0 );
277 case 2:
278 return sign * RealGradient( 0.125*(-1.0-eta), 0.0, 0.0 );
279 case 3:
280 return sign * RealGradient( 0.0, -0.125*(-1.0+xi), 0.0 );
281 case 8:
282 return sign * RealGradient( -0.125*(1.0-eta), 0.0, 0.0 );
283 case 9:
284 return sign * RealGradient( 0.0, -0.125*(1.0+xi), 0.0 );
285 case 10:
286 return sign * RealGradient( 0.125*(1.0+eta), 0.0, 0.0 );
287 case 11:
288 return sign * RealGradient( 0.0, -0.125*(1.0-xi), 0.0 );
289
290 default:
291 libmesh_error_msg("Invalid i = " << i);
292 } // switch(i)
293
294 } // j = 2
295
296 default:
297 libmesh_error_msg("Invalid j = " << j);
298 }
299 }
300
301 case TET10:
302 case TET14:
303 {
304 switch (j)
305 {
306 // d()/dxi
307 case 0:
308 {
309 switch(i)
310 {
311 case 0:
312 return sign * RealGradient( 0.0, -1.0, -1.0 );
313 case 1:
314 return sign * RealGradient( 0.0, -1.0, 0.0 );
315 case 2:
316 return sign * RealGradient( 0.0, 1.0, 0.0 );
317 case 3:
318 return sign * RealGradient( 0.0, 0.0, 1.0 );
319 case 4:
320 return sign * RealGradient( 0.0, 0.0, -1.0 );
321 case 5:
322 return RealGradient();
323
324 default:
325 libmesh_error_msg("Invalid i = " << i);
326 } // switch(i)
327
328 } // j = 0
329
330 // d()/deta
331 case 1:
332 {
333 switch(i)
334 {
335 case 0:
336 case 1:
337 return sign * RealGradient( 1.0, 0.0, 0.0 );
338 case 2:
339 return sign * RealGradient( -1.0, 0.0, -1.0 );
340 case 3:
341 return sign * RealGradient( 0.0, 0.0, 1.0 );
342 case 4:
343 return RealGradient();
344 case 5:
345 return sign * RealGradient( 0.0, 0.0, -1.0 );
346
347 default:
348 libmesh_error_msg("Invalid i = " << i);
349 } // switch(i)
350
351 } // j = 1
352
353 // d()/dzeta
354 case 2:
355 {
356 switch(i)
357 {
358 case 0:
359 return sign * RealGradient( 1.0, 0.0, 0.0 );
360 case 1:
361 return RealGradient();
362 case 2:
363 case 5:
364 return sign * RealGradient( 0.0, 1.0, 0.0 );
365 case 3:
366 return sign * RealGradient( -1.0, -1.0, 0.0 );
367 case 4:
368 return sign * RealGradient( 1.0, 0.0, 0.0 );
369
370 default:
371 libmesh_error_msg("Invalid i = " << i);
372 } // switch(i)
373
374 } // j = 2
375
376 default:
377 libmesh_error_msg("Invalid j = " << j);
378 }
379 }
380
381 default:
382 libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(elem->type()));
383 }
384 }
385 // unsupported order
386 default:
387 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
397template <>
398RealGradient FE<3,NEDELEC_ONE>::shape_deriv(const ElemType,
399 const Order,
400 const unsigned int,
401 const unsigned int,
402 const Point &)
403{
404 libmesh_error_msg("Nedelec elements require the element type \nbecause edge orientation is needed.");
405 return RealGradient();
406}
407
408
409template <>
410RealGradient 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 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
423template <>
424RealGradient 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 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 libmesh_assert_less (j, 6);
442
443 const Order totalorder = order + add_p_level*elem->p_level();
444 libmesh_assert_less(i, n_dofs(elem->type(), totalorder));
445
446 const char sign = elem->positive_edge_orientation(i) ? 1 : -1;
447
448 switch (totalorder)
449 {
450 // linear Nedelec (first kind) shape function second derivatives
451 case FIRST:
452 {
453 switch (elem->type())
454 {
455 case HEX20:
456 case HEX27:
457 {
458 switch (j)
459 {
460 // d^2()/dxi^2, d^2()/deta^2, d^2()/dzeta^2
461 case 0:
462 case 2:
463 case 5:
464 return RealGradient();
465
466 // d^2()/dxideta
467 case 1:
468 {
469 switch(i)
470 {
471 case 0:
472 case 1:
473 case 2:
474 case 3:
475 case 8:
476 case 9:
477 case 10:
478 case 11:
479 return RealGradient();
480 case 4:
481 case 6:
482 return sign * RealGradient( 0.0, 0.0, -0.125 );
483 case 5:
484 case 7:
485 return sign * RealGradient( 0.0, 0.0, 0.125 );
486
487 default:
488 libmesh_error_msg("Invalid i = " << i);
489 } // switch(i)
490
491 } // j = 1
492
493 // d^2()/dxidzeta
494 case 3:
495 {
496 switch(i)
497 {
498 case 0:
499 case 2:
500 case 4:
501 case 5:
502 case 6:
503 case 7:
504 case 8:
505 case 10:
506 return RealGradient();
507 case 1:
508 case 3:
509 case 11:
510 return sign * RealGradient( 0.0, 0.125 );
511 case 9:
512 return sign * RealGradient( 0.0, -0.125, 0.0 );
513
514 default:
515 libmesh_error_msg("Invalid i = " << i);
516 } // switch(i)
517
518 } // j = 3
519
520 // d^2()/detadzeta
521 case 4:
522 {
523 switch(i)
524 {
525 case 0:
526 return sign * RealGradient( -0.125, 0.0, 0.0 );
527 case 1:
528 case 3:
529 case 4:
530 case 5:
531 case 6:
532 case 7:
533 case 9:
534 case 11:
535 return RealGradient();
536 case 2:
537 case 8:
538 case 10:
539 return sign * RealGradient( 0.125, 0.0, 0.0 );
540
541 default:
542 libmesh_error_msg("Invalid i = " << i);
543 } // switch(i)
544
545 } // j = 4
546
547 default:
548 libmesh_error_msg("Invalid j = " << j);
549 }
550 }
551
552 // All second derivatives for linear tets are zero.
553 case TET10:
554 case TET14:
555 return RealGradient();
556
557 default:
558 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 default:
566 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
577template <>
578RealGradient 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 libmesh_error_msg("Nedelec elements require the element type \nbecause edge orientation is needed.");
585 return RealGradient();
586}
587
588
589template <>
590RealGradient 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 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
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition elem.h:96
libMesh::Elem::positive_edge_orientation
bool positive_edge_orientation(const unsigned int i) const
Definition elem.C:3634
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Elem::p_level
unsigned int p_level() const
Definition elem.h:3125
libMesh::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition fe_type.h:197
libMesh::FEType::order
OrderWrapper order
The approximation order of the element (at 0 p-refinement level).
Definition fe_type.h:203
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::FE::shape_second_deriv
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::FE::shape_deriv
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:40
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition enum_elem_type.h:33
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition libmesh_common.h:569
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::TOLERANCE
static constexpr Real TOLERANCE
Definition libmesh_common.h:151
libMesh::RealGradient
RealVectorValue RealGradient
Definition hp_coarsentest.h:49
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition libmesh_common.h:144
Generated on Fri Aug 21 2026 13:32:28 for libMesh by doxygen 1.9.8