fe_szabab_shape_1D.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 
20 // C++ includes
21 
22 // Local includes
23 #include "libmesh/libmesh_config.h"
24 
25 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
26 
27 #include "libmesh/fe.h"
28 #include "libmesh/elem.h"
29 
30 namespace libMesh
31 {
32 
33 
34 template <>
36  const Order libmesh_dbg_var(order),
37  const unsigned int i,
38  const Point& p)
39 {
40  const Real xi = p(0);
41  const Real xi2 = xi*xi;
42 
43 
44  // Use this libmesh_assert rather than a switch with a single entry...
45  // It will go away in optimized mode, essentially has the same effect.
46  libmesh_assert_less_equal (order, SEVENTH);
47 
48 // switch (order)
49 // {
50 // case FIRST:
51 // case SECOND:
52 // case THIRD:
53 // case FOURTH:
54 // case FIFTH:
55 // case SIXTH:
56 // case SEVENTH:
57 
58  switch(i)
59  {
60  //nodal shape functions
61  case 0: return 1./2.-1./2.*xi;
62  case 1: return 1./2.+1./2.*xi;
63  case 2: return 1./4. *2.4494897427831780982*(xi2-1.);
64  case 3: return 1./4. *3.1622776601683793320*(xi2-1.)*xi;
65  case 4: return 1./16. *3.7416573867739413856*((5.*xi2-6.)*xi2+1.);
66  case 5: return 3./16. *1.4142135623730950488*(3.+(-10.+7.*xi2)*xi2)*xi;
67  case 6: return 1./32. *4.6904157598234295546*(-1.+(15.+(-35.+21.*xi2)*xi2)*xi2);
68  case 7: return 1./32. *5.0990195135927848300*(-5.+(35.+(-63.+33.*xi2)*xi2)*xi2)*xi;
69  case 8: return 1./256.*5.4772255750516611346*(5.+(-140.+(630.+(-924.+429.*xi2)*xi2)*xi2)*xi2);
70 
71  default:
72  libMesh::err << "Invalid shape function index!" << std::endl;
73  libmesh_error();
74  }
75 
76 // default:
77 // {
78 // libMesh::err << "ERROR: Unsupported polynomial order!" << std::endl;
79 // libmesh_error();
80 // }
81 // }
82 
83  libmesh_error();
84  return 0.;
85 }
86 
87 
88 
89 template <>
91  const Order order,
92  const unsigned int i,
93  const Point& p)
94 {
95  libmesh_assert(elem);
96 
97  return FE<1,SZABAB>::shape(elem->type(), static_cast<Order>(order + elem->p_level()), i, p);
98 }
99 
100 
101 
102 template <>
104  const Order libmesh_dbg_var(order),
105  const unsigned int i,
106  const unsigned int libmesh_dbg_var(j),
107  const Point& p)
108 {
109  // only d()/dxi in 1D!
110  libmesh_assert_equal_to (j, 0);
111 
112  const Real xi = p(0);
113  const Real xi2 = xi*xi;
114 
115  // Use this libmesh_assert rather than a switch with a single entry...
116  // It will go away in optimized mode, essentially has the same effect.
117  libmesh_assert_less_equal (order, SEVENTH);
118 
119 // switch (order)
120 // {
121 // case FIRST:
122 // case SECOND:
123 // case THIRD:
124 // case FOURTH:
125 // case FIFTH:
126 // case SIXTH:
127 // case SEVENTH:
128 
129  switch(i)
130  {
131  case 0: return -1./2.;
132  case 1: return 1./2.;
133  case 2: return 1./2.*2.4494897427831780982*xi;
134  case 3: return -1./4.*3.1622776601683793320+3./4.*3.1622776601683793320*xi2;
135  case 4: return 1./16.*3.7416573867739413856*(-12.+20*xi2)*xi;
136  case 5: return 9./16.*1.4142135623730950488+(-45./8.*1.4142135623730950488+105./16.*1.4142135623730950488*xi2)*xi2;
137  case 6: return 1./32.*4.6904157598234295546*(30.+(-140.+126.*xi2)*xi2)*xi;
138  case 7: return -5./32.*5.0990195135927848300+(105./32.*5.0990195135927848300+(-315./32.*5.0990195135927848300+231./32.*5.0990195135927848300*xi2)*xi2)*xi2;
139  case 8: return 1./256.*5.4772255750516611346*(-280.+(2520.+(-5544.+3432.*xi2)*xi2)*xi2)*xi;
140 
141  default:
142  libMesh::err << "Invalid shape function index!" << std::endl;
143  libmesh_error();
144  }
145 
146 // default:
147 // {
148 // libMesh::err << "ERROR: Unsupported polynomial order!" << std::endl;
149 // libmesh_error();
150 // }
151 // }
152 
153  libmesh_error();
154  return 0.;
155 }
156 
157 
158 
159 template <>
161  const Order order,
162  const unsigned int i,
163  const unsigned int j,
164  const Point& p)
165 {
166  libmesh_assert(elem);
167 
168  return FE<1,SZABAB>::shape_deriv(elem->type(),
169  static_cast<Order>(order + elem->p_level()), i, j, p);
170 }
171 
172 
173 
174 template <>
176  const Order,
177  const unsigned int,
178  const unsigned int,
179  const Point&)
180 {
181  static bool warning_given = false;
182 
183  if (!warning_given)
184  libMesh::err << "Second derivatives for Szabab elements "
185  << " are not yet implemented!"
186  << std::endl;
187 
188  warning_given = true;
189  return 0.;
190 }
191 
192 
193 
194 template <>
196  const Order,
197  const unsigned int,
198  const unsigned int,
199  const Point&)
200 {
201  static bool warning_given = false;
202 
203  if (!warning_given)
204  libMesh::err << "Second derivatives for Szabab elements "
205  << " are not yet implemented!"
206  << std::endl;
207 
208  warning_given = true;
209  return 0.;
210 }
211 
212 } // namespace libMesh
213 
214 
215 #endif //LIBMESH_ENABLE_HIGHER_ORDER_SHAPES

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

Hosted By:
SourceForge.net Logo