fe_l2_hierarchic_shape_2D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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 #include "libmesh/number_lookups.h"
25 #include "libmesh/utility.h"
26 
27 
28 namespace libMesh
29 {
30 
31 template <>
32 Real FE<2,L2_HIERARCHIC>::shape(const ElemType,
33  const Order,
34  const unsigned int,
35  const Point &)
36 {
37  libmesh_error_msg("Hierarchic polynomials require the element type \nbecause edge orientation is needed.");
38  return 0.;
39 }
40 
41 
42 
43 template <>
44 Real FE<2,L2_HIERARCHIC>::shape(const Elem * elem,
45  const Order order,
46  const unsigned int i,
47  const Point & p)
48 {
49  libmesh_assert(elem);
50 
51  const Order totalorder = static_cast<Order>(order+elem->p_level());
52  libmesh_assert_greater (totalorder, 0);
53 
54  switch (elem->type())
55  {
56  case TRI3:
57  case TRISHELL3:
58  case TRI6:
59  {
60  const Real zeta1 = p(0);
61  const Real zeta2 = p(1);
62  const Real zeta0 = 1. - zeta1 - zeta2;
63 
64  libmesh_assert_less (i, (totalorder+1u)*(totalorder+2u)/2);
65  libmesh_assert (elem->type() == TRI6 || totalorder < 2);
66 
67  // Vertex DoFs
68  if (i == 0)
69  return zeta0;
70  else if (i == 1)
71  return zeta1;
72  else if (i == 2)
73  return zeta2;
74  // Edge DoFs
75  else if (i < totalorder + 2u)
76  {
77  // Avoid returning NaN on vertices!
78  if (zeta0 + zeta1 == 0.)
79  return 0.;
80 
81  const unsigned int basisorder = i - 1;
82  // Get factors to account for edge-flipping
83  Real f0 = 1;
84  if (basisorder%2 && (elem->point(0) > elem->point(1)))
85  f0 = -1.;
86 
87  Real edgeval = (zeta1 - zeta0) / (zeta1 + zeta0);
88  Real crossfunc = zeta0 + zeta1;
89  for (unsigned int n=1; n != basisorder; ++n)
90  crossfunc *= (zeta0 + zeta1);
91 
92  return f0 * crossfunc *
93  FE<1,L2_HIERARCHIC>::shape(EDGE3, totalorder,
94  basisorder, edgeval);
95  }
96  else if (i < 2u*totalorder + 1)
97  {
98  // Avoid returning NaN on vertices!
99  if (zeta1 + zeta2 == 0.)
100  return 0.;
101 
102  const unsigned int basisorder = i - totalorder;
103  // Get factors to account for edge-flipping
104  Real f1 = 1;
105  if (basisorder%2 && (elem->point(1) > elem->point(2)))
106  f1 = -1.;
107 
108  Real edgeval = (zeta2 - zeta1) / (zeta2 + zeta1);
109  Real crossfunc = zeta2 + zeta1;
110  for (unsigned int n=1; n != basisorder; ++n)
111  crossfunc *= (zeta2 + zeta1);
112 
113  return f1 * crossfunc *
114  FE<1,L2_HIERARCHIC>::shape(EDGE3, totalorder,
115  basisorder, edgeval);
116  }
117  else if (i < 3u*totalorder)
118  {
119  // Avoid returning NaN on vertices!
120  if (zeta0 + zeta2 == 0.)
121  return 0.;
122 
123  const unsigned int basisorder = i - (2u*totalorder) + 1;
124  // Get factors to account for edge-flipping
125  Real f2 = 1;
126  if (basisorder%2 && (elem->point(2) > elem->point(0)))
127  f2 = -1.;
128 
129  Real edgeval = (zeta0 - zeta2) / (zeta0 + zeta2);
130  Real crossfunc = zeta0 + zeta2;
131  for (unsigned int n=1; n != basisorder; ++n)
132  crossfunc *= (zeta0 + zeta2);
133 
134  return f2 * crossfunc *
135  FE<1,L2_HIERARCHIC>::shape(EDGE3, totalorder,
136  basisorder, edgeval);
137  }
138  // Interior DoFs
139  else
140  {
141  const unsigned int basisnum = i - (3u*totalorder);
142  unsigned int exp0 = triangular_number_column[basisnum] + 1;
143  unsigned int exp1 = triangular_number_row[basisnum] + 1 -
144  triangular_number_column[basisnum];
145 
146  Real returnval = 1;
147  for (unsigned int n = 0; n != exp0; ++n)
148  returnval *= zeta0;
149  for (unsigned int n = 0; n != exp1; ++n)
150  returnval *= zeta1;
151  returnval *= zeta2;
152  return returnval;
153  }
154  }
155 
156  // Hierarchic shape functions on the quadrilateral.
157  case QUAD4:
158  case QUADSHELL4:
159  libmesh_assert_less (totalorder, 2);
160  libmesh_fallthrough();
161  case QUAD8:
162  case QUADSHELL8:
163  case QUAD9:
164  {
165  // Compute quad shape functions as a tensor-product
166  const Real xi = p(0);
167  const Real eta = p(1);
168 
169  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
170 
171  // Example i, i0, i1 values for totalorder = 5:
172  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
173  // static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 0, 0, 0, 0, 2, 3, 3, 2, 4, 4, 4, 3, 2, 5, 5, 5, 5, 4, 3, 2};
174  // static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 2, 2, 3, 3, 2, 3, 4, 4, 4, 2, 3, 4, 5, 5, 5, 5};
175 
176  unsigned int i0, i1;
177 
178  // Vertex DoFs
179  if (i == 0)
180  { i0 = 0; i1 = 0; }
181  else if (i == 1)
182  { i0 = 1; i1 = 0; }
183  else if (i == 2)
184  { i0 = 1; i1 = 1; }
185  else if (i == 3)
186  { i0 = 0; i1 = 1; }
187  // Edge DoFs
188  else if (i < totalorder + 3u)
189  { i0 = i - 2; i1 = 0; }
190  else if (i < 2u*totalorder + 2)
191  { i0 = 1; i1 = i - totalorder - 1; }
192  else if (i < 3u*totalorder + 1)
193  { i0 = i - 2u*totalorder; i1 = 1; }
194  else if (i < 4u*totalorder)
195  { i0 = 0; i1 = i - 3u*totalorder + 1; }
196  // Interior DoFs
197  else
198  {
199  unsigned int basisnum = i - 4*totalorder;
200  i0 = square_number_column[basisnum] + 2;
201  i1 = square_number_row[basisnum] + 2;
202  }
203 
204  // Flip odd degree of freedom values if necessary
205  // to keep continuity on sides
206  Real f = 1.;
207 
208  if ((i0%2) && (i0 > 2) && (i1 == 0))
209  f = (elem->point(0) > elem->point(1))?-1.:1.;
210  else if ((i0%2) && (i0>2) && (i1 == 1))
211  f = (elem->point(3) > elem->point(2))?-1.:1.;
212  else if ((i0 == 0) && (i1%2) && (i1>2))
213  f = (elem->point(0) > elem->point(3))?-1.:1.;
214  else if ((i0 == 1) && (i1%2) && (i1>2))
215  f = (elem->point(1) > elem->point(2))?-1.:1.;
216 
217  return f*(FE<1,L2_HIERARCHIC>::shape(EDGE3, totalorder, i0, xi)*
218  FE<1,L2_HIERARCHIC>::shape(EDGE3, totalorder, i1, eta));
219  }
220 
221  default:
222  libmesh_error_msg("ERROR: Unsupported element type = " << elem->type());
223  }
224 
225  return 0.;
226 }
227 
228 
229 
230 template <>
231 Real FE<2,L2_HIERARCHIC>::shape_deriv(const ElemType,
232  const Order,
233  const unsigned int,
234  const unsigned int,
235  const Point &)
236 {
237  libmesh_error_msg("Hierarchic polynomials require the element type \nbecause edge orientation is needed.");
238  return 0.;
239 }
240 
241 
242 
243 template <>
244 Real FE<2,L2_HIERARCHIC>::shape_deriv(const Elem * elem,
245  const Order order,
246  const unsigned int i,
247  const unsigned int j,
248  const Point & p)
249 {
250  libmesh_assert(elem);
251 
252  const ElemType type = elem->type();
253 
254  const Order totalorder = static_cast<Order>(order+elem->p_level());
255 
256  libmesh_assert_greater (totalorder, 0);
257 
258  switch (type)
259  {
260  // 1st & 2nd-order Hierarchics.
261  case TRI3:
262  case TRISHELL3:
263  case TRI6:
264  {
265  const Real eps = 1.e-6;
266 
267  libmesh_assert_less (j, 2);
268 
269  switch (j)
270  {
271  // d()/dxi
272  case 0:
273  {
274  const Point pp(p(0)+eps, p(1));
275  const Point pm(p(0)-eps, p(1));
276 
277  return (FE<2,L2_HIERARCHIC>::shape(elem, order, i, pp) -
278  FE<2,L2_HIERARCHIC>::shape(elem, order, i, pm))/2./eps;
279  }
280 
281  // d()/deta
282  case 1:
283  {
284  const Point pp(p(0), p(1)+eps);
285  const Point pm(p(0), p(1)-eps);
286 
287  return (FE<2,L2_HIERARCHIC>::shape(elem, order, i, pp) -
288  FE<2,L2_HIERARCHIC>::shape(elem, order, i, pm))/2./eps;
289  }
290 
291 
292  default:
293  libmesh_error_msg("Invalid derivative index j = " << j);
294  }
295  }
296 
297  case QUAD4:
298  case QUADSHELL4:
299  libmesh_assert_less (totalorder, 2);
300  libmesh_fallthrough();
301  case QUAD8:
302  case QUADSHELL8:
303  case QUAD9:
304  {
305  // Compute quad shape functions as a tensor-product
306  const Real xi = p(0);
307  const Real eta = p(1);
308 
309  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
310 
311  // Example i, i0, i1 values for totalorder = 5:
312  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
313  // static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 0, 0, 0, 0, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5};
314  // static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5};
315 
316  unsigned int i0, i1;
317 
318  // Vertex DoFs
319  if (i == 0)
320  { i0 = 0; i1 = 0; }
321  else if (i == 1)
322  { i0 = 1; i1 = 0; }
323  else if (i == 2)
324  { i0 = 1; i1 = 1; }
325  else if (i == 3)
326  { i0 = 0; i1 = 1; }
327  // Edge DoFs
328  else if (i < totalorder + 3u)
329  { i0 = i - 2; i1 = 0; }
330  else if (i < 2u*totalorder + 2)
331  { i0 = 1; i1 = i - totalorder - 1; }
332  else if (i < 3u*totalorder + 1u)
333  { i0 = i - 2u*totalorder; i1 = 1; }
334  else if (i < 4u*totalorder)
335  { i0 = 0; i1 = i - 3u*totalorder + 1; }
336  // Interior DoFs
337  else
338  {
339  unsigned int basisnum = i - 4*totalorder;
340  i0 = square_number_column[basisnum] + 2;
341  i1 = square_number_row[basisnum] + 2;
342  }
343 
344  // Flip odd degree of freedom values if necessary
345  // to keep continuity on sides
346  Real f = 1.;
347 
348  if ((i0%2) && (i0 > 2) && (i1 == 0))
349  f = (elem->point(0) > elem->point(1))?-1.:1.;
350  else if ((i0%2) && (i0>2) && (i1 == 1))
351  f = (elem->point(3) > elem->point(2))?-1.:1.;
352  else if ((i0 == 0) && (i1%2) && (i1>2))
353  f = (elem->point(0) > elem->point(3))?-1.:1.;
354  else if ((i0 == 1) && (i1%2) && (i1>2))
355  f = (elem->point(1) > elem->point(2))?-1.:1.;
356 
357  switch (j)
358  {
359  // d()/dxi
360  case 0:
361  return f*(FE<1,L2_HIERARCHIC>::shape_deriv(EDGE3, totalorder, i0, 0, xi)*
362  FE<1,L2_HIERARCHIC>::shape (EDGE3, totalorder, i1, eta));
363 
364  // d()/deta
365  case 1:
366  return f*(FE<1,L2_HIERARCHIC>::shape (EDGE3, totalorder, i0, xi)*
367  FE<1,L2_HIERARCHIC>::shape_deriv(EDGE3, totalorder, i1, 0, eta));
368 
369  default:
370  libmesh_error_msg("Invalid derivative index j = " << j);
371  }
372  }
373 
374  default:
375  libmesh_error_msg("ERROR: Unsupported element type = " << type);
376  }
377 
378  return 0.;
379 }
380 
381 
382 
383 template <>
384 Real FE<2,L2_HIERARCHIC>::shape_second_deriv(const ElemType,
385  const Order,
386  const unsigned int,
387  const unsigned int,
388  const Point &)
389 {
390  libmesh_error_msg("Hierarchic polynomials require the element type \nbecause edge orientation is needed.");
391  return 0.;
392 }
393 
394 
395 
396 template <>
397 Real FE<2,L2_HIERARCHIC>::shape_second_deriv(const Elem * elem,
398  const Order order,
399  const unsigned int i,
400  const unsigned int j,
401  const Point & p)
402 {
403  libmesh_assert(elem);
404 
405  // I have been lazy here and am using finite differences
406  // to compute the derivatives!
407  const Real eps = 1.e-6;
408  Point pp, pm;
409  unsigned int prevj = libMesh::invalid_uint;
410 
411  switch (j)
412  {
413  // d^2()/dxi^2
414  case 0:
415  {
416  pp = Point(p(0)+eps, p(1));
417  pm = Point(p(0)-eps, p(1));
418  prevj = 0;
419  break;
420  }
421 
422  // d^2()/dxideta
423  case 1:
424  {
425  pp = Point(p(0), p(1)+eps);
426  pm = Point(p(0), p(1)-eps);
427  prevj = 0;
428  break;
429  }
430 
431  // d^2()/deta^2
432  case 2:
433  {
434  pp = Point(p(0), p(1)+eps);
435  pm = Point(p(0), p(1)-eps);
436  prevj = 1;
437  break;
438  }
439  default:
440  libmesh_error_msg("Invalid derivative index j = " << j);
441  }
442  return (FE<2,L2_HIERARCHIC>::shape_deriv(elem, order, i, prevj, pp) -
443  FE<2,L2_HIERARCHIC>::shape_deriv(elem, order, i, prevj, pm)
444  )/2./eps;
445 }
446 
447 } // namespace libMesh
libMesh::invalid_uint
const unsigned int invalid_uint
Definition: libmesh.h:245
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::Elem
The base class for all geometric element types.
Definition: elem.h:100
libMesh::triangular_number_row
const unsigned char triangular_number_row[]
Definition: number_lookups.C:28
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:2555
libMesh::square_number_column
const unsigned char square_number_column[]
Definition: number_lookups.C:56
libMesh::FE
Template class which generates the different FE families and orders.
Definition: fe.h:89
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131
libMesh::triangular_number_column
const unsigned char triangular_number_column[]
Definition: number_lookups.C:41
libMesh::square_number_row
const unsigned char square_number_row[]
Definition: number_lookups.C:69
number_lookups.h
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Point
A geometric point in (x,y,z) space.
Definition: point.h:38
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition: elem.h:1892
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 by doxygen