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 : // libmesh includes
20 : #include "libmesh/fe.h"
21 : #include "libmesh/elem.h"
22 : #include "libmesh/fe_interface.h"
23 : #include "libmesh/enum_to_string.h"
24 :
25 : // Anonymous namespace for persistent variables.
26 : // This allows us to cache the global-to-local mapping transformation
27 : // This should also screw up multithreading royally
28 : namespace
29 : {
30 : using namespace libMesh;
31 :
32 : // Keep track of which element was most recently used to generate
33 : // cached data
34 : static dof_id_type old_elem_id = DofObject::invalid_id;
35 : static const Elem * old_elem_ptr = nullptr;
36 :
37 : // Coefficient naming: d(1)d(2n) is the coefficient of the
38 : // global shape function corresponding to value 1 in terms of the
39 : // local shape function corresponding to normal derivative 2
40 : static Real d1xd1x, d2xd2x;
41 :
42 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
43 :
44 : Real clough_raw_shape_second_deriv(const unsigned int basis_num,
45 : const unsigned int deriv_type,
46 : const Point & p);
47 :
48 : #endif
49 :
50 : Real clough_raw_shape_deriv(const unsigned int basis_num,
51 : const unsigned int deriv_type,
52 : const Point & p);
53 : Real clough_raw_shape(const unsigned int basis_num,
54 : const Point & p);
55 :
56 :
57 : // Compute the static coefficients for an element
58 0 : void clough_compute_coefs(const Elem * elem)
59 : {
60 : // Using static globals for old_elem_id, etc. will fail
61 : // horribly with more than one thread.
62 0 : libmesh_assert_equal_to (libMesh::n_threads(), 1);
63 :
64 : // Coefficients are cached from old elements; we rely on that cache
65 : // except in dbg mode
66 : #ifndef DEBUG
67 0 : if (elem->id() == old_elem_id &&
68 0 : elem == old_elem_ptr)
69 0 : return;
70 : #endif
71 :
72 0 : const FEFamily mapping_family = FEMap::map_fe_type(*elem);
73 0 : const FEType map_fe_type(elem->default_order(), mapping_family);
74 :
75 : // Note: we explicitly don't consider the elem->p_level() when
76 : // computing the number of mapping shape functions.
77 : const int n_mapping_shape_functions =
78 0 : FEInterface::n_shape_functions(map_fe_type, /*extra_order=*/0, elem);
79 :
80 : // Degrees of freedom are at vertices and edge midpoints
81 0 : std::vector<Point> dofpt;
82 0 : dofpt.push_back(Point(0));
83 0 : dofpt.push_back(Point(1));
84 :
85 : // Mapping functions - first derivatives at each dofpt
86 0 : std::vector<Real> dxdxi(2);
87 0 : std::vector<Real> dxidx(2);
88 :
89 : FEInterface::shape_deriv_ptr shape_deriv_ptr =
90 0 : FEInterface::shape_deriv_function(map_fe_type, elem);
91 :
92 0 : for (int p = 0; p != 2; ++p)
93 : {
94 0 : for (int i = 0; i != n_mapping_shape_functions; ++i)
95 : {
96 : const Real ddxi = shape_deriv_ptr
97 0 : (map_fe_type, elem, i, 0, dofpt[p], /*add_p_level=*/false);
98 0 : dxdxi[p] += dofpt[p](0) * ddxi;
99 : }
100 : }
101 :
102 : // Calculate derivative scaling factors
103 :
104 : #ifdef DEBUG
105 : // The cached factors should equal our calculations
106 0 : if (elem->id() == old_elem_id &&
107 0 : elem == old_elem_ptr)
108 : {
109 0 : libmesh_assert_equal_to(d1xd1x, dxdxi[0]);
110 0 : libmesh_assert_equal_to(d2xd2x, dxdxi[1]);
111 : }
112 : #endif
113 :
114 0 : old_elem_id = elem->id();
115 0 : old_elem_ptr = elem;
116 :
117 0 : d1xd1x = dxdxi[0];
118 0 : d2xd2x = dxdxi[1];
119 0 : }
120 :
121 :
122 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
123 :
124 : // Return shape function second derivatives on the unit interval
125 0 : Real clough_raw_shape_second_deriv(const unsigned int basis_num,
126 : const unsigned int deriv_type,
127 : const Point & p)
128 : {
129 0 : Real xi = p(0);
130 :
131 0 : switch (deriv_type)
132 : {
133 :
134 : // second derivative in xi-xi direction
135 0 : case 0:
136 0 : switch (basis_num)
137 : {
138 0 : case 0:
139 0 : return -6 + 12*xi;
140 0 : case 1:
141 0 : return 6 - 12*xi;
142 0 : case 2:
143 0 : return -4 + 6*xi;
144 0 : case 3:
145 0 : return -2 + 6*xi;
146 :
147 0 : default:
148 0 : libmesh_error_msg("Invalid shape function index i = " <<
149 : basis_num);
150 : }
151 :
152 0 : default:
153 0 : libmesh_error_msg("Invalid shape function derivative j = " <<
154 : deriv_type);
155 : }
156 : }
157 :
158 : #endif
159 :
160 :
161 0 : Real clough_raw_shape_deriv(const unsigned int basis_num,
162 : const unsigned int deriv_type,
163 : const Point & p)
164 : {
165 0 : Real xi = p(0);
166 :
167 0 : switch (deriv_type)
168 : {
169 0 : case 0:
170 0 : switch (basis_num)
171 : {
172 0 : case 0:
173 0 : return -6*xi + 6*xi*xi;
174 0 : case 1:
175 0 : return 6*xi - 6*xi*xi;
176 0 : case 2:
177 0 : return 1 - 4*xi + 3*xi*xi;
178 0 : case 3:
179 0 : return -2*xi + 3*xi*xi;
180 :
181 0 : default:
182 0 : libmesh_error_msg("Invalid shape function index i = " <<
183 : basis_num);
184 : }
185 :
186 0 : default:
187 0 : libmesh_error_msg("Invalid shape function derivative j = " <<
188 : deriv_type);
189 : }
190 : }
191 :
192 0 : Real clough_raw_shape(const unsigned int basis_num,
193 : const Point & p)
194 : {
195 0 : Real xi = p(0);
196 :
197 0 : switch (basis_num)
198 : {
199 0 : case 0:
200 0 : return 1 - 3*xi*xi + 2*xi*xi*xi;
201 0 : case 1:
202 0 : return 3*xi*xi - 2*xi*xi*xi;
203 0 : case 2:
204 0 : return xi - 2*xi*xi + xi*xi*xi;
205 0 : case 3:
206 0 : return -xi*xi + xi*xi*xi;
207 :
208 0 : default:
209 0 : libmesh_error_msg("Invalid shape function index i = " <<
210 : basis_num);
211 : }
212 : }
213 :
214 :
215 : } // end anonymous namespace
216 :
217 :
218 : namespace libMesh
219 : {
220 :
221 :
222 0 : LIBMESH_DEFAULT_VECTORIZED_FE(1,CLOUGH)
223 :
224 :
225 : template <>
226 0 : Real FE<1,CLOUGH>::shape(const Elem * elem,
227 : const Order order,
228 : const unsigned int i,
229 : const Point & p,
230 : const bool add_p_level)
231 : {
232 0 : libmesh_assert(elem);
233 :
234 0 : clough_compute_coefs(elem);
235 :
236 0 : const ElemType type = elem->type();
237 :
238 : const Order totalorder =
239 0 : order + add_p_level*elem->p_level();
240 :
241 0 : switch (totalorder)
242 : {
243 : // 3rd-order C1 cubic element
244 0 : case THIRD:
245 : {
246 0 : switch (type)
247 : {
248 : // C1 functions on the C1 cubic edge
249 0 : case EDGE2:
250 : case EDGE3:
251 : {
252 0 : libmesh_assert_less (i, 4);
253 :
254 0 : switch (i)
255 : {
256 0 : case 0:
257 0 : return clough_raw_shape(0, p);
258 0 : case 1:
259 0 : return clough_raw_shape(1, p);
260 0 : case 2:
261 0 : return d1xd1x * clough_raw_shape(2, p);
262 0 : case 3:
263 0 : return d2xd2x * clough_raw_shape(3, p);
264 0 : default:
265 0 : libmesh_error_msg("Invalid shape function index i = " << i);
266 : }
267 : }
268 0 : default:
269 0 : libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
270 : }
271 : }
272 : // by default throw an error
273 0 : default:
274 0 : libmesh_error_msg("ERROR: Unsupported polynomial order = " << totalorder);
275 : }
276 : }
277 :
278 :
279 :
280 : template <>
281 0 : Real FE<1,CLOUGH>::shape(const ElemType,
282 : const Order,
283 : const unsigned int,
284 : const Point &)
285 : {
286 0 : libmesh_error_msg("Clough-Tocher elements require the real element \nto construct gradient-based degrees of freedom.");
287 : return 0.;
288 : }
289 :
290 : template <>
291 0 : Real FE<1,CLOUGH>::shape(const FEType fet,
292 : const Elem * elem,
293 : const unsigned int i,
294 : const Point & p,
295 : const bool add_p_level)
296 : {
297 0 : return FE<1,CLOUGH>::shape(elem, fet.order, i, p, add_p_level);
298 : }
299 :
300 :
301 :
302 :
303 : template <>
304 0 : Real FE<1,CLOUGH>::shape_deriv(const ElemType,
305 : const Order,
306 : const unsigned int,
307 : const unsigned int,
308 : const Point &)
309 : {
310 0 : libmesh_error_msg("Clough-Tocher elements require the real element \nto construct gradient-based degrees of freedom.");
311 : return 0.;
312 : }
313 :
314 :
315 :
316 : template <>
317 0 : Real FE<1,CLOUGH>::shape_deriv(const Elem * elem,
318 : const Order order,
319 : const unsigned int i,
320 : const unsigned int j,
321 : const Point & p,
322 : const bool add_p_level)
323 : {
324 0 : libmesh_assert(elem);
325 :
326 0 : clough_compute_coefs(elem);
327 :
328 0 : const ElemType type = elem->type();
329 :
330 : const Order totalorder =
331 0 : order + add_p_level*elem->p_level();
332 :
333 0 : switch (totalorder)
334 : {
335 : // 3rd-order C1 cubic element
336 0 : case THIRD:
337 : {
338 0 : switch (type)
339 : {
340 : // C1 functions on the C1 cubic edge
341 0 : case EDGE2:
342 : case EDGE3:
343 : {
344 0 : switch (i)
345 : {
346 0 : case 0:
347 0 : return clough_raw_shape_deriv(0, j, p);
348 0 : case 1:
349 0 : return clough_raw_shape_deriv(1, j, p);
350 0 : case 2:
351 0 : return d1xd1x * clough_raw_shape_deriv(2, j, p);
352 0 : case 3:
353 0 : return d2xd2x * clough_raw_shape_deriv(3, j, p);
354 0 : default:
355 0 : libmesh_error_msg("Invalid shape function index i = " << i);
356 : }
357 : }
358 0 : default:
359 0 : libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
360 : }
361 : }
362 : // by default throw an error
363 0 : default:
364 0 : libmesh_error_msg("ERROR: Unsupported polynomial order = " << totalorder);
365 : }
366 : }
367 :
368 :
369 : template <>
370 0 : Real FE<1,CLOUGH>::shape_deriv(const FEType fet,
371 : const Elem * elem,
372 : const unsigned int i,
373 : const unsigned int j,
374 : const Point & p,
375 : const bool add_p_level)
376 : {
377 0 : return FE<1,CLOUGH>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
378 : }
379 :
380 :
381 :
382 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
383 :
384 : template <>
385 0 : Real FE<1,CLOUGH>::shape_second_deriv(const Elem * elem,
386 : const Order order,
387 : const unsigned int i,
388 : const unsigned int j,
389 : const Point & p,
390 : const bool add_p_level)
391 : {
392 0 : libmesh_assert(elem);
393 :
394 0 : clough_compute_coefs(elem);
395 :
396 0 : const ElemType type = elem->type();
397 :
398 : const Order totalorder =
399 0 : order + add_p_level*elem->p_level();
400 :
401 0 : switch (totalorder)
402 : {
403 : // 3rd-order C1 cubic element
404 0 : case THIRD:
405 : {
406 0 : switch (type)
407 : {
408 : // C1 functions on the C1 cubic edge
409 0 : case EDGE2:
410 : case EDGE3:
411 : {
412 0 : switch (i)
413 : {
414 0 : case 0:
415 0 : return clough_raw_shape_second_deriv(0, j, p);
416 0 : case 1:
417 0 : return clough_raw_shape_second_deriv(1, j, p);
418 0 : case 2:
419 0 : return d1xd1x * clough_raw_shape_second_deriv(2, j, p);
420 0 : case 3:
421 0 : return d2xd2x * clough_raw_shape_second_deriv(3, j, p);
422 0 : default:
423 0 : libmesh_error_msg("Invalid shape function index i = " << i);
424 : }
425 : }
426 0 : default:
427 0 : libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
428 : }
429 : }
430 : // by default throw an error
431 0 : default:
432 0 : libmesh_error_msg("ERROR: Unsupported polynomial order = " << totalorder);
433 : }
434 : }
435 :
436 :
437 : template <>
438 0 : Real FE<1,CLOUGH>::shape_second_deriv(const ElemType,
439 : const Order,
440 : const unsigned int,
441 : const unsigned int,
442 : const Point &)
443 : {
444 0 : libmesh_error_msg("Clough-Tocher elements require the real element \nto construct gradient-based degrees of freedom.");
445 : return 0.;
446 : }
447 :
448 : template <>
449 0 : Real FE<1,CLOUGH>::shape_second_deriv(const FEType fet,
450 : const Elem * elem,
451 : const unsigned int i,
452 : const unsigned int j,
453 : const Point & p,
454 : const bool add_p_level)
455 : {
456 0 : return FE<1,CLOUGH>::shape_second_deriv(elem, fet.order, i, j, p, add_p_level);
457 : }
458 :
459 :
460 :
461 : #endif
462 :
463 : } // namespace libMesh
|
