Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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 :
31 : namespace libMesh
32 : {
33 :
34 :
35 16092 : LIBMESH_DEFAULT_VECTORIZED_FE(1,SZABAB)
36 :
37 :
38 : template <>
39 1209384 : Real FE<1,SZABAB>::shape(const ElemType,
40 : const Order libmesh_dbg_var(order),
41 : const unsigned int i,
42 : const Point & p)
43 : {
44 1209384 : const Real xi = p(0);
45 1209384 : const Real xi2 = xi*xi;
46 :
47 :
48 : // Use this libmesh_assert rather than a switch with a single entry...
49 : // It will go away in optimized mode, essentially has the same effect.
50 100782 : libmesh_assert_less_equal (order, SEVENTH);
51 :
52 : // switch (order)
53 : // {
54 : // case FIRST:
55 : // case SECOND:
56 : // case THIRD:
57 : // case FOURTH:
58 : // case FIFTH:
59 : // case SIXTH:
60 : // case SEVENTH:
61 :
62 1209384 : switch(i)
63 : {
64 : //nodal shape functions
65 296208 : case 0: return 1./2.-1./2.*xi;
66 296208 : case 1: return 1./2.+1./2.*xi;
67 270792 : case 2: return 1./4. *2.4494897427831780982*(xi2-1.);
68 216396 : case 3: return 1./4. *3.1622776601683793320*(xi2-1.)*xi;
69 129780 : case 4: return 1./16. *3.7416573867739413856*((5.*xi2-6.)*xi2+1.);
70 0 : case 5: return 3./16. *1.4142135623730950488*(3.+(-10.+7.*xi2)*xi2)*xi;
71 0 : case 6: return 1./32. *4.6904157598234295546*(-1.+(15.+(-35.+21.*xi2)*xi2)*xi2);
72 0 : case 7: return 1./32. *5.0990195135927848300*(-5.+(35.+(-63.+33.*xi2)*xi2)*xi2)*xi;
73 0 : case 8: return 1./256.*5.4772255750516611346*(5.+(-140.+(630.+(-924.+429.*xi2)*xi2)*xi2)*xi2);
74 :
75 0 : default:
76 0 : libmesh_error_msg("Invalid shape function index!");
77 : }
78 : }
79 :
80 :
81 :
82 : template <>
83 13032 : Real FE<1,SZABAB>::shape(const Elem * elem,
84 : const Order order,
85 : const unsigned int i,
86 : const Point & p,
87 : const bool add_p_level)
88 : {
89 1086 : libmesh_assert(elem);
90 :
91 15204 : return FE<1,SZABAB>::shape(elem->type(), order + add_p_level*elem->p_level(), i, p);
92 : }
93 :
94 :
95 : template <>
96 0 : Real FE<1,SZABAB>::shape(const FEType fet,
97 : const Elem * elem,
98 : const unsigned int i,
99 : const Point & p,
100 : const bool add_p_level)
101 : {
102 0 : libmesh_assert(elem);
103 :
104 0 : return FE<1,SZABAB>::shape(elem->type(), fet.order + add_p_level*elem->p_level(), i, p);
105 : }
106 :
107 :
108 :
109 : template <>
110 396216 : Real FE<1,SZABAB>::shape_deriv(const ElemType,
111 : const Order libmesh_dbg_var(order),
112 : const unsigned int i,
113 : const unsigned int libmesh_dbg_var(j),
114 : const Point & p)
115 : {
116 : // only d()/dxi in 1D!
117 33018 : libmesh_assert_equal_to (j, 0);
118 :
119 396216 : const Real xi = p(0);
120 396216 : const Real xi2 = xi*xi;
121 :
122 : // Use this libmesh_assert rather than a switch with a single entry...
123 : // It will go away in optimized mode, essentially has the same effect.
124 33018 : libmesh_assert_less_equal (order, SEVENTH);
125 :
126 : // switch (order)
127 : // {
128 : // case FIRST:
129 : // case SECOND:
130 : // case THIRD:
131 : // case FOURTH:
132 : // case FIFTH:
133 : // case SIXTH:
134 : // case SEVENTH:
135 :
136 396216 : switch(i)
137 : {
138 8394 : case 0:return -1./2.;
139 100728 : case 1:return 1./2.;
140 87912 : case 2:return 1./2.*2.4494897427831780982*xi;
141 67788 : case 3:return -1./4.*3.1622776601683793320+3./4.*3.1622776601683793320*xi2;
142 39060 : case 4:return 1./16.*3.7416573867739413856*(-12.+20*xi2)*xi;
143 0 : case 5:return 9./16.*1.4142135623730950488+(-45./8.*1.4142135623730950488+105./16.*1.4142135623730950488*xi2)*xi2;
144 0 : case 6:return 1./32.*4.6904157598234295546*(30.+(-140.+126.*xi2)*xi2)*xi;
145 0 : case 7:return -5./32.*5.0990195135927848300+(105./32.*5.0990195135927848300+(-315./32.*5.0990195135927848300+231./32.*5.0990195135927848300*xi2)*xi2)*xi2;
146 0 : case 8:return 1./256.*5.4772255750516611346*(-280.+(2520.+(-5544.+3432.*xi2)*xi2)*xi2)*xi;
147 :
148 0 : default:
149 0 : libmesh_error_msg("Invalid shape function index!");
150 : }
151 : }
152 :
153 :
154 :
155 : template <>
156 6984 : Real FE<1,SZABAB>::shape_deriv(const Elem * elem,
157 : const Order order,
158 : const unsigned int i,
159 : const unsigned int j,
160 : const Point & p,
161 : const bool add_p_level)
162 : {
163 582 : libmesh_assert(elem);
164 :
165 7566 : return FE<1,SZABAB>::shape_deriv(elem->type(),
166 7566 : order + add_p_level*elem->p_level(), i, j, p);
167 : }
168 :
169 :
170 :
171 : template <>
172 0 : Real FE<1,SZABAB>::shape_deriv(const FEType fet,
173 : const Elem * elem,
174 : const unsigned int i,
175 : const unsigned int j,
176 : const Point & p,
177 : const bool add_p_level)
178 : {
179 0 : libmesh_assert(elem);
180 :
181 0 : return FE<1,SZABAB>::shape_deriv(elem->type(),
182 0 : fet.order + add_p_level*elem->p_level(),
183 : i,
184 : j,
185 0 : p);
186 : }
187 :
188 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
189 :
190 : template <>
191 0 : Real FE<1,SZABAB>::shape_second_deriv(const ElemType,
192 : const Order,
193 : const unsigned int,
194 : const unsigned int,
195 : const Point &)
196 : {
197 : static bool warning_given = false;
198 :
199 0 : if (!warning_given)
200 0 : libMesh::err << "Second derivatives for Szabab elements "
201 0 : << " are not yet implemented!"
202 0 : << std::endl;
203 :
204 0 : warning_given = true;
205 0 : return 0.;
206 : }
207 :
208 :
209 :
210 : template <>
211 0 : Real FE<1,SZABAB>::shape_second_deriv(const Elem *,
212 : const Order,
213 : const unsigned int,
214 : const unsigned int,
215 : const Point &,
216 : const bool)
217 : {
218 : static bool warning_given = false;
219 :
220 0 : if (!warning_given)
221 0 : libMesh::err << "Second derivatives for Szabab elements "
222 0 : << " are not yet implemented!"
223 0 : << std::endl;
224 :
225 0 : warning_given = true;
226 0 : return 0.;
227 : }
228 :
229 :
230 : template <>
231 0 : Real FE<1,SZABAB>::shape_second_deriv(const FEType fet,
232 : const Elem * elem,
233 : const unsigned int i,
234 : const unsigned int j,
235 : const Point & p,
236 : const bool add_p_level)
237 : {
238 0 : libmesh_assert(elem);
239 :
240 0 : return FE<1,SZABAB>::shape_second_deriv(elem->type(),
241 0 : fet.order + add_p_level*elem->p_level(),
242 : i,
243 : j,
244 0 : p);
245 : }
246 :
247 :
248 : #endif
249 :
250 : } // namespace libMesh
251 :
252 :
253 : #endif //LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
|
