fe_hermite_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 namespace
27 {
28 using namespace libMesh;
29 
30 // Compute the static coefficients for an element
31 void hermite_compute_coefs(const Elem* elem, Real & d1xd1x, Real & d2xd2x)
32 {
33  const Order mapping_order (elem->default_order());
34  const ElemType mapping_elem_type (elem->type());
35  const int n_mapping_shape_functions =
36  FE<1,LAGRANGE>::n_shape_functions(mapping_elem_type,
37  mapping_order);
38 
39  // Degrees of freedom are at vertices and edge midpoints
40  std::vector<Point> dofpt;
41  dofpt.push_back(Point(-1));
42  dofpt.push_back(Point(1));
43 
44  // Mapping functions - first derivatives at each dofpt
45  std::vector<Real> dxdxi(2);
46  std::vector<Real> dxidx(2);
47 
48  for (int p = 0; p != 2; ++p)
49  {
50  dxdxi[p] = 0;
51  for (int i = 0; i != n_mapping_shape_functions; ++i)
52  {
54  (mapping_elem_type, mapping_order, i, 0, dofpt[p]);
55  dxdxi[p] += elem->point(i)(0) * ddxi;
56  }
57  }
58 
59  // Calculate derivative scaling factors
60 
61  d1xd1x = dxdxi[0];
62  d2xd2x = dxdxi[1];
63 }
64 
65 
66 } // end anonymous namespace
67 
68 
69 namespace libMesh
70 {
71 
72 
73 template<>
75  (const unsigned int i, const Real xi)
76 {
77  using Utility::pow;
78 
79  switch (i)
80  {
81  case 0:
82  return 1.5 * xi;
83  case 1:
84  return -1.5 * xi;
85  case 2:
86  return 0.5 * (-1. + 3.*xi);
87  case 3:
88  return 0.5 * (1. + 3.*xi);
89  case 4:
90  return (8.*xi*xi + 4.*(xi*xi-1.))/24.;
91  case 5:
92  return (8.*xi*xi*xi + 12.*xi*(xi*xi-1.))/120.;
93 // case 6:
94 // return (8.*pow<4>(xi) + 20.*xi*xi*(xi*xi-1.) +
95 // 2.*(xi*xi-1)*(xi*xi-1))/720.;
96  default:
97  Real denominator = 720., xipower = 1.;
98  for (unsigned n=6; n != i; ++n)
99  {
100  xipower *= xi;
101  denominator *= (n+1);
102  }
103  return (8.*pow<4>(xi)*xipower +
104  (8.*(i-4)+4.)*xi*xi*xipower*(xi*xi-1.) +
105  (i-4)*(i-5)*xipower*(xi*xi-1.)*(xi*xi-1.))/denominator;
106  }
107 
108  libmesh_error();
109  return 0.;
110 }
111 
112 
113 
114 template<>
116  (const unsigned int i, const Real xi)
117 {
118  switch (i)
119  {
120  case 0:
121  return 0.75 * (-1. + xi*xi);
122  case 1:
123  return 0.75 * (1. - xi*xi);
124  case 2:
125  return 0.25 * (-1. - 2.*xi + 3.*xi*xi);
126  case 3:
127  return 0.25 * (-1. + 2.*xi + 3.*xi*xi);
128  case 4:
129  return 4.*xi * (xi*xi-1.)/24.;
130  case 5:
131  return (4*xi*xi*(xi*xi-1.) + (xi*xi-1.)*(xi*xi-1.))/120.;
132 // case 6:
133 // return (4*xi*xi*xi*(xi*xi-1.) + 2*xi*(xi*xi-1.)*(xi*xi-1.))/720.;
134  default:
135  Real denominator = 720., xipower = 1.;
136  for (unsigned n=6; n != i; ++n)
137  {
138  xipower *= xi;
139  denominator *= (n+1);
140  }
141  return (4*xi*xi*xi*xipower*(xi*xi-1.) +
142  (i-4)*xi*xipower*(xi*xi-1.)*(xi*xi-1.))/denominator;
143  }
144 
145  libmesh_error();
146  return 0.;
147 }
148 
149 template<>
151  (const unsigned int i, const Real xi)
152 {
153  switch (i)
154  {
155  case 0:
156  return 0.25 * (2. - 3.*xi + xi*xi*xi);
157  case 1:
158  return 0.25 * (2. + 3.*xi - xi*xi*xi);
159  case 2:
160  return 0.25 * (1. - xi - xi*xi + xi*xi*xi);
161  case 3:
162  return 0.25 * (-1. - xi + xi*xi + xi*xi*xi);
163  // All high order terms have the form x^(p-4)(x^2-1)^2/p!
164  case 4:
165  return (xi*xi-1.) * (xi*xi-1.)/24.;
166  case 5:
167  return xi * (xi*xi-1.) * (xi*xi-1.)/120.;
168 // case 6:
169 // return xi*xi * (xi*xi-1.) * (xi*xi-1.)/720.;
170  default:
171  Real denominator = 720., xipower = 1.;
172  for (unsigned n=6; n != i; ++n)
173  {
174  xipower *= xi;
175  denominator *= (n+1);
176  }
177  return (xi*xi*xipower*(xi*xi-1.)*(xi*xi-1.))/denominator;
178 
179  }
180 
181  libmesh_error();
182  return 0.;
183 }
184 
185 
186 template <>
188  const Order,
189  const unsigned int,
190  const Point&)
191 {
192  libMesh::err << "Hermite elements require the real element\n"
193  << "to construct gradient-based degrees of freedom."
194  << std::endl;
195 
196  libmesh_error();
197  return 0.;
198 }
199 
200 
201 
202 template <>
204  const Order order,
205  const unsigned int i,
206  const Point& p)
207 {
208  libmesh_assert(elem);
209 
210  // Coefficient naming: d(1)d(2n) is the coefficient of the
211  // global shape function corresponding to value 1 in terms of the
212  // local shape function corresponding to normal derivative 2
213  Real d1xd1x, d2xd2x;
214 
215  hermite_compute_coefs(elem, d1xd1x, d2xd2x);
216 
217  const ElemType type = elem->type();
218 
219  const Order totalorder = static_cast<Order>(order + elem->p_level());
220 
221  switch (totalorder)
222  {
223  // Hermite cubic shape functions
224  case THIRD:
225  {
226  switch (type)
227  {
228  // C1 functions on the C1 cubic edge
229  case EDGE2:
230  case EDGE3:
231  {
232  libmesh_assert_less (i, 4);
233 
234  switch (i)
235  {
236  case 0:
237  return FEHermite<1>::hermite_raw_shape(0, p(0));
238  case 1:
239  return d1xd1x * FEHermite<1>::hermite_raw_shape(2, p(0));
240  case 2:
241  return FEHermite<1>::hermite_raw_shape(1, p(0));
242  case 3:
243  return d2xd2x * FEHermite<1>::hermite_raw_shape(3, p(0));
244  default:
245  return FEHermite<1>::hermite_raw_shape(i, p(0));
246  }
247  }
248  default:
249  libMesh::err << "ERROR: Unsupported element type!" << std::endl;
250  libmesh_error();
251  }
252  }
253  // by default throw an error
254  default:
255  libMesh::err << "ERROR: Unsupported polynomial order!" << std::endl;
256  libmesh_error();
257  }
258 
259  libmesh_error();
260  return 0.;
261 }
262 
263 
264 
265 template <>
267  const Order,
268  const unsigned int,
269  const unsigned int,
270  const Point&)
271 {
272  libMesh::err << "Hermite elements require the real element\n"
273  << "to construct gradient-based degrees of freedom."
274  << std::endl;
275 
276  libmesh_error();
277  return 0.;
278 }
279 
280 
281 
282 template <>
284  const Order order,
285  const unsigned int i,
286  const unsigned int,
287  const Point& p)
288 {
289  libmesh_assert(elem);
290 
291  // Coefficient naming: d(1)d(2n) is the coefficient of the
292  // global shape function corresponding to value 1 in terms of the
293  // local shape function corresponding to normal derivative 2
294  Real d1xd1x, d2xd2x;
295 
296  hermite_compute_coefs(elem, d1xd1x, d2xd2x);
297 
298  const ElemType type = elem->type();
299 
300  const Order totalorder = static_cast<Order>(order + elem->p_level());
301 
302  switch (totalorder)
303  {
304  // Hermite cubic shape functions
305  case THIRD:
306  {
307  switch (type)
308  {
309  // C1 functions on the C1 cubic edge
310  case EDGE2:
311  case EDGE3:
312  {
313  switch (i)
314  {
315  case 0:
317  case 1:
318  return d1xd1x * FEHermite<1>::hermite_raw_shape_deriv(2, p(0));
319  case 2:
321  case 3:
322  return d2xd2x * FEHermite<1>::hermite_raw_shape_deriv(3, p(0));
323  default:
325  }
326  }
327  default:
328  libMesh::err << "ERROR: Unsupported element type!" << std::endl;
329  libmesh_error();
330  }
331  }
332  // by default throw an error
333  default:
334  libMesh::err << "ERROR: Unsupported polynomial order!" << std::endl;
335  libmesh_error();
336  }
337 
338  libmesh_error();
339  return 0.;
340 }
341 
342 
343 
344 template <>
346  const Order order,
347  const unsigned int i,
348  const unsigned int,
349  const Point& p)
350 {
351  libmesh_assert(elem);
352 
353  // Coefficient naming: d(1)d(2n) is the coefficient of the
354  // global shape function corresponding to value 1 in terms of the
355  // local shape function corresponding to normal derivative 2
356  Real d1xd1x, d2xd2x;
357 
358  hermite_compute_coefs(elem, d1xd1x, d2xd2x);
359 
360  const ElemType type = elem->type();
361 
362  const Order totalorder = static_cast<Order>(order + elem->p_level());
363 
364  switch (totalorder)
365  {
366  // Hermite cubic shape functions
367  case THIRD:
368  {
369  switch (type)
370  {
371  // C1 functions on the C1 cubic edge
372  case EDGE2:
373  case EDGE3:
374  {
375  switch (i)
376  {
377  case 0:
379  case 1:
380  return d1xd1x * FEHermite<1>::hermite_raw_shape_second_deriv(2, p(0));
381  case 2:
383  case 3:
384  return d2xd2x * FEHermite<1>::hermite_raw_shape_second_deriv(3, p(0));
385  default:
387  }
388  }
389  default:
390  libMesh::err << "ERROR: Unsupported element type!" << std::endl;
391  libmesh_error();
392  }
393  }
394  // by default throw an error
395  default:
396  libMesh::err << "ERROR: Unsupported polynomial order!" << std::endl;
397  libmesh_error();
398  }
399 
400  libmesh_error();
401  return 0.;
402 }
403 
404 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo