fe_l2_hierarchic_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 // C++ inlcludes
20 
21 // Local includes
22 #include "libmesh/fe.h"
23 #include "libmesh/elem.h"
24 #include "libmesh/utility.h"
25 
26 
27 namespace libMesh
28 {
29 
30 
31 
32 template <>
34  const Order order,
35  const unsigned int i,
36  const Point& p)
37 {
38  libmesh_assert_less (i, order+1u);
39 
40  // Declare that we are using our own special power function
41  // from the Utility namespace. This saves typing later.
42  using Utility::pow;
43 
44  const Real xi = p(0);
45 
46  Real returnval = 1.;
47 
48  switch (i)
49  {
50  case 0:
51  returnval = .5*(1. - xi);
52  break;
53  case 1:
54  returnval = .5*(1. + xi);
55  break;
56  // All even-terms have the same form.
57  // (xi^p - 1.)/p!
58  case 2:
59  returnval = (xi*xi - 1.)/2.;
60  break;
61  case 4:
62  returnval = (pow<4>(xi) - 1.)/24.;
63  break;
64  case 6:
65  returnval = (pow<6>(xi) - 1.)/720.;
66  break;
67 
68  // All odd-terms have the same form.
69  // (xi^p - xi)/p!
70  case 3:
71  returnval = (xi*xi*xi - xi)/6.;
72  break;
73  case 5:
74  returnval = (pow<5>(xi) - xi)/120.;
75  break;
76  case 7:
77  returnval = (pow<7>(xi) - xi)/5040.;
78  break;
79  default:
80  Real denominator = 1.;
81  for (unsigned int n=1; n <= i; ++n)
82  {
83  returnval *= xi;
84  denominator *= n;
85  }
86  // Odd:
87  if (i % 2)
88  returnval = (returnval - xi)/denominator;
89  // Even:
90  else
91  returnval = (returnval - 1.)/denominator;
92  break;
93  }
94 
95  return returnval;
96 }
97 
98 
99 
100 template <>
102  const Order order,
103  const unsigned int i,
104  const Point& p)
105 {
106  libmesh_assert(elem);
107 
108  return FE<1,L2_HIERARCHIC>::shape(elem->type(), static_cast<Order>(order + elem->p_level()), i, p);
109 }
110 
111 
112 
113 template <>
115  const Order order,
116  const unsigned int i,
117  const unsigned int j,
118  const Point& p)
119 {
120  // only d()/dxi in 1D!
121 
122  libmesh_assert_equal_to (j, 0);
123  libmesh_assert_less (i, order+1u);
124 
125  // Declare that we are using our own special power function
126  // from the Utility namespace. This saves typing later.
127  using Utility::pow;
128 
129  const Real xi = p(0);
130 
131  Real returnval = 1.;
132 
133  switch (i)
134  {
135  case 0:
136  returnval = -.5;
137  break;
138  case 1:
139  returnval = .5;
140  break;
141  // All even-terms have the same form.
142  // xi^(p-1)/(p-1)!
143  case 2:
144  returnval = xi;
145  break;
146  case 4:
147  returnval = pow<3>(xi)/6.;
148  break;
149  case 6:
150  returnval = pow<5>(xi)/120.;
151  break;
152  // All odd-terms have the same form.
153  // (p*xi^(p-1) - 1.)/p!
154  case 3:
155  returnval = (3*xi*xi - 1.)/6.;
156  break;
157  case 5:
158  returnval = (5.*pow<4>(xi) - 1.)/120.;
159  break;
160  case 7:
161  returnval = (7.*pow<6>(xi) - 1.)/5040.;
162  break;
163  default:
164  Real denominator = 1.;
165  for (unsigned int n=1; n != i; ++n)
166  {
167  returnval *= xi;
168  denominator *= n;
169  }
170  // Odd:
171  if (i % 2)
172  returnval = (i * returnval - 1.)/denominator/i;
173  // Even:
174  else
175  returnval = returnval/denominator;
176  break;
177  }
178 
179  return returnval;
180 }
181 
182 
183 
184 template <>
186  const Order order,
187  const unsigned int i,
188  const unsigned int j,
189  const Point& p)
190 {
191  libmesh_assert(elem);
192 
193  return FE<1,L2_HIERARCHIC>::shape_deriv(elem->type(),
194  static_cast<Order>(order + elem->p_level()), i, j, p);
195 }
196 
197 
198 
199 template <>
201  const Order order,
202  const unsigned int i,
203  const unsigned int j,
204  const Point& p)
205 {
206  // only d2()/d2xi in 1D!
207 
208  libmesh_assert_equal_to (j, 0);
209  libmesh_assert_less (i, order+1u);
210 
211  // Declare that we are using our own special power function
212  // from the Utility namespace. This saves typing later.
213  using Utility::pow;
214 
215  const Real xi = p(0);
216 
217  Real returnval = 1.;
218 
219  switch (i)
220  {
221  case 0:
222  case 1:
223  returnval = 0;
224  break;
225  // All terms have the same form.
226  // xi^(p-2)/(p-2)!
227  case 2:
228  returnval = 1;
229  break;
230  case 3:
231  returnval = xi;
232  break;
233  case 4:
234  returnval = pow<2>(xi)/2.;
235  break;
236  case 5:
237  returnval = pow<3>(xi)/6.;
238  break;
239  case 6:
240  returnval = pow<4>(xi)/24.;
241  break;
242  case 7:
243  returnval = pow<5>(xi)/120.;
244  break;
245 
246  default:
247  Real denominator = 1.;
248  for (unsigned int n=1; n != i; ++n)
249  {
250  returnval *= xi;
251  denominator *= n;
252  }
253  // Odd:
254  if (i % 2)
255  returnval = (i * returnval - 1.)/denominator/i;
256  // Even:
257  else
258  returnval = returnval/denominator;
259  break;
260  }
261 
262  return returnval;
263 }
264 
265 
266 
267 template <>
269  const Order order,
270  const unsigned int i,
271  const unsigned int j,
272  const Point& p)
273 {
274  libmesh_assert(elem);
275 
277  static_cast<Order>(order + elem->p_level()), i, j, p);
278 }
279 
280 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo