fe_nedelec_one_shape_3D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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 <>
29 RealGradient FE<3,NEDELEC_ONE>::shape(const ElemType,
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 <>
41 RealGradient FE<3,NEDELEC_ONE>::shape(const Elem * elem,
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 #else
186  return RealGradient();
187 #endif
188 }
189 
190 
191 
192 
193 template <>
194 RealGradient FE<3,NEDELEC_ONE>::shape_deriv(const ElemType,
195  const Order,
196  const unsigned int,
197  const unsigned int,
198  const Point &)
199 {
200  libmesh_error_msg("Nedelec elements require the element type \nbecause edge orientation is needed.");
201  return RealGradient();
202 }
203 
204 template <>
205 RealGradient FE<3,NEDELEC_ONE>::shape_deriv(const Elem * elem,
206  const Order order,
207  const unsigned int i,
208  const unsigned int j,
209  const Point & p)
210 {
211 #if LIBMESH_DIM == 3
212  libmesh_assert(elem);
213  libmesh_assert_less (j, 3);
214 
215  const Order totalorder = static_cast<Order>(order + elem->p_level());
216 
217  switch (totalorder)
218  {
219  case FIRST:
220  {
221  switch (elem->type())
222  {
223  case HEX20:
224  case HEX27:
225  {
226  libmesh_assert_less (i, 12);
227 
228  const Real xi = p(0);
229  const Real eta = p(1);
230  const Real zeta = p(2);
231 
232  // Even with a loose inverse_map tolerance we ought to
233  // be nearly on the element interior in master
234  // coordinates
235  libmesh_assert_less_equal ( std::fabs(xi), 1.0+TOLERANCE );
236  libmesh_assert_less_equal ( std::fabs(eta), 1.0+TOLERANCE );
237  libmesh_assert_less_equal ( std::fabs(zeta), 1.0+TOLERANCE );
238 
239  switch (j)
240  {
241  // d()/dxi
242  case 0:
243  {
244  switch(i)
245  {
246  case 0:
247  case 2:
248  case 8:
249  case 10:
250  return RealGradient();
251  case 1:
252  {
253  if (elem->point(1) > elem->point(2))
254  return RealGradient( 0.0, -0.125*(1.0-zeta) );
255  else
256  return RealGradient( 0.0, 0.125*(1.0-zeta) );
257  }
258  case 3:
259  {
260  if (elem->point(3) > elem->point(0))
261  return RealGradient( 0.0, 0.125*(-1.0+zeta) );
262  else
263  return RealGradient( 0.0, -0.125*(-1.0+zeta) );
264  }
265  case 4:
266  {
267  if (elem->point(0) > elem->point(4))
268  return RealGradient( 0.0, 0.0, -0.125*(-1.0+eta) );
269  else
270  return RealGradient( 0.0, 0.0, 0.125*(-1.0+eta) );
271  }
272  case 5:
273  {
274  if (elem->point(1) > elem->point(5))
275  return RealGradient( 0.0, 0.0, -0.125*(1.0-eta) );
276  else
277  return RealGradient( 0.0, 0.0, 0.125*(1.0-eta) );
278  }
279  case 6:
280  {
281  if (elem->point(2) > elem->point(6))
282  return RealGradient( 0.0, 0.0, -0.125*(1.0+eta) );
283  else
284  return RealGradient( 0.0, 0.0, 0.125*(1.0+eta) );
285  }
286  case 7:
287  {
288  if (elem->point(3) > elem->point(7))
289  return RealGradient( 0.0, 0.0, -0.125*(-1.0-eta) );
290  else
291  return RealGradient( 0.0, 0.0, 0.125*(-1.0-eta) );
292  }
293  case 9:
294  {
295  if (elem->point(5) > elem->point(6))
296  return RealGradient( 0.0, -0.125*(1.0+zeta), 0.0 );
297  else
298  return RealGradient( 0.0, 0.125*(1.0+zeta), 0.0 );
299  }
300  case 11:
301  {
302  if (elem->point(4) > elem->point(7))
303  return RealGradient( 0.0, -0.125*(-1.0-zeta), 0.0 );
304  else
305  return RealGradient( 0.0, 0.125*(-1.0-zeta), 0.0 );
306  }
307  default:
308  libmesh_error_msg("Invalid i = " << i);
309  } // switch(i)
310 
311  } // j=0
312 
313  // d()/deta
314  case 1:
315  {
316  switch(i)
317  {
318  case 1:
319  case 3:
320  case 9:
321  case 11:
322  return RealGradient();
323  case 0:
324  {
325  if (elem->point(0) > elem->point(1))
326  return RealGradient( -0.125*(-1.0+zeta), 0.0, 0.0 );
327  else
328  return RealGradient( 0.125*(-1.0+zeta), 0.0, 0.0 );
329  }
330  case 2:
331  {
332  if (elem->point(2) > elem->point(3))
333  return RealGradient( 0.125*(1.0-zeta), 0.0, 0.0 );
334  else
335  return RealGradient( -0.125*(1.0-zeta), 0.0, 0.0 );
336  }
337  case 4:
338  {
339  if (elem->point(0) > elem->point(4))
340  return RealGradient( 0.0, 0.0, -0.125*(-1.0+xi) );
341  else
342  return RealGradient( 0.0, 0.0, 0.125*(-1.0+xi) );
343  }
344  case 5:
345  {
346  if (elem->point(1) > elem->point(5))
347  return RealGradient( 0.0, 0.0, -0.125*(-1.0-xi) );
348  else
349  return RealGradient( 0.0, 0.0, 0.125*(-1.0-xi) );
350  }
351  case 6:
352  {
353  if (elem->point(2) > elem->point(6))
354  return RealGradient( 0.0, 0.0, -0.125*(1.0+xi) );
355  else
356  return RealGradient( 0.0, 0.0, 0.125*(1.0+xi) );
357  }
358  case 7:
359  {
360  if (elem->point(3) > elem->point(7))
361  return RealGradient( 0.0, 0.0, -0.125*(1.0-xi) );
362  else
363  return RealGradient( 0.0, 0.0, 0.125*(1.0-xi) );
364  }
365  case 8:
366  {
367  if (elem->point(4) > elem->point(5))
368  return RealGradient( -0.125*(-1.0-zeta), 0.0, 0.0 );
369  else
370  return RealGradient( 0.125*(-1.0-zeta), 0.0, 0.0 );
371  }
372  case 10:
373  {
374  if (elem->point(7) > elem->point(6))
375  return RealGradient( -0.125*(1.0+zeta), 0.0, 0.0 );
376  else
377  return RealGradient( 0.125*(1.0+zeta), 0.0, 0.0 );
378  }
379  default:
380  libmesh_error_msg("Invalid i = " << i);
381  } // switch(i)
382 
383  } // j=1
384 
385  // d()/dzeta
386  case 2:
387  {
388  switch(i)
389  {
390  case 4:
391  case 5:
392  case 6:
393  case 7:
394  return RealGradient();
395 
396  case 0:
397  {
398  if (elem->point(0) > elem->point(1))
399  return RealGradient( -0.125*(-1.0+eta), 0.0, 0.0 );
400  else
401  return RealGradient( 0.125*(-1.0+eta), 0.0, 0.0 );
402  }
403  case 1:
404  {
405  if (elem->point(1) > elem->point(2))
406  return RealGradient( 0.0, -0.125*(-1.0-xi), 0.0 );
407  else
408  return RealGradient( 0.0, 0.125*(-1.0-xi), 0.0 );
409  }
410  case 2:
411  {
412  if (elem->point(2) > elem->point(3))
413  return RealGradient( 0.125*(-1.0-eta), 0.0, 0.0 );
414  else
415  return RealGradient( -0.125*(-1.0-eta), 0.0, 0.0 );
416  }
417  case 3:
418  {
419  if (elem->point(3) > elem->point(0))
420  return RealGradient( 0.0, 0.125*(-1.0+xi), 0.0 );
421  else
422  return RealGradient( 0.0, -0.125*(-1.0+xi), 0.0 );
423  }
424  case 8:
425  {
426  if (elem->point(4) > elem->point(5))
427  return RealGradient( -0.125*(1.0-eta), 0.0, 0.0 );
428  else
429  return RealGradient( 0.125*(1.0-eta), 0.0, 0.0 );
430  }
431  case 9:
432  {
433  if (elem->point(5) > elem->point(6))
434  return RealGradient( 0.0, -0.125*(1.0+xi), 0.0 );
435  else
436  return RealGradient( 0.0, 0.125*(1.0+xi), 0.0 );
437  }
438  case 10:
439  {
440  if (elem->point(7) > elem->point(6))
441  return RealGradient( -0.125*(1.0+eta), 0.0, 0.0 );
442  else
443  return RealGradient( 0.125*(1.0+eta), 0.0, 0.0 );
444  }
445  case 11:
446  {
447  if (elem->point(4) > elem->point(7))
448  return RealGradient( 0.0, -0.125*(1.0-xi), 0.0 );
449  else
450  return RealGradient( 0.0, 0.125*(1.0-xi), 0.0 );
451  }
452  default:
453  libmesh_error_msg("Invalid i = " << i);
454  } // switch(i)
455 
456  } // j = 2
457 
458  default:
459  libmesh_error_msg("Invalid j = " << j);
460  }
461 
462  return RealGradient();
463  }
464 
465  case TET10:
466  {
467  libmesh_assert_less (i, 6);
468 
469  libmesh_not_implemented();
470  return RealGradient();
471  }
472 
473  default:
474  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << elem->type());
475  }
476  }
477 
478  // unsupported order
479  default:
480  libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << totalorder);
481  }
482 
483 #else
484  return RealGradient();
485 #endif
486 }
487 
488 
489 
490 template <>
491 RealGradient FE<3,NEDELEC_ONE>::shape_second_deriv(const ElemType,
492  const Order,
493  const unsigned int,
494  const unsigned int,
495  const Point &)
496 {
497  libmesh_error_msg("Nedelec elements require the element type \nbecause edge orientation is needed.");
498  return RealGradient();
499 }
500 
501 
502 
503 template <>
504 RealGradient FE<3,NEDELEC_ONE>::shape_second_deriv(const Elem * elem,
505  const Order order,
506  const unsigned int i,
507  const unsigned int j,
508  const Point & libmesh_dbg_var(p))
509 {
510 #if LIBMESH_DIM == 3
511 
512  libmesh_assert(elem);
513 
514  // j = 0 ==> d^2 phi / dxi^2
515  // j = 1 ==> d^2 phi / dxi deta
516  // j = 2 ==> d^2 phi / deta^2
517  // j = 3 ==> d^2 phi / dxi dzeta
518  // j = 4 ==> d^2 phi / deta dzeta
519  // j = 5 ==> d^2 phi / dzeta^2
520  libmesh_assert_less (j, 6);
521 
522  const Order totalorder = static_cast<Order>(order + elem->p_level());
523 
524  switch (totalorder)
525  {
526  // linear Lagrange shape functions
527  case FIRST:
528  {
529  switch (elem->type())
530  {
531  case HEX20:
532  case HEX27:
533  {
534  libmesh_assert_less (i, 12);
535 
536 #ifndef NDEBUG
537  const Real xi = p(0);
538  const Real eta = p(1);
539  const Real zeta = p(2);
540 #endif
541 
542  libmesh_assert_less_equal ( std::fabs(xi), 1.0+TOLERANCE );
543  libmesh_assert_less_equal ( std::fabs(eta), 1.0+TOLERANCE );
544  libmesh_assert_less_equal ( std::fabs(zeta), 1.0+TOLERANCE );
545 
546  switch (j)
547  {
548  // d^2()/dxi^2
549  case 0:
550  {
551  // All d^2()/dxi^2 derivatives for linear hexes are zero.
552  return RealGradient();
553  } // j=0
554 
555  // d^2()/dxideta
556  case 1:
557  {
558  switch(i)
559  {
560  case 0:
561  case 1:
562  case 2:
563  case 3:
564  case 8:
565  case 9:
566  case 10:
567  case 11:
568  return RealGradient();
569  case 4:
570  {
571  if (elem->point(0) > elem->point(4))
572  return RealGradient( 0.0, 0.0, -0.125 );
573  else
574  return RealGradient( 0.0, 0.0, 0.125 );
575  }
576  case 5:
577  {
578  if (elem->point(1) > elem->point(5))
579  return RealGradient( 0.0, 0.0, 0.125 );
580  else
581  return RealGradient( 0.0, 0.0, -0.125 );
582  }
583  case 6:
584  {
585  if (elem->point(2) > elem->point(6))
586  return RealGradient( 0.0, 0.0, -0.125 );
587  else
588  return RealGradient( 0.0, 0.0, 0.125 );
589  }
590  case 7:
591  {
592  if (elem->point(3) > elem->point(7))
593  return RealGradient( 0.0, 0.0, 0.125 );
594  else
595  return RealGradient( 0.0, 0.0, -0.125 );
596  }
597  default:
598  libmesh_error_msg("Invalid i = " << i);
599  } // switch(i)
600 
601  } // j=1
602 
603  // d^2()/deta^2
604  case 2:
605  {
606  // All d^2()/deta^2 derivatives for linear hexes are zero.
607  return RealGradient();
608  } // j = 2
609 
610  // d^2()/dxidzeta
611  case 3:
612  {
613  switch(i)
614  {
615  case 0:
616  case 2:
617  case 4:
618  case 5:
619  case 6:
620  case 7:
621  case 8:
622  case 10:
623  return RealGradient();
624 
625  case 1:
626  {
627  if (elem->point(1) > elem->point(2))
628  return RealGradient( 0.0, 0.125 );
629  else
630  return RealGradient( 0.0, -0.125 );
631  }
632  case 3:
633  {
634  if (elem->point(3) > elem->point(0))
635  return RealGradient( 0.0, -0.125 );
636  else
637  return RealGradient( 0.0, 0.125 );
638  }
639  case 9:
640  {
641  if (elem->point(5) > elem->point(6))
642  return RealGradient( 0.0, -0.125, 0.0 );
643  else
644  return RealGradient( 0.0, 0.125, 0.0 );
645  }
646  case 11:
647  {
648  if (elem->point(4) > elem->point(7))
649  return RealGradient( 0.0, 0.125, 0.0 );
650  else
651  return RealGradient( 0.0, -0.125, 0.0 );
652  }
653  default:
654  libmesh_error_msg("Invalid i = " << i);
655  } // switch(i)
656 
657  } // j = 3
658 
659  // d^2()/detadzeta
660  case 4:
661  {
662  switch(i)
663  {
664  case 1:
665  case 3:
666  case 4:
667  case 5:
668  case 6:
669  case 7:
670  case 9:
671  case 11:
672  return RealGradient();
673 
674  case 0:
675  {
676  if (elem->point(0) > elem->point(1))
677  return RealGradient( -0.125, 0.0, 0.0 );
678  else
679  return RealGradient( 0.125, 0.0, 0.0 );
680  }
681  case 2:
682  {
683  if (elem->point(2) > elem->point(3))
684  return RealGradient( 0.125, 0.0, 0.0 );
685  else
686  return RealGradient( -0.125, 0.0, 0.0 );
687  }
688  case 8:
689  {
690  if (elem->point(4) > elem->point(5))
691  return RealGradient( 0.125, 0.0, 0.0 );
692  else
693  return RealGradient( -0.125, 0.0, 0.0 );
694  }
695  case 10:
696  {
697  if (elem->point(7) > elem->point(6))
698  return RealGradient( -0.125, 0.0, 0.0 );
699  else
700  return RealGradient( 0.125, 0.0, 0.0 );
701  }
702  default:
703  libmesh_error_msg("Invalid i = " << i);
704  } // switch(i)
705 
706  } // j = 4
707 
708  // d^2()/dzeta^2
709  case 5:
710  {
711  // All d^2()/dzeta^2 derivatives for linear hexes are zero.
712  return RealGradient();
713  } // j = 5
714 
715  default:
716  libmesh_error_msg("Invalid j = " << j);
717  }
718 
719  return RealGradient();
720  }
721 
722  case TET10:
723  {
724  libmesh_assert_less (i, 6);
725 
726  libmesh_not_implemented();
727  return RealGradient();
728  }
729 
730  default:
731  libmesh_error_msg("ERROR: Unsupported 3D element type!: " << elem->type());
732 
733  } //switch(type)
734 
735  } // case FIRST:
736  // unsupported order
737  default:
738  libmesh_error_msg("ERROR: Unsupported 3D FE order!: " << totalorder);
739  }
740 
741 #else
742  return RealGradient();
743 #endif
744 }
745 
746 } // namespace libMesh
libMesh::RealGradient
RealVectorValue RealGradient
Definition: hp_coarsentest.h:50
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::Elem
The base class for all geometric element types.
Definition: elem.h:100
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:2555
libMesh::TOLERANCE
static const Real TOLERANCE
Definition: libmesh_common.h:139
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Point
A geometric point in (x,y,z) space.
Definition: point.h:38
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition: elem.h:1892
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 by doxygen