Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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 :
23 :
24 : // Anonymous namespace for functions shared by HIERARCHIC and
25 : // L2_HIERARCHIC implementations. Implementations appear at the bottom
26 : // of this file.
27 : namespace
28 : {
29 : using namespace libMesh;
30 :
31 : Real fe_hierarchic_1D_shape(const ElemType,
32 : const Order libmesh_dbg_var(order),
33 : const unsigned int i,
34 : const Point & p);
35 :
36 : Real fe_hierarchic_1D_shape_deriv(const ElemType,
37 : const Order libmesh_dbg_var(order),
38 : const unsigned int i,
39 : const unsigned int libmesh_dbg_var(j),
40 : const Point & p);
41 :
42 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
43 :
44 : Real fe_hierarchic_1D_shape_second_deriv(const ElemType,
45 : const Order libmesh_dbg_var(order),
46 : const unsigned int i,
47 : const unsigned int libmesh_dbg_var(j),
48 : const Point & p);
49 :
50 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
51 :
52 : } // anonymous namespace
53 :
54 :
55 :
56 : namespace libMesh
57 : {
58 :
59 :
60 126172 : LIBMESH_DEFAULT_VECTORIZED_FE(1,HIERARCHIC)
61 30888 : LIBMESH_DEFAULT_VECTORIZED_FE(1,L2_HIERARCHIC)
62 7200 : LIBMESH_DEFAULT_VECTORIZED_FE(1,SIDE_HIERARCHIC)
63 :
64 :
65 : template <>
66 3190779726 : Real FE<1,HIERARCHIC>::shape(const ElemType elem_type,
67 : const Order order,
68 : const unsigned int i,
69 : const Point & p)
70 : {
71 3190779726 : return fe_hierarchic_1D_shape(elem_type, order, i, p);
72 : }
73 :
74 :
75 :
76 : template <>
77 2718657436 : Real FE<1,L2_HIERARCHIC>::shape(const ElemType elem_type,
78 : const Order order,
79 : const unsigned int i,
80 : const Point & p)
81 : {
82 2718657436 : return fe_hierarchic_1D_shape(elem_type, order, i, p);
83 : }
84 :
85 :
86 :
87 : template <>
88 0 : Real FE<1,SIDE_HIERARCHIC>::shape(const ElemType,
89 : const Order,
90 : const unsigned int i,
91 : const Point & p)
92 : {
93 0 : unsigned int right_side = p(0) > 0; // 0 false, 1 true
94 0 : return (right_side == i);
95 : }
96 :
97 :
98 :
99 : template <>
100 125815 : Real FE<1,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 9623 : libmesh_assert(elem);
107 :
108 145061 : return fe_hierarchic_1D_shape(elem->type(), order + add_p_level*elem->p_level(), i, p);
109 : }
110 :
111 :
112 :
113 : template <>
114 0 : Real FE<1,HIERARCHIC>::shape(const FEType fet,
115 : const Elem * elem,
116 : const unsigned int i,
117 : const Point & p,
118 : const bool add_p_level)
119 : {
120 0 : libmesh_assert(elem);
121 0 : return fe_hierarchic_1D_shape(elem->type(), fet.order + add_p_level*elem->p_level(), i, p);
122 : }
123 :
124 :
125 :
126 :
127 :
128 : template <>
129 20784 : Real FE<1,L2_HIERARCHIC>::shape(const Elem * elem,
130 : const Order order,
131 : const unsigned int i,
132 : const Point & p,
133 : const bool add_p_level)
134 : {
135 1732 : libmesh_assert(elem);
136 :
137 24248 : return fe_hierarchic_1D_shape(elem->type(), order + add_p_level*elem->p_level(), i, p);
138 : }
139 :
140 : template <>
141 0 : Real FE<1,L2_HIERARCHIC>::shape(const FEType fet,
142 : const Elem * elem,
143 : const unsigned int i,
144 : const Point & p,
145 : const bool add_p_level)
146 : {
147 0 : libmesh_assert(elem);
148 0 : return fe_hierarchic_1D_shape(elem->type(), fet.order + add_p_level*elem->p_level(), i, p);
149 : }
150 :
151 :
152 :
153 : template <>
154 2400 : Real FE<1,SIDE_HIERARCHIC>::shape(const Elem *,
155 : const Order,
156 : const unsigned int i,
157 : const Point & p,
158 : const bool)
159 : {
160 2400 : unsigned int right_side = p(0) > 0; // 0 false, 1 true
161 2400 : return (right_side == i);
162 : }
163 :
164 : template <>
165 0 : Real FE<1,SIDE_HIERARCHIC>::shape(const FEType,
166 : const Elem *,
167 : const unsigned int i,
168 : const Point & p,
169 : const bool)
170 : {
171 0 : unsigned int right_side = p(0) > 0; // 0 false, 1 true
172 0 : return (right_side == i);
173 : }
174 :
175 :
176 : template <>
177 91593542 : Real FE<1,HIERARCHIC>::shape_deriv(const ElemType elem_type,
178 : const Order order,
179 : const unsigned int i,
180 : const unsigned int j,
181 : const Point & p)
182 : {
183 91593542 : return fe_hierarchic_1D_shape_deriv(elem_type, order, i, j, p);
184 : }
185 :
186 :
187 :
188 : template <>
189 5062176 : Real FE<1,L2_HIERARCHIC>::shape_deriv(const ElemType elem_type,
190 : const Order order,
191 : const unsigned int i,
192 : const unsigned int j,
193 : const Point & p)
194 : {
195 5062176 : return fe_hierarchic_1D_shape_deriv(elem_type, order, i, j, p);
196 : }
197 :
198 :
199 :
200 : template <>
201 0 : Real FE<1,SIDE_HIERARCHIC>::shape_deriv(const ElemType,
202 : const Order,
203 : const unsigned int,
204 : const unsigned int,
205 : const Point &)
206 : {
207 0 : return 0;
208 : }
209 :
210 :
211 :
212 : template <>
213 59372 : Real FE<1,HIERARCHIC>::shape_deriv(const Elem * elem,
214 : const Order order,
215 : const unsigned int i,
216 : const unsigned int j,
217 : const Point & p,
218 : const bool add_p_level)
219 : {
220 4574 : libmesh_assert(elem);
221 :
222 63946 : return fe_hierarchic_1D_shape_deriv(elem->type(),
223 63946 : order + add_p_level*elem->p_level(), i, j, p);
224 : }
225 :
226 :
227 :
228 : template <>
229 0 : Real FE<1,HIERARCHIC>::shape_deriv(const FEType fet,
230 : const Elem * elem,
231 : const unsigned int i,
232 : const unsigned int j,
233 : const Point & p,
234 : const bool add_p_level)
235 : {
236 0 : libmesh_assert(elem);
237 0 : return fe_hierarchic_1D_shape_deriv(elem->type(), fet.order + add_p_level*elem->p_level(), i, j, p);
238 : }
239 :
240 :
241 :
242 :
243 : template <>
244 12468 : Real FE<1,L2_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 : const bool add_p_level)
250 : {
251 1039 : libmesh_assert(elem);
252 :
253 13507 : return fe_hierarchic_1D_shape_deriv(elem->type(),
254 13507 : order + add_p_level*elem->p_level(), i, j, p);
255 : }
256 :
257 :
258 :
259 : template <>
260 0 : Real FE<1,L2_HIERARCHIC>::shape_deriv(const FEType fet,
261 : const Elem * elem,
262 : const unsigned int i,
263 : const unsigned int j,
264 : const Point & p,
265 : const bool add_p_level)
266 : {
267 0 : libmesh_assert(elem);
268 0 : return fe_hierarchic_1D_shape_deriv(elem->type(), fet.order + add_p_level*elem->p_level(), i, j, p);
269 : }
270 :
271 :
272 :
273 : template <>
274 2400 : Real FE<1,SIDE_HIERARCHIC>::shape_deriv(const Elem *,
275 : const Order,
276 : const unsigned int,
277 : const unsigned int,
278 : const Point &,
279 : const bool)
280 : {
281 2400 : return 0;
282 : }
283 :
284 :
285 :
286 : template <>
287 0 : Real FE<1,SIDE_HIERARCHIC>::shape_deriv(const FEType,
288 : const Elem *,
289 : const unsigned int,
290 : const unsigned int,
291 : const Point &,
292 : const bool)
293 : {
294 0 : return 0;
295 : }
296 :
297 :
298 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
299 :
300 : template <>
301 103680 : Real FE<1,HIERARCHIC>::shape_second_deriv(const ElemType elem_type,
302 : const Order order,
303 : const unsigned int i,
304 : const unsigned int j,
305 : const Point & p)
306 : {
307 103680 : return fe_hierarchic_1D_shape_second_deriv(elem_type, order, i, j, p);
308 : }
309 :
310 :
311 :
312 :
313 : template <>
314 0 : Real FE<1,L2_HIERARCHIC>::shape_second_deriv(const ElemType elem_type,
315 : const Order order,
316 : const unsigned int i,
317 : const unsigned int j,
318 : const Point & p)
319 : {
320 0 : return fe_hierarchic_1D_shape_second_deriv(elem_type, order, i, j, p);
321 : }
322 :
323 :
324 :
325 : template <>
326 0 : Real FE<1,SIDE_HIERARCHIC>::shape_second_deriv(const ElemType,
327 : const Order,
328 : const unsigned int,
329 : const unsigned int,
330 : const Point &)
331 : {
332 0 : return 0;
333 : }
334 :
335 :
336 :
337 : template <>
338 27130 : Real FE<1,HIERARCHIC>::shape_second_deriv(const Elem * elem,
339 : const Order order,
340 : const unsigned int i,
341 : const unsigned int j,
342 : const Point & p,
343 : const bool add_p_level)
344 : {
345 2101 : libmesh_assert(elem);
346 :
347 29231 : return fe_hierarchic_1D_shape_second_deriv(elem->type(),
348 29231 : order + add_p_level*elem->p_level(), i, j, p);
349 : }
350 :
351 :
352 :
353 : template <>
354 0 : Real FE<1,HIERARCHIC>::shape_second_deriv(const FEType fet,
355 : const Elem * elem,
356 : const unsigned int i,
357 : const unsigned int j,
358 : const Point & p,
359 : const bool add_p_level)
360 : {
361 0 : libmesh_assert(elem);
362 0 : return fe_hierarchic_1D_shape_second_deriv(elem->type(),
363 0 : fet.order + add_p_level*elem->p_level(), i, j, p);
364 : }
365 :
366 :
367 : template <>
368 12468 : Real FE<1,L2_HIERARCHIC>::shape_second_deriv(const Elem * elem,
369 : const Order order,
370 : const unsigned int i,
371 : const unsigned int j,
372 : const Point & p,
373 : const bool add_p_level)
374 : {
375 1039 : libmesh_assert(elem);
376 :
377 13507 : return fe_hierarchic_1D_shape_second_deriv(elem->type(),
378 13507 : order + add_p_level*elem->p_level(), i, j, p);
379 : }
380 :
381 :
382 : template <>
383 0 : Real FE<1,L2_HIERARCHIC>::shape_second_deriv(const FEType fet,
384 : const Elem * elem,
385 : const unsigned int i,
386 : const unsigned int j,
387 : const Point & p,
388 : const bool add_p_level)
389 : {
390 0 : libmesh_assert(elem);
391 0 : return fe_hierarchic_1D_shape_second_deriv(elem->type(),
392 0 : fet.order + add_p_level*elem->p_level(), i, j, p);
393 : }
394 :
395 :
396 : template <>
397 2400 : Real FE<1,SIDE_HIERARCHIC>::shape_second_deriv(const Elem *,
398 : const Order,
399 : const unsigned int,
400 : const unsigned int,
401 : const Point &,
402 : const bool)
403 : {
404 2400 : return 0.;
405 : }
406 :
407 :
408 : template <>
409 0 : Real FE<1,SIDE_HIERARCHIC>::shape_second_deriv(const FEType,
410 : const Elem *,
411 : const unsigned int,
412 : const unsigned int,
413 : const Point &,
414 : const bool)
415 : {
416 0 : return 0.;
417 : }
418 :
419 : #endif
420 :
421 : } // namespace libMesh
422 :
423 :
424 :
425 : namespace
426 : {
427 : using namespace libMesh;
428 :
429 5422116432 : Real fe_hierarchic_1D_shape(const ElemType,
430 : const Order order,
431 : const unsigned int i,
432 : const Point & p)
433 : {
434 487467329 : libmesh_assert_less (i, order+1u);
435 :
436 : // If we were to define p=0 here, it wouldn't be hierarchic
437 5422116432 : libmesh_error_msg_if (order <= 0,
438 : "HIERARCHIC FE families do not support p=0");
439 :
440 5422116432 : const Real xi = p(0);
441 :
442 487467329 : Real returnval = 1.;
443 :
444 5422116432 : switch (i)
445 : {
446 : // The vertex functions, left unscaled so that their coefficients remain the values of the
447 : // finite element solution at the two vertices
448 997638657 : case 0:
449 997638657 : returnval = .5*(1. - xi);
450 997638657 : break;
451 997620737 : case 1:
452 997620737 : returnval = .5*(1. + xi);
453 997620737 : break;
454 :
455 : // The bubble functions, xi^p - 1 for even p and xi^p - xi for odd, each vanishing at both
456 : // vertices and scaled to unit H1 seminorm over the interval
457 308689493 : default:
458 12643593820 : for (unsigned int n=1; n <= i; ++n)
459 9525426275 : returnval *= xi;
460 :
461 3426857038 : returnval = (returnval - ((i % 2) ? xi : 1.)) *
462 3426857038 : fe_hierarchic_bubble_scaling(i);
463 3426857038 : break;
464 : }
465 :
466 5422116432 : return returnval;
467 : }
468 :
469 :
470 :
471 89491527 : Real fe_hierarchic_1D_shape_deriv(const ElemType,
472 : const Order order,
473 : const unsigned int i,
474 : const unsigned int libmesh_dbg_var(j),
475 : const Point & p)
476 : {
477 : // only d()/dxi in 1D!
478 7236031 : libmesh_assert_equal_to (j, 0);
479 7236031 : libmesh_assert_less (i, order+1u);
480 :
481 : // If we were to define p=0 here, it wouldn't be hierarchic
482 89491527 : libmesh_error_msg_if (order <= 0,
483 : "HIERARCHIC FE families do not support p=0");
484 :
485 89491527 : const Real xi = p(0);
486 :
487 7236031 : Real returnval = 1.;
488 :
489 89491527 : switch (i)
490 : {
491 1893539 : case 0:
492 1893539 : returnval = -.5;
493 1893539 : break;
494 23348015 : case 1:
495 1893539 : returnval = .5;
496 23348015 : break;
497 :
498 : // The bubbles differentiate to p xi^(p-1), less the one an odd bubble's linear term
499 : // contributes, under the same scaling their shape functions carry
500 3448953 : default:
501 107702112 : for (unsigned int n=1; n != i; ++n)
502 68355568 : returnval *= xi;
503 :
504 42795497 : returnval = (Real(i) * returnval - ((i % 2) ? 1. : 0.)) *
505 42795497 : fe_hierarchic_bubble_scaling(i);
506 42795497 : break;
507 : }
508 :
509 89491527 : return returnval;
510 : }
511 :
512 :
513 :
514 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
515 :
516 131498 : Real fe_hierarchic_1D_shape_second_deriv(const ElemType,
517 : const Order order,
518 : const unsigned int i,
519 : const unsigned int libmesh_dbg_var(j),
520 : const Point & p)
521 : {
522 : // only d2()/d2xi in 1D!
523 11780 : libmesh_assert_equal_to (j, 0);
524 11780 : libmesh_assert_less (i, order+1u);
525 :
526 : // If we were to define p=0 here, it wouldn't be hierarchic
527 131498 : libmesh_error_msg_if (order <= 0,
528 : "HIERARCHIC FE families do not support p=0");
529 :
530 131498 : const Real xi = p(0);
531 :
532 11780 : Real returnval = 1.;
533 :
534 131498 : switch (i)
535 : {
536 6164 : case 0:
537 : case 1:
538 6164 : returnval = 0;
539 6164 : break;
540 :
541 : // Both parities differentiate twice to p (p-1) xi^(p-2), the linear term of an odd bubble
542 : // dropping out, under the same scaling their shape functions carry
543 5616 : default:
544 111846 : for (unsigned int n=2; n != i; ++n)
545 54884 : returnval *= xi;
546 :
547 125156 : returnval = Real(i) * (Real(i) - 1.) * returnval *
548 62578 : fe_hierarchic_bubble_scaling(i);
549 62578 : break;
550 : }
551 :
552 131498 : return returnval;
553 : }
554 :
555 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
556 :
557 : } // anonymous namespace
|
