libMesh
fe_clough_shape_1D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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++ includes
20 
21 // Local includes
22 #include "libmesh/fe.h"
23 #include "libmesh/elem.h"
24 
25 
26 // Anonymous namespace for persistent 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 // Keep track of which element was most recently used to generate
34 // cached data
35 static dof_id_type old_elem_id = DofObject::invalid_id;
36 static const Elem * old_elem_ptr = libmesh_nullptr;
37 
38 // Coefficient naming: d(1)d(2n) is the coefficient of the
39 // global shape function corresponding to value 1 in terms of the
40 // local shape function corresponding to normal derivative 2
41 static Real d1xd1x, d2xd2x;
42 
43 Real clough_raw_shape_second_deriv(const unsigned int basis_num,
44  const unsigned int deriv_type,
45  const Point & p);
46 Real clough_raw_shape_deriv(const unsigned int basis_num,
47  const unsigned int deriv_type,
48  const Point & p);
49 Real clough_raw_shape(const unsigned int basis_num,
50  const Point & p);
51 
52 
53 // Compute the static coefficients for an element
54 void clough_compute_coefs(const Elem * elem)
55 {
56  // Using static globals for old_elem_id, etc. will fail
57  // horribly with more than one thread.
58  libmesh_assert_equal_to (libMesh::n_threads(), 1);
59 
60  // Coefficients are cached from old elements; we rely on that cache
61  // except in dbg mode
62 #ifndef DEBUG
63  if (elem->id() == old_elem_id &&
64  elem == old_elem_ptr)
65  return;
66 #endif
67 
68  const Order mapping_order (elem->default_order());
69  const ElemType mapping_elem_type (elem->type());
70  const int n_mapping_shape_functions =
71  FE<1,LAGRANGE>::n_shape_functions(mapping_elem_type,
72  mapping_order);
73 
74  // Degrees of freedom are at vertices and edge midpoints
75  std::vector<Point> dofpt;
76  dofpt.push_back(Point(0));
77  dofpt.push_back(Point(1));
78 
79  // Mapping functions - first derivatives at each dofpt
80  std::vector<Real> dxdxi(2);
81  std::vector<Real> dxidx(2);
82 
83  for (int p = 0; p != 2; ++p)
84  {
85  for (int i = 0; i != n_mapping_shape_functions; ++i)
86  {
88  (mapping_elem_type, mapping_order, i, 0, dofpt[p]);
89  dxdxi[p] += dofpt[p](0) * ddxi;
90  }
91  }
92 
93  // Calculate derivative scaling factors
94 
95 #ifdef DEBUG
96  // The cached factors should equal our calculations
97  if (elem->id() == old_elem_id &&
98  elem == old_elem_ptr)
99  {
100  libmesh_assert_equal_to(d1xd1x, dxdxi[0]);
101  libmesh_assert_equal_to(d2xd2x, dxdxi[1]);
102  }
103 #endif
104 
105  old_elem_id = elem->id();
106  old_elem_ptr = elem;
107 
108  d1xd1x = dxdxi[0];
109  d2xd2x = dxdxi[1];
110 }
111 
112 
113 // Return shape function second derivatives on the unit interval
114 Real clough_raw_shape_second_deriv(const unsigned int basis_num,
115  const unsigned int deriv_type,
116  const Point & p)
117 {
118  Real xi = p(0);
119 
120  switch (deriv_type)
121  {
122 
123  // second derivative in xi-xi direction
124  case 0:
125  switch (basis_num)
126  {
127  case 0:
128  return -6 + 12*xi;
129  case 1:
130  return 6 - 12*xi;
131  case 2:
132  return -4 + 6*xi;
133  case 3:
134  return -2 + 6*xi;
135 
136  default:
137  libmesh_error_msg("Invalid shape function index i = " <<
138  basis_num);
139  }
140 
141  default:
142  libmesh_error_msg("Invalid shape function derivative j = " <<
143  deriv_type);
144  }
145 
146  libmesh_error_msg("We'll never get here!");
147  return 0.;
148 }
149 
150 
151 
152 Real clough_raw_shape_deriv(const unsigned int basis_num,
153  const unsigned int deriv_type,
154  const Point & p)
155 {
156  Real xi = p(0);
157 
158  switch (deriv_type)
159  {
160  case 0:
161  switch (basis_num)
162  {
163  case 0:
164  return -6*xi + 6*xi*xi;
165  case 1:
166  return 6*xi - 6*xi*xi;
167  case 2:
168  return 1 - 4*xi + 3*xi*xi;
169  case 3:
170  return -2*xi + 3*xi*xi;
171 
172  default:
173  libmesh_error_msg("Invalid shape function index i = " <<
174  basis_num);
175  }
176 
177  default:
178  libmesh_error_msg("Invalid shape function derivative j = " <<
179  deriv_type);
180  }
181 
182  libmesh_error_msg("We'll never get here!");
183  return 0.;
184 }
185 
186 Real clough_raw_shape(const unsigned int basis_num,
187  const Point & p)
188 {
189  Real xi = p(0);
190 
191  switch (basis_num)
192  {
193  case 0:
194  return 1 - 3*xi*xi + 2*xi*xi*xi;
195  case 1:
196  return 3*xi*xi - 2*xi*xi*xi;
197  case 2:
198  return xi - 2*xi*xi + xi*xi*xi;
199  case 3:
200  return -xi*xi + xi*xi*xi;
201 
202  default:
203  libmesh_error_msg("Invalid shape function index i = " <<
204  basis_num);
205  }
206 
207  libmesh_error_msg("We'll never get here!");
208  return 0.;
209 }
210 
211 
212 } // end anonymous namespace
213 
214 
215 namespace libMesh
216 {
217 
218 
219 template <>
221  const Order,
222  const unsigned int,
223  const Point &)
224 {
225  libmesh_error_msg("Clough-Tocher elements require the real element \nto construct gradient-based degrees of freedom.");
226  return 0.;
227 }
228 
229 
230 
231 template <>
233  const Order order,
234  const unsigned int i,
235  const Point & p)
236 {
237  libmesh_assert(elem);
238 
239  clough_compute_coefs(elem);
240 
241  const ElemType type = elem->type();
242 
243  const Order totalorder = static_cast<Order>(order + elem->p_level());
244 
245  switch (totalorder)
246  {
247  // 3rd-order C1 cubic element
248  case THIRD:
249  {
250  switch (type)
251  {
252  // C1 functions on the C1 cubic edge
253  case EDGE2:
254  case EDGE3:
255  {
256  libmesh_assert_less (i, 4);
257 
258  switch (i)
259  {
260  case 0:
261  return clough_raw_shape(0, p);
262  case 1:
263  return clough_raw_shape(1, p);
264  case 2:
265  return d1xd1x * clough_raw_shape(2, p);
266  case 3:
267  return d2xd2x * clough_raw_shape(3, p);
268  default:
269  libmesh_error_msg("Invalid shape function index i = " << i);
270  }
271  }
272  default:
273  libmesh_error_msg("ERROR: Unsupported element type = " << type);
274  }
275  }
276  // by default throw an error
277  default:
278  libmesh_error_msg("ERROR: Unsupported polynomial order = " << totalorder);
279  }
280 
281  libmesh_error_msg("We'll never get here!");
282  return 0.;
283 }
284 
285 
286 
287 template <>
289  const Order,
290  const unsigned int,
291  const unsigned int,
292  const Point &)
293 {
294  libmesh_error_msg("Clough-Tocher elements require the real element \nto construct gradient-based degrees of freedom.");
295  return 0.;
296 }
297 
298 
299 
300 template <>
302  const Order order,
303  const unsigned int i,
304  const unsigned int j,
305  const Point & p)
306 {
307  libmesh_assert(elem);
308 
309  clough_compute_coefs(elem);
310 
311  const ElemType type = elem->type();
312 
313  const Order totalorder = static_cast<Order>(order + elem->p_level());
314 
315  switch (totalorder)
316  {
317  // 3rd-order C1 cubic element
318  case THIRD:
319  {
320  switch (type)
321  {
322  // C1 functions on the C1 cubic edge
323  case EDGE2:
324  case EDGE3:
325  {
326  switch (i)
327  {
328  case 0:
329  return clough_raw_shape_deriv(0, j, p);
330  case 1:
331  return clough_raw_shape_deriv(1, j, p);
332  case 2:
333  return d1xd1x * clough_raw_shape_deriv(2, j, p);
334  case 3:
335  return d2xd2x * clough_raw_shape_deriv(3, j, p);
336  default:
337  libmesh_error_msg("Invalid shape function index i = " << i);
338  }
339  }
340  default:
341  libmesh_error_msg("ERROR: Unsupported element type = " << type);
342  }
343  }
344  // by default throw an error
345  default:
346  libmesh_error_msg("ERROR: Unsupported polynomial order = " << totalorder);
347  }
348 
349  libmesh_error_msg("We'll never get here!");
350  return 0.;
351 }
352 
353 
354 
355 template <>
357  const Order order,
358  const unsigned int i,
359  const unsigned int j,
360  const Point & p)
361 {
362  libmesh_assert(elem);
363 
364  clough_compute_coefs(elem);
365 
366  const ElemType type = elem->type();
367 
368  const Order totalorder = static_cast<Order>(order + elem->p_level());
369 
370  switch (totalorder)
371  {
372  // 3rd-order C1 cubic element
373  case THIRD:
374  {
375  switch (type)
376  {
377  // C1 functions on the C1 cubic edge
378  case EDGE2:
379  case EDGE3:
380  {
381  switch (i)
382  {
383  case 0:
384  return clough_raw_shape_second_deriv(0, j, p);
385  case 1:
386  return clough_raw_shape_second_deriv(1, j, p);
387  case 2:
388  return d1xd1x * clough_raw_shape_second_deriv(2, j, p);
389  case 3:
390  return d2xd2x * clough_raw_shape_second_deriv(3, j, p);
391  default:
392  libmesh_error_msg("Invalid shape function index i = " << i);
393  }
394  }
395  default:
396  libmesh_error_msg("ERROR: Unsupported element type = " << type);
397  }
398  }
399  // by default throw an error
400  default:
401  libmesh_error_msg("ERROR: Unsupported polynomial order = " << totalorder);
402  }
403 
404  libmesh_error_msg("We'll never get here!");
405  return 0.;
406 }
407 
408 } // namespace libMesh
unsigned int n_threads()
Definition: libmesh_base.h:125
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
virtual ElemType type() const =0
unsigned int p_level() const
Definition: elem.h:2422
ElemType
Defines an enum for geometric element types.
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
const class libmesh_nullptr_t libmesh_nullptr
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
The libMesh namespace provides an interface to certain functionality in the library.
libmesh_assert(j)
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
Definition: dof_object.h:324
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual unsigned int n_shape_functions() const libmesh_override
Definition: fe.C:36
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
dof_id_type id() const
Definition: dof_object.h:632
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
uint8_t dof_id_type
Definition: id_types.h:64