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 : // Local includes
19 : #include "libmesh/fe.h"
20 : #include "libmesh/elem.h"
21 : #include "libmesh/fe_interface.h"
22 : #include "libmesh/utility.h"
23 : #include "libmesh/enum_to_string.h"
24 :
25 : namespace
26 : {
27 : using namespace libMesh;
28 :
29 : // Compute the static coefficients for an element
30 70047714 : void hermite_compute_coefs(const Elem * elem, Real & d1xd1x, Real & d2xd2x)
31 : {
32 70047714 : const FEFamily mapping_family = FEMap::map_fe_type(*elem);
33 70047714 : const FEType map_fe_type(elem->default_order(), mapping_family);
34 :
35 : // Note: we explicitly don't consider the elem->p_level() when
36 : // computing the number of mapping shape functions.
37 : const int n_mapping_shape_functions =
38 70047714 : FEInterface::n_shape_functions (map_fe_type, /*extra_order=*/0, elem);
39 :
40 : // Degrees of freedom are at vertices and edge midpoints
41 70047714 : static const Point dofpt[2] = {Point(-1), Point(1)};
42 :
43 : // Mapping functions - first derivatives at each dofpt
44 75458015 : std::vector<Real> dxdxi(2);
45 70047714 : std::vector<Real> dxidx(2);
46 :
47 : FEInterface::shape_deriv_ptr shape_deriv_ptr =
48 70047714 : FEInterface::shape_deriv_function(map_fe_type, elem);
49 :
50 210143142 : for (int p = 0; p != 2; ++p)
51 : {
52 150879278 : dxdxi[p] = 0;
53 423146832 : for (int i = 0; i != n_mapping_shape_functions; ++i)
54 : {
55 : const Real ddxi = shape_deriv_ptr
56 283051404 : (map_fe_type, elem, i, 0, dofpt[p], /*add_p_level=*/false);
57 326651016 : dxdxi[p] += elem->point(i)(0) * ddxi;
58 : }
59 : }
60 :
61 : // Calculate derivative scaling factors
62 70047714 : d1xd1x = dxdxi[0];
63 70047714 : d2xd2x = dxdxi[1];
64 70047714 : }
65 :
66 :
67 : } // end anonymous namespace
68 :
69 :
70 : namespace libMesh
71 : {
72 :
73 :
74 15023404 : LIBMESH_DEFAULT_VECTORIZED_FE(1,HERMITE)
75 :
76 :
77 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
78 :
79 : template<>
80 25815432 : Real FEHermite<1>::hermite_raw_shape_second_deriv (const unsigned int i, const Real xi)
81 : {
82 : using Utility::pow;
83 :
84 25815432 : switch (i)
85 : {
86 6424143 : case 0:
87 6424143 : return 1.5 * xi;
88 6424143 : case 1:
89 6424143 : return -1.5 * xi;
90 6424143 : case 2:
91 6424143 : return 0.5 * (-1. + 3.*xi);
92 6424143 : case 3:
93 6424143 : return 0.5 * (1. + 3.*xi);
94 92852 : case 4:
95 92852 : return (8.*xi*xi + 4.*(xi*xi-1.))/24.;
96 11534 : case 5:
97 11534 : return (8.*xi*xi*xi + 12.*xi*(xi*xi-1.))/120.;
98 : // case 6:
99 : // return (8.*pow<4>(xi) + 20.*xi*xi*(xi*xi-1.) +
100 : // 2.*(xi*xi-1)*(xi*xi-1))/720.;
101 1116 : default:
102 1116 : Real denominator = 720., xipower = 1.;
103 23452 : for (unsigned n=6; n != i; ++n)
104 : {
105 8978 : xipower *= xi;
106 8978 : denominator *= (n+1);
107 2252 : }
108 15610 : return (8.*pow<4>(xi)*xipower +
109 16726 : (8.*(i-4)+4.)*xi*xi*xipower*(xi*xi-1.) +
110 14474 : (i-4)*(i-5)*xipower*(xi*xi-1.)*(xi*xi-1.))/denominator;
111 : }
112 : }
113 :
114 : #endif
115 :
116 :
117 : template<>
118 63094384 : Real FEHermite<1>::hermite_raw_shape_deriv(const unsigned int i, const Real xi)
119 : {
120 63094384 : switch (i)
121 : {
122 15703346 : case 0:
123 15703346 : return 0.75 * (-1. + xi*xi);
124 15703346 : case 1:
125 15703346 : return 0.75 * (1. - xi*xi);
126 15703346 : case 2:
127 15703346 : return 0.25 * (-1. - 2.*xi + 3.*xi*xi);
128 15703346 : case 3:
129 15703346 : return 0.25 * (-1. + 2.*xi + 3.*xi*xi);
130 228105 : case 4:
131 228105 : return 4.*xi * (xi*xi-1.)/24.;
132 23570 : case 5:
133 23570 : return (4*xi*xi*(xi*xi-1.) + (xi*xi-1.)*(xi*xi-1.))/120.;
134 : // case 6:
135 : // return (4*xi*xi*xi*(xi*xi-1.) + 2*xi*(xi*xi-1.)*(xi*xi-1.))/720.;
136 2444 : default:
137 2444 : Real denominator = 720., xipower = 1.;
138 48029 : for (unsigned n=6; n != i; ++n)
139 : {
140 18704 : xipower *= xi;
141 18704 : denominator *= (n+1);
142 4896 : }
143 34221 : return (4*xi*xi*xi*xipower*(xi*xi-1.) +
144 29325 : (i-4)*xi*xipower*(xi*xi-1.)*(xi*xi-1.))/denominator;
145 : }
146 : }
147 :
148 : template<>
149 168105690 : Real FEHermite<1>::hermite_raw_shape(const unsigned int i, const Real xi)
150 : {
151 168105690 : switch (i)
152 : {
153 41919202 : case 0:
154 41919202 : return 0.25 * (2. - 3.*xi + xi*xi*xi);
155 41919202 : case 1:
156 41919202 : return 0.25 * (2. + 3.*xi - xi*xi*xi);
157 41919202 : case 2:
158 41919202 : return 0.25 * (1. - xi - xi*xi + xi*xi*xi);
159 41919202 : case 3:
160 41919202 : return 0.25 * (-1. - xi + xi*xi + xi*xi*xi);
161 : // All high order terms have the form x^(p-4)(x^2-1)^2/p!
162 359071 : case 4:
163 359071 : return (xi*xi-1.) * (xi*xi-1.)/24.;
164 30857 : case 5:
165 30857 : return xi * (xi*xi-1.) * (xi*xi-1.)/120.;
166 : // case 6:
167 : // return xi*xi * (xi*xi-1.) * (xi*xi-1.)/720.;
168 3284 : default:
169 3284 : Real denominator = 720., xipower = 1.;
170 63856 : for (unsigned n=6; n != i; ++n)
171 : {
172 24902 : xipower *= xi;
173 24902 : denominator *= (n+1);
174 6553 : }
175 38954 : return (xi*xi*xipower*(xi*xi-1.)*(xi*xi-1.))/denominator;
176 :
177 : }
178 : }
179 :
180 :
181 : template <>
182 60547162 : Real FE<1,HERMITE>::shape(const Elem * elem,
183 : const Order libmesh_dbg_var(order),
184 : const unsigned int i,
185 : const Point & p,
186 : const bool libmesh_dbg_var(add_p_level))
187 : {
188 5182605 : libmesh_assert(elem);
189 :
190 : // Coefficient naming: d(1)d(2n) is the coefficient of the
191 : // global shape function corresponding to value 1 in terms of the
192 : // local shape function corresponding to normal derivative 2
193 : Real d1xd1x, d2xd2x;
194 :
195 60547162 : hermite_compute_coefs(elem, d1xd1x, d2xd2x);
196 :
197 60547162 : const ElemType type = elem->type();
198 :
199 : #ifndef NDEBUG
200 : const unsigned int totalorder =
201 5182605 : order + add_p_level * elem->p_level();
202 : #endif
203 :
204 60547162 : switch (type)
205 : {
206 : // C1 functions on the C1 cubic edge
207 60547162 : case EDGE2:
208 : case EDGE3:
209 : {
210 5182605 : libmesh_assert_less (i, totalorder+1);
211 :
212 60547162 : switch (i)
213 : {
214 15106370 : case 0:
215 15106370 : return FEHermite<1>::hermite_raw_shape(0, p(0));
216 15106370 : case 1:
217 15106370 : return d1xd1x * FEHermite<1>::hermite_raw_shape(2, p(0));
218 15106370 : case 2:
219 15106370 : return FEHermite<1>::hermite_raw_shape(1, p(0));
220 15106370 : case 3:
221 15106370 : return d2xd2x * FEHermite<1>::hermite_raw_shape(3, p(0));
222 121682 : default:
223 121682 : return FEHermite<1>::hermite_raw_shape(i, p(0));
224 : }
225 : }
226 0 : default:
227 0 : libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
228 : }
229 : }
230 :
231 :
232 :
233 : template <>
234 0 : Real FE<1,HERMITE>::shape(const ElemType,
235 : const Order,
236 : const unsigned int,
237 : const Point &)
238 : {
239 0 : libmesh_error_msg("Hermite elements require the real element \nto construct gradient-based degrees of freedom.");
240 : return 0.;
241 : }
242 :
243 :
244 : template <>
245 0 : Real FE<1,HERMITE>::shape(const FEType fet,
246 : const Elem * elem,
247 : const unsigned int i,
248 : const Point & p,
249 : const bool add_p_level)
250 : {
251 0 : return FE<1,HERMITE>::shape(elem, fet.order, i, p, add_p_level);
252 : }
253 :
254 :
255 : template <>
256 4891944 : Real FE<1,HERMITE>::shape_deriv(const Elem * elem,
257 : const Order libmesh_dbg_var(order),
258 : const unsigned int i,
259 : const unsigned int,
260 : const Point & p,
261 : const bool libmesh_dbg_var(add_p_level))
262 : {
263 125389 : libmesh_assert(elem);
264 :
265 : // Coefficient naming: d(1)d(2n) is the coefficient of the
266 : // global shape function corresponding to value 1 in terms of the
267 : // local shape function corresponding to normal derivative 2
268 : Real d1xd1x, d2xd2x;
269 :
270 4891944 : hermite_compute_coefs(elem, d1xd1x, d2xd2x);
271 :
272 4891944 : const ElemType type = elem->type();
273 :
274 : #ifndef NDEBUG
275 : const unsigned int totalorder =
276 125389 : order + add_p_level * elem->p_level();
277 : #endif
278 :
279 4891944 : switch (type)
280 : {
281 : // C1 functions on the C1 cubic edge
282 4891944 : case EDGE2:
283 : case EDGE3:
284 : {
285 125389 : libmesh_assert_less (i, totalorder+1);
286 :
287 4891944 : switch (i)
288 : {
289 1199386 : case 0:
290 1199386 : return FEHermite<1>::hermite_raw_shape_deriv(0, p(0));
291 1199386 : case 1:
292 1199386 : return d1xd1x * FEHermite<1>::hermite_raw_shape_deriv(2, p(0));
293 1199386 : case 2:
294 1199386 : return FEHermite<1>::hermite_raw_shape_deriv(1, p(0));
295 1199386 : case 3:
296 1199386 : return d2xd2x * FEHermite<1>::hermite_raw_shape_deriv(3, p(0));
297 94400 : default:
298 94400 : return FEHermite<1>::hermite_raw_shape_deriv(i, p(0));
299 : }
300 : }
301 0 : default:
302 0 : libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
303 : }
304 : }
305 :
306 :
307 : template <>
308 0 : Real FE<1,HERMITE>::shape_deriv(const ElemType,
309 : const Order,
310 : const unsigned int,
311 : const unsigned int,
312 : const Point &)
313 : {
314 0 : libmesh_error_msg("Hermite elements require the real element \nto construct gradient-based degrees of freedom.");
315 : return 0.;
316 : }
317 :
318 : template <>
319 0 : Real FE<1,HERMITE>::shape_deriv(const FEType fet,
320 : const Elem * elem,
321 : const unsigned int i,
322 : const unsigned int j,
323 : const Point & p,
324 : const bool add_p_level)
325 : {
326 0 : return FE<1,HERMITE>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
327 : }
328 :
329 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
330 :
331 :
332 : template <>
333 4608608 : Real FE<1,HERMITE>::shape_second_deriv(const Elem * elem,
334 : const Order libmesh_dbg_var(order),
335 : const unsigned int i,
336 : const unsigned int,
337 : const Point & p,
338 : const bool libmesh_dbg_var(add_p_level))
339 : {
340 102307 : libmesh_assert(elem);
341 :
342 : // Coefficient naming: d(1)d(2n) is the coefficient of the
343 : // global shape function corresponding to value 1 in terms of the
344 : // local shape function corresponding to normal derivative 2
345 : Real d1xd1x, d2xd2x;
346 :
347 4608608 : hermite_compute_coefs(elem, d1xd1x, d2xd2x);
348 :
349 4608608 : const ElemType type = elem->type();
350 :
351 : #ifndef NDEBUG
352 : const unsigned int totalorder =
353 102307 : order + add_p_level * elem->p_level();
354 : #endif
355 :
356 4608608 : switch (type)
357 : {
358 : // C1 functions on the C1 cubic edge
359 4608608 : case EDGE2:
360 : case EDGE3:
361 : {
362 102307 : libmesh_assert_less (i, totalorder+1);
363 :
364 4608608 : switch (i)
365 : {
366 1139687 : case 0:
367 1139687 : return FEHermite<1>::hermite_raw_shape_second_deriv(0, p(0));
368 1139687 : case 1:
369 1139687 : return d1xd1x * FEHermite<1>::hermite_raw_shape_second_deriv(2, p(0));
370 1139687 : case 2:
371 1139687 : return FEHermite<1>::hermite_raw_shape_second_deriv(1, p(0));
372 1139687 : case 3:
373 1139687 : return d2xd2x * FEHermite<1>::hermite_raw_shape_second_deriv(3, p(0));
374 49860 : default:
375 49860 : return FEHermite<1>::hermite_raw_shape_second_deriv(i, p(0));
376 : }
377 : }
378 0 : default:
379 0 : libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
380 : }
381 : }
382 :
383 :
384 : template <>
385 0 : Real FE<1,HERMITE>::shape_second_deriv(const ElemType,
386 : const Order,
387 : const unsigned int,
388 : const unsigned int,
389 : const Point &)
390 : {
391 0 : libmesh_error_msg("Hermite elements require the real element \nto construct gradient-based degrees of freedom.");
392 : return 0.;
393 : }
394 :
395 : template <>
396 0 : Real FE<1,HERMITE>::shape_second_deriv(const FEType fet,
397 : const Elem * elem,
398 : const unsigned int i,
399 : const unsigned int j,
400 : const Point & p,
401 : const bool add_p_level)
402 : {
403 0 : return FE<1,HERMITE>::shape_second_deriv(elem, fet.order, i, j, p, add_p_level);
404 : }
405 :
406 :
407 : #endif
408 :
409 : } // namespace libMesh
|
