libMesh
fe_nedelec_one_shape_3D.C
Go to the documentation of this file.
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 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  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 
135 template <>
136 RealGradient 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 
146 template <>
147 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  return FE<3,NEDELEC_ONE>::shape(elem, fet.order, i, p, add_p_level);
154 }
155 
156 
157 template <>
158 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  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 
397 template <>
398 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  libmesh_error_msg("Nedelec elements require the element type \nbecause edge orientation is needed.");
405  return RealGradient();
406 }
407 
408 
409 template <>
410 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  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 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  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 
577 template <>
578 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  libmesh_error_msg("Nedelec elements require the element type \nbecause edge orientation is needed.");
585  return RealGradient();
586 }
587 
588 
589 template <>
590 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  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::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
libMesh::RealGradient
RealVectorValue RealGradient
Definition: hp_coarsentest.h:49
libMesh::Order
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::TOLERANCE
static constexpr Real TOLERANCE
Definition: libmesh_common.h:151
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
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::Elem::p_level
unsigned int p_level() const
Definition: elem.h:3108
libMesh::FEType::order
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:215
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition: libmesh_common.h:515
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::Elem::positive_edge_orientation
bool positive_edge_orientation(const unsigned int i) const
Definition: elem.C:3589
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
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)
Generated on Thu Jul 17 2025 01:29:07 for libMesh by   doxygen 1.8.14