libMesh
fe_nedelec_one_shape_3D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 
25 namespace libMesh
26 {
27 
28 template <>
30  const Order,
31  const unsigned int,
32  const Point &)
33 {
34  libmesh_error_msg("Nedelec elements require the element type \nbecause edge orientation is needed.");
35  return RealGradient();
36 }
37 
38 
39 
40 template <>
42  const Order order,
43  const unsigned int i,
44  const Point & p)
45 {
46 #if LIBMESH_DIM == 3
47  libmesh_assert(elem);
48 
49  const Order totalorder = static_cast<Order>(order + elem->p_level());
50 
51  switch (totalorder)
52  {
53  // linear Lagrange shape functions
54  case FIRST:
55  {
56  switch (elem->type())
57  {
58  case HEX20:
59  case HEX27:
60  {
61  libmesh_assert_less (i, 12);
62 
63  const Real xi = p(0);
64  const Real eta = p(1);
65  const Real zeta = p(2);
66 
67  // Even with a loose inverse_map tolerance we ought to
68  // be nearly on the element interior in master
69  // coordinates
70  libmesh_assert_less_equal ( std::fabs(xi), 1.0+10*TOLERANCE );
71  libmesh_assert_less_equal ( std::fabs(eta), 1.0+10*TOLERANCE );
72  libmesh_assert_less_equal ( std::fabs(zeta), 1.0+10*TOLERANCE );
73 
74  switch(i)
75  {
76  case 0:
77  {
78  if (elem->point(0) > elem->point(1))
79  return RealGradient( -0.125*(1.0-eta-zeta+eta*zeta), 0.0, 0.0 );
80  else
81  return RealGradient( 0.125*(1.0-eta-zeta+eta*zeta), 0.0, 0.0 );
82  }
83  case 1:
84  {
85  if (elem->point(1) > elem->point(2))
86  return RealGradient( 0.0, -0.125*(1.0+xi-zeta-xi*zeta), 0.0 );
87  else
88  return RealGradient( 0.0, 0.125*(1.0+xi-zeta-xi*zeta), 0.0 );
89  }
90  case 2:
91  {
92  if (elem->point(2) > elem->point(3))
93  return RealGradient( 0.125*(1.0+eta-zeta-eta*zeta), 0.0, 0.0 );
94  else
95  return RealGradient( -0.125*(1.0+eta-zeta-eta*zeta), 0.0, 0.0 );
96  }
97  case 3:
98  {
99  if (elem->point(3) > elem->point(0))
100  return RealGradient( 0.0, 0.125*(1.0-xi-zeta+xi*zeta), 0.0 );
101  else
102  return RealGradient( 0.0, -0.125*(1.0-xi-zeta+xi*zeta), 0.0 );
103  }
104  case 4:
105  {
106  if (elem->point(0) > elem->point(4))
107  return RealGradient( 0.0, 0.0, -0.125*(1.0-xi-eta+xi*eta) );
108  else
109  return RealGradient( 0.0, 0.0, 0.125*(1.0-xi-eta+xi*eta) );
110  }
111  case 5:
112  {
113  if (elem->point(1) > elem->point(5))
114  return RealGradient( 0.0, 0.0, -0.125*(1.0+xi-eta-xi*eta) );
115  else
116  return RealGradient( 0.0, 0.0, 0.125*(1.0+xi-eta-xi*eta) );
117  }
118  case 6:
119  {
120  if (elem->point(2) > elem->point(6))
121  return RealGradient( 0.0, 0.0, -0.125*(1.0+xi+eta+xi*eta) );
122  else
123  return RealGradient( 0.0, 0.0, 0.125*(1.0+xi+eta+xi*eta) );
124  }
125  case 7:
126  {
127  if (elem->point(3) > elem->point(7))
128  return RealGradient( 0.0, 0.0, -0.125*(1.0-xi+eta-xi*eta) );
129  else
130  return RealGradient( 0.0, 0.0, 0.125*(1.0-xi+eta-xi*eta) );
131  }
132  case 8:
133  {
134  if (elem->point(4) > elem->point(5))
135  return RealGradient( -0.125*(1.0-eta+zeta-eta*zeta), 0.0, 0.0 );
136  else
137  return RealGradient( 0.125*(1.0-eta+zeta-eta*zeta), 0.0, 0.0 );
138  }
139  case 9:
140  {
141  if (elem->point(5) > elem->point(6))
142  return RealGradient( 0.0, -0.125*(1.0+xi+zeta+xi*zeta), 0.0 );
143  else
144  return RealGradient( 0.0, 0.125*(1.0+xi+zeta+xi*zeta), 0.0 );
145  }
146  case 10:
147  {
148  if (elem->point(7) > elem->point(6))
149  return RealGradient( -0.125*(1.0+eta+zeta+eta*zeta), 0.0, 0.0 );
150  else
151  return RealGradient( 0.125*(1.0+eta+zeta+eta*zeta), 0.0, 0.0 );
152  }
153  case 11:
154  {
155  if (elem->point(4) > elem->point(7))
156  return RealGradient( 0.0, -0.125*(1.0-xi+zeta-xi*zeta), 0.0 );
157  else
158  return RealGradient( 0.0, 0.125*(1.0-xi+zeta-xi*zeta), 0.0 );
159  }
160 
161  default:
162  libmesh_error_msg("Invalid i = " << i);
163  }
164 
165  return RealGradient();
166  }
167 
168  case TET10:
169  {
170  libmesh_assert_less (i, 6);
171 
172  libmesh_not_implemented();
173  return RealGradient();
174  }
175 
176  default:
177  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << elem->type());
178  }
179  }
180 
181  // unsupported order
182  default:
183  libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << totalorder);
184  }
185 #endif
186 
187  libmesh_error_msg("We'll never get here!");
188  return RealGradient();
189 }
190 
191 
192 
193 
194 template <>
196  const Order,
197  const unsigned int,
198  const unsigned int,
199  const Point &)
200 {
201  libmesh_error_msg("Nedelec elements require the element type \nbecause edge orientation is needed.");
202  return RealGradient();
203 }
204 
205 template <>
207  const Order order,
208  const unsigned int i,
209  const unsigned int j,
210  const Point & p)
211 {
212 #if LIBMESH_DIM == 3
213  libmesh_assert(elem);
214  libmesh_assert_less (j, 3);
215 
216  const Order totalorder = static_cast<Order>(order + elem->p_level());
217 
218  switch (totalorder)
219  {
220  case FIRST:
221  {
222  switch (elem->type())
223  {
224  case HEX20:
225  case HEX27:
226  {
227  libmesh_assert_less (i, 12);
228 
229  const Real xi = p(0);
230  const Real eta = p(1);
231  const Real zeta = p(2);
232 
233  // Even with a loose inverse_map tolerance we ought to
234  // be nearly on the element interior in master
235  // coordinates
236  libmesh_assert_less_equal ( std::fabs(xi), 1.0+TOLERANCE );
237  libmesh_assert_less_equal ( std::fabs(eta), 1.0+TOLERANCE );
238  libmesh_assert_less_equal ( std::fabs(zeta), 1.0+TOLERANCE );
239 
240  switch (j)
241  {
242  // d()/dxi
243  case 0:
244  {
245  switch(i)
246  {
247  case 0:
248  case 2:
249  case 8:
250  case 10:
251  return RealGradient();
252  case 1:
253  {
254  if (elem->point(1) > elem->point(2))
255  return RealGradient( 0.0, -0.125*(1.0-zeta) );
256  else
257  return RealGradient( 0.0, 0.125*(1.0-zeta) );
258  }
259  case 3:
260  {
261  if (elem->point(3) > elem->point(0))
262  return RealGradient( 0.0, 0.125*(-1.0+zeta) );
263  else
264  return RealGradient( 0.0, -0.125*(-1.0+zeta) );
265  }
266  case 4:
267  {
268  if (elem->point(0) > elem->point(4))
269  return RealGradient( 0.0, 0.0, -0.125*(-1.0+eta) );
270  else
271  return RealGradient( 0.0, 0.0, 0.125*(-1.0+eta) );
272  }
273  case 5:
274  {
275  if (elem->point(1) > elem->point(5))
276  return RealGradient( 0.0, 0.0, -0.125*(1.0-eta) );
277  else
278  return RealGradient( 0.0, 0.0, 0.125*(1.0-eta) );
279  }
280  case 6:
281  {
282  if (elem->point(2) > elem->point(6))
283  return RealGradient( 0.0, 0.0, -0.125*(1.0+eta) );
284  else
285  return RealGradient( 0.0, 0.0, 0.125*(1.0+eta) );
286  }
287  case 7:
288  {
289  if (elem->point(3) > elem->point(7))
290  return RealGradient( 0.0, 0.0, -0.125*(-1.0-eta) );
291  else
292  return RealGradient( 0.0, 0.0, 0.125*(-1.0-eta) );
293  }
294  case 9:
295  {
296  if (elem->point(5) > elem->point(6))
297  return RealGradient( 0.0, -0.125*(1.0+zeta), 0.0 );
298  else
299  return RealGradient( 0.0, 0.125*(1.0+zeta), 0.0 );
300  }
301  case 11:
302  {
303  if (elem->point(4) > elem->point(7))
304  return RealGradient( 0.0, -0.125*(-1.0-zeta), 0.0 );
305  else
306  return RealGradient( 0.0, 0.125*(-1.0-zeta), 0.0 );
307  }
308  default:
309  libmesh_error_msg("Invalid i = " << i);
310  } // switch(i)
311 
312  } // j=0
313 
314  // d()/deta
315  case 1:
316  {
317  switch(i)
318  {
319  case 1:
320  case 3:
321  case 9:
322  case 11:
323  return RealGradient();
324  case 0:
325  {
326  if (elem->point(0) > elem->point(1))
327  return RealGradient( -0.125*(-1.0+zeta), 0.0, 0.0 );
328  else
329  return RealGradient( 0.125*(-1.0+zeta), 0.0, 0.0 );
330  }
331  case 2:
332  {
333  if (elem->point(2) > elem->point(3))
334  return RealGradient( 0.125*(1.0-zeta), 0.0, 0.0 );
335  else
336  return RealGradient( -0.125*(1.0-zeta), 0.0, 0.0 );
337  }
338  case 4:
339  {
340  if (elem->point(0) > elem->point(4))
341  return RealGradient( 0.0, 0.0, -0.125*(-1.0+xi) );
342  else
343  return RealGradient( 0.0, 0.0, 0.125*(-1.0+xi) );
344  }
345  case 5:
346  {
347  if (elem->point(1) > elem->point(5))
348  return RealGradient( 0.0, 0.0, -0.125*(-1.0-xi) );
349  else
350  return RealGradient( 0.0, 0.0, 0.125*(-1.0-xi) );
351  }
352  case 6:
353  {
354  if (elem->point(2) > elem->point(6))
355  return RealGradient( 0.0, 0.0, -0.125*(1.0+xi) );
356  else
357  return RealGradient( 0.0, 0.0, 0.125*(1.0+xi) );
358  }
359  case 7:
360  {
361  if (elem->point(3) > elem->point(7))
362  return RealGradient( 0.0, 0.0, -0.125*(1.0-xi) );
363  else
364  return RealGradient( 0.0, 0.0, 0.125*(1.0-xi) );
365  }
366  case 8:
367  {
368  if (elem->point(4) > elem->point(5))
369  return RealGradient( -0.125*(-1.0-zeta), 0.0, 0.0 );
370  else
371  return RealGradient( 0.125*(-1.0-zeta), 0.0, 0.0 );
372  }
373  case 10:
374  {
375  if (elem->point(7) > elem->point(6))
376  return RealGradient( -0.125*(1.0+zeta), 0.0, 0.0 );
377  else
378  return RealGradient( 0.125*(1.0+zeta), 0.0, 0.0 );
379  }
380  default:
381  libmesh_error_msg("Invalid i = " << i);
382  } // switch(i)
383 
384  } // j=1
385 
386  // d()/dzeta
387  case 2:
388  {
389  switch(i)
390  {
391  case 4:
392  case 5:
393  case 6:
394  case 7:
395  return RealGradient();
396 
397  case 0:
398  {
399  if (elem->point(0) > elem->point(1))
400  return RealGradient( -0.125*(-1.0+eta), 0.0, 0.0 );
401  else
402  return RealGradient( 0.125*(-1.0+eta), 0.0, 0.0 );
403  }
404  case 1:
405  {
406  if (elem->point(1) > elem->point(2))
407  return RealGradient( 0.0, -0.125*(-1.0-xi), 0.0 );
408  else
409  return RealGradient( 0.0, 0.125*(-1.0-xi), 0.0 );
410  }
411  case 2:
412  {
413  if (elem->point(2) > elem->point(3))
414  return RealGradient( 0.125*(-1.0-eta), 0.0, 0.0 );
415  else
416  return RealGradient( -0.125*(-1.0-eta), 0.0, 0.0 );
417  }
418  case 3:
419  {
420  if (elem->point(3) > elem->point(0))
421  return RealGradient( 0.0, 0.125*(-1.0+xi), 0.0 );
422  else
423  return RealGradient( 0.0, -0.125*(-1.0+xi), 0.0 );
424  }
425  case 8:
426  {
427  if (elem->point(4) > elem->point(5))
428  return RealGradient( -0.125*(1.0-eta), 0.0, 0.0 );
429  else
430  return RealGradient( 0.125*(1.0-eta), 0.0, 0.0 );
431  }
432  case 9:
433  {
434  if (elem->point(5) > elem->point(6))
435  return RealGradient( 0.0, -0.125*(1.0+xi), 0.0 );
436  else
437  return RealGradient( 0.0, 0.125*(1.0+xi), 0.0 );
438  }
439  case 10:
440  {
441  if (elem->point(7) > elem->point(6))
442  return RealGradient( -0.125*(1.0+eta), 0.0, 0.0 );
443  else
444  return RealGradient( 0.125*(1.0+eta), 0.0, 0.0 );
445  }
446  case 11:
447  {
448  if (elem->point(4) > elem->point(7))
449  return RealGradient( 0.0, -0.125*(1.0-xi), 0.0 );
450  else
451  return RealGradient( 0.0, 0.125*(1.0-xi), 0.0 );
452  }
453  default:
454  libmesh_error_msg("Invalid i = " << i);
455  } // switch(i)
456 
457  } // j = 2
458 
459  default:
460  libmesh_error_msg("Invalid j = " << j);
461  }
462 
463  return RealGradient();
464  }
465 
466  case TET10:
467  {
468  libmesh_assert_less (i, 6);
469 
470  libmesh_not_implemented();
471  return RealGradient();
472  }
473 
474  default:
475  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << elem->type());
476  }
477  }
478 
479  // unsupported order
480  default:
481  libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << totalorder);
482  }
483 
484 #endif
485 
486  libmesh_error_msg("We'll never get here!");
487  return RealGradient();
488 }
489 
490 
491 
492 template <>
494  const Order,
495  const unsigned int,
496  const unsigned int,
497  const Point &)
498 {
499  libmesh_error_msg("Nedelec elements require the element type \nbecause edge orientation is needed.");
500  return RealGradient();
501 }
502 
503 
504 
505 template <>
507  const Order order,
508  const unsigned int i,
509  const unsigned int j,
510  const Point & libmesh_dbg_var(p))
511 {
512 #if LIBMESH_DIM == 3
513 
514  libmesh_assert(elem);
515 
516  // j = 0 ==> d^2 phi / dxi^2
517  // j = 1 ==> d^2 phi / dxi deta
518  // j = 2 ==> d^2 phi / deta^2
519  // j = 3 ==> d^2 phi / dxi dzeta
520  // j = 4 ==> d^2 phi / deta dzeta
521  // j = 5 ==> d^2 phi / dzeta^2
522  libmesh_assert_less (j, 6);
523 
524  const Order totalorder = static_cast<Order>(order + elem->p_level());
525 
526  switch (totalorder)
527  {
528  // linear Lagrange shape functions
529  case FIRST:
530  {
531  switch (elem->type())
532  {
533  case HEX20:
534  case HEX27:
535  {
536  libmesh_assert_less (i, 12);
537 
538 #ifndef NDEBUG
539  const Real xi = p(0);
540  const Real eta = p(1);
541  const Real zeta = p(2);
542 #endif
543 
544  libmesh_assert_less_equal ( std::fabs(xi), 1.0+TOLERANCE );
545  libmesh_assert_less_equal ( std::fabs(eta), 1.0+TOLERANCE );
546  libmesh_assert_less_equal ( std::fabs(zeta), 1.0+TOLERANCE );
547 
548  switch (j)
549  {
550  // d^2()/dxi^2
551  case 0:
552  {
553  // All d^2()/dxi^2 derivatives for linear hexes are zero.
554  return RealGradient();
555  } // j=0
556 
557  // d^2()/dxideta
558  case 1:
559  {
560  switch(i)
561  {
562  case 0:
563  case 1:
564  case 2:
565  case 3:
566  case 8:
567  case 9:
568  case 10:
569  case 11:
570  return RealGradient();
571  case 4:
572  {
573  if (elem->point(0) > elem->point(4))
574  return RealGradient( 0.0, 0.0, -0.125 );
575  else
576  return RealGradient( 0.0, 0.0, 0.125 );
577  }
578  case 5:
579  {
580  if (elem->point(1) > elem->point(5))
581  return RealGradient( 0.0, 0.0, 0.125 );
582  else
583  return RealGradient( 0.0, 0.0, -0.125 );
584  }
585  case 6:
586  {
587  if (elem->point(2) > elem->point(6))
588  return RealGradient( 0.0, 0.0, -0.125 );
589  else
590  return RealGradient( 0.0, 0.0, 0.125 );
591  }
592  case 7:
593  {
594  if (elem->point(3) > elem->point(7))
595  return RealGradient( 0.0, 0.0, 0.125 );
596  else
597  return RealGradient( 0.0, 0.0, -0.125 );
598  }
599  default:
600  libmesh_error_msg("Invalid i = " << i);
601  } // switch(i)
602 
603  } // j=1
604 
605  // d^2()/deta^2
606  case 2:
607  {
608  // All d^2()/deta^2 derivatives for linear hexes are zero.
609  return RealGradient();
610  } // j = 2
611 
612  // d^2()/dxidzeta
613  case 3:
614  {
615  switch(i)
616  {
617  case 0:
618  case 2:
619  case 4:
620  case 5:
621  case 6:
622  case 7:
623  case 8:
624  case 10:
625  return RealGradient();
626 
627  case 1:
628  {
629  if (elem->point(1) > elem->point(2))
630  return RealGradient( 0.0, 0.125 );
631  else
632  return RealGradient( 0.0, -0.125 );
633  }
634  case 3:
635  {
636  if (elem->point(3) > elem->point(0))
637  return RealGradient( 0.0, -0.125 );
638  else
639  return RealGradient( 0.0, 0.125 );
640  }
641  case 9:
642  {
643  if (elem->point(5) > elem->point(6))
644  return RealGradient( 0.0, -0.125, 0.0 );
645  else
646  return RealGradient( 0.0, 0.125, 0.0 );
647  }
648  case 11:
649  {
650  if (elem->point(4) > elem->point(7))
651  return RealGradient( 0.0, 0.125, 0.0 );
652  else
653  return RealGradient( 0.0, -0.125, 0.0 );
654  }
655  default:
656  libmesh_error_msg("Invalid i = " << i);
657  } // switch(i)
658 
659  } // j = 3
660 
661  // d^2()/detadzeta
662  case 4:
663  {
664  switch(i)
665  {
666  case 1:
667  case 3:
668  case 4:
669  case 5:
670  case 6:
671  case 7:
672  case 9:
673  case 11:
674  return RealGradient();
675 
676  case 0:
677  {
678  if (elem->point(0) > elem->point(1))
679  return RealGradient( -0.125, 0.0, 0.0 );
680  else
681  return RealGradient( 0.125, 0.0, 0.0 );
682  }
683  case 2:
684  {
685  if (elem->point(2) > elem->point(3))
686  return RealGradient( 0.125, 0.0, 0.0 );
687  else
688  return RealGradient( -0.125, 0.0, 0.0 );
689  }
690  case 8:
691  {
692  if (elem->point(4) > elem->point(5))
693  return RealGradient( 0.125, 0.0, 0.0 );
694  else
695  return RealGradient( -0.125, 0.0, 0.0 );
696  }
697  case 10:
698  {
699  if (elem->point(7) > elem->point(6))
700  return RealGradient( -0.125, 0.0, 0.0 );
701  else
702  return RealGradient( 0.125, 0.0, 0.0 );
703  }
704  default:
705  libmesh_error_msg("Invalid i = " << i);
706  } // switch(i)
707 
708  } // j = 4
709 
710  // d^2()/dzeta^2
711  case 5:
712  {
713  // All d^2()/dzeta^2 derivatives for linear hexes are zero.
714  return RealGradient();
715  } // j = 5
716 
717  default:
718  libmesh_error_msg("Invalid j = " << j);
719  }
720 
721  return RealGradient();
722  }
723 
724  case TET10:
725  {
726  libmesh_assert_less (i, 6);
727 
728  libmesh_not_implemented();
729  return RealGradient();
730  }
731 
732  default:
733  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << elem->type());
734 
735  } //switch(type)
736 
737  } // case FIRST:
738  // unsupported order
739  default:
740  libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << totalorder);
741  }
742 
743 #endif
744 
745  libmesh_error_msg("We'll never get here!");
746  return RealGradient();
747 }
748 
749 } // namespace libMesh
RealVectorValue RealGradient
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
virtual ElemType type() const =0
unsigned int p_level() const
Definition: elem.h:2422
ElemType
Defines an enum for geometric element types.
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
static const Real TOLERANCE
The libMesh namespace provides an interface to certain functionality in the library.
libmesh_assert(j)
PetscErrorCode Vec Mat libmesh_dbg_var(j)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Point & point(const unsigned int i) const
Definition: elem.h:1809
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)