libMesh
fe_szabab_shape_1D.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 
20 // C++ includes
21 
22 // Local includes
23 #include "libmesh/libmesh_config.h"
24 
25 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
26 
27 #include "libmesh/fe.h"
28 #include "libmesh/elem.h"
29 
30 namespace libMesh
31 {
32 
33 
34 template <>
36  const Order libmesh_dbg_var(order),
37  const unsigned int i,
38  const Point & p)
39 {
40  const Real xi = p(0);
41  const Real xi2 = xi*xi;
42 
43 
44  // Use this libmesh_assert rather than a switch with a single entry...
45  // It will go away in optimized mode, essentially has the same effect.
46  libmesh_assert_less_equal (order, SEVENTH);
47 
48  // switch (order)
49  // {
50  // case FIRST:
51  // case SECOND:
52  // case THIRD:
53  // case FOURTH:
54  // case FIFTH:
55  // case SIXTH:
56  // case SEVENTH:
57 
58  switch(i)
59  {
60  //nodal shape functions
61  case 0: return 1./2.-1./2.*xi;
62  case 1: return 1./2.+1./2.*xi;
63  case 2: return 1./4. *2.4494897427831780982*(xi2-1.);
64  case 3: return 1./4. *3.1622776601683793320*(xi2-1.)*xi;
65  case 4: return 1./16. *3.7416573867739413856*((5.*xi2-6.)*xi2+1.);
66  case 5: return 3./16. *1.4142135623730950488*(3.+(-10.+7.*xi2)*xi2)*xi;
67  case 6: return 1./32. *4.6904157598234295546*(-1.+(15.+(-35.+21.*xi2)*xi2)*xi2);
68  case 7: return 1./32. *5.0990195135927848300*(-5.+(35.+(-63.+33.*xi2)*xi2)*xi2)*xi;
69  case 8: return 1./256.*5.4772255750516611346*(5.+(-140.+(630.+(-924.+429.*xi2)*xi2)*xi2)*xi2);
70 
71  default:
72  libmesh_error_msg("Invalid shape function index!");
73  }
74 
75  libmesh_error_msg("We'll never get here!");
76  return 0.;
77 }
78 
79 
80 
81 template <>
83  const Order order,
84  const unsigned int i,
85  const Point & p)
86 {
87  libmesh_assert(elem);
88 
89  return FE<1,SZABAB>::shape(elem->type(), static_cast<Order>(order + elem->p_level()), i, p);
90 }
91 
92 
93 
94 template <>
96  const Order libmesh_dbg_var(order),
97  const unsigned int i,
98  const unsigned int libmesh_dbg_var(j),
99  const Point & p)
100 {
101  // only d()/dxi in 1D!
102  libmesh_assert_equal_to (j, 0);
103 
104  const Real xi = p(0);
105  const Real xi2 = xi*xi;
106 
107  // Use this libmesh_assert rather than a switch with a single entry...
108  // It will go away in optimized mode, essentially has the same effect.
109  libmesh_assert_less_equal (order, SEVENTH);
110 
111  // switch (order)
112  // {
113  // case FIRST:
114  // case SECOND:
115  // case THIRD:
116  // case FOURTH:
117  // case FIFTH:
118  // case SIXTH:
119  // case SEVENTH:
120 
121  switch(i)
122  {
123  case 0:return -1./2.;
124  case 1:return 1./2.;
125  case 2:return 1./2.*2.4494897427831780982*xi;
126  case 3:return -1./4.*3.1622776601683793320+3./4.*3.1622776601683793320*xi2;
127  case 4:return 1./16.*3.7416573867739413856*(-12.+20*xi2)*xi;
128  case 5:return 9./16.*1.4142135623730950488+(-45./8.*1.4142135623730950488+105./16.*1.4142135623730950488*xi2)*xi2;
129  case 6:return 1./32.*4.6904157598234295546*(30.+(-140.+126.*xi2)*xi2)*xi;
130  case 7:return -5./32.*5.0990195135927848300+(105./32.*5.0990195135927848300+(-315./32.*5.0990195135927848300+231./32.*5.0990195135927848300*xi2)*xi2)*xi2;
131  case 8:return 1./256.*5.4772255750516611346*(-280.+(2520.+(-5544.+3432.*xi2)*xi2)*xi2)*xi;
132 
133  default:
134  libmesh_error_msg("Invalid shape function index!");
135  }
136 
137  libmesh_error_msg("We'll never get here!");
138  return 0.;
139 }
140 
141 
142 
143 template <>
145  const Order order,
146  const unsigned int i,
147  const unsigned int j,
148  const Point & p)
149 {
150  libmesh_assert(elem);
151 
152  return FE<1,SZABAB>::shape_deriv(elem->type(),
153  static_cast<Order>(order + elem->p_level()), i, j, p);
154 }
155 
156 
157 
158 template <>
160  const Order,
161  const unsigned int,
162  const unsigned int,
163  const Point &)
164 {
165  static bool warning_given = false;
166 
167  if (!warning_given)
168  libMesh::err << "Second derivatives for Szabab elements "
169  << " are not yet implemented!"
170  << std::endl;
171 
172  warning_given = true;
173  return 0.;
174 }
175 
176 
177 
178 template <>
180  const Order,
181  const unsigned int,
182  const unsigned int,
183  const Point &)
184 {
185  static bool warning_given = false;
186 
187  if (!warning_given)
188  libMesh::err << "Second derivatives for Szabab elements "
189  << " are not yet implemented!"
190  << std::endl;
191 
192  warning_given = true;
193  return 0.;
194 }
195 
196 } // namespace libMesh
197 
198 
199 #endif //LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
OStreamProxy err
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)
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
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)