fe_clough_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 
25 
26 // Anonymous namespace for persistant variables.
27 // This allows us to cache the global-to-local mapping transformation
28 // This should also screw up multithreading royally
29 namespace
30 {
31  using namespace libMesh;
32 
33  static dof_id_type old_elem_id = DofObject::invalid_id;
34  // Coefficient naming: d(1)d(2n) is the coefficient of the
35  // global shape function corresponding to value 1 in terms of the
36  // local shape function corresponding to normal derivative 2
37  static Real d1xd1x, d2xd2x;
38 
39 Real clough_raw_shape_second_deriv(const unsigned int basis_num,
40  const unsigned int deriv_type,
41  const Point& p);
42 Real clough_raw_shape_deriv(const unsigned int basis_num,
43  const unsigned int deriv_type,
44  const Point& p);
45 Real clough_raw_shape(const unsigned int basis_num,
46  const Point& p);
47 
48 
49 // Compute the static coefficients for an element
50 void clough_compute_coefs(const Elem* elem)
51 {
52  // Using static globals for old_elem_id, etc. will fail
53  // horribly with more than one thread.
54  libmesh_assert_equal_to (libMesh::n_threads(), 1);
55 
56  // Coefficients are cached from old elements
57  if (elem->id() == old_elem_id)
58  return;
59 
60  old_elem_id = elem->id();
61 
62  const Order mapping_order (elem->default_order());
63  const ElemType mapping_elem_type (elem->type());
64  const int n_mapping_shape_functions =
65  FE<1,LAGRANGE>::n_shape_functions(mapping_elem_type,
66  mapping_order);
67 
68  // Degrees of freedom are at vertices and edge midpoints
69  std::vector<Point> dofpt;
70  dofpt.push_back(Point(0));
71  dofpt.push_back(Point(1));
72 
73  // Mapping functions - first derivatives at each dofpt
74  std::vector<Real> dxdxi(2);
75  std::vector<Real> dxidx(2);
76 
77  for (int p = 0; p != 2; ++p)
78  {
79  for (int i = 0; i != n_mapping_shape_functions; ++i)
80  {
82  (mapping_elem_type, mapping_order, i, 0, dofpt[p]);
83  dxdxi[p] += dofpt[p](0) * ddxi;
84  }
85  }
86 
87  // Calculate derivative scaling factors
88 
89  d1xd1x = dxdxi[0];
90  d2xd2x = dxdxi[1];
91 }
92 
93 
94  // Return shape function second derivatives on the unit interval
95 Real clough_raw_shape_second_deriv(const unsigned int basis_num,
96  const unsigned int deriv_type,
97  const Point& p)
98 {
99  Real xi = p(0);
100 
101  switch (deriv_type)
102  {
103 
104  // second derivative in xi-xi direction
105  case 0:
106  switch (basis_num)
107  {
108  case 0:
109  return -6 + 12*xi;
110  case 1:
111  return 6 - 12*xi;
112  case 2:
113  return -4 + 6*xi;
114  case 3:
115  return -2 + 6*xi;
116  }
117  }
118 
119  libmesh_error();
120  return 0.;
121 }
122 
123 
124 
125 Real clough_raw_shape_deriv(const unsigned int basis_num,
126  const unsigned int deriv_type,
127  const Point& p)
128 {
129  Real xi = p(0);
130 
131  switch (deriv_type)
132  {
133  case 0:
134  switch (basis_num)
135  {
136  case 0:
137  return -6*xi + 6*xi*xi;
138  case 1:
139  return 6*xi - 6*xi*xi;
140  case 2:
141  return 1 - 4*xi + 3*xi*xi;
142  case 3:
143  return -2*xi + 3*xi*xi;
144  }
145  }
146 
147  libmesh_error();
148  return 0.;
149 }
150 
151 Real clough_raw_shape(const unsigned int basis_num,
152  const Point& p)
153 {
154  Real xi = p(0);
155 
156  switch (basis_num)
157  {
158  case 0:
159  return 1 - 3*xi*xi + 2*xi*xi*xi;
160  case 1:
161  return 3*xi*xi - 2*xi*xi*xi;
162  case 2:
163  return xi - 2*xi*xi + xi*xi*xi;
164  case 3:
165  return -xi*xi + xi*xi*xi;
166  }
167 
168  libmesh_error();
169  return 0.;
170 }
171 
172 
173 } // end anonymous namespace
174 
175 
176 namespace libMesh
177 {
178 
179 
180 template <>
182  const Order,
183  const unsigned int,
184  const Point&)
185 {
186  libMesh::err << "Clough-Tocher elements require the real element\n"
187  << "to construct gradient-based degrees of freedom."
188  << std::endl;
189 
190  libmesh_error();
191  return 0.;
192 }
193 
194 
195 
196 template <>
198  const Order order,
199  const unsigned int i,
200  const Point& p)
201 {
202  libmesh_assert(elem);
203 
204  clough_compute_coefs(elem);
205 
206  const ElemType type = elem->type();
207 
208  const Order totalorder = static_cast<Order>(order + elem->p_level());
209 
210  switch (totalorder)
211  {
212  // 3rd-order C1 cubic element
213  case THIRD:
214  {
215  switch (type)
216  {
217  // C1 functions on the C1 cubic edge
218  case EDGE2:
219  case EDGE3:
220  {
221  libmesh_assert_less (i, 4);
222 
223  switch (i)
224  {
225  case 0:
226  return clough_raw_shape(0, p);
227  case 1:
228  return clough_raw_shape(1, p);
229  case 2:
230  return d1xd1x * clough_raw_shape(2, p);
231  case 3:
232  return d2xd2x * clough_raw_shape(3, p);
233  default:
234  libmesh_error();
235  }
236  }
237  default:
238  libMesh::err << "ERROR: Unsupported element type!" << std::endl;
239  libmesh_error();
240  }
241  }
242  // by default throw an error
243  default:
244  libMesh::err << "ERROR: Unsupported polynomial order!" << std::endl;
245  libmesh_error();
246  }
247 
248  libmesh_error();
249  return 0.;
250 }
251 
252 
253 
254 template <>
256  const Order,
257  const unsigned int,
258  const unsigned int,
259  const Point&)
260 {
261  libMesh::err << "Clough-Tocher elements require the real element\n"
262  << "to construct gradient-based degrees of freedom."
263  << std::endl;
264 
265  libmesh_error();
266  return 0.;
267 }
268 
269 
270 
271 template <>
273  const Order order,
274  const unsigned int i,
275  const unsigned int j,
276  const Point& p)
277 {
278  libmesh_assert(elem);
279 
280  clough_compute_coefs(elem);
281 
282  const ElemType type = elem->type();
283 
284  const Order totalorder = static_cast<Order>(order + elem->p_level());
285 
286  switch (totalorder)
287  {
288  // 3rd-order C1 cubic element
289  case THIRD:
290  {
291  switch (type)
292  {
293  // C1 functions on the C1 cubic edge
294  case EDGE2:
295  case EDGE3:
296  {
297  switch (i)
298  {
299  case 0:
300  return clough_raw_shape_deriv(0, j, p);
301  case 1:
302  return clough_raw_shape_deriv(1, j, p);
303  case 2:
304  return d1xd1x * clough_raw_shape_deriv(2, j, p);
305  case 3:
306  return d2xd2x * clough_raw_shape_deriv(3, j, p);
307  default:
308  libmesh_error();
309  }
310  }
311  default:
312  libMesh::err << "ERROR: Unsupported element type!" << std::endl;
313  libmesh_error();
314  }
315  }
316  // by default throw an error
317  default:
318  libMesh::err << "ERROR: Unsupported polynomial order!" << std::endl;
319  libmesh_error();
320  }
321 
322  libmesh_error();
323  return 0.;
324 }
325 
326 
327 
328 template <>
330  const Order order,
331  const unsigned int i,
332  const unsigned int j,
333  const Point& p)
334 {
335  libmesh_assert(elem);
336 
337  clough_compute_coefs(elem);
338 
339  const ElemType type = elem->type();
340 
341  const Order totalorder = static_cast<Order>(order + elem->p_level());
342 
343  switch (totalorder)
344  {
345  // 3rd-order C1 cubic element
346  case THIRD:
347  {
348  switch (type)
349  {
350  // C1 functions on the C1 cubic edge
351  case EDGE2:
352  case EDGE3:
353  {
354  switch (i)
355  {
356  case 0:
357  return clough_raw_shape_second_deriv(0, j, p);
358  case 1:
359  return clough_raw_shape_second_deriv(1, j, p);
360  case 2:
361  return d1xd1x * clough_raw_shape_second_deriv(2, j, p);
362  case 3:
363  return d2xd2x * clough_raw_shape_second_deriv(3, j, p);
364  default:
365  libmesh_error();
366  }
367  }
368  default:
369  libMesh::err << "ERROR: Unsupported element type!" << std::endl;
370  libmesh_error();
371  }
372  }
373  // by default throw an error
374  default:
375  libMesh::err << "ERROR: Unsupported polynomial order!" << std::endl;
376  libmesh_error();
377  }
378 
379  libmesh_error();
380  return 0.;
381 }
382 
383 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo