fe_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,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,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,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,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,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,HIERARCHIC>::shape(EDGE3, totalorder, i0, xi)*
218  FE<1,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,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,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,HIERARCHIC>::shape(elem, order, i, pp) -
278  FE<2,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,HIERARCHIC>::shape(elem, order, i, pp) -
288  FE<2,HIERARCHIC>::shape(elem, order, i, pm))/2./eps;
289  }
290 
291  default:
292  libmesh_error_msg("Invalid derivative index j = " << j);
293  }
294  }
295 
296  case QUAD4:
297  case QUADSHELL4:
298  libmesh_assert_less (totalorder, 2);
299  libmesh_fallthrough();
300  case QUAD8:
301  case QUADSHELL8:
302  case QUAD9:
303  {
304  // Compute quad shape functions as a tensor-product
305  const Real xi = p(0);
306  const Real eta = p(1);
307 
308  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
309 
310  // Example i, i0, i1 values for totalorder = 5:
311  // 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
312  // 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};
313  // 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};
314 
315  unsigned int i0, i1;
316 
317  // Vertex DoFs
318  if (i == 0)
319  { i0 = 0; i1 = 0; }
320  else if (i == 1)
321  { i0 = 1; i1 = 0; }
322  else if (i == 2)
323  { i0 = 1; i1 = 1; }
324  else if (i == 3)
325  { i0 = 0; i1 = 1; }
326  // Edge DoFs
327  else if (i < totalorder + 3u)
328  { i0 = i - 2; i1 = 0; }
329  else if (i < 2u*totalorder + 2)
330  { i0 = 1; i1 = i - totalorder - 1; }
331  else if (i < 3u*totalorder + 1u)
332  { i0 = i - 2u*totalorder; i1 = 1; }
333  else if (i < 4u*totalorder)
334  { i0 = 0; i1 = i - 3u*totalorder + 1; }
335  // Interior DoFs
336  else
337  {
338  unsigned int basisnum = i - 4*totalorder;
339  i0 = square_number_column[basisnum] + 2;
340  i1 = square_number_row[basisnum] + 2;
341  }
342 
343  // Flip odd degree of freedom values if necessary
344  // to keep continuity on sides
345  Real f = 1.;
346 
347  if ((i0%2) && (i0 > 2) && (i1 == 0))
348  f = (elem->point(0) > elem->point(1))?-1.:1.;
349  else if ((i0%2) && (i0>2) && (i1 == 1))
350  f = (elem->point(3) > elem->point(2))?-1.:1.;
351  else if ((i0 == 0) && (i1%2) && (i1>2))
352  f = (elem->point(0) > elem->point(3))?-1.:1.;
353  else if ((i0 == 1) && (i1%2) && (i1>2))
354  f = (elem->point(1) > elem->point(2))?-1.:1.;
355 
356  switch (j)
357  {
358  // d()/dxi
359  case 0:
360  return f*(FE<1,HIERARCHIC>::shape_deriv(EDGE3, totalorder, i0, 0, xi)*
361  FE<1,HIERARCHIC>::shape (EDGE3, totalorder, i1, eta));
362 
363  // d()/deta
364  case 1:
365  return f*(FE<1,HIERARCHIC>::shape (EDGE3, totalorder, i0, xi)*
366  FE<1,HIERARCHIC>::shape_deriv(EDGE3, totalorder, i1, 0, eta));
367 
368  default:
369  libmesh_error_msg("Invalid derivative index j = " << j);
370  }
371  }
372 
373  default:
374  libmesh_error_msg("ERROR: Unsupported element type = " << type);
375  }
376 
377  return 0.;
378 }
379 
380 
381 
382 template <>
383 Real FE<2,HIERARCHIC>::shape_second_deriv(const ElemType,
384  const Order,
385  const unsigned int,
386  const unsigned int,
387  const Point &)
388 {
389  libmesh_error_msg("Hierarchic polynomials require the element type \nbecause edge orientation is needed.");
390  return 0.;
391 }
392 
393 
394 
395 template <>
396 Real FE<2,HIERARCHIC>::shape_second_deriv(const Elem * elem,
397  const Order order,
398  const unsigned int i,
399  const unsigned int j,
400  const Point & p)
401 {
402  libmesh_assert(elem);
403 
404  // I have been lazy here and am using finite differences
405  // to compute the derivatives!
406  const Real eps = 1.e-6;
407  Point pp, pm;
408  unsigned int prevj = libMesh::invalid_uint;
409 
410  switch (j)
411  {
412  // d^2()/dxi^2
413  case 0:
414  {
415  pp = Point(p(0)+eps, p(1));
416  pm = Point(p(0)-eps, p(1));
417  prevj = 0;
418  break;
419  }
420 
421  // d^2()/dxideta
422  case 1:
423  {
424  pp = Point(p(0), p(1)+eps);
425  pm = Point(p(0), p(1)-eps);
426  prevj = 0;
427  break;
428  }
429 
430  // d^2()/deta^2
431  case 2:
432  {
433  pp = Point(p(0), p(1)+eps);
434  pm = Point(p(0), p(1)-eps);
435  prevj = 1;
436  break;
437  }
438  default:
439  libmesh_error_msg("Invalid derivative index j = " << j);
440  }
441 
442  return (FE<2,HIERARCHIC>::shape_deriv(elem, order, i, prevj, pp) -
443  FE<2,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