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