fe_nedelec_one_shape_2D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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++ inlcludes
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 #if LIBMESH_DIM > 1
35  libMesh::err << "Nedelec elements require the element type\n"
36  << "because edge orientation is needed."
37  << std::endl;
38  libmesh_error();
39 #endif // LIBMESH_DIM > 1
40 
41  libmesh_error();
42  return RealGradient();
43 }
44 
45 
46 // An excellent discussion of Nedelec shape functions is given in
47 // http://www.dealii.org/developer/reports/nedelec/nedelec.pdf
48 template <>
50  const Order order,
51  const unsigned int i,
52  const Point& p)
53 {
54 #if LIBMESH_DIM > 1
55  libmesh_assert(elem);
56 
57  const Order total_order = static_cast<Order>(order + elem->p_level());
58 
59  switch (total_order)
60  {
61  case FIRST:
62  {
63  switch (elem->type())
64  {
65  case QUAD8:
66  case QUAD9:
67  {
68  libmesh_assert_less (i, 4);
69 
70  const Real xi = p(0);
71  const Real eta = p(1);
72 
73  // Even with a loose inverse_map tolerance we ought to
74  // be nearly on the element interior in master
75  // coordinates
76  libmesh_assert_less_equal ( std::fabs(xi), 1.0+10*TOLERANCE );
77  libmesh_assert_less_equal ( std::fabs(eta), 1.0+10*TOLERANCE );
78 
79  switch(i)
80  {
81  case 0:
82  {
83  if( elem->point(0) > elem->point(1) )
84  return RealGradient( -0.25*(1.0-eta), 0.0 );
85  else
86  return RealGradient( 0.25*(1.0-eta), 0.0 );
87  }
88  case 1:
89  {
90  if( elem->point(1) > elem->point(2) )
91  return RealGradient( 0.0, -0.25*(1.0+xi) );
92  else
93  return RealGradient( 0.0, 0.25*(1.0+xi) );
94  }
95 
96  case 2:
97  {
98  if( elem->point(2) > elem->point(3) )
99  return RealGradient( 0.25*(1.0+eta), 0.0 );
100  else
101  return RealGradient( -0.25*(1.0+eta), 0.0 );
102  }
103  case 3:
104  {
105  if( elem->point(3) > elem->point(0) )
106  return RealGradient( 0.0, -0.25*(xi-1.0) );
107  else
108  return RealGradient( 0.0, 0.25*(xi-1.0) );
109  }
110 
111  default:
112  libmesh_error();
113 
114  }
115 
116  return RealGradient();
117  }
118 
119  case TRI6:
120  {
121  const Real xi = p(0);
122  const Real eta = p(1);
123 
124  libmesh_assert_less (i, 3);
125 
126  switch(i)
127  {
128  case 0:
129  {
130  if( elem->point(0) > elem->point(1) )
131  return RealGradient( -1.0+eta, -xi );
132  else
133  return RealGradient( 1.0-eta, xi );
134  }
135  case 1:
136  {
137  if( elem->point(1) > elem->point(2) )
138  return RealGradient( eta, -xi );
139  else
140  return RealGradient( -eta, xi );
141  }
142 
143  case 2:
144  {
145  if( elem->point(2) > elem->point(0) )
146  return RealGradient( eta, -xi+1.0 );
147  else
148  return RealGradient( -eta, xi-1.0 );
149  }
150 
151  default:
152  libmesh_error();
153 
154  }
155  }
156 
157  default:
158  {
159  libMesh::err << "ERROR: Unsupported 2D element type!: " << elem->type()
160  << std::endl;
161  libmesh_error();
162  }
163  }
164  }
165 
166  // unsupported order
167  default:
168  {
169  libMesh::err << "ERROR: Unsupported 2D FE order!: " << total_order
170  << std::endl;
171  libmesh_error();
172  }
173  }
174 #endif // LIBMESH_DIM > 1
175 
176  libmesh_error();
177  return RealGradient();
178 }
179 
180 
181 
182 template <>
184  const Order,
185  const unsigned int,
186  const unsigned int,
187  const Point&)
188 {
189 #if LIBMESH_DIM > 1
190  libMesh::err << "Nedelec elements require the element type\n"
191  << "because edge orientation is needed."
192  << std::endl;
193  libmesh_error();
194 #endif // LIBMESH_DIM > 1
195 
196  libmesh_error();
197  return RealGradient();
198 }
199 
200 
201 
202 template <>
204  const Order order,
205  const unsigned int i,
206  const unsigned int j,
207  const Point&)
208 {
209 #if LIBMESH_DIM > 1
210  libmesh_assert(elem);
211  libmesh_assert_less (j, 2);
212 
213  const Order total_order = static_cast<Order>(order + elem->p_level());
214 
215  switch (total_order)
216  {
217  // linear Lagrange shape functions
218  case FIRST:
219  {
220  switch (elem->type())
221  {
222  case QUAD8:
223  case QUAD9:
224  {
225  libmesh_assert_less (i, 4);
226 
227  switch (j)
228  {
229  // d()/dxi
230  case 0:
231  {
232  switch(i)
233  {
234  case 0:
235  case 2:
236  return RealGradient();
237  case 1:
238  {
239  if( elem->point(1) > elem->point(2) )
240  return RealGradient( 0.0, -0.25 );
241  else
242  return RealGradient( 0.0, 0.25 );
243  }
244  case 3:
245  {
246  if( elem->point(3) > elem->point(0) )
247  return RealGradient( 0.0, -0.25 );
248  else
249  return RealGradient( 0.0, 0.25 );
250  }
251  default:
252  libmesh_error();
253  }
254  } // j=0
255 
256  // d()/deta
257  case 1:
258  {
259  switch(i)
260  {
261  case 1:
262  case 3:
263  return RealGradient();
264  case 0:
265  {
266  if( elem->point(0) > elem->point(1) )
267  return RealGradient( 0.25 );
268  else
269  return RealGradient( -0.25 );
270  }
271  case 2:
272  {
273  if( elem->point(2) > elem->point(3) )
274  return RealGradient( 0.25 );
275  else
276  return RealGradient( -0.25 );
277  }
278  default:
279  libmesh_error();
280  }
281  } // j=1
282 
283  default:
284  libmesh_error();
285  }
286 
287  return RealGradient();
288  }
289 
290  case TRI6:
291  {
292  libmesh_assert_less (i, 3);
293 
294  // Account for edge flipping
295  Real f = 1.0;
296 
297  switch(i)
298  {
299  case 0:
300  {
301  if( elem->point(0) > elem->point(1) )
302  f = -1.0;
303  break;
304  }
305  case 1:
306  {
307  if( elem->point(1) > elem->point(2) )
308  f = -1.0;
309  break;
310  }
311  case 2:
312  {
313  if( elem->point(2) > elem->point(0) )
314  f = -1.0;
315  break;
316  }
317  default:
318  libmesh_error();
319  }
320 
321  switch (j)
322  {
323  // d()/dxi
324  case 0:
325  {
326  return RealGradient( 0.0, f*1.0);
327  }
328  // d()/deta
329  case 1:
330  {
331  return RealGradient( f*(-1.0) );
332  }
333  default:
334  libmesh_error();
335  }
336  }
337 
338  default:
339  {
340  libMesh::err << "ERROR: Unsupported 2D element type!: " << elem->type()
341  << std::endl;
342  libmesh_error();
343  }
344  }
345  }
346  // unsupported order
347  default:
348  {
349  libMesh::err << "ERROR: Unsupported 2D FE order!: " << total_order
350  << std::endl;
351  libmesh_error();
352  }
353  }
354 #endif // LIBMESH_DIM > 1
355 
356  libmesh_error();
357  return RealGradient();
358 }
359 
360 
361 
362 
363 template <>
365  const Order,
366  const unsigned int,
367  const unsigned int,
368  const Point&)
369 {
370 #if LIBMESH_DIM > 1
371  libMesh::err << "Nedelec elements require the element type\n"
372  << "because edge orientation is needed."
373  << std::endl;
374  libmesh_error();
375 #endif // LIBMESH_DIM > 1
376 
377  libmesh_error();
378  return RealGradient();
379 }
380 
381 
382 
383 template <>
385  const Order order,
386  const unsigned int i,
387  const unsigned int j,
388  const Point&)
389 {
390 #if LIBMESH_DIM > 1
391  libmesh_assert(elem);
392 
393  // j = 0 ==> d^2 phi / dxi^2
394  // j = 1 ==> d^2 phi / dxi deta
395  // j = 2 ==> d^2 phi / deta^2
396  libmesh_assert_less (j, 3);
397 
398  const Order total_order = static_cast<Order>(order + elem->p_level());
399 
400  switch (total_order)
401  {
402  // linear Lagrange shape functions
403  case FIRST:
404  {
405  switch (elem->type())
406  {
407  case QUAD8:
408  case QUAD9:
409  {
410  libmesh_assert_less (i, 4);
411  // All second derivatives for linear quads are zero.
412  return RealGradient();
413  }
414 
415  case TRI6:
416  {
417  libmesh_assert_less (i, 3);
418  // All second derivatives for linear triangles are zero.
419  return RealGradient();
420  }
421 
422  default:
423  {
424  libMesh::err << "ERROR: Unsupported 2D element type!: " << elem->type()
425  << std::endl;
426  libmesh_error();
427  }
428 
429  } // end switch (type)
430  } // end case FIRST
431 
432  // unsupported order
433  default:
434  {
435  libMesh::err << "ERROR: Unsupported 2D FE order!: " << total_order
436  << std::endl;
437  libmesh_error();
438  }
439 
440  } // end switch (order)
441 
442 #endif // LIBMESH_DIM > 1
443 
444  libmesh_error();
445  return RealGradient();
446 }
447 
448 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo