libMesh
fe_nedelec_one_shape_2D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 // C++ includes
20 
21 // Local includes
22 #include "libmesh/fe.h"
23 #include "libmesh/elem.h"
24 
25 namespace libMesh
26 {
27 
28 template <>
30  const Order,
31  const unsigned int,
32  const Point &)
33 {
34  libmesh_error_msg("Nedelec elements require the element type \nbecause edge orientation is needed.");
35  return RealGradient();
36 }
37 
38 
39 // An excellent discussion of Nedelec shape functions is given in
40 // http://www.dealii.org/developer/reports/nedelec/nedelec.pdf
41 template <>
43  const Order order,
44  const unsigned int i,
45  const Point & p)
46 {
47 #if LIBMESH_DIM > 1
48  libmesh_assert(elem);
49 
50  const Order total_order = static_cast<Order>(order + elem->p_level());
51 
52  switch (total_order)
53  {
54  case FIRST:
55  {
56  switch (elem->type())
57  {
58  case QUAD8:
59  case QUAD9:
60  {
61  libmesh_assert_less (i, 4);
62 
63  const Real xi = p(0);
64  const Real eta = p(1);
65 
66  // Even with a loose inverse_map tolerance we ought to
67  // be nearly on the element interior in master
68  // coordinates
69  libmesh_assert_less_equal ( std::fabs(xi), 1.0+10*TOLERANCE );
70  libmesh_assert_less_equal ( std::fabs(eta), 1.0+10*TOLERANCE );
71 
72  switch(i)
73  {
74  case 0:
75  {
76  if (elem->point(0) > elem->point(1))
77  return RealGradient( -0.25*(1.0-eta), 0.0 );
78  else
79  return RealGradient( 0.25*(1.0-eta), 0.0 );
80  }
81  case 1:
82  {
83  if (elem->point(1) > elem->point(2))
84  return RealGradient( 0.0, -0.25*(1.0+xi) );
85  else
86  return RealGradient( 0.0, 0.25*(1.0+xi) );
87  }
88 
89  case 2:
90  {
91  if (elem->point(2) > elem->point(3))
92  return RealGradient( 0.25*(1.0+eta), 0.0 );
93  else
94  return RealGradient( -0.25*(1.0+eta), 0.0 );
95  }
96  case 3:
97  {
98  if (elem->point(3) > elem->point(0))
99  return RealGradient( 0.0, -0.25*(xi-1.0) );
100  else
101  return RealGradient( 0.0, 0.25*(xi-1.0) );
102  }
103 
104  default:
105  libmesh_error_msg("Invalid i = " << i);
106  }
107 
108  return RealGradient();
109  }
110 
111  case TRI6:
112  {
113  const Real xi = p(0);
114  const Real eta = p(1);
115 
116  libmesh_assert_less (i, 3);
117 
118  switch(i)
119  {
120  case 0:
121  {
122  if (elem->point(0) > elem->point(1))
123  return RealGradient( -1.0+eta, -xi );
124  else
125  return RealGradient( 1.0-eta, xi );
126  }
127  case 1:
128  {
129  if (elem->point(1) > elem->point(2))
130  return RealGradient( eta, -xi );
131  else
132  return RealGradient( -eta, xi );
133  }
134 
135  case 2:
136  {
137  if (elem->point(2) > elem->point(0))
138  return RealGradient( eta, -xi+1.0 );
139  else
140  return RealGradient( -eta, xi-1.0 );
141  }
142 
143  default:
144  libmesh_error_msg("Invalid i = " << i);
145  }
146  }
147 
148  default:
149  libmesh_error_msg("ERROR: Unsupported 2D element type!: " << elem->type());
150  }
151  }
152 
153  // unsupported order
154  default:
155  libmesh_error_msg("ERROR: Unsupported 2D FE order!: " << total_order);
156  }
157 #endif // LIBMESH_DIM > 1
158 
159  libmesh_error_msg("We'll never get here!");
160  return RealGradient();
161 }
162 
163 
164 
165 template <>
167  const Order,
168  const unsigned int,
169  const unsigned int,
170  const Point &)
171 {
172  libmesh_error_msg("Nedelec elements require the element type \nbecause edge orientation is needed.");
173  return RealGradient();
174 }
175 
176 
177 
178 template <>
180  const Order order,
181  const unsigned int i,
182  const unsigned int j,
183  const Point &)
184 {
185 #if LIBMESH_DIM > 1
186  libmesh_assert(elem);
187  libmesh_assert_less (j, 2);
188 
189  const Order total_order = static_cast<Order>(order + elem->p_level());
190 
191  switch (total_order)
192  {
193  // linear Lagrange shape functions
194  case FIRST:
195  {
196  switch (elem->type())
197  {
198  case QUAD8:
199  case QUAD9:
200  {
201  libmesh_assert_less (i, 4);
202 
203  switch (j)
204  {
205  // d()/dxi
206  case 0:
207  {
208  switch(i)
209  {
210  case 0:
211  case 2:
212  return RealGradient();
213  case 1:
214  {
215  if (elem->point(1) > elem->point(2))
216  return RealGradient( 0.0, -0.25 );
217  else
218  return RealGradient( 0.0, 0.25 );
219  }
220  case 3:
221  {
222  if (elem->point(3) > elem->point(0))
223  return RealGradient( 0.0, -0.25 );
224  else
225  return RealGradient( 0.0, 0.25 );
226  }
227  default:
228  libmesh_error_msg("Invalid i = " << i);
229  }
230  } // j=0
231 
232  // d()/deta
233  case 1:
234  {
235  switch(i)
236  {
237  case 1:
238  case 3:
239  return RealGradient();
240  case 0:
241  {
242  if (elem->point(0) > elem->point(1))
243  return RealGradient( 0.25 );
244  else
245  return RealGradient( -0.25 );
246  }
247  case 2:
248  {
249  if (elem->point(2) > elem->point(3))
250  return RealGradient( 0.25 );
251  else
252  return RealGradient( -0.25 );
253  }
254  default:
255  libmesh_error_msg("Invalid i = " << i);
256  }
257  } // j=1
258 
259  default:
260  libmesh_error_msg("Invalid j = " << j);
261  }
262 
263  return RealGradient();
264  }
265 
266  case TRI6:
267  {
268  libmesh_assert_less (i, 3);
269 
270  // Account for edge flipping
271  Real f = 1.0;
272 
273  switch(i)
274  {
275  case 0:
276  {
277  if (elem->point(0) > elem->point(1))
278  f = -1.0;
279  break;
280  }
281  case 1:
282  {
283  if (elem->point(1) > elem->point(2))
284  f = -1.0;
285  break;
286  }
287  case 2:
288  {
289  if (elem->point(2) > elem->point(0))
290  f = -1.0;
291  break;
292  }
293  default:
294  libmesh_error_msg("Invalid i = " << i);
295  }
296 
297  switch (j)
298  {
299  // d()/dxi
300  case 0:
301  {
302  return RealGradient( 0.0, f*1.0);
303  }
304  // d()/deta
305  case 1:
306  {
307  return RealGradient( f*(-1.0) );
308  }
309  default:
310  libmesh_error_msg("Invalid j = " << j);
311  }
312  }
313 
314  default:
315  libmesh_error_msg("ERROR: Unsupported 2D element type!: " << elem->type());
316  }
317  }
318  // unsupported order
319  default:
320  libmesh_error_msg("ERROR: Unsupported 2D FE order!: " << total_order);
321  }
322 #endif // LIBMESH_DIM > 1
323 
324  libmesh_error_msg("We'll never get here!");
325  return RealGradient();
326 }
327 
328 
329 
330 
331 template <>
333  const Order,
334  const unsigned int,
335  const unsigned int,
336  const Point &)
337 {
338  libmesh_error_msg("Nedelec elements require the element type \nbecause edge orientation is needed.");
339  return RealGradient();
340 }
341 
342 
343 
344 template <>
346  const Order order,
347  const unsigned int libmesh_dbg_var(i),
348  const unsigned int libmesh_dbg_var(j),
349  const Point &)
350 {
351 #if LIBMESH_DIM > 1
352  libmesh_assert(elem);
353 
354  // j = 0 ==> d^2 phi / dxi^2
355  // j = 1 ==> d^2 phi / dxi deta
356  // j = 2 ==> d^2 phi / deta^2
357  libmesh_assert_less (j, 3);
358 
359  const Order total_order = static_cast<Order>(order + elem->p_level());
360 
361  switch (total_order)
362  {
363  // linear Lagrange shape functions
364  case FIRST:
365  {
366  switch (elem->type())
367  {
368  case QUAD8:
369  case QUAD9:
370  {
371  libmesh_assert_less (i, 4);
372  // All second derivatives for linear quads are zero.
373  return RealGradient();
374  }
375 
376  case TRI6:
377  {
378  libmesh_assert_less (i, 3);
379  // All second derivatives for linear triangles are zero.
380  return RealGradient();
381  }
382 
383  default:
384  libmesh_error_msg("ERROR: Unsupported 2D element type!: " << elem->type());
385 
386  } // end switch (type)
387  } // end case FIRST
388 
389  // unsupported order
390  default:
391  libmesh_error_msg("ERROR: Unsupported 2D FE order!: " << total_order);
392 
393  } // end switch (order)
394 
395 #endif // LIBMESH_DIM > 1
396 
397  libmesh_error_msg("We'll never get here!");
398  return RealGradient();
399 }
400 
401 } // namespace libMesh
RealVectorValue RealGradient
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
virtual ElemType type() const =0
unsigned int p_level() const
Definition: elem.h:2422
ElemType
Defines an enum for geometric element types.
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
static const Real TOLERANCE
The libMesh namespace provides an interface to certain functionality in the library.
libmesh_assert(j)
PetscErrorCode Vec Mat libmesh_dbg_var(j)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Point & point(const unsigned int i) const
Definition: elem.h:1809
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)