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

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo