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/enum_to_string.h"
23 :
24 : namespace libMesh
25 : {
26 :
27 : template <>
28 82028382 : RealGradient FE<3,RAVIART_THOMAS>::shape(const Elem * elem,
29 : const Order order,
30 : const unsigned int i,
31 : const Point & p,
32 : const bool add_p_level)
33 : {
34 : #if LIBMESH_DIM == 3
35 6836142 : libmesh_assert(elem);
36 :
37 88864524 : const Order totalorder = order + add_p_level*elem->p_level();
38 6836142 : libmesh_assert_less(i, n_dofs(elem->type(), totalorder));
39 :
40 82028382 : const char sign = elem->positive_face_orientation(i) ? 1 : -1;
41 :
42 82028382 : const Real xi = p(0);
43 82028382 : const Real eta = p(1);
44 82028382 : const Real zeta = p(2);
45 :
46 82028382 : switch (totalorder)
47 : {
48 : // linear Raviart-Thomas shape functions
49 82028382 : case FIRST:
50 : {
51 82028382 : switch (elem->type())
52 : {
53 46575342 : case HEX27:
54 : {
55 : // Even with a loose inverse_map tolerance we ought to
56 : // be nearly on the element interior in master
57 : // coordinates
58 3881646 : libmesh_assert_less_equal ( std::fabs(xi), 1.0+10*TOLERANCE );
59 3881646 : libmesh_assert_less_equal ( std::fabs(eta), 1.0+10*TOLERANCE );
60 3881646 : libmesh_assert_less_equal ( std::fabs(zeta), 1.0+10*TOLERANCE );
61 :
62 46575342 : switch(i)
63 : {
64 7762557 : case 0:
65 7762557 : return sign * RealGradient( 0.0, 0.0, 0.125*(zeta-1.0) );
66 7762557 : case 1:
67 7762557 : return sign * RealGradient( 0.0, 0.125*(eta-1.0), 0.0 );
68 7762557 : case 2:
69 7762557 : return sign * RealGradient( 0.125*(xi+1.0), 0.0, 0.0 );
70 7762557 : case 3:
71 7762557 : return sign * RealGradient( 0.0, 0.125*(1.0+eta), 0.0 );
72 7762557 : case 4:
73 7762557 : return sign * RealGradient( 0.125*(xi-1.0), 0.0, 0.0 );
74 7762557 : case 5:
75 7762557 : return sign * RealGradient( 0.0, 0.0, 0.125*(1.0+zeta) );
76 :
77 0 : default:
78 0 : libmesh_error_msg("Invalid i = " << i);
79 : }
80 : }
81 :
82 35453040 : case TET14:
83 : {
84 35453040 : switch(i)
85 : {
86 8863260 : case 0:
87 8863260 : return sign * RealGradient( 2.0*xi, 2.0*eta, 2.0*zeta-2.0 );
88 8863260 : case 1:
89 8863260 : return sign * RealGradient( 2.0*xi, 2.0*eta-2.0, 2.0*zeta );
90 8863260 : case 2:
91 8863260 : return sign * RealGradient( 2.0*xi, 2.0*eta, 2.0*zeta );
92 8863260 : case 3:
93 8863260 : return sign * RealGradient( 2.0*xi-2.0, 2.0*eta, 2.0*zeta );
94 :
95 0 : default:
96 0 : libmesh_error_msg("Invalid i = " << i);
97 : }
98 : }
99 :
100 0 : default:
101 0 : libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(elem->type()));
102 : }
103 : }
104 :
105 : // unsupported order
106 0 : default:
107 0 : libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << totalorder);
108 : }
109 :
110 : #else // LIBMESH_DIM != 3
111 : libmesh_ignore(elem, order, i, p, add_p_level);
112 : libmesh_not_implemented();
113 : #endif
114 : }
115 :
116 :
117 : template <>
118 10363464 : RealGradient FE<3,L2_RAVIART_THOMAS>::shape(const Elem * elem,
119 : const Order order,
120 : const unsigned int i,
121 : const Point & p,
122 : const bool add_p_level)
123 : {
124 10363464 : return FE<3,RAVIART_THOMAS>::shape(elem, order, i, p, add_p_level);
125 : }
126 :
127 :
128 : template <>
129 0 : RealGradient FE<3,RAVIART_THOMAS>::shape(const ElemType,
130 : const Order,
131 : const unsigned int,
132 : const Point &)
133 : {
134 0 : libmesh_error_msg("Raviart-Thomas elements require the element type \nbecause face orientation is needed.");
135 : return RealGradient();
136 : }
137 :
138 :
139 : template <>
140 0 : RealGradient FE<3,L2_RAVIART_THOMAS>::shape(const ElemType,
141 : const Order,
142 : const unsigned int,
143 : const Point &)
144 : {
145 0 : libmesh_error_msg("Raviart-Thomas elements require the element type \nbecause face orientation is needed.");
146 : return RealGradient();
147 : }
148 :
149 :
150 : template <>
151 0 : RealGradient FE<3,RAVIART_THOMAS>::shape(const FEType fet,
152 : const Elem * elem,
153 : const unsigned int i,
154 : const Point & p,
155 : const bool add_p_level)
156 : {
157 0 : return FE<3,RAVIART_THOMAS>::shape(elem, fet.order, i, p, add_p_level);
158 : }
159 :
160 :
161 : template <>
162 0 : RealGradient FE<3,L2_RAVIART_THOMAS>::shape(const FEType fet,
163 : const Elem * elem,
164 : const unsigned int i,
165 : const Point & p,
166 : const bool add_p_level)
167 : {
168 0 : return FE<3,L2_RAVIART_THOMAS>::shape(elem, fet.order, i, p, add_p_level);
169 : }
170 :
171 :
172 : template <>
173 126464562 : RealGradient FE<3,RAVIART_THOMAS>::shape_deriv(const Elem * elem,
174 : const Order order,
175 : const unsigned int i,
176 : const unsigned int j,
177 : const Point &,
178 : const bool add_p_level)
179 : {
180 : #if LIBMESH_DIM == 3
181 10534860 : libmesh_assert(elem);
182 10534860 : libmesh_assert_less (j, 3);
183 :
184 136999422 : const Order totalorder = order + add_p_level*elem->p_level();
185 10534860 : libmesh_assert_less(i, n_dofs(elem->type(), totalorder));
186 :
187 126464562 : const char sign = elem->positive_face_orientation(i) ? 1 : -1;
188 :
189 126464562 : switch (totalorder)
190 : {
191 : // linear Raviart-Thomas shape function first derivatives
192 126464562 : case FIRST:
193 : {
194 126464562 : switch (elem->type())
195 : {
196 74627298 : case HEX27:
197 : {
198 74627298 : switch (j)
199 : {
200 : // d()/dxi
201 24875766 : case 0:
202 : {
203 22804146 : switch(i)
204 : {
205 1381080 : case 0:
206 : case 1:
207 : case 3:
208 : case 5:
209 1381080 : return RealGradient();
210 8291922 : case 2:
211 : case 4:
212 690540 : return sign * RealGradient( 0.125, 0.0, 0.0 );
213 :
214 0 : default:
215 0 : libmesh_error_msg("Invalid i = " << i);
216 : } // switch(i)
217 :
218 : } // j = 0
219 :
220 : // d()/deta
221 24875766 : case 1:
222 : {
223 22804146 : switch(i)
224 : {
225 1381080 : case 0:
226 : case 2:
227 : case 4:
228 : case 5:
229 1381080 : return RealGradient();
230 8291922 : case 1:
231 : case 3:
232 690540 : return sign * RealGradient( 0.0, 0.125, 0.0 );
233 :
234 0 : default:
235 0 : libmesh_error_msg("Invalid i = " << i);
236 : } // switch(i)
237 :
238 : } // j = 1
239 :
240 : // d()/dzeta
241 24875766 : case 2:
242 : {
243 24875766 : switch(i)
244 : {
245 1381080 : case 1:
246 : case 2:
247 : case 3:
248 : case 4:
249 1381080 : return RealGradient();
250 8291922 : case 0:
251 : case 5:
252 690540 : return sign * RealGradient( 0.0, 0.0, 0.125 );
253 :
254 0 : default:
255 0 : libmesh_error_msg("Invalid i = " << i);
256 : } // switch(i)
257 :
258 : } // j = 2
259 :
260 0 : default:
261 0 : libmesh_error_msg("Invalid j = " << j);
262 : }
263 : }
264 :
265 51837264 : case TET14:
266 : {
267 51837264 : switch (j)
268 : {
269 : // d()/dxi
270 17279088 : case 0:
271 : {
272 17279088 : switch(i)
273 : {
274 17279088 : case 0:
275 : case 1:
276 : case 2:
277 : case 3:
278 1440000 : return sign * RealGradient( 2.0, 0.0, 0.0 );
279 :
280 0 : default:
281 0 : libmesh_error_msg("Invalid i = " << i);
282 : } // switch(i)
283 :
284 : } // j = 0
285 :
286 : // d()/deta
287 17279088 : case 1:
288 : {
289 17279088 : switch(i)
290 : {
291 17279088 : case 0:
292 : case 1:
293 : case 2:
294 : case 3:
295 1440000 : return sign * RealGradient( 0.0, 2.0, 0.0 );
296 :
297 0 : default:
298 0 : libmesh_error_msg("Invalid i = " << i);
299 : } // switch(i)
300 :
301 : } // j = 1
302 :
303 : // d()/dzeta
304 17279088 : case 2:
305 : {
306 17279088 : switch(i)
307 : {
308 17279088 : case 0:
309 : case 1:
310 : case 2:
311 : case 3:
312 1440000 : return sign * RealGradient( 0.0, 0.0, 2.0 );
313 :
314 0 : default:
315 0 : libmesh_error_msg("Invalid i = " << i);
316 : } // switch(i)
317 :
318 : } // j = 2
319 :
320 0 : default:
321 0 : libmesh_error_msg("Invalid j = " << j);
322 : }
323 : }
324 :
325 0 : default:
326 0 : libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(elem->type()));
327 : }
328 : }
329 : // unsupported order
330 0 : default:
331 0 : libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << totalorder);
332 : }
333 :
334 : #else // LIBMESH_DIM != 3
335 : libmesh_ignore(elem, order, i, j, p, add_p_level);
336 : libmesh_not_implemented();
337 : #endif
338 : }
339 :
340 :
341 : template <>
342 10333656 : RealGradient FE<3,L2_RAVIART_THOMAS>::shape_deriv(const Elem * elem,
343 : const Order order,
344 : const unsigned int i,
345 : const unsigned int j,
346 : const Point & p,
347 : const bool add_p_level)
348 : {
349 10333656 : return FE<3,RAVIART_THOMAS>::shape_deriv(elem, order, i, j, p, add_p_level);
350 : }
351 :
352 :
353 : template <>
354 0 : RealGradient FE<3,RAVIART_THOMAS>::shape_deriv(const ElemType,
355 : const Order,
356 : const unsigned int,
357 : const unsigned int,
358 : const Point &)
359 : {
360 0 : libmesh_error_msg("Raviart-Thomas elements require the element type \nbecause face orientation is needed.");
361 : return RealGradient();
362 : }
363 :
364 :
365 : template <>
366 0 : RealGradient FE<3,L2_RAVIART_THOMAS>::shape_deriv(const ElemType,
367 : const Order,
368 : const unsigned int,
369 : const unsigned int,
370 : const Point &)
371 : {
372 0 : libmesh_error_msg("Raviart-Thomas elements require the element type \nbecause face orientation is needed.");
373 : return RealGradient();
374 : }
375 :
376 :
377 : template <>
378 0 : RealGradient FE<3,RAVIART_THOMAS>::shape_deriv(const FEType fet,
379 : const Elem * elem,
380 : const unsigned int i,
381 : const unsigned int j,
382 : const Point & p,
383 : const bool add_p_level)
384 : {
385 0 : return FE<3,RAVIART_THOMAS>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
386 : }
387 :
388 :
389 : template <>
390 0 : RealGradient FE<3,L2_RAVIART_THOMAS>::shape_deriv(const FEType fet,
391 : const Elem * elem,
392 : const unsigned int i,
393 : const unsigned int j,
394 : const Point & p,
395 : const bool add_p_level)
396 : {
397 0 : return FE<3,L2_RAVIART_THOMAS>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
398 : }
399 :
400 :
401 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
402 :
403 : template <>
404 0 : RealGradient FE<3,RAVIART_THOMAS>::shape_second_deriv(const Elem * elem,
405 : const Order order,
406 : const unsigned int libmesh_dbg_var(i),
407 : const unsigned int libmesh_dbg_var(j),
408 : const Point &,
409 : const bool add_p_level)
410 : {
411 : #if LIBMESH_DIM == 3
412 :
413 0 : libmesh_assert(elem);
414 :
415 : // j = 0 ==> d^2 phi / dxi^2
416 : // j = 1 ==> d^2 phi / dxi deta
417 : // j = 2 ==> d^2 phi / deta^2
418 : // j = 3 ==> d^2 phi / dxi dzeta
419 : // j = 4 ==> d^2 phi / deta dzeta
420 : // j = 5 ==> d^2 phi / dzeta^2
421 0 : libmesh_assert_less (j, 6);
422 :
423 0 : const Order totalorder = order + add_p_level*elem->p_level();
424 0 : libmesh_assert_less(i, n_dofs(elem->type(), totalorder));
425 :
426 0 : switch (totalorder)
427 : {
428 : // linear Raviart-Thomas shape function second derivatives
429 0 : case FIRST:
430 : {
431 0 : switch (elem->type())
432 : {
433 : // All second derivatives for linear hexes and tets are zero.
434 0 : case HEX27:
435 : case TET14:
436 0 : return RealGradient();
437 :
438 0 : default:
439 0 : libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(elem->type()));
440 :
441 : } //switch(type)
442 : }
443 :
444 : // unsupported order
445 0 : default:
446 0 : libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << totalorder);
447 : }
448 :
449 : #else // LIBMESH_DIM != 3
450 : libmesh_assert(true || p(0));
451 : libmesh_ignore(elem, order, i, j, add_p_level);
452 : libmesh_not_implemented();
453 : #endif
454 : }
455 :
456 :
457 : template <>
458 0 : RealGradient FE<3,L2_RAVIART_THOMAS>::shape_second_deriv(const Elem * elem,
459 : const Order order,
460 : const unsigned int i,
461 : const unsigned int j,
462 : const Point & p,
463 : const bool add_p_level)
464 : {
465 0 : return FE<3,RAVIART_THOMAS>::shape_second_deriv(elem, order, i, j, p, add_p_level);
466 : }
467 :
468 :
469 : template <>
470 0 : RealGradient FE<3,RAVIART_THOMAS>::shape_second_deriv(const ElemType,
471 : const Order,
472 : const unsigned int,
473 : const unsigned int,
474 : const Point &)
475 : {
476 0 : libmesh_error_msg("Raviart-Thomas elements require the element type \nbecause face orientation is needed.");
477 : return RealGradient();
478 : }
479 :
480 :
481 : template <>
482 0 : RealGradient FE<3,L2_RAVIART_THOMAS>::shape_second_deriv(const ElemType,
483 : const Order,
484 : const unsigned int,
485 : const unsigned int,
486 : const Point &)
487 : {
488 0 : libmesh_error_msg("Raviart-Thomas elements require the element type \nbecause face orientation is needed.");
489 : return RealGradient();
490 : }
491 :
492 :
493 : template <>
494 0 : RealGradient FE<3,RAVIART_THOMAS>::shape_second_deriv(const FEType fet,
495 : const Elem * elem,
496 : const unsigned int i,
497 : const unsigned int j,
498 : const Point & p,
499 : const bool add_p_level)
500 : {
501 0 : return FE<3,RAVIART_THOMAS>::shape_second_deriv(elem, fet.order, i, j, p, add_p_level);
502 : }
503 :
504 :
505 : template <>
506 0 : RealGradient FE<3,L2_RAVIART_THOMAS>::shape_second_deriv(const FEType fet,
507 : const Elem * elem,
508 : const unsigned int i,
509 : const unsigned int j,
510 : const Point & p,
511 : const bool add_p_level)
512 : {
513 0 : return FE<3,L2_RAVIART_THOMAS>::shape_second_deriv(elem, fet.order, i, j, p, add_p_level);
514 : }
515 :
516 : #endif
517 :
518 : } // namespace libMesh
|
