libMesh
fe_hierarchic_shape_2D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2019 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/number_lookups.h"
23 
24 
25 // Anonymous namespace for functions shared by HIERARCHIC and
26 // L2_HIERARCHIC implementations. Implementations appear at the bottom
27 // of this file.
28 namespace
29 {
30 using namespace libMesh;
31 
32 Real fe_hierarchic_2D_shape(const Elem * elem,
33  const Order order,
34  const unsigned int i,
35  const Point & p,
36  const bool add_p_level);
37 
38 Real fe_hierarchic_2D_shape_deriv(const Elem * elem,
39  const Order order,
40  const unsigned int i,
41  const unsigned int j,
42  const Point & p,
43  const bool add_p_level);
44 
45 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
46 
47 Real fe_hierarchic_2D_shape_second_deriv(const Elem * elem,
48  const Order order,
49  const unsigned int i,
50  const unsigned int j,
51  const Point & p,
52  const bool add_p_level);
53 
54 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
55 
56 } // anonymous namespace
57 
58 
59 
60 namespace libMesh
61 {
62 
63 template <>
64 Real FE<2,HIERARCHIC>::shape(const ElemType,
65  const Order,
66  const unsigned int,
67  const Point &)
68 {
69  libmesh_error_msg("Hierarchic shape functions require an Elem for edge orientation.");
70  return 0.;
71 }
72 
73 
74 
75 template <>
76 Real FE<2,L2_HIERARCHIC>::shape(const ElemType,
77  const Order,
78  const unsigned int,
79  const Point &)
80 {
81  libmesh_error_msg("Hierarchic shape functions require an Elem for edge orientation.");
82  return 0.;
83 }
84 
85 
86 
87 template <>
88 Real FE<2,HIERARCHIC>::shape(const Elem * elem,
89  const Order order,
90  const unsigned int i,
91  const Point & p,
92  const bool add_p_level)
93 {
94  return fe_hierarchic_2D_shape(elem, order, i, p, add_p_level);
95 }
96 
97 
98 
99 template <>
100 Real FE<2,L2_HIERARCHIC>::shape(const Elem * elem,
101  const Order order,
102  const unsigned int i,
103  const Point & p,
104  const bool add_p_level)
105 {
106  return fe_hierarchic_2D_shape(elem, order, i, p, add_p_level);
107 }
108 
109 
110 
111 template <>
112 Real FE<2,HIERARCHIC>::shape_deriv(const ElemType,
113  const Order,
114  const unsigned int,
115  const unsigned int,
116  const Point &)
117 {
118  libmesh_error_msg("Hierarchic shape functions require an Elem for edge orientation.");
119  return 0.;
120 }
121 
122 
123 
124 template <>
125 Real FE<2,L2_HIERARCHIC>::shape_deriv(const ElemType,
126  const Order,
127  const unsigned int,
128  const unsigned int,
129  const Point &)
130 {
131  libmesh_error_msg("Hierarchic shape functions require an Elem for edge orientation.");
132  return 0.;
133 }
134 
135 
136 
137 template <>
138 Real FE<2,HIERARCHIC>::shape_deriv(const Elem * elem,
139  const Order order,
140  const unsigned int i,
141  const unsigned int j,
142  const Point & p,
143  const bool add_p_level)
144 {
145  return fe_hierarchic_2D_shape_deriv(elem, order, i, j, p, add_p_level);
146 }
147 
148 
149 
150 template <>
151 Real FE<2,L2_HIERARCHIC>::shape_deriv(const Elem * elem,
152  const Order order,
153  const unsigned int i,
154  const unsigned int j,
155  const Point & p,
156  const bool add_p_level)
157 {
158  return fe_hierarchic_2D_shape_deriv(elem, order, i, j, p, add_p_level);
159 }
160 
161 
162 
163 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
164 
165 template <>
166 Real FE<2,HIERARCHIC>::shape_second_deriv(const ElemType,
167  const Order,
168  const unsigned int,
169  const unsigned int,
170  const Point &)
171 {
172  libmesh_error_msg("Hierarchic shape functions require an Elem for edge orientation.");
173  return 0.;
174 }
175 
176 
177 
178 template <>
179 Real FE<2,L2_HIERARCHIC>::shape_second_deriv(const ElemType,
180  const Order,
181  const unsigned int,
182  const unsigned int,
183  const Point &)
184 {
185  libmesh_error_msg("Hierarchic shape functions require an Elem for edge orientation.");
186  return 0.;
187 }
188 
189 
190 
191 template <>
192 Real FE<2,HIERARCHIC>::shape_second_deriv(const Elem * elem,
193  const Order order,
194  const unsigned int i,
195  const unsigned int j,
196  const Point & p,
197  const bool add_p_level)
198 {
199  return fe_hierarchic_2D_shape_second_deriv(elem, order, i, j, p, add_p_level);
200 }
201 
202 
203 
204 template <>
205 Real FE<2,L2_HIERARCHIC>::shape_second_deriv(const Elem * elem,
206  const Order order,
207  const unsigned int i,
208  const unsigned int j,
209  const Point & p,
210  const bool add_p_level)
211 {
212  return fe_hierarchic_2D_shape_second_deriv(elem, order, i, j, p, add_p_level);
213 }
214 
215 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
216 
217 } // namespace libMesh
218 
219 
220 
221 namespace
222 {
223 using namespace libMesh;
224 
225 Real fe_hierarchic_2D_shape(const Elem * elem,
226  const Order order,
227  const unsigned int i,
228  const Point & p,
229  const bool add_p_level)
230 {
231  libmesh_assert(elem);
232 
233  const Order totalorder =
234  static_cast<Order>(order+add_p_level*elem->p_level());
235  libmesh_assert_greater (totalorder, 0);
236 
237  switch (elem->type())
238  {
239  case TRI3:
240  case TRISHELL3:
241  case TRI6:
242  {
243  const Real zeta1 = p(0);
244  const Real zeta2 = p(1);
245  const Real zeta0 = 1. - zeta1 - zeta2;
246 
247  libmesh_assert_less (i, (totalorder+1u)*(totalorder+2u)/2);
248  libmesh_assert (elem->type() == TRI6 || totalorder < 2);
249 
250  // Vertex DoFs
251  if (i == 0)
252  return zeta0;
253  else if (i == 1)
254  return zeta1;
255  else if (i == 2)
256  return zeta2;
257  // Edge DoFs
258  else if (i < totalorder + 2u)
259  {
260  // Avoid returning NaN on vertices!
261  if (zeta0 + zeta1 == 0.)
262  return 0.;
263 
264  const unsigned int basisorder = i - 1;
265  // Get factors to account for edge-flipping
266  Real f0 = 1;
267  if (basisorder%2 && (elem->point(0) > elem->point(1)))
268  f0 = -1.;
269 
270  Real edgeval = (zeta1 - zeta0) / (zeta1 + zeta0);
271  Real crossfunc = zeta0 + zeta1;
272  for (unsigned int n=1; n != basisorder; ++n)
273  crossfunc *= (zeta0 + zeta1);
274 
275  return f0 * crossfunc *
276  FE<1,HIERARCHIC>::shape(EDGE3, totalorder,
277  basisorder, edgeval);
278  }
279  else if (i < 2u*totalorder + 1)
280  {
281  // Avoid returning NaN on vertices!
282  if (zeta1 + zeta2 == 0.)
283  return 0.;
284 
285  const unsigned int basisorder = i - totalorder;
286  // Get factors to account for edge-flipping
287  Real f1 = 1;
288  if (basisorder%2 && (elem->point(1) > elem->point(2)))
289  f1 = -1.;
290 
291  Real edgeval = (zeta2 - zeta1) / (zeta2 + zeta1);
292  Real crossfunc = zeta2 + zeta1;
293  for (unsigned int n=1; n != basisorder; ++n)
294  crossfunc *= (zeta2 + zeta1);
295 
296  return f1 * crossfunc *
297  FE<1,HIERARCHIC>::shape(EDGE3, totalorder,
298  basisorder, edgeval);
299  }
300  else if (i < 3u*totalorder)
301  {
302  // Avoid returning NaN on vertices!
303  if (zeta0 + zeta2 == 0.)
304  return 0.;
305 
306  const unsigned int basisorder = i - (2u*totalorder) + 1;
307  // Get factors to account for edge-flipping
308  Real f2 = 1;
309  if (basisorder%2 && (elem->point(2) > elem->point(0)))
310  f2 = -1.;
311 
312  Real edgeval = (zeta0 - zeta2) / (zeta0 + zeta2);
313  Real crossfunc = zeta0 + zeta2;
314  for (unsigned int n=1; n != basisorder; ++n)
315  crossfunc *= (zeta0 + zeta2);
316 
317  return f2 * crossfunc *
318  FE<1,HIERARCHIC>::shape(EDGE3, totalorder,
319  basisorder, edgeval);
320  }
321  // Interior DoFs
322  else
323  {
324  const unsigned int basisnum = i - (3u*totalorder);
325  unsigned int exp0 = triangular_number_column[basisnum] + 1;
326  unsigned int exp1 = triangular_number_row[basisnum] + 1 -
327  triangular_number_column[basisnum];
328 
329  Real returnval = 1;
330  for (unsigned int n = 0; n != exp0; ++n)
331  returnval *= zeta0;
332  for (unsigned int n = 0; n != exp1; ++n)
333  returnval *= zeta1;
334  returnval *= zeta2;
335  return returnval;
336  }
337  }
338 
339  // Hierarchic shape functions on the quadrilateral.
340  case QUAD4:
341  case QUADSHELL4:
342  libmesh_assert_less (totalorder, 2);
343  libmesh_fallthrough();
344  case QUAD8:
345  case QUADSHELL8:
346  case QUAD9:
347  {
348  // Compute quad shape functions as a tensor-product
349  const Real xi = p(0);
350  const Real eta = p(1);
351 
352  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
353 
354  // Example i, i0, i1 values for totalorder = 5:
355  // 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
356  // 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};
357  // 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};
358 
359  unsigned int i0, i1;
360 
361  // Vertex DoFs
362  if (i == 0)
363  { i0 = 0; i1 = 0; }
364  else if (i == 1)
365  { i0 = 1; i1 = 0; }
366  else if (i == 2)
367  { i0 = 1; i1 = 1; }
368  else if (i == 3)
369  { i0 = 0; i1 = 1; }
370  // Edge DoFs
371  else if (i < totalorder + 3u)
372  { i0 = i - 2; i1 = 0; }
373  else if (i < 2u*totalorder + 2)
374  { i0 = 1; i1 = i - totalorder - 1; }
375  else if (i < 3u*totalorder + 1)
376  { i0 = i - 2u*totalorder; i1 = 1; }
377  else if (i < 4u*totalorder)
378  { i0 = 0; i1 = i - 3u*totalorder + 1; }
379  // Interior DoFs
380  else
381  {
382  unsigned int basisnum = i - 4*totalorder;
383  i0 = square_number_column[basisnum] + 2;
384  i1 = square_number_row[basisnum] + 2;
385  }
386 
387  // Flip odd degree of freedom values if necessary
388  // to keep continuity on sides
389  Real f = 1.;
390 
391  if ((i0%2) && (i0 > 2) && (i1 == 0))
392  f = (elem->point(0) > elem->point(1))?-1.:1.;
393  else if ((i0%2) && (i0>2) && (i1 == 1))
394  f = (elem->point(3) > elem->point(2))?-1.:1.;
395  else if ((i0 == 0) && (i1%2) && (i1>2))
396  f = (elem->point(0) > elem->point(3))?-1.:1.;
397  else if ((i0 == 1) && (i1%2) && (i1>2))
398  f = (elem->point(1) > elem->point(2))?-1.:1.;
399 
400  return f*(FE<1,HIERARCHIC>::shape(EDGE3, totalorder, i0, xi)*
401  FE<1,HIERARCHIC>::shape(EDGE3, totalorder, i1, eta));
402  }
403 
404  default:
405  libmesh_error_msg("ERROR: Unsupported element type = " << elem->type());
406  }
407 
408  return 0.;
409 }
410 
411 
412 
413 Real fe_hierarchic_2D_shape_deriv(const Elem * elem,
414  const Order order,
415  const unsigned int i,
416  const unsigned int j,
417  const Point & p,
418  const bool add_p_level)
419 {
420  libmesh_assert(elem);
421 
422  const ElemType type = elem->type();
423 
424  const Order totalorder =
425  static_cast<Order>(order+add_p_level*elem->p_level());
426 
427  libmesh_assert_greater (totalorder, 0);
428 
429  switch (type)
430  {
431  // 1st & 2nd-order Hierarchics.
432  case TRI3:
433  case TRISHELL3:
434  case TRI6:
435  {
436  const Real eps = 1.e-6;
437 
438  libmesh_assert_less (j, 2);
439 
440  switch (j)
441  {
442  // d()/dxi
443  case 0:
444  {
445  const Point pp(p(0)+eps, p(1));
446  const Point pm(p(0)-eps, p(1));
447 
448  return (FE<2,HIERARCHIC>::shape(elem, order, i, pp) -
449  FE<2,HIERARCHIC>::shape(elem, order, i, pm))/2./eps;
450  }
451 
452  // d()/deta
453  case 1:
454  {
455  const Point pp(p(0), p(1)+eps);
456  const Point pm(p(0), p(1)-eps);
457 
458  return (FE<2,HIERARCHIC>::shape(elem, order, i, pp) -
459  FE<2,HIERARCHIC>::shape(elem, order, i, pm))/2./eps;
460  }
461 
462  default:
463  libmesh_error_msg("Invalid derivative index j = " << j);
464  }
465  }
466 
467  case QUAD4:
468  case QUADSHELL4:
469  libmesh_assert_less (totalorder, 2);
470  libmesh_fallthrough();
471  case QUAD8:
472  case QUADSHELL8:
473  case QUAD9:
474  {
475  // Compute quad shape functions as a tensor-product
476  const Real xi = p(0);
477  const Real eta = p(1);
478 
479  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
480 
481  // Example i, i0, i1 values for totalorder = 5:
482  // 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
483  // 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};
484  // 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};
485 
486  unsigned int i0, i1;
487 
488  // Vertex DoFs
489  if (i == 0)
490  { i0 = 0; i1 = 0; }
491  else if (i == 1)
492  { i0 = 1; i1 = 0; }
493  else if (i == 2)
494  { i0 = 1; i1 = 1; }
495  else if (i == 3)
496  { i0 = 0; i1 = 1; }
497  // Edge DoFs
498  else if (i < totalorder + 3u)
499  { i0 = i - 2; i1 = 0; }
500  else if (i < 2u*totalorder + 2)
501  { i0 = 1; i1 = i - totalorder - 1; }
502  else if (i < 3u*totalorder + 1u)
503  { i0 = i - 2u*totalorder; i1 = 1; }
504  else if (i < 4u*totalorder)
505  { i0 = 0; i1 = i - 3u*totalorder + 1; }
506  // Interior DoFs
507  else
508  {
509  unsigned int basisnum = i - 4*totalorder;
510  i0 = square_number_column[basisnum] + 2;
511  i1 = square_number_row[basisnum] + 2;
512  }
513 
514  // Flip odd degree of freedom values if necessary
515  // to keep continuity on sides
516  Real f = 1.;
517 
518  if ((i0%2) && (i0 > 2) && (i1 == 0))
519  f = (elem->point(0) > elem->point(1))?-1.:1.;
520  else if ((i0%2) && (i0>2) && (i1 == 1))
521  f = (elem->point(3) > elem->point(2))?-1.:1.;
522  else if ((i0 == 0) && (i1%2) && (i1>2))
523  f = (elem->point(0) > elem->point(3))?-1.:1.;
524  else if ((i0 == 1) && (i1%2) && (i1>2))
525  f = (elem->point(1) > elem->point(2))?-1.:1.;
526 
527  switch (j)
528  {
529  // d()/dxi
530  case 0:
531  return f*(FE<1,HIERARCHIC>::shape_deriv(EDGE3, totalorder, i0, 0, xi)*
532  FE<1,HIERARCHIC>::shape (EDGE3, totalorder, i1, eta));
533 
534  // d()/deta
535  case 1:
536  return f*(FE<1,HIERARCHIC>::shape (EDGE3, totalorder, i0, xi)*
537  FE<1,HIERARCHIC>::shape_deriv(EDGE3, totalorder, i1, 0, eta));
538 
539  default:
540  libmesh_error_msg("Invalid derivative index j = " << j);
541  }
542  }
543 
544  default:
545  libmesh_error_msg("ERROR: Unsupported element type = " << type);
546  }
547 
548  return 0.;
549 }
550 
551 
552 
553 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
554 
555 Real fe_hierarchic_2D_shape_second_deriv(const Elem * elem,
556  const Order order,
557  const unsigned int i,
558  const unsigned int j,
559  const Point & p,
560  const bool add_p_level)
561 {
562  libmesh_assert(elem);
563 
564  // I have been lazy here and am using finite differences
565  // to compute the derivatives!
566  const Real eps = 1.e-6;
567  Point pp, pm;
568  unsigned int prevj = libMesh::invalid_uint;
569 
570  switch (j)
571  {
572  // d^2()/dxi^2
573  case 0:
574  {
575  pp = Point(p(0)+eps, p(1));
576  pm = Point(p(0)-eps, p(1));
577  prevj = 0;
578  break;
579  }
580 
581  // d^2()/dxideta
582  case 1:
583  {
584  pp = Point(p(0), p(1)+eps);
585  pm = Point(p(0), p(1)-eps);
586  prevj = 0;
587  break;
588  }
589 
590  // d^2()/deta^2
591  case 2:
592  {
593  pp = Point(p(0), p(1)+eps);
594  pm = Point(p(0), p(1)-eps);
595  prevj = 1;
596  break;
597  }
598  default:
599  libmesh_error_msg("Invalid derivative index j = " << j);
600  }
601 
602  return (FE<2,HIERARCHIC>::shape_deriv(elem, order, i, prevj, pp, add_p_level) -
603  FE<2,HIERARCHIC>::shape_deriv(elem, order, i, prevj, pm, add_p_level)
604  )/2./eps;
605 }
606 
607 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
608 
609 } // anonymous namespace
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)
libMesh::invalid_uint
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value.
Definition: libmesh.h:249
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
Definition: factoryfunction.C:55
libMesh::Order
Order
Definition: enum_order.h:40
libMesh::Elem::p_level
unsigned int p_level() const
Definition: elem.h:2512
libMesh::triangular_number_row
const unsigned char triangular_number_row[]
Definition: number_lookups.C:28
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition: elem.h:1955
libMesh::square_number_column
const unsigned char square_number_column[]
Definition: number_lookups.C:56
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
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::FE
A specific instantiation of the FEBase class.
Definition: fe.h:95
libMesh::QUAD4
Definition: enum_elem_type.h:41
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
libMesh::TRI3
Definition: enum_elem_type.h:39
libMesh::triangular_number_column
const unsigned char triangular_number_column[]
Definition: number_lookups.C:41
libMesh::TRI6
Definition: enum_elem_type.h:40
libMesh::QUADSHELL8
Definition: enum_elem_type.h:73
libMesh::EDGE3
Definition: enum_elem_type.h:36
libMesh::square_number_row
const unsigned char square_number_row[]
Definition: number_lookups.C:69
libMesh::QUADSHELL4
Definition: enum_elem_type.h:72
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:100
libMesh::QUAD9
Definition: enum_elem_type.h:43
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:121
libMesh::TRISHELL3
Definition: enum_elem_type.h:71
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::QUAD8
Definition: enum_elem_type.h:42
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
Generated on Sat Jan 25 2020 12:06:57 for libMesh by   doxygen 1.8.16