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 : // Local includes
20 : #include "libmesh/fe.h"
21 : #include "libmesh/elem.h"
22 : #include "libmesh/fe_interface.h"
23 : #include "libmesh/number_lookups.h"
24 : #include "libmesh/enum_to_string.h"
25 :
26 : namespace
27 : {
28 : using namespace libMesh;
29 :
30 : // Compute the static coefficients for an element
31 33598328 : void hermite_compute_coefs(const Elem * elem, std::vector<std::vector<Real>> & dxdxi
32 : #ifdef DEBUG
33 : , std::vector<Real> & dxdeta, std::vector<Real> & dydxi
34 : #endif
35 : )
36 : {
37 33598328 : const FEFamily mapping_family = FEMap::map_fe_type(*elem);
38 33598328 : const FEType map_fe_type(elem->default_order(), mapping_family);
39 :
40 : // Note: we explicitly don't consider the elem->p_level() when
41 : // computing the number of mapping shape functions.
42 : const int n_mapping_shape_functions =
43 33598328 : FEInterface::n_shape_functions (map_fe_type, /*extra_order=*/0, elem);
44 :
45 33598328 : static const Point dofpt[2] = {Point(-1,-1), Point(1,1)};
46 :
47 : FEInterface::shape_deriv_ptr shape_deriv_ptr =
48 33598328 : FEInterface::shape_deriv_function(map_fe_type, elem);
49 :
50 100794984 : for (int p = 0; p != 2; ++p)
51 : {
52 67196656 : dxdxi[0][p] = 0;
53 71278356 : dxdxi[1][p] = 0;
54 : #ifdef DEBUG
55 4081700 : dxdeta[p] = 0;
56 4081700 : dydxi[p] = 0;
57 : #endif
58 671966560 : for (int i = 0; i != n_mapping_shape_functions; ++i)
59 : {
60 : const Real ddxi = shape_deriv_ptr
61 604769904 : (map_fe_type, elem, i, 0, dofpt[p], /*add_p_level=*/false);
62 : const Real ddeta = shape_deriv_ptr
63 604769904 : (map_fe_type, elem, i, 1, dofpt[p], /*add_p_level=*/false);
64 :
65 641505204 : dxdxi[0][p] += elem->point(i)(0) * ddxi;
66 641505204 : dxdxi[1][p] += elem->point(i)(1) * ddeta;
67 : // dxdeta and dydxi should be 0!
68 : #ifdef DEBUG
69 36735300 : dxdeta[p] += elem->point(i)(0) * ddeta;
70 36735300 : dydxi[p] += elem->point(i)(1) * ddxi;
71 : #endif
72 : }
73 : // No singular elements!
74 4081700 : libmesh_assert(dxdxi[0][p]);
75 4081700 : libmesh_assert(dxdxi[1][p]);
76 : // No non-rectilinear or non-axis-aligned elements!
77 : #ifdef DEBUG
78 4081700 : libmesh_assert_less (std::abs(dxdeta[p]), 1e-9);
79 4081700 : libmesh_assert_less (std::abs(dydxi[p]), 1e-9);
80 : #endif
81 : }
82 33598328 : }
83 :
84 :
85 :
86 33598328 : Real hermite_bases_2D (std::vector<unsigned int> & bases1D,
87 : const std::vector<std::vector<Real>> & dxdxi,
88 : const Order & o,
89 : unsigned int i)
90 : {
91 2040850 : bases1D.clear();
92 33598328 : bases1D.resize(2,0);
93 2040850 : Real coef = 1.0;
94 :
95 33598328 : unsigned int e = o-3;
96 :
97 : // Nodes
98 33598328 : if (i < 16)
99 : {
100 33091808 : switch (i / 4)
101 : {
102 499660 : case 0:
103 499660 : break;
104 999320 : case 1:
105 8272952 : bases1D[0] = 1;
106 8272952 : break;
107 999320 : case 2:
108 8272952 : bases1D[0] = 1;
109 8272952 : bases1D[1] = 1;
110 8272952 : break;
111 999320 : case 3:
112 8272952 : bases1D[1] = 1;
113 8272952 : break;
114 0 : default:
115 0 : libmesh_error_msg("Invalid basis node = " << i/4);
116 : }
117 :
118 33091808 : unsigned int basisnum = i%4;
119 33091808 : switch (basisnum)
120 : {
121 499660 : case 0: // DoF = value at node
122 499660 : coef = 1.0;
123 499660 : break;
124 999320 : case 1: // DoF = x derivative at node
125 8272952 : coef = dxdxi[0][bases1D[0]];
126 8272952 : bases1D[0] += 2; break;
127 999320 : case 2: // DoF = y derivative at node
128 8272952 : coef = dxdxi[1][bases1D[1]];
129 8272952 : bases1D[1] += 2; break;
130 999320 : case 3: // DoF = xy derivative at node
131 8772612 : coef = dxdxi[0][bases1D[0]] * dxdxi[1][bases1D[1]];
132 8272952 : bases1D[0] += 2; bases1D[1] += 2; break;
133 0 : default:
134 0 : libmesh_error_msg("Invalid basisnum = " << basisnum);
135 : }
136 : }
137 : // Edges
138 506520 : else if (i < 16 + 4*2*e)
139 : {
140 450240 : unsigned int basisnum = (i - 16) % (2*e);
141 450240 : switch ((i - 16) / (2*e))
142 : {
143 112560 : case 0:
144 112560 : bases1D[0] = basisnum/2 + 4;
145 112560 : bases1D[1] = basisnum%2*2;
146 112560 : if (basisnum%2)
147 56280 : coef = dxdxi[1][0];
148 9380 : break;
149 112560 : case 1:
150 112560 : bases1D[0] = basisnum%2*2 + 1;
151 112560 : bases1D[1] = basisnum/2 + 4;
152 112560 : if (basisnum%2)
153 56280 : coef = dxdxi[0][1];
154 9380 : break;
155 112560 : case 2:
156 112560 : bases1D[0] = basisnum/2 + 4;
157 112560 : bases1D[1] = basisnum%2*2 + 1;
158 112560 : if (basisnum%2)
159 56280 : coef = dxdxi[1][1];
160 9380 : break;
161 112560 : case 3:
162 112560 : bases1D[0] = basisnum%2*2;
163 112560 : bases1D[1] = basisnum/2 + 4;
164 112560 : if (basisnum%2)
165 56280 : coef = dxdxi[0][0];
166 9380 : break;
167 0 : default:
168 0 : libmesh_error_msg("Invalid basisnum = " << basisnum);
169 : }
170 : }
171 : // Interior
172 : else
173 : {
174 56280 : unsigned int basisnum = i - 16 - 8*e;
175 56280 : bases1D[0] = square_number_row[basisnum]+4;
176 56280 : bases1D[1] = square_number_column[basisnum]+4;
177 : }
178 :
179 : // No singular elements
180 2040850 : libmesh_assert(coef);
181 33598328 : return coef;
182 : }
183 :
184 : } // end anonymous namespace
185 :
186 :
187 : namespace libMesh
188 : {
189 :
190 :
191 5979471 : LIBMESH_DEFAULT_VECTORIZED_FE(2,HERMITE)
192 :
193 :
194 : template <>
195 6660764 : Real FE<2,HERMITE>::shape(const Elem * elem,
196 : const Order order,
197 : const unsigned int i,
198 : const Point & p,
199 : const bool add_p_level)
200 : {
201 429109 : libmesh_assert(elem);
202 :
203 7948091 : std::vector<std::vector<Real>> dxdxi(2, std::vector<Real>(2, 0));
204 :
205 : #ifdef DEBUG
206 1287327 : std::vector<Real> dxdeta(2), dydxi(2);
207 : #endif
208 :
209 6660764 : hermite_compute_coefs(elem,dxdxi
210 : #ifdef DEBUG
211 : ,dxdeta,dydxi
212 : #endif
213 : );
214 :
215 6660764 : const ElemType type = elem->type();
216 :
217 : const Order totalorder =
218 7518982 : order + add_p_level*elem->p_level();
219 :
220 6231655 : switch (type)
221 : {
222 0 : case QUAD4:
223 : case QUADSHELL4:
224 0 : libmesh_assert_less (totalorder, 4);
225 : libmesh_fallthrough();
226 : case QUAD8:
227 : case QUADSHELL8:
228 : case QUAD9:
229 : case QUADSHELL9:
230 : {
231 429109 : libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
232 :
233 429109 : std::vector<unsigned int> bases1D;
234 :
235 6660764 : Real coef = hermite_bases_2D(bases1D, dxdxi, totalorder, i);
236 :
237 6660764 : return coef * FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
238 13321528 : FEHermite<1>::hermite_raw_shape(bases1D[1],p(1));
239 : }
240 0 : default:
241 0 : libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
242 : }
243 5802546 : }
244 :
245 :
246 :
247 : template <>
248 0 : Real FE<2,HERMITE>::shape(const ElemType,
249 : const Order,
250 : const unsigned int,
251 : const Point &)
252 : {
253 0 : libmesh_error_msg("Hermite elements require the real element \nto construct gradient-based degrees of freedom.");
254 : return 0.;
255 : }
256 :
257 :
258 : template <>
259 0 : Real FE<2,HERMITE>::shape(const FEType fet,
260 : const Elem * elem,
261 : const unsigned int i,
262 : const Point & p,
263 : const bool add_p_level)
264 : {
265 0 : return FE<2,HERMITE>::shape(elem, fet.order, i, p, add_p_level);
266 : }
267 :
268 :
269 : template <>
270 13070880 : Real FE<2,HERMITE>::shape_deriv(const Elem * elem,
271 : const Order order,
272 : const unsigned int i,
273 : const unsigned int j,
274 : const Point & p,
275 : const bool add_p_level)
276 : {
277 835176 : libmesh_assert(elem);
278 835176 : libmesh_assert (j == 0 || j == 1);
279 :
280 15576408 : std::vector<std::vector<Real>> dxdxi(2, std::vector<Real>(2, 0));
281 :
282 : #ifdef DEBUG
283 2505528 : std::vector<Real> dxdeta(2), dydxi(2);
284 : #endif
285 :
286 13070880 : hermite_compute_coefs(elem,dxdxi
287 : #ifdef DEBUG
288 : ,dxdeta,dydxi
289 : #endif
290 : );
291 :
292 13070880 : const ElemType type = elem->type();
293 :
294 : const Order totalorder =
295 14741232 : order + add_p_level*elem->p_level();
296 :
297 12235704 : switch (type)
298 : {
299 0 : case QUAD4:
300 : case QUADSHELL4:
301 0 : libmesh_assert_less (totalorder, 4);
302 : libmesh_fallthrough();
303 : case QUAD8:
304 : case QUADSHELL8:
305 : case QUAD9:
306 : case QUADSHELL9:
307 : {
308 835176 : libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
309 :
310 1670352 : std::vector<unsigned int> bases1D;
311 :
312 13070880 : Real coef = hermite_bases_2D(bases1D, dxdxi, totalorder, i);
313 :
314 13070880 : switch (j)
315 : {
316 6535440 : case 0:
317 6535440 : return coef *
318 6535440 : FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) *
319 6535440 : FEHermite<1>::hermite_raw_shape(bases1D[1],p(1));
320 6535440 : case 1:
321 6535440 : return coef *
322 6535440 : FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
323 6535440 : FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1));
324 0 : default:
325 0 : libmesh_error_msg("Invalid derivative index j = " << j);
326 : }
327 : }
328 0 : default:
329 0 : libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
330 : }
331 11400528 : }
332 :
333 :
334 : template <>
335 0 : Real FE<2,HERMITE>::shape_deriv(const ElemType,
336 : const Order,
337 : const unsigned int,
338 : const unsigned int,
339 : const Point &)
340 : {
341 0 : libmesh_error_msg("Hermite elements require the real element \nto construct gradient-based degrees of freedom.");
342 : return 0.;
343 : }
344 :
345 :
346 : template <>
347 0 : Real FE<2,HERMITE>::shape_deriv(const FEType fet,
348 : const Elem * elem,
349 : const unsigned int i,
350 : const unsigned int j,
351 : const Point & p,
352 : const bool add_p_level)
353 : {
354 0 : return FE<2,HERMITE>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
355 : }
356 :
357 :
358 :
359 :
360 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
361 :
362 :
363 : template <>
364 13866684 : Real FE<2,HERMITE>::shape_second_deriv(const Elem * elem,
365 : const Order order,
366 : const unsigned int i,
367 : const unsigned int j,
368 : const Point & p,
369 : const bool add_p_level)
370 : {
371 776565 : libmesh_assert(elem);
372 776565 : libmesh_assert (j == 0 || j == 1 || j == 2);
373 :
374 16196379 : std::vector<std::vector<Real>> dxdxi(2, std::vector<Real>(2, 0));
375 :
376 : #ifdef DEBUG
377 2329695 : std::vector<Real> dxdeta(2), dydxi(2);
378 : #endif
379 :
380 13866684 : hermite_compute_coefs(elem,dxdxi
381 : #ifdef DEBUG
382 : ,dxdeta,dydxi
383 : #endif
384 : );
385 :
386 13866684 : const ElemType type = elem->type();
387 :
388 : const Order totalorder =
389 15419814 : order + add_p_level*elem->p_level();
390 :
391 13090119 : switch (type)
392 : {
393 0 : case QUAD4:
394 : case QUADSHELL4:
395 0 : libmesh_assert_less (totalorder, 4);
396 : libmesh_fallthrough();
397 : case QUAD8:
398 : case QUADSHELL8:
399 : case QUAD9:
400 : case QUADSHELL9:
401 : {
402 776565 : libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
403 :
404 1553130 : std::vector<unsigned int> bases1D;
405 :
406 13866684 : Real coef = hermite_bases_2D(bases1D, dxdxi, totalorder, i);
407 :
408 13866684 : switch (j)
409 : {
410 4622228 : case 0:
411 4622228 : return coef *
412 4622228 : FEHermite<1>::hermite_raw_shape_second_deriv(bases1D[0],p(0)) *
413 4622228 : FEHermite<1>::hermite_raw_shape(bases1D[1],p(1));
414 4622228 : case 1:
415 4622228 : return coef *
416 4622228 : FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) *
417 4622228 : FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1));
418 4622228 : case 2:
419 4622228 : return coef *
420 4622228 : FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
421 4622228 : FEHermite<1>::hermite_raw_shape_second_deriv(bases1D[1],p(1));
422 0 : default:
423 0 : libmesh_error_msg("Invalid derivative index j = " << j);
424 : }
425 : }
426 0 : default:
427 0 : libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
428 : }
429 12313554 : }
430 :
431 :
432 : template <>
433 0 : Real FE<2,HERMITE>::shape_second_deriv(const ElemType,
434 : const Order,
435 : const unsigned int,
436 : const unsigned int,
437 : const Point &)
438 : {
439 0 : libmesh_error_msg("Hermite elements require the real element \nto construct gradient-based degrees of freedom.");
440 : return 0.;
441 : }
442 :
443 :
444 : template <>
445 0 : Real FE<2,HERMITE>::shape_second_deriv(const FEType fet,
446 : const Elem * elem,
447 : const unsigned int i,
448 : const unsigned int j,
449 : const Point & p,
450 : const bool add_p_level)
451 : {
452 0 : return FE<2,HERMITE>::shape_second_deriv(elem, fet.order, i, j, p, add_p_level);
453 : }
454 :
455 :
456 : #endif
457 :
458 : } // namespace libMesh
|
