libMesh
fe_nedelec_one_shape_3D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2024 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 RealGradient FE<3,NEDELEC_ONE>::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  libmesh_assert(elem);
36 
37  const Order totalorder = static_cast<Order>(order + add_p_level * elem->p_level());
38 
39  switch (totalorder)
40  {
41  // linear Nedelec (first kind) shape functions
42  case FIRST:
43  {
44  switch (elem->type())
45  {
46  case HEX20:
47  case HEX27:
48  {
49  libmesh_assert_less (i, 12);
50 
51  const Real xi = p(0);
52  const Real eta = p(1);
53  const Real zeta = p(2);
54 
55  // Even with a loose inverse_map tolerance we ought to
56  // be nearly on the element interior in master
57  // coordinates
58  libmesh_assert_less_equal ( std::fabs(xi), 1.0+10*TOLERANCE );
59  libmesh_assert_less_equal ( std::fabs(eta), 1.0+10*TOLERANCE );
60  libmesh_assert_less_equal ( std::fabs(zeta), 1.0+10*TOLERANCE );
61 
62  switch(i)
63  {
64  case 0:
65  {
66  if (elem->point(0) > elem->point(1))
67  return RealGradient( -0.125*(1.0-eta-zeta+eta*zeta), 0.0, 0.0 );
68  else
69  return RealGradient( 0.125*(1.0-eta-zeta+eta*zeta), 0.0, 0.0 );
70  }
71  case 1:
72  {
73  if (elem->point(1) > elem->point(2))
74  return RealGradient( 0.0, -0.125*(1.0+xi-zeta-xi*zeta), 0.0 );
75  else
76  return RealGradient( 0.0, 0.125*(1.0+xi-zeta-xi*zeta), 0.0 );
77  }
78  case 2:
79  {
80  if (elem->point(2) > elem->point(3))
81  return RealGradient( 0.125*(1.0+eta-zeta-eta*zeta), 0.0, 0.0 );
82  else
83  return RealGradient( -0.125*(1.0+eta-zeta-eta*zeta), 0.0, 0.0 );
84  }
85  case 3:
86  {
87  if (elem->point(3) > elem->point(0))
88  return RealGradient( 0.0, 0.125*(1.0-xi-zeta+xi*zeta), 0.0 );
89  else
90  return RealGradient( 0.0, -0.125*(1.0-xi-zeta+xi*zeta), 0.0 );
91  }
92  case 4:
93  {
94  if (elem->point(0) > elem->point(4))
95  return RealGradient( 0.0, 0.0, -0.125*(1.0-xi-eta+xi*eta) );
96  else
97  return RealGradient( 0.0, 0.0, 0.125*(1.0-xi-eta+xi*eta) );
98  }
99  case 5:
100  {
101  if (elem->point(1) > elem->point(5))
102  return RealGradient( 0.0, 0.0, -0.125*(1.0+xi-eta-xi*eta) );
103  else
104  return RealGradient( 0.0, 0.0, 0.125*(1.0+xi-eta-xi*eta) );
105  }
106  case 6:
107  {
108  if (elem->point(2) > elem->point(6))
109  return RealGradient( 0.0, 0.0, -0.125*(1.0+xi+eta+xi*eta) );
110  else
111  return RealGradient( 0.0, 0.0, 0.125*(1.0+xi+eta+xi*eta) );
112  }
113  case 7:
114  {
115  if (elem->point(3) > elem->point(7))
116  return RealGradient( 0.0, 0.0, -0.125*(1.0-xi+eta-xi*eta) );
117  else
118  return RealGradient( 0.0, 0.0, 0.125*(1.0-xi+eta-xi*eta) );
119  }
120  case 8:
121  {
122  if (elem->point(4) > elem->point(5))
123  return RealGradient( -0.125*(1.0-eta+zeta-eta*zeta), 0.0, 0.0 );
124  else
125  return RealGradient( 0.125*(1.0-eta+zeta-eta*zeta), 0.0, 0.0 );
126  }
127  case 9:
128  {
129  if (elem->point(5) > elem->point(6))
130  return RealGradient( 0.0, -0.125*(1.0+xi+zeta+xi*zeta), 0.0 );
131  else
132  return RealGradient( 0.0, 0.125*(1.0+xi+zeta+xi*zeta), 0.0 );
133  }
134  case 10:
135  {
136  if (elem->point(7) > elem->point(6))
137  return RealGradient( -0.125*(1.0+eta+zeta+eta*zeta), 0.0, 0.0 );
138  else
139  return RealGradient( 0.125*(1.0+eta+zeta+eta*zeta), 0.0, 0.0 );
140  }
141  case 11:
142  {
143  if (elem->point(4) > elem->point(7))
144  return RealGradient( 0.0, -0.125*(1.0-xi+zeta-xi*zeta), 0.0 );
145  else
146  return RealGradient( 0.0, 0.125*(1.0-xi+zeta-xi*zeta), 0.0 );
147  }
148 
149  default:
150  libmesh_error_msg("Invalid i = " << i);
151  }
152 
153  return RealGradient();
154  }
155 
156  case TET10:
157  case TET14:
158  {
159  libmesh_assert_less (i, 6);
160 
161  const Real xi = p(0);
162  const Real eta = p(1);
163  const Real zeta = p(2);
164 
165  switch(i)
166  {
167  case 0:
168  {
169  if (elem->point(0) > elem->point(1))
170  return RealGradient( -1.0+eta+zeta, -xi, -xi );
171  else
172  return RealGradient( 1.0-eta-zeta, xi, xi );
173  }
174  case 1:
175  {
176  if (elem->point(1) > elem->point(2))
177  return RealGradient( eta, -xi, 0.0 );
178  else
179  return RealGradient( -eta, xi, 0.0 );
180  }
181  case 2:
182  {
183  if (elem->point(0) > elem->point(2))
184  return RealGradient( -eta, -1.0+xi+zeta, -eta );
185  else
186  return RealGradient( eta, 1.0-xi-zeta, eta );
187  }
188  case 3:
189  {
190  if (elem->point(0) > elem->point(3))
191  return RealGradient( -zeta, -zeta, -1.0+xi+eta );
192  else
193  return RealGradient( zeta, zeta, 1.0-xi-eta );
194  }
195  case 4:
196  {
197  if (elem->point(1) > elem->point(3))
198  return RealGradient( zeta, 0.0, -xi );
199  else
200  return RealGradient( -zeta, 0.0, xi );
201  }
202  case 5:
203  {
204  if (elem->point(2) > elem->point(3))
205  return RealGradient( 0.0, zeta, -eta );
206  else
207  return RealGradient( 0.0, -zeta, eta );
208  }
209  default:
210  libmesh_error_msg("Invalid i = " << i);
211  }
212  return RealGradient();
213  }
214 
215  default:
216  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(elem->type()));
217  }
218  }
219 
220  // unsupported order
221  default:
222  libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << totalorder);
223  }
224 #else // LIBMESH_DIM != 3
225  libmesh_ignore(elem, order, i, p, add_p_level);
226  libmesh_not_implemented();
227 #endif
228 }
229 
230 
231 
232 template <>
233 RealGradient FE<3,NEDELEC_ONE>::shape(const ElemType,
234  const Order,
235  const unsigned int,
236  const Point &)
237 {
238  libmesh_error_msg("Nedelec elements require the element type \nbecause edge orientation is needed.");
239  return RealGradient();
240 }
241 
242 
243 
244 template <>
245 RealGradient FE<3,NEDELEC_ONE>::shape(const FEType fet,
246  const Elem * elem,
247  const unsigned int i,
248  const Point & p,
249  const bool add_p_level)
250 {
251  return FE<3,NEDELEC_ONE>::shape(elem, fet.order, i, p, add_p_level);
252 }
253 
254 
255 template <>
256 RealGradient FE<3,NEDELEC_ONE>::shape_deriv(const Elem * elem,
257  const Order order,
258  const unsigned int i,
259  const unsigned int j,
260  const Point & p,
261  const bool add_p_level)
262 {
263 #if LIBMESH_DIM == 3
264  libmesh_assert(elem);
265  libmesh_assert_less (j, 3);
266 
267  const Order totalorder = static_cast<Order>(order + add_p_level * elem->p_level());
268 
269  switch (totalorder)
270  {
271  // linear Nedelec (first kind) shape function first derivatives
272  case FIRST:
273  {
274  switch (elem->type())
275  {
276  case HEX20:
277  case HEX27:
278  {
279  libmesh_assert_less (i, 12);
280 
281  const Real xi = p(0);
282  const Real eta = p(1);
283  const Real zeta = p(2);
284 
285  // Even with a loose inverse_map tolerance we ought to
286  // be nearly on the element interior in master
287  // coordinates
288  libmesh_assert_less_equal ( std::fabs(xi), 1.0+TOLERANCE );
289  libmesh_assert_less_equal ( std::fabs(eta), 1.0+TOLERANCE );
290  libmesh_assert_less_equal ( std::fabs(zeta), 1.0+TOLERANCE );
291 
292  switch (j)
293  {
294  // d()/dxi
295  case 0:
296  {
297  switch(i)
298  {
299  case 0:
300  case 2:
301  case 8:
302  case 10:
303  return RealGradient();
304  case 1:
305  {
306  if (elem->point(1) > elem->point(2))
307  return RealGradient( 0.0, -0.125*(1.0-zeta) );
308  else
309  return RealGradient( 0.0, 0.125*(1.0-zeta) );
310  }
311  case 3:
312  {
313  if (elem->point(3) > elem->point(0))
314  return RealGradient( 0.0, 0.125*(-1.0+zeta) );
315  else
316  return RealGradient( 0.0, -0.125*(-1.0+zeta) );
317  }
318  case 4:
319  {
320  if (elem->point(0) > elem->point(4))
321  return RealGradient( 0.0, 0.0, -0.125*(-1.0+eta) );
322  else
323  return RealGradient( 0.0, 0.0, 0.125*(-1.0+eta) );
324  }
325  case 5:
326  {
327  if (elem->point(1) > elem->point(5))
328  return RealGradient( 0.0, 0.0, -0.125*(1.0-eta) );
329  else
330  return RealGradient( 0.0, 0.0, 0.125*(1.0-eta) );
331  }
332  case 6:
333  {
334  if (elem->point(2) > elem->point(6))
335  return RealGradient( 0.0, 0.0, -0.125*(1.0+eta) );
336  else
337  return RealGradient( 0.0, 0.0, 0.125*(1.0+eta) );
338  }
339  case 7:
340  {
341  if (elem->point(3) > elem->point(7))
342  return RealGradient( 0.0, 0.0, -0.125*(-1.0-eta) );
343  else
344  return RealGradient( 0.0, 0.0, 0.125*(-1.0-eta) );
345  }
346  case 9:
347  {
348  if (elem->point(5) > elem->point(6))
349  return RealGradient( 0.0, -0.125*(1.0+zeta), 0.0 );
350  else
351  return RealGradient( 0.0, 0.125*(1.0+zeta), 0.0 );
352  }
353  case 11:
354  {
355  if (elem->point(4) > elem->point(7))
356  return RealGradient( 0.0, -0.125*(-1.0-zeta), 0.0 );
357  else
358  return RealGradient( 0.0, 0.125*(-1.0-zeta), 0.0 );
359  }
360  default:
361  libmesh_error_msg("Invalid i = " << i);
362  } // switch(i)
363 
364  } // j = 0
365 
366  // d()/deta
367  case 1:
368  {
369  switch(i)
370  {
371  case 1:
372  case 3:
373  case 9:
374  case 11:
375  return RealGradient();
376  case 0:
377  {
378  if (elem->point(0) > elem->point(1))
379  return RealGradient( -0.125*(-1.0+zeta), 0.0, 0.0 );
380  else
381  return RealGradient( 0.125*(-1.0+zeta), 0.0, 0.0 );
382  }
383  case 2:
384  {
385  if (elem->point(2) > elem->point(3))
386  return RealGradient( 0.125*(1.0-zeta), 0.0, 0.0 );
387  else
388  return RealGradient( -0.125*(1.0-zeta), 0.0, 0.0 );
389  }
390  case 4:
391  {
392  if (elem->point(0) > elem->point(4))
393  return RealGradient( 0.0, 0.0, -0.125*(-1.0+xi) );
394  else
395  return RealGradient( 0.0, 0.0, 0.125*(-1.0+xi) );
396  }
397  case 5:
398  {
399  if (elem->point(1) > elem->point(5))
400  return RealGradient( 0.0, 0.0, -0.125*(-1.0-xi) );
401  else
402  return RealGradient( 0.0, 0.0, 0.125*(-1.0-xi) );
403  }
404  case 6:
405  {
406  if (elem->point(2) > elem->point(6))
407  return RealGradient( 0.0, 0.0, -0.125*(1.0+xi) );
408  else
409  return RealGradient( 0.0, 0.0, 0.125*(1.0+xi) );
410  }
411  case 7:
412  {
413  if (elem->point(3) > elem->point(7))
414  return RealGradient( 0.0, 0.0, -0.125*(1.0-xi) );
415  else
416  return RealGradient( 0.0, 0.0, 0.125*(1.0-xi) );
417  }
418  case 8:
419  {
420  if (elem->point(4) > elem->point(5))
421  return RealGradient( -0.125*(-1.0-zeta), 0.0, 0.0 );
422  else
423  return RealGradient( 0.125*(-1.0-zeta), 0.0, 0.0 );
424  }
425  case 10:
426  {
427  if (elem->point(7) > elem->point(6))
428  return RealGradient( -0.125*(1.0+zeta), 0.0, 0.0 );
429  else
430  return RealGradient( 0.125*(1.0+zeta), 0.0, 0.0 );
431  }
432  default:
433  libmesh_error_msg("Invalid i = " << i);
434  } // switch(i)
435 
436  } // j = 1
437 
438  // d()/dzeta
439  case 2:
440  {
441  switch(i)
442  {
443  case 4:
444  case 5:
445  case 6:
446  case 7:
447  return RealGradient();
448 
449  case 0:
450  {
451  if (elem->point(0) > elem->point(1))
452  return RealGradient( -0.125*(-1.0+eta), 0.0, 0.0 );
453  else
454  return RealGradient( 0.125*(-1.0+eta), 0.0, 0.0 );
455  }
456  case 1:
457  {
458  if (elem->point(1) > elem->point(2))
459  return RealGradient( 0.0, -0.125*(-1.0-xi), 0.0 );
460  else
461  return RealGradient( 0.0, 0.125*(-1.0-xi), 0.0 );
462  }
463  case 2:
464  {
465  if (elem->point(2) > elem->point(3))
466  return RealGradient( 0.125*(-1.0-eta), 0.0, 0.0 );
467  else
468  return RealGradient( -0.125*(-1.0-eta), 0.0, 0.0 );
469  }
470  case 3:
471  {
472  if (elem->point(3) > elem->point(0))
473  return RealGradient( 0.0, 0.125*(-1.0+xi), 0.0 );
474  else
475  return RealGradient( 0.0, -0.125*(-1.0+xi), 0.0 );
476  }
477  case 8:
478  {
479  if (elem->point(4) > elem->point(5))
480  return RealGradient( -0.125*(1.0-eta), 0.0, 0.0 );
481  else
482  return RealGradient( 0.125*(1.0-eta), 0.0, 0.0 );
483  }
484  case 9:
485  {
486  if (elem->point(5) > elem->point(6))
487  return RealGradient( 0.0, -0.125*(1.0+xi), 0.0 );
488  else
489  return RealGradient( 0.0, 0.125*(1.0+xi), 0.0 );
490  }
491  case 10:
492  {
493  if (elem->point(7) > elem->point(6))
494  return RealGradient( -0.125*(1.0+eta), 0.0, 0.0 );
495  else
496  return RealGradient( 0.125*(1.0+eta), 0.0, 0.0 );
497  }
498  case 11:
499  {
500  if (elem->point(4) > elem->point(7))
501  return RealGradient( 0.0, -0.125*(1.0-xi), 0.0 );
502  else
503  return RealGradient( 0.0, 0.125*(1.0-xi), 0.0 );
504  }
505  default:
506  libmesh_error_msg("Invalid i = " << i);
507  } // switch(i)
508 
509  } // j = 2
510 
511  default:
512  libmesh_error_msg("Invalid j = " << j);
513  }
514 
515  return RealGradient();
516  }
517 
518  case TET10:
519  case TET14:
520  {
521  libmesh_assert_less (i, 6);
522 
523 
524  switch (j)
525  {
526  // d()/dxi
527  case 0:
528  {
529  switch(i)
530  {
531  case 0:
532  {
533  if (elem->point(0) > elem->point(1))
534  return RealGradient( 0.0, -1.0, -1.0 );
535  else
536  return RealGradient( 0.0, 1.0, 1.0 );
537  }
538  case 1:
539  {
540  if (elem->point(1) > elem->point(2))
541  return RealGradient( 0.0, -1.0, 0.0 );
542  else
543  return RealGradient( 0.0, 1.0, 0.0 );
544  }
545  case 2:
546  {
547  if (elem->point(0) > elem->point(2))
548  return RealGradient( 0.0, 1.0, 0.0 );
549  else
550  return RealGradient( 0.0, -1.0, 0.0 );
551  }
552  case 3:
553  {
554  if (elem->point(0) > elem->point(3))
555  return RealGradient( 0.0, 0.0, 1.0 );
556  else
557  return RealGradient( 0.0, 0.0, -1.0 );
558  }
559  case 4:
560  {
561  if (elem->point(1) > elem->point(3))
562  return RealGradient( 0.0, 0.0, -1.0 );
563  else
564  return RealGradient( 0.0, 0.0, 1.0 );
565  }
566  case 5:
567  {
568  return RealGradient();
569  }
570 
571  default:
572  libmesh_error_msg("Invalid i = " << i);
573  } // switch(i)
574 
575  } // j = 0
576 
577  // d()/deta
578  case 1:
579  {
580  switch(i)
581  {
582  case 0:
583  {
584  if (elem->point(0) > elem->point(1))
585  return RealGradient( 1.0, 0.0, 0.0 );
586  else
587  return RealGradient( -1.0, 0.0, 0.0 );
588  }
589  case 1:
590  {
591  if (elem->point(1) > elem->point(2))
592  return RealGradient( 1.0, 0.0, 0.0 );
593  else
594  return RealGradient( -1.0, 0.0, 0.0 );
595  }
596  case 2:
597  {
598  if (elem->point(0) > elem->point(2))
599  return RealGradient( -1.0, 0.0, -1.0 );
600  else
601  return RealGradient( 1.0, 0.0, 1.0 );
602  }
603  case 3:
604  {
605  if (elem->point(0) > elem->point(3))
606  return RealGradient( 0.0, 0.0, 1.0 );
607  else
608  return RealGradient( 0.0, 0.0, -1.0 );
609  }
610  case 4:
611  {
612  return RealGradient();
613  }
614  case 5:
615  {
616  if (elem->point(2) > elem->point(3))
617  return RealGradient( 0.0, 0.0, -1.0 );
618  else
619  return RealGradient( 0.0, 0.0, 1.0 );
620  }
621  default:
622  libmesh_error_msg("Invalid i = " << i);
623  } // switch(i)
624 
625  } // j = 1
626 
627  // d()/dzeta
628  case 2:
629  {
630  switch(i)
631  {
632  case 0:
633  {
634  if (elem->point(0) > elem->point(1))
635  return RealGradient( 1.0, 0.0, 0.0 );
636  else
637  return RealGradient( -1.0, 0.0, 0.0 );
638  }
639  case 1:
640  {
641  return RealGradient();
642  }
643  case 2:
644  {
645  if (elem->point(0) > elem->point(2))
646  return RealGradient( 0.0, 1.0, 0.0 );
647  else
648  return RealGradient( 0.0, -1.0, 0.0 );
649  }
650  case 3:
651  {
652  if (elem->point(0) > elem->point(3))
653  return RealGradient( -1.0, -1.0, 0.0 );
654  else
655  return RealGradient( 1.0, 1.0, 0.0 );
656  }
657  case 4:
658  {
659  if (elem->point(1) > elem->point(3))
660  return RealGradient( 1.0, 0.0, 0.0 );
661  else
662  return RealGradient( -1.0, 0.0, 0.0 );
663  }
664  case 5:
665  {
666  if (elem->point(2) > elem->point(3))
667  return RealGradient( 0.0, 1.0, 0.0 );
668  else
669  return RealGradient( 0.0, -1.0, 0.0 );
670  }
671  default:
672  libmesh_error_msg("Invalid i = " << i);
673  } // switch(i)
674  } // j = 2
675 
676  default:
677  libmesh_error_msg("Invalid j = " << j);
678  }
679  return RealGradient();
680  }
681 
682  default:
683  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(elem->type()));
684  }
685  }
686  // unsupported order
687  default:
688  libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << totalorder);
689  }
690 
691 #else // LIBMESH_DIM != 3
692  libmesh_ignore(elem, order, i, j, p, add_p_level);
693  libmesh_not_implemented();
694 #endif
695 }
696 
697 
698 
699 
700 template <>
701 RealGradient FE<3,NEDELEC_ONE>::shape_deriv(const ElemType,
702  const Order,
703  const unsigned int,
704  const unsigned int,
705  const Point &)
706 {
707  libmesh_error_msg("Nedelec elements require the element type \nbecause edge orientation is needed.");
708  return RealGradient();
709 }
710 
711 
712 template <>
713 RealGradient FE<3,NEDELEC_ONE>::shape_deriv(const FEType fet,
714  const Elem * elem,
715  const unsigned int i,
716  const unsigned int j,
717  const Point & p,
718  const bool add_p_level)
719 {
720  return FE<3,NEDELEC_ONE>::shape_deriv(elem, fet.order, i, j, p, add_p_level);
721 }
722 
723 
724 
725 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
726 
727 template <>
728 RealGradient FE<3,NEDELEC_ONE>::shape_second_deriv(const Elem * elem,
729  const Order order,
730  const unsigned int i,
731  const unsigned int j,
732  const Point & libmesh_dbg_var(p),
733  const bool add_p_level)
734 {
735 #if LIBMESH_DIM == 3
736 
737  libmesh_assert(elem);
738 
739  // j = 0 ==> d^2 phi / dxi^2
740  // j = 1 ==> d^2 phi / dxi deta
741  // j = 2 ==> d^2 phi / deta^2
742  // j = 3 ==> d^2 phi / dxi dzeta
743  // j = 4 ==> d^2 phi / deta dzeta
744  // j = 5 ==> d^2 phi / dzeta^2
745  libmesh_assert_less (j, 6);
746 
747  const Order totalorder = static_cast<Order>(order + add_p_level * elem->p_level());
748 
749  switch (totalorder)
750  {
751  // linear Nedelec (first kind) shape function second derivatives
752  case FIRST:
753  {
754  switch (elem->type())
755  {
756  case HEX20:
757  case HEX27:
758  {
759  libmesh_assert_less (i, 12);
760 
761 #ifndef NDEBUG
762  const Real xi = p(0);
763  const Real eta = p(1);
764  const Real zeta = p(2);
765 #endif
766 
767  libmesh_assert_less_equal ( std::fabs(xi), 1.0+TOLERANCE );
768  libmesh_assert_less_equal ( std::fabs(eta), 1.0+TOLERANCE );
769  libmesh_assert_less_equal ( std::fabs(zeta), 1.0+TOLERANCE );
770 
771  switch (j)
772  {
773  // d^2()/dxi^2
774  case 0:
775  {
776  // All d^2()/dxi^2 derivatives for linear hexes are zero.
777  return RealGradient();
778  } // j = 0
779 
780  // d^2()/dxideta
781  case 1:
782  {
783  switch(i)
784  {
785  case 0:
786  case 1:
787  case 2:
788  case 3:
789  case 8:
790  case 9:
791  case 10:
792  case 11:
793  return RealGradient();
794  case 4:
795  {
796  if (elem->point(0) > elem->point(4))
797  return RealGradient( 0.0, 0.0, -0.125 );
798  else
799  return RealGradient( 0.0, 0.0, 0.125 );
800  }
801  case 5:
802  {
803  if (elem->point(1) > elem->point(5))
804  return RealGradient( 0.0, 0.0, 0.125 );
805  else
806  return RealGradient( 0.0, 0.0, -0.125 );
807  }
808  case 6:
809  {
810  if (elem->point(2) > elem->point(6))
811  return RealGradient( 0.0, 0.0, -0.125 );
812  else
813  return RealGradient( 0.0, 0.0, 0.125 );
814  }
815  case 7:
816  {
817  if (elem->point(3) > elem->point(7))
818  return RealGradient( 0.0, 0.0, 0.125 );
819  else
820  return RealGradient( 0.0, 0.0, -0.125 );
821  }
822  default:
823  libmesh_error_msg("Invalid i = " << i);
824  } // switch(i)
825 
826  } // j = 1
827 
828  // d^2()/deta^2
829  case 2:
830  {
831  // All d^2()/deta^2 derivatives for linear hexes are zero.
832  return RealGradient();
833  } // j = 2
834 
835  // d^2()/dxidzeta
836  case 3:
837  {
838  switch(i)
839  {
840  case 0:
841  case 2:
842  case 4:
843  case 5:
844  case 6:
845  case 7:
846  case 8:
847  case 10:
848  return RealGradient();
849 
850  case 1:
851  {
852  if (elem->point(1) > elem->point(2))
853  return RealGradient( 0.0, 0.125 );
854  else
855  return RealGradient( 0.0, -0.125 );
856  }
857  case 3:
858  {
859  if (elem->point(3) > elem->point(0))
860  return RealGradient( 0.0, -0.125 );
861  else
862  return RealGradient( 0.0, 0.125 );
863  }
864  case 9:
865  {
866  if (elem->point(5) > elem->point(6))
867  return RealGradient( 0.0, -0.125, 0.0 );
868  else
869  return RealGradient( 0.0, 0.125, 0.0 );
870  }
871  case 11:
872  {
873  if (elem->point(4) > elem->point(7))
874  return RealGradient( 0.0, 0.125, 0.0 );
875  else
876  return RealGradient( 0.0, -0.125, 0.0 );
877  }
878  default:
879  libmesh_error_msg("Invalid i = " << i);
880  } // switch(i)
881 
882  } // j = 3
883 
884  // d^2()/detadzeta
885  case 4:
886  {
887  switch(i)
888  {
889  case 1:
890  case 3:
891  case 4:
892  case 5:
893  case 6:
894  case 7:
895  case 9:
896  case 11:
897  return RealGradient();
898 
899  case 0:
900  {
901  if (elem->point(0) > elem->point(1))
902  return RealGradient( -0.125, 0.0, 0.0 );
903  else
904  return RealGradient( 0.125, 0.0, 0.0 );
905  }
906  case 2:
907  {
908  if (elem->point(2) > elem->point(3))
909  return RealGradient( 0.125, 0.0, 0.0 );
910  else
911  return RealGradient( -0.125, 0.0, 0.0 );
912  }
913  case 8:
914  {
915  if (elem->point(4) > elem->point(5))
916  return RealGradient( 0.125, 0.0, 0.0 );
917  else
918  return RealGradient( -0.125, 0.0, 0.0 );
919  }
920  case 10:
921  {
922  if (elem->point(7) > elem->point(6))
923  return RealGradient( -0.125, 0.0, 0.0 );
924  else
925  return RealGradient( 0.125, 0.0, 0.0 );
926  }
927  default:
928  libmesh_error_msg("Invalid i = " << i);
929  } // switch(i)
930 
931  } // j = 4
932 
933  // d^2()/dzeta^2
934  case 5:
935  {
936  // All d^2()/dzeta^2 derivatives for linear hexes are zero.
937  return RealGradient();
938  } // j = 5
939 
940  default:
941  libmesh_error_msg("Invalid j = " << j);
942  }
943 
944  return RealGradient();
945  }
946 
947  case TET10:
948  case TET14:
949  {
950  libmesh_assert_less (i, 6);
951  // All second derivatives for linear tets are zero.
952  return RealGradient();
953  }
954 
955  default:
956  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(elem->type()));
957 
958  } //switch(type)
959 
960  } // case FIRST:
961  // unsupported order
962  default:
963  libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << totalorder);
964  }
965 
966 #else // LIBMESH_DIM != 3
967  libmesh_assert(true || p(0));
968  libmesh_ignore(elem, order, i, j, add_p_level);
969  libmesh_not_implemented();
970 #endif
971 }
972 
973 
974 
975 template <>
976 RealGradient FE<3,NEDELEC_ONE>::shape_second_deriv(const ElemType,
977  const Order,
978  const unsigned int,
979  const unsigned int,
980  const Point &)
981 {
982  libmesh_error_msg("Nedelec elements require the element type \nbecause edge orientation is needed.");
983  return RealGradient();
984 }
985 
986 
987 
988 template <>
989 RealGradient FE<3,NEDELEC_ONE>::shape_second_deriv(const FEType fet,
990  const Elem * elem,
991  const unsigned int i,
992  const unsigned int j,
993  const Point & p,
994  const bool add_p_level)
995 {
996  return FE<3,NEDELEC_ONE>::shape_second_deriv(elem, fet.order, i, j, p, add_p_level);
997 }
998 
999 
1000 
1001 #endif
1002 
1003 } // namespace libMesh
libMesh::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:182
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
libMesh::RealGradient
RealVectorValue RealGradient
Definition: hp_coarsentest.h:49
libMesh::Order
Order
defines an enum for polynomial orders.
Definition: enum_order.h:40
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::TOLERANCE
static constexpr Real TOLERANCE
Definition: libmesh_common.h:150
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
libMesh::FE::shape_deriv
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::Elem::p_level
unsigned int p_level() const
Definition: elem.h:2945
libMesh::FEType::order
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:201
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition: libmesh_common.h:516
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:143
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition: elem.h:2277
libMesh::FE::shape_second_deriv
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
Generated on Mon Apr 15 2024 22:00:34 for libMesh by   doxygen 1.8.14