libMesh
fe_nedelec_one_shape_2D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2024 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 // An excellent discussion of Nedelec shape functions is given in
28 // https://www.dealii.org/reports/nedelec/nedelec.pdf
29 template <>
31  const Order order,
32  const unsigned int i,
33  const Point & p,
34  const bool add_p_level)
35 {
36 #if LIBMESH_DIM > 1
37  libmesh_assert(elem);
38 
39  const Order total_order = static_cast<Order>(order + add_p_level * elem->p_level());
40 
41  switch (total_order)
42  {
43  // linear Nedelec (first kind) shape functions
44  case FIRST:
45  {
46  switch (elem->type())
47  {
48  case QUAD8:
49  case QUAD9:
50  {
51  libmesh_assert_less (i, 4);
52 
53  const Real xi = p(0);
54  const Real eta = p(1);
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 
62  switch(i)
63  {
64  case 0:
65  {
66  if (elem->point(0) > elem->point(1))
67  return RealGradient( -0.25*(1.0-eta), 0.0 );
68  else
69  return RealGradient( 0.25*(1.0-eta), 0.0 );
70  }
71  case 1:
72  {
73  if (elem->point(1) > elem->point(2))
74  return RealGradient( 0.0, -0.25*(1.0+xi) );
75  else
76  return RealGradient( 0.0, 0.25*(1.0+xi) );
77  }
78 
79  case 2:
80  {
81  if (elem->point(2) > elem->point(3))
82  return RealGradient( 0.25*(1.0+eta), 0.0 );
83  else
84  return RealGradient( -0.25*(1.0+eta), 0.0 );
85  }
86  case 3:
87  {
88  if (elem->point(3) > elem->point(0))
89  return RealGradient( 0.0, -0.25*(xi-1.0) );
90  else
91  return RealGradient( 0.0, 0.25*(xi-1.0) );
92  }
93 
94  default:
95  libmesh_error_msg("Invalid i = " << i);
96  }
97 
98  return RealGradient();
99  }
100 
101  case TRI6:
102  case TRI7:
103  {
104  const Real xi = p(0);
105  const Real eta = p(1);
106 
107  libmesh_assert_less (i, 3);
108 
109  switch(i)
110  {
111  case 0:
112  {
113  if (elem->point(0) > elem->point(1))
114  return RealGradient( -1.0+eta, -xi );
115  else
116  return RealGradient( 1.0-eta, xi );
117  }
118  case 1:
119  {
120  if (elem->point(1) > elem->point(2))
121  return RealGradient( eta, -xi );
122  else
123  return RealGradient( -eta, xi );
124  }
125 
126  case 2:
127  {
128  if (elem->point(2) > elem->point(0))
129  return RealGradient( eta, -xi+1.0 );
130  else
131  return RealGradient( -eta, xi-1.0 );
132  }
133 
134  default:
135  libmesh_error_msg("Invalid i = " << i);
136  }
137  }
138 
139  default:
140  libmesh_error_msg("ERROR: Unsupported 2D element type!: " << Utility::enum_to_string(elem->type()));
141  }
142  }
143 
144  // unsupported order
145  default:
146  libmesh_error_msg("ERROR: Unsupported 2D FE order!: " << total_order);
147  }
148 #else // LIBMESH_DIM > 1
149  libmesh_ignore(elem, order, i, p, add_p_level);
150  libmesh_not_implemented();
151 #endif
152 }
153 
154 
155 
156 
157 template <>
159  const Order,
160  const unsigned int,
161  const Point &)
162 {
163  libmesh_error_msg("Nedelec elements require the element type \nbecause edge orientation is needed.");
164  return RealGradient();
165 }
166 
167 
168 template <>
170  const Elem * elem,
171  const unsigned int i,
172  const Point & p,
173  const bool add_p_level)
174 {
175  return FE<2,NEDELEC_ONE>::shape(elem, fet.order, i, p, add_p_level);
176 }
177 
178 template <>
180  const Order order,
181  const unsigned int i,
182  const unsigned int j,
183  const Point &,
184  const bool add_p_level)
185 {
186 #if LIBMESH_DIM > 1
187  libmesh_assert(elem);
188  libmesh_assert_less (j, 2);
189 
190  const Order total_order = static_cast<Order>(order + add_p_level * elem->p_level());
191 
192  switch (total_order)
193  {
194  // linear Nedelec (first kind) shape function first derivatives
195  case FIRST:
196  {
197  switch (elem->type())
198  {
199  case QUAD8:
200  case QUAD9:
201  {
202  libmesh_assert_less (i, 4);
203 
204  switch (j)
205  {
206  // d()/dxi
207  case 0:
208  {
209  switch(i)
210  {
211  case 0:
212  case 2:
213  return RealGradient();
214  case 1:
215  {
216  if (elem->point(1) > elem->point(2))
217  return RealGradient( 0.0, -0.25 );
218  else
219  return RealGradient( 0.0, 0.25 );
220  }
221  case 3:
222  {
223  if (elem->point(3) > elem->point(0))
224  return RealGradient( 0.0, -0.25 );
225  else
226  return RealGradient( 0.0, 0.25 );
227  }
228  default:
229  libmesh_error_msg("Invalid i = " << i);
230  }
231  } // j = 0
232 
233  // d()/deta
234  case 1:
235  {
236  switch(i)
237  {
238  case 1:
239  case 3:
240  return RealGradient();
241  case 0:
242  {
243  if (elem->point(0) > elem->point(1))
244  return RealGradient( 0.25 );
245  else
246  return RealGradient( -0.25 );
247  }
248  case 2:
249  {
250  if (elem->point(2) > elem->point(3))
251  return RealGradient( 0.25 );
252  else
253  return RealGradient( -0.25 );
254  }
255  default:
256  libmesh_error_msg("Invalid i = " << i);
257  }
258  } // j = 1
259 
260  default:
261  libmesh_error_msg("Invalid j = " << j);
262  }
263 
264  return RealGradient();
265  }
266 
267  case TRI6:
268  case TRI7:
269  {
270  libmesh_assert_less (i, 3);
271 
272  // Account for edge flipping
273  Real f = 1.0;
274 
275  switch(i)
276  {
277  case 0:
278  {
279  if (elem->point(0) > elem->point(1))
280  f = -1.0;
281  break;
282  }
283  case 1:
284  {
285  if (elem->point(1) > elem->point(2))
286  f = -1.0;
287  break;
288  }
289  case 2:
290  {
291  if (elem->point(2) > elem->point(0))
292  f = -1.0;
293  break;
294  }
295  default:
296  libmesh_error_msg("Invalid i = " << i);
297  }
298 
299  switch (j)
300  {
301  // d()/dxi
302  case 0:
303  {
304  return RealGradient( 0.0, f*1.0);
305  }
306  // d()/deta
307  case 1:
308  {
309  return RealGradient( f*(-1.0) );
310  }
311  default:
312  libmesh_error_msg("Invalid j = " << j);
313  }
314  }
315 
316  default:
317  libmesh_error_msg("ERROR: Unsupported 2D element type!: " << Utility::enum_to_string(elem->type()));
318  }
319  }
320  // unsupported order
321  default:
322  libmesh_error_msg("ERROR: Unsupported 2D FE order!: " << total_order);
323  }
324 #else // LIBMESH_DIM > 1
325  libmesh_ignore(elem, order, i, j, add_p_level);
326  libmesh_not_implemented();
327 #endif
328 }
329 
330 
331 
332 template <>
334  const Order,
335  const unsigned int,
336  const unsigned int,
337  const Point &)
338 {
339  libmesh_error_msg("Nedelec elements require the element type \nbecause edge orientation is needed.");
340  return RealGradient();
341 }
342 
343 template <>
345  const Elem * elem,
346  const unsigned int i,
347  const unsigned int j,
348  const Point & p,
349  const bool add_p_level)
350 {
351  return FE<2,NEDELEC_ONE>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
352 }
353 
354 
355 
356 
357 
358 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
359 
360 template <>
362  const Order order,
363  const unsigned int libmesh_dbg_var(i),
364  const unsigned int libmesh_dbg_var(j),
365  const Point &,
366  const bool add_p_level)
367 {
368 #if LIBMESH_DIM > 1
369  libmesh_assert(elem);
370 
371  // j = 0 ==> d^2 phi / dxi^2
372  // j = 1 ==> d^2 phi / dxi deta
373  // j = 2 ==> d^2 phi / deta^2
374  libmesh_assert_less (j, 3);
375 
376  const Order total_order = static_cast<Order>(order + add_p_level * elem->p_level());
377 
378  switch (total_order)
379  {
380  // linear Nedelec (first kind) shape function second derivatives
381  case FIRST:
382  {
383  switch (elem->type())
384  {
385  case QUAD8:
386  case QUAD9:
387  {
388  libmesh_assert_less (i, 4);
389  // All second derivatives for linear quads are zero.
390  return RealGradient();
391  }
392 
393  case TRI6:
394  case TRI7:
395  {
396  libmesh_assert_less (i, 3);
397  // All second derivatives for linear triangles are zero.
398  return RealGradient();
399  }
400 
401  default:
402  libmesh_error_msg("ERROR: Unsupported 2D element type!: " << Utility::enum_to_string(elem->type()));
403 
404  } // end switch (type)
405  } // end case FIRST
406 
407  // unsupported order
408  default:
409  libmesh_error_msg("ERROR: Unsupported 2D FE order!: " << total_order);
410 
411  } // end switch (order)
412 
413 #else // LIBMESH_DIM > 1
414  libmesh_assert(true || i || j);
415  libmesh_ignore(elem, order, add_p_level);
416  libmesh_not_implemented();
417 #endif
418 }
419 
420 
421 
422 template <>
424  const Order,
425  const unsigned int,
426  const unsigned int,
427  const Point &)
428 {
429  libmesh_error_msg("Nedelec elements require the element type \nbecause edge orientation is needed.");
430  return RealGradient();
431 }
432 
433 
434 template <>
436  const Elem * elem,
437  const unsigned int i,
438  const unsigned int j,
439  const Point & p,
440  const bool add_p_level)
441 {
442  return FE<2,NEDELEC_ONE>::shape_second_deriv(elem, fet.order, i, j, p, add_p_level);
443 }
444 
445 
446 
447 #endif
448 
449 } // namespace libMesh
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:182
ElemType
Defines an enum for geometric element types.
RealVectorValue RealGradient
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
static constexpr Real TOLERANCE
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
unsigned int p_level() const
Definition: elem.h:2945
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:201
The libMesh namespace provides an interface to certain functionality in the library.
void libmesh_ignore(const Args &...)
libmesh_assert(ctx)
std::string enum_to_string(const T e)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual ElemType type() const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
const Point & point(const unsigned int i) const
Definition: elem.h:2277
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)