libMesh
fe_clough_shape_2D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2019 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 #include "libmesh/fe_interface.h"
25 
26 
27 // Anonymous namespace for persistent variables.
28 // This allows us to cache the global-to-local mapping transformation
29 // FIXME: This should also screw up multithreading royally
30 namespace
31 {
32 using namespace libMesh;
33 
34 // Keep track of which element was most recently used to generate
35 // cached data
36 static dof_id_type old_elem_id = DofObject::invalid_id;
37 static const Elem * old_elem_ptr = nullptr;
38 
39 // Coefficient naming: d(1)d(2n) is the coefficient of the
40 // global shape function corresponding to value 1 in terms of the
41 // local shape function corresponding to normal derivative 2
42 static Real d1d2n, d1d3n, d2d3n, d2d1n, d3d1n, d3d2n;
43 static Real d1xd1x, d1xd1y, d1xd2n, d1xd3n;
44 static Real d1yd1x, d1yd1y, d1yd2n, d1yd3n;
45 static Real d2xd2x, d2xd2y, d2xd3n, d2xd1n;
46 static Real d2yd2x, d2yd2y, d2yd3n, d2yd1n;
47 static Real d3xd3x, d3xd3y, d3xd1n, d3xd2n;
48 static Real d3yd3x, d3yd3y, d3yd1n, d3yd2n;
49 static Real d1nd1n, d2nd2n, d3nd3n;
50 // Normal vector naming: N01x is the x component of the
51 // unit vector at point 0 normal to (possibly curved) side 01
52 static Real N01x, N01y, N10x, N10y;
53 static Real N02x, N02y, N20x, N20y;
54 static Real N21x, N21y, N12x, N12y;
55 
56 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
57 
58 Real clough_raw_shape_second_deriv(const unsigned int basis_num,
59  const unsigned int deriv_type,
60  const Point & p);
61 #endif
62 
63 Real clough_raw_shape_deriv(const unsigned int basis_num,
64  const unsigned int deriv_type,
65  const Point & p);
66 Real clough_raw_shape(const unsigned int basis_num,
67  const Point & p);
68 unsigned char subtriangle_lookup(const Point & p);
69 
70 
71 // Compute the static coefficients for an element
72 void clough_compute_coefs(const Elem * elem)
73 {
74  // Using static globals for old_elem_id, etc. will fail
75  // horribly with more than one thread.
76  libmesh_assert_equal_to (libMesh::n_threads(), 1);
77 
78  // Coefficients are cached from old elements; we rely on that cache
79  // except in dbg mode
80 #ifndef DEBUG
81  if (elem->id() == old_elem_id &&
82  elem == old_elem_ptr)
83  return;
84 #endif
85 
86  const FEFamily mapping_family = FEMap::map_fe_type(*elem);
87  const Order mapping_order (elem->default_order());
88  const ElemType mapping_elem_type (elem->type());
89 
90  const FEType map_fe_type(mapping_order, mapping_family);
91 
92  const int n_mapping_shape_functions =
93  FEInterface::n_shape_functions(2, map_fe_type, mapping_elem_type);
94 
95  // Degrees of freedom are at vertices and edge midpoints
96  std::vector<Point> dofpt;
97  dofpt.push_back(Point(0,0));
98  dofpt.push_back(Point(1,0));
99  dofpt.push_back(Point(0,1));
100  dofpt.push_back(Point(1/2.,1/2.));
101  dofpt.push_back(Point(0,1/2.));
102  dofpt.push_back(Point(1/2.,0));
103 
104  // Mapping functions - first derivatives at each dofpt
105  std::vector<Real> dxdxi(6), dxdeta(6), dydxi(6), dydeta(6);
106  std::vector<Real> dxidx(6), detadx(6), dxidy(6), detady(6);
107 
108  FEInterface::shape_deriv_ptr shape_deriv_ptr =
109  FEInterface::shape_deriv_function(2, map_fe_type);
110 
111  for (int p = 0; p != 6; ++p)
112  {
113  // libMesh::err << p << ' ' << dofpt[p];
114  for (int i = 0; i != n_mapping_shape_functions; ++i)
115  {
116  const Real ddxi = shape_deriv_ptr
117  (elem, mapping_order, i, 0, dofpt[p], false);
118  const Real ddeta = shape_deriv_ptr
119  (elem, mapping_order, i, 1, dofpt[p], false);
120 
121  // libMesh::err << ddxi << ' ';
122  // libMesh::err << ddeta << std::endl;
123 
124  dxdxi[p] += elem->point(i)(0) * ddxi;
125  dydxi[p] += elem->point(i)(1) * ddxi;
126  dxdeta[p] += elem->point(i)(0) * ddeta;
127  dydeta[p] += elem->point(i)(1) * ddeta;
128  }
129 
130  // for (int i = 0; i != 12; ++i)
131  // libMesh::err << i << ' ' << clough_raw_shape(i, dofpt[p]) << std::endl;
132 
133  // libMesh::err << elem->point(p)(0) << ' ';
134  // libMesh::err << elem->point(p)(1) << ' ';
135  // libMesh::err << dxdxi[p] << ' ';
136  // libMesh::err << dydxi[p] << ' ';
137  // libMesh::err << dxdeta[p] << ' ';
138  // libMesh::err << dydeta[p] << std::endl << std::endl;
139 
140  const Real inv_jac = 1. / (dxdxi[p]*dydeta[p] -
141  dxdeta[p]*dydxi[p]);
142  dxidx[p] = dydeta[p] * inv_jac;
143  dxidy[p] = - dxdeta[p] * inv_jac;
144  detadx[p] = - dydxi[p] * inv_jac;
145  detady[p] = dxdxi[p] * inv_jac;
146  }
147 
148  // Calculate midpoint normal vectors
149  Real N1x = dydeta[3] - dydxi[3];
150  Real N1y = dxdxi[3] - dxdeta[3];
151  Real Nlength = std::sqrt(static_cast<Real>(N1x*N1x + N1y*N1y));
152  N1x /= Nlength; N1y /= Nlength;
153 
154  Real N2x = - dydeta[4];
155  Real N2y = dxdeta[4];
156  Nlength = std::sqrt(static_cast<Real>(N2x*N2x + N2y*N2y));
157  N2x /= Nlength; N2y /= Nlength;
158 
159  Real N3x = dydxi[5];
160  Real N3y = - dxdxi[5];
161  Nlength = std::sqrt(static_cast<Real>(N3x*N3x + N3y*N3y));
162  N3x /= Nlength; N3y /= Nlength;
163 
164  // Calculate corner normal vectors (used for reduced element)
165  N01x = dydxi[0];
166  N01y = - dxdxi[0];
167  Nlength = std::sqrt(static_cast<Real>(N01x*N01x + N01y*N01y));
168  N01x /= Nlength; N01y /= Nlength;
169 
170  N10x = dydxi[1];
171  N10y = - dxdxi[1];
172  Nlength = std::sqrt(static_cast<Real>(N10x*N10x + N10y*N10y));
173  N10x /= Nlength; N10y /= Nlength;
174 
175  N02x = - dydeta[0];
176  N02y = dxdeta[0];
177  Nlength = std::sqrt(static_cast<Real>(N02x*N02x + N02y*N02y));
178  N02x /= Nlength; N02y /= Nlength;
179 
180  N20x = - dydeta[2];
181  N20y = dxdeta[2];
182  Nlength = std::sqrt(static_cast<Real>(N20x*N20x + N20y*N20y));
183  N20x /= Nlength; N20y /= Nlength;
184 
185  N12x = dydeta[1] - dydxi[1];
186  N12y = dxdxi[1] - dxdeta[1];
187  Nlength = std::sqrt(static_cast<Real>(N12x*N12x + N12y*N12y));
188  N12x /= Nlength; N12y /= Nlength;
189 
190  N21x = dydeta[1] - dydxi[1];
191  N21y = dxdxi[1] - dxdeta[1];
192  Nlength = std::sqrt(static_cast<Real>(N21x*N21x + N21y*N21y));
193  N21x /= Nlength; N21y /= Nlength;
194 
195  // for (int i=0; i != 6; ++i) {
196  // libMesh::err << elem->node_id(i) << ' ';
197  // }
198  // libMesh::err << std::endl;
199 
200  // for (int i=0; i != 6; ++i) {
201  // libMesh::err << elem->point(i)(0) << ' ';
202  // libMesh::err << elem->point(i)(1) << ' ';
203  // }
204  // libMesh::err << std::endl;
205 
206 
207  // give normal vectors a globally unique orientation
208 
209  if (elem->point(2) < elem->point(1))
210  {
211  // libMesh::err << "Flipping nodes " << elem->node_id(2);
212  // libMesh::err << " and " << elem->node_id(1);
213  // libMesh::err << " around node " << elem->node_id(4);
214  // libMesh::err << std::endl;
215  N1x = -N1x; N1y = -N1y;
216  N12x = -N12x; N12y = -N12y;
217  N21x = -N21x; N21y = -N21y;
218  }
219  else
220  {
221  // libMesh::err << "Not flipping nodes " << elem->node_id(2);
222  // libMesh::err << " and " << elem->node_id(1);
223  // libMesh::err << " around node " << elem->node_id(4);
224  // libMesh::err << std::endl;
225  }
226  if (elem->point(0) < elem->point(2))
227  {
228  // libMesh::err << "Flipping nodes " << elem->node_id(0);
229  // libMesh::err << " and " << elem->node_id(2);
230  // libMesh::err << " around node " << elem->node_id(5);
231  // libMesh::err << std::endl;
232  // libMesh::err << N2x << ' ' << N2y << std::endl;
233  N2x = -N2x; N2y = -N2y;
234  N02x = -N02x; N02y = -N02y;
235  N20x = -N20x; N20y = -N20y;
236  // libMesh::err << N2x << ' ' << N2y << std::endl;
237  }
238  else
239  {
240  // libMesh::err << "Not flipping nodes " << elem->node_id(0);
241  // libMesh::err << " and " << elem->node_id(2);
242  // libMesh::err << " around node " << elem->node_id(5);
243  // libMesh::err << std::endl;
244  }
245  if (elem->point(1) < elem->point(0))
246  {
247  // libMesh::err << "Flipping nodes " << elem->node_id(1);
248  // libMesh::err << " and " << elem->node_id(0);
249  // libMesh::err << " around node " << elem->node_id(3);
250  // libMesh::err << std::endl;
251  N3x = -N3x;
252  N3y = -N3y;
253  N01x = -N01x; N01y = -N01y;
254  N10x = -N10x; N10y = -N10y;
255  }
256  else
257  {
258  // libMesh::err << "Not flipping nodes " << elem->node_id(1);
259  // libMesh::err << " and " << elem->node_id(0);
260  // libMesh::err << " around node " << elem->node_id(3);
261  // libMesh::err << std::endl;
262  }
263 
264  // libMesh::err << N2x << ' ' << N2y << std::endl;
265 
266  // Cache basis function gradients
267  // FIXME: the raw_shape calls shouldn't be done on every element!
268  // FIXME: I should probably be looping, too...
269  // Gradient naming: d(1)d(2n)d(xi) is the xi component of the
270  // gradient of the
271  // local basis function corresponding to value 1 at the node
272  // corresponding to normal vector 2
273 
274  Real d1d2ndxi = clough_raw_shape_deriv(0, 0, dofpt[4]);
275  Real d1d2ndeta = clough_raw_shape_deriv(0, 1, dofpt[4]);
276  Real d1d2ndx = d1d2ndxi * dxidx[4] + d1d2ndeta * detadx[4];
277  Real d1d2ndy = d1d2ndxi * dxidy[4] + d1d2ndeta * detady[4];
278  Real d1d3ndxi = clough_raw_shape_deriv(0, 0, dofpt[5]);
279  Real d1d3ndeta = clough_raw_shape_deriv(0, 1, dofpt[5]);
280  Real d1d3ndx = d1d3ndxi * dxidx[5] + d1d3ndeta * detadx[5];
281  Real d1d3ndy = d1d3ndxi * dxidy[5] + d1d3ndeta * detady[5];
282  Real d2d3ndxi = clough_raw_shape_deriv(1, 0, dofpt[5]);
283  Real d2d3ndeta = clough_raw_shape_deriv(1, 1, dofpt[5]);
284  Real d2d3ndx = d2d3ndxi * dxidx[5] + d2d3ndeta * detadx[5];
285  Real d2d3ndy = d2d3ndxi * dxidy[5] + d2d3ndeta * detady[5];
286  Real d2d1ndxi = clough_raw_shape_deriv(1, 0, dofpt[3]);
287  Real d2d1ndeta = clough_raw_shape_deriv(1, 1, dofpt[3]);
288  Real d2d1ndx = d2d1ndxi * dxidx[3] + d2d1ndeta * detadx[3];
289  Real d2d1ndy = d2d1ndxi * dxidy[3] + d2d1ndeta * detady[3];
290  Real d3d1ndxi = clough_raw_shape_deriv(2, 0, dofpt[3]);
291  Real d3d1ndeta = clough_raw_shape_deriv(2, 1, dofpt[3]);
292  Real d3d1ndx = d3d1ndxi * dxidx[3] + d3d1ndeta * detadx[3];
293  Real d3d1ndy = d3d1ndxi * dxidy[3] + d3d1ndeta * detady[3];
294  Real d3d2ndxi = clough_raw_shape_deriv(2, 0, dofpt[4]);
295  Real d3d2ndeta = clough_raw_shape_deriv(2, 1, dofpt[4]);
296  Real d3d2ndx = d3d2ndxi * dxidx[4] + d3d2ndeta * detadx[4];
297  Real d3d2ndy = d3d2ndxi * dxidy[4] + d3d2ndeta * detady[4];
298  Real d1xd2ndxi = clough_raw_shape_deriv(3, 0, dofpt[4]);
299  Real d1xd2ndeta = clough_raw_shape_deriv(3, 1, dofpt[4]);
300  Real d1xd2ndx = d1xd2ndxi * dxidx[4] + d1xd2ndeta * detadx[4];
301  Real d1xd2ndy = d1xd2ndxi * dxidy[4] + d1xd2ndeta * detady[4];
302  Real d1xd3ndxi = clough_raw_shape_deriv(3, 0, dofpt[5]);
303  Real d1xd3ndeta = clough_raw_shape_deriv(3, 1, dofpt[5]);
304  Real d1xd3ndx = d1xd3ndxi * dxidx[5] + d1xd3ndeta * detadx[5];
305  Real d1xd3ndy = d1xd3ndxi * dxidy[5] + d1xd3ndeta * detady[5];
306  Real d1yd2ndxi = clough_raw_shape_deriv(4, 0, dofpt[4]);
307  Real d1yd2ndeta = clough_raw_shape_deriv(4, 1, dofpt[4]);
308  Real d1yd2ndx = d1yd2ndxi * dxidx[4] + d1yd2ndeta * detadx[4];
309  Real d1yd2ndy = d1yd2ndxi * dxidy[4] + d1yd2ndeta * detady[4];
310  Real d1yd3ndxi = clough_raw_shape_deriv(4, 0, dofpt[5]);
311  Real d1yd3ndeta = clough_raw_shape_deriv(4, 1, dofpt[5]);
312  Real d1yd3ndx = d1yd3ndxi * dxidx[5] + d1yd3ndeta * detadx[5];
313  Real d1yd3ndy = d1yd3ndxi * dxidy[5] + d1yd3ndeta * detady[5];
314  Real d2xd3ndxi = clough_raw_shape_deriv(5, 0, dofpt[5]);
315  Real d2xd3ndeta = clough_raw_shape_deriv(5, 1, dofpt[5]);
316  Real d2xd3ndx = d2xd3ndxi * dxidx[5] + d2xd3ndeta * detadx[5];
317  Real d2xd3ndy = d2xd3ndxi * dxidy[5] + d2xd3ndeta * detady[5];
318  Real d2xd1ndxi = clough_raw_shape_deriv(5, 0, dofpt[3]);
319  Real d2xd1ndeta = clough_raw_shape_deriv(5, 1, dofpt[3]);
320  Real d2xd1ndx = d2xd1ndxi * dxidx[3] + d2xd1ndeta * detadx[3];
321  Real d2xd1ndy = d2xd1ndxi * dxidy[3] + d2xd1ndeta * detady[3];
322  Real d2yd3ndxi = clough_raw_shape_deriv(6, 0, dofpt[5]);
323  Real d2yd3ndeta = clough_raw_shape_deriv(6, 1, dofpt[5]);
324  Real d2yd3ndx = d2yd3ndxi * dxidx[5] + d2yd3ndeta * detadx[5];
325  Real d2yd3ndy = d2yd3ndxi * dxidy[5] + d2yd3ndeta * detady[5];
326  Real d2yd1ndxi = clough_raw_shape_deriv(6, 0, dofpt[3]);
327  Real d2yd1ndeta = clough_raw_shape_deriv(6, 1, dofpt[3]);
328  Real d2yd1ndx = d2yd1ndxi * dxidx[3] + d2yd1ndeta * detadx[3];
329  Real d2yd1ndy = d2yd1ndxi * dxidy[3] + d2yd1ndeta * detady[3];
330  Real d3xd1ndxi = clough_raw_shape_deriv(7, 0, dofpt[3]);
331  Real d3xd1ndeta = clough_raw_shape_deriv(7, 1, dofpt[3]);
332  Real d3xd1ndx = d3xd1ndxi * dxidx[3] + d3xd1ndeta * detadx[3];
333  Real d3xd1ndy = d3xd1ndxi * dxidy[3] + d3xd1ndeta * detady[3];
334  Real d3xd2ndxi = clough_raw_shape_deriv(7, 0, dofpt[4]);
335  Real d3xd2ndeta = clough_raw_shape_deriv(7, 1, dofpt[4]);
336  Real d3xd2ndx = d3xd2ndxi * dxidx[4] + d3xd2ndeta * detadx[4];
337  Real d3xd2ndy = d3xd2ndxi * dxidy[4] + d3xd2ndeta * detady[4];
338  Real d3yd1ndxi = clough_raw_shape_deriv(8, 0, dofpt[3]);
339  Real d3yd1ndeta = clough_raw_shape_deriv(8, 1, dofpt[3]);
340  Real d3yd1ndx = d3yd1ndxi * dxidx[3] + d3yd1ndeta * detadx[3];
341  Real d3yd1ndy = d3yd1ndxi * dxidy[3] + d3yd1ndeta * detady[3];
342  Real d3yd2ndxi = clough_raw_shape_deriv(8, 0, dofpt[4]);
343  Real d3yd2ndeta = clough_raw_shape_deriv(8, 1, dofpt[4]);
344  Real d3yd2ndx = d3yd2ndxi * dxidx[4] + d3yd2ndeta * detadx[4];
345  Real d3yd2ndy = d3yd2ndxi * dxidy[4] + d3yd2ndeta * detady[4];
346  Real d1nd1ndxi = clough_raw_shape_deriv(9, 0, dofpt[3]);
347  Real d1nd1ndeta = clough_raw_shape_deriv(9, 1, dofpt[3]);
348  Real d1nd1ndx = d1nd1ndxi * dxidx[3] + d1nd1ndeta * detadx[3];
349  Real d1nd1ndy = d1nd1ndxi * dxidy[3] + d1nd1ndeta * detady[3];
350  Real d2nd2ndxi = clough_raw_shape_deriv(10, 0, dofpt[4]);
351  Real d2nd2ndeta = clough_raw_shape_deriv(10, 1, dofpt[4]);
352  Real d2nd2ndx = d2nd2ndxi * dxidx[4] + d2nd2ndeta * detadx[4];
353  Real d2nd2ndy = d2nd2ndxi * dxidy[4] + d2nd2ndeta * detady[4];
354  Real d3nd3ndxi = clough_raw_shape_deriv(11, 0, dofpt[5]);
355  Real d3nd3ndeta = clough_raw_shape_deriv(11, 1, dofpt[5]);
356  Real d3nd3ndx = d3nd3ndxi * dxidx[3] + d3nd3ndeta * detadx[3];
357  Real d3nd3ndy = d3nd3ndxi * dxidy[3] + d3nd3ndeta * detady[3];
358 
359  Real d1xd1dxi = clough_raw_shape_deriv(3, 0, dofpt[0]);
360  Real d1xd1deta = clough_raw_shape_deriv(3, 1, dofpt[0]);
361  Real d1xd1dx = d1xd1dxi * dxidx[0] + d1xd1deta * detadx[0];
362  Real d1xd1dy = d1xd1dxi * dxidy[0] + d1xd1deta * detady[0];
363  Real d1yd1dxi = clough_raw_shape_deriv(4, 0, dofpt[0]);
364  Real d1yd1deta = clough_raw_shape_deriv(4, 1, dofpt[0]);
365  Real d1yd1dx = d1yd1dxi * dxidx[0] + d1yd1deta * detadx[0];
366  Real d1yd1dy = d1yd1dxi * dxidy[0] + d1yd1deta * detady[0];
367  Real d2xd2dxi = clough_raw_shape_deriv(5, 0, dofpt[1]);
368  Real d2xd2deta = clough_raw_shape_deriv(5, 1, dofpt[1]);
369  Real d2xd2dx = d2xd2dxi * dxidx[1] + d2xd2deta * detadx[1];
370  Real d2xd2dy = d2xd2dxi * dxidy[1] + d2xd2deta * detady[1];
371 
372  // libMesh::err << dofpt[4](0) << ' ';
373  // libMesh::err << dofpt[4](1) << ' ';
374  // libMesh::err << (int)subtriangle_lookup(dofpt[5]) << ' ';
375  // libMesh::err << dxdxi[4] << ' ';
376  // libMesh::err << dxdeta[4] << ' ';
377  // libMesh::err << dydxi[4] << ' ';
378  // libMesh::err << dydeta[4] << ' ';
379  // libMesh::err << dxidx[4] << ' ';
380  // libMesh::err << dxidy[4] << ' ';
381  // libMesh::err << detadx[4] << ' ';
382  // libMesh::err << detady[4] << ' ';
383  // libMesh::err << N1x << ' ';
384  // libMesh::err << N1y << ' ';
385  // libMesh::err << d2yd1ndxi << ' ';
386  // libMesh::err << d2yd1ndeta << ' ';
387  // libMesh::err << d2yd1ndx << ' ';
388  // libMesh::err << d2yd1ndy << std::endl;
389 
390  Real d2yd2dxi = clough_raw_shape_deriv(6, 0, dofpt[1]);
391  Real d2yd2deta = clough_raw_shape_deriv(6, 1, dofpt[1]);
392  Real d2yd2dx = d2yd2dxi * dxidx[1] + d2yd2deta * detadx[1];
393  Real d2yd2dy = d2yd2dxi * dxidy[1] + d2yd2deta * detady[1];
394  Real d3xd3dxi = clough_raw_shape_deriv(7, 0, dofpt[2]);
395  Real d3xd3deta = clough_raw_shape_deriv(7, 1, dofpt[2]);
396  Real d3xd3dx = d3xd3dxi * dxidx[1] + d3xd3deta * detadx[1];
397  Real d3xd3dy = d3xd3dxi * dxidy[1] + d3xd3deta * detady[1];
398  Real d3yd3dxi = clough_raw_shape_deriv(8, 0, dofpt[2]);
399  Real d3yd3deta = clough_raw_shape_deriv(8, 1, dofpt[2]);
400  Real d3yd3dx = d3yd3dxi * dxidx[1] + d3yd3deta * detadx[1];
401  Real d3yd3dy = d3yd3dxi * dxidy[1] + d3yd3deta * detady[1];
402 
403  // Calculate normal derivatives
404 
405  Real d1nd1ndn = d1nd1ndx * N1x + d1nd1ndy * N1y;
406  Real d1xd2ndn = d1xd2ndx * N2x + d1xd2ndy * N2y;
407  Real d1xd3ndn = d1xd3ndx * N3x + d1xd3ndy * N3y;
408  Real d1yd2ndn = d1yd2ndx * N2x + d1yd2ndy * N2y;
409  Real d1yd3ndn = d1yd3ndx * N3x + d1yd3ndy * N3y;
410 
411  Real d2nd2ndn = d2nd2ndx * N2x + d2nd2ndy * N2y;
412  Real d2xd3ndn = d2xd3ndx * N3x + d2xd3ndy * N3y;
413  Real d2xd1ndn = d2xd1ndx * N1x + d2xd1ndy * N1y;
414  Real d2yd3ndn = d2yd3ndx * N3x + d2yd3ndy * N3y;
415  Real d2yd1ndn = d2yd1ndx * N1x + d2yd1ndy * N1y;
416 
417  Real d3nd3ndn = d3nd3ndx * N3x + d3nd3ndy * N3y;
418  Real d3xd1ndn = d3xd1ndx * N1x + d3xd1ndy * N1y;
419  Real d3xd2ndn = d3xd2ndx * N2x + d3xd2ndy * N2y;
420  Real d3yd1ndn = d3yd1ndx * N1x + d3yd1ndy * N1y;
421  Real d3yd2ndn = d3yd2ndx * N2x + d3yd2ndy * N2y;
422 
423  // Calculate midpoint scaling factors
424 
425  d1nd1n = 1. / d1nd1ndn;
426  d2nd2n = 1. / d2nd2ndn;
427  d3nd3n = 1. / d3nd3ndn;
428 
429  // Calculate midpoint derivative adjustments to nodal value
430  // interpolant functions
431 
432 #ifdef DEBUG
433  // The cached factors should equal our calculations if we're still
434  // operating on the cached element
435  if (elem->id() == old_elem_id &&
436  elem == old_elem_ptr)
437  {
438  libmesh_assert_equal_to(d1d2n, -(d1d2ndx * N2x + d1d2ndy * N2y) / d2nd2ndn);
439  libmesh_assert_equal_to(d1d3n, -(d1d3ndx * N3x + d1d3ndy * N3y) / d3nd3ndn);
440  libmesh_assert_equal_to(d2d3n, -(d2d3ndx * N3x + d2d3ndy * N3y) / d3nd3ndn);
441  libmesh_assert_equal_to(d2d1n, -(d2d1ndx * N1x + d2d1ndy * N1y) / d1nd1ndn);
442  libmesh_assert_equal_to(d3d1n, -(d3d1ndx * N1x + d3d1ndy * N1y) / d1nd1ndn);
443  libmesh_assert_equal_to(d3d2n, -(d3d2ndx * N2x + d3d2ndy * N2y) / d2nd2ndn);
444  libmesh_assert_equal_to(d1xd1x, 1. / (d1xd1dx - d1xd1dy * d1yd1dx / d1yd1dy));
445  libmesh_assert_equal_to(d1xd1y, 1. / (d1yd1dx - d1xd1dx * d1yd1dy / d1xd1dy));
446  libmesh_assert_equal_to(d1yd1y, 1. / (d1yd1dy - d1yd1dx * d1xd1dy / d1xd1dx));
447  libmesh_assert_equal_to(d1yd1x, 1. / (d1xd1dy - d1yd1dy * d1xd1dx / d1yd1dx));
448  libmesh_assert_equal_to(d2xd2x, 1. / (d2xd2dx - d2xd2dy * d2yd2dx / d2yd2dy));
449  libmesh_assert_equal_to(d2xd2y, 1. / (d2yd2dx - d2xd2dx * d2yd2dy / d2xd2dy));
450  libmesh_assert_equal_to(d2yd2y, 1. / (d2yd2dy - d2yd2dx * d2xd2dy / d2xd2dx));
451  libmesh_assert_equal_to(d2yd2x, 1. / (d2xd2dy - d2yd2dy * d2xd2dx / d2yd2dx));
452  libmesh_assert_equal_to(d3xd3x, 1. / (d3xd3dx - d3xd3dy * d3yd3dx / d3yd3dy));
453  libmesh_assert_equal_to(d3xd3y, 1. / (d3yd3dx - d3xd3dx * d3yd3dy / d3xd3dy));
454  libmesh_assert_equal_to(d3yd3y, 1. / (d3yd3dy - d3yd3dx * d3xd3dy / d3xd3dx));
455  libmesh_assert_equal_to(d3yd3x, 1. / (d3xd3dy - d3yd3dy * d3xd3dx / d3yd3dx));
456  libmesh_assert_equal_to(d1xd2n, -(d1xd1x * d1xd2ndn + d1xd1y * d1yd2ndn) / d2nd2ndn);
457  libmesh_assert_equal_to(d1yd2n, -(d1yd1y * d1yd2ndn + d1yd1x * d1xd2ndn) / d2nd2ndn);
458  libmesh_assert_equal_to(d1xd3n, -(d1xd1x * d1xd3ndn + d1xd1y * d1yd3ndn) / d3nd3ndn);
459  libmesh_assert_equal_to(d1yd3n, -(d1yd1y * d1yd3ndn + d1yd1x * d1xd3ndn) / d3nd3ndn);
460  libmesh_assert_equal_to(d2xd3n, -(d2xd2x * d2xd3ndn + d2xd2y * d2yd3ndn) / d3nd3ndn);
461  libmesh_assert_equal_to(d2yd3n, -(d2yd2y * d2yd3ndn + d2yd2x * d2xd3ndn) / d3nd3ndn);
462  libmesh_assert_equal_to(d2xd1n, -(d2xd2x * d2xd1ndn + d2xd2y * d2yd1ndn) / d1nd1ndn);
463  libmesh_assert_equal_to(d2yd1n, -(d2yd2y * d2yd1ndn + d2yd2x * d2xd1ndn) / d1nd1ndn);
464  libmesh_assert_equal_to(d3xd1n, -(d3xd3x * d3xd1ndn + d3xd3y * d3yd1ndn) / d1nd1ndn);
465  libmesh_assert_equal_to(d3yd1n, -(d3yd3y * d3yd1ndn + d3yd3x * d3xd1ndn) / d1nd1ndn);
466  libmesh_assert_equal_to(d3xd2n, -(d3xd3x * d3xd2ndn + d3xd3y * d3yd2ndn) / d2nd2ndn);
467  libmesh_assert_equal_to(d3yd2n, -(d3yd3y * d3yd2ndn + d3yd3x * d3xd2ndn) / d2nd2ndn);
468  }
469 #endif
470 
471  d1d2n = -(d1d2ndx * N2x + d1d2ndy * N2y) / d2nd2ndn;
472  d1d3n = -(d1d3ndx * N3x + d1d3ndy * N3y) / d3nd3ndn;
473  d2d3n = -(d2d3ndx * N3x + d2d3ndy * N3y) / d3nd3ndn;
474  d2d1n = -(d2d1ndx * N1x + d2d1ndy * N1y) / d1nd1ndn;
475  d3d1n = -(d3d1ndx * N1x + d3d1ndy * N1y) / d1nd1ndn;
476  d3d2n = -(d3d2ndx * N2x + d3d2ndy * N2y) / d2nd2ndn;
477 
478  // Calculate nodal derivative scaling factors
479 
480  d1xd1x = 1. / (d1xd1dx - d1xd1dy * d1yd1dx / d1yd1dy);
481  d1xd1y = 1. / (d1yd1dx - d1xd1dx * d1yd1dy / d1xd1dy);
482  // d1xd1y = - d1xd1x * (d1xd1dy / d1yd1dy);
483  d1yd1y = 1. / (d1yd1dy - d1yd1dx * d1xd1dy / d1xd1dx);
484  d1yd1x = 1. / (d1xd1dy - d1yd1dy * d1xd1dx / d1yd1dx);
485  // d1yd1x = - d1yd1y * (d1yd1dx / d1xd1dx);
486  d2xd2x = 1. / (d2xd2dx - d2xd2dy * d2yd2dx / d2yd2dy);
487  d2xd2y = 1. / (d2yd2dx - d2xd2dx * d2yd2dy / d2xd2dy);
488  // d2xd2y = - d2xd2x * (d2xd2dy / d2yd2dy);
489  d2yd2y = 1. / (d2yd2dy - d2yd2dx * d2xd2dy / d2xd2dx);
490  d2yd2x = 1. / (d2xd2dy - d2yd2dy * d2xd2dx / d2yd2dx);
491  // d2yd2x = - d2yd2y * (d2yd2dx / d2xd2dx);
492  d3xd3x = 1. / (d3xd3dx - d3xd3dy * d3yd3dx / d3yd3dy);
493  d3xd3y = 1. / (d3yd3dx - d3xd3dx * d3yd3dy / d3xd3dy);
494  // d3xd3y = - d3xd3x * (d3xd3dy / d3yd3dy);
495  d3yd3y = 1. / (d3yd3dy - d3yd3dx * d3xd3dy / d3xd3dx);
496  d3yd3x = 1. / (d3xd3dy - d3yd3dy * d3xd3dx / d3yd3dx);
497  // d3yd3x = - d3yd3y * (d3yd3dx / d3xd3dx);
498 
499  // libMesh::err << d1xd1dx << ' ';
500  // libMesh::err << d1xd1dy << ' ';
501  // libMesh::err << d1yd1dx << ' ';
502  // libMesh::err << d1yd1dy << ' ';
503  // libMesh::err << d2xd2dx << ' ';
504  // libMesh::err << d2xd2dy << ' ';
505  // libMesh::err << d2yd2dx << ' ';
506  // libMesh::err << d2yd2dy << ' ';
507  // libMesh::err << d3xd3dx << ' ';
508  // libMesh::err << d3xd3dy << ' ';
509  // libMesh::err << d3yd3dx << ' ';
510  // libMesh::err << d3yd3dy << std::endl;
511 
512  // Calculate midpoint derivative adjustments to nodal derivative
513  // interpolant functions
514 
515  d1xd2n = -(d1xd1x * d1xd2ndn + d1xd1y * d1yd2ndn) / d2nd2ndn;
516  d1yd2n = -(d1yd1y * d1yd2ndn + d1yd1x * d1xd2ndn) / d2nd2ndn;
517  d1xd3n = -(d1xd1x * d1xd3ndn + d1xd1y * d1yd3ndn) / d3nd3ndn;
518  d1yd3n = -(d1yd1y * d1yd3ndn + d1yd1x * d1xd3ndn) / d3nd3ndn;
519  d2xd3n = -(d2xd2x * d2xd3ndn + d2xd2y * d2yd3ndn) / d3nd3ndn;
520  d2yd3n = -(d2yd2y * d2yd3ndn + d2yd2x * d2xd3ndn) / d3nd3ndn;
521  d2xd1n = -(d2xd2x * d2xd1ndn + d2xd2y * d2yd1ndn) / d1nd1ndn;
522  d2yd1n = -(d2yd2y * d2yd1ndn + d2yd2x * d2xd1ndn) / d1nd1ndn;
523  d3xd1n = -(d3xd3x * d3xd1ndn + d3xd3y * d3yd1ndn) / d1nd1ndn;
524  d3yd1n = -(d3yd3y * d3yd1ndn + d3yd3x * d3xd1ndn) / d1nd1ndn;
525  d3xd2n = -(d3xd3x * d3xd2ndn + d3xd3y * d3yd2ndn) / d2nd2ndn;
526  d3yd2n = -(d3yd3y * d3yd2ndn + d3yd3x * d3xd2ndn) / d2nd2ndn;
527 
528  old_elem_id = elem->id();
529  old_elem_ptr = elem;
530 
531  // Cross your fingers
532  // libMesh::err << d1nd1ndn << ' ';
533  // libMesh::err << d2xd1ndn << ' ';
534  // libMesh::err << d2yd1ndn << ' ';
535  // libMesh::err << std::endl;
536 
537  // libMesh::err << "Transform variables: ";
538  // libMesh::err << d1nd1n << ' ';
539  // libMesh::err << d2nd2n << ' ';
540  // libMesh::err << d3nd3n << ' ';
541  // libMesh::err << d1d2n << ' ';
542  // libMesh::err << d1d3n << ' ';
543  // libMesh::err << d2d3n << ' ';
544  // libMesh::err << d2d1n << ' ';
545  // libMesh::err << d3d1n << ' ';
546  // libMesh::err << d3d2n << std::endl;
547  // libMesh::err << d1xd1x << ' ';
548  // libMesh::err << d1xd1y << ' ';
549  // libMesh::err << d1yd1x << ' ';
550  // libMesh::err << d1yd1y << ' ';
551  // libMesh::err << d2xd2x << ' ';
552  // libMesh::err << d2xd2y << ' ';
553  // libMesh::err << d2yd2x << ' ';
554  // libMesh::err << d2yd2y << ' ';
555  // libMesh::err << d3xd3x << ' ';
556  // libMesh::err << d3xd3y << ' ';
557  // libMesh::err << d3yd3x << ' ';
558  // libMesh::err << d3yd3y << std::endl;
559  // libMesh::err << d1xd2n << ' ';
560  // libMesh::err << d1yd2n << ' ';
561  // libMesh::err << d1xd3n << ' ';
562  // libMesh::err << d1yd3n << ' ';
563  // libMesh::err << d2xd3n << ' ';
564  // libMesh::err << d2yd3n << ' ';
565  // libMesh::err << d2xd1n << ' ';
566  // libMesh::err << d2yd1n << ' ';
567  // libMesh::err << d3xd1n << ' ';
568  // libMesh::err << d3yd1n << ' ';
569  // libMesh::err << d3xd2n << ' ';
570  // libMesh::err << d3yd2n << ' ';
571  // libMesh::err << std::endl;
572  // libMesh::err << std::endl;
573 }
574 
575 
576 unsigned char subtriangle_lookup(const Point & p)
577 {
578  if ((p(0) >= p(1)) && (p(0) + 2 * p(1) <= 1))
579  return 0;
580  if ((p(0) <= p(1)) && (p(1) + 2 * p(0) <= 1))
581  return 2;
582  return 1;
583 }
584 
585 
586 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
587 // Return shape function second derivatives on the unit right
588 // triangle
589 Real clough_raw_shape_second_deriv(const unsigned int basis_num,
590  const unsigned int deriv_type,
591  const Point & p)
592 {
593  Real xi = p(0), eta = p(1);
594 
595  switch (deriv_type)
596  {
597 
598  // second derivative in xi-xi direction
599  case 0:
600  switch (basis_num)
601  {
602  case 0:
603  switch (subtriangle_lookup(p))
604  {
605  case 0:
606  return -6 + 12*xi;
607  case 1:
608  return -30 + 42*xi + 42*eta;
609  case 2:
610  return -6 + 18*xi - 6*eta;
611 
612  default:
613  libmesh_error_msg("Invalid subtriangle lookup = " <<
614  subtriangle_lookup(p));
615  }
616  case 1:
617  switch (subtriangle_lookup(p))
618  {
619  case 0:
620  return 6 - 12*xi;
621  case 1:
622  return 18 - 27*xi - 21*eta;
623  case 2:
624  return 6 - 15*xi + 3*eta;
625 
626  default:
627  libmesh_error_msg("Invalid subtriangle lookup = " <<
628  subtriangle_lookup(p));
629  }
630  case 2:
631  switch (subtriangle_lookup(p))
632  {
633  case 0:
634  return 0;
635  case 1:
636  return 12 - 15*xi - 21*eta;
637  case 2:
638  return -3*xi + 3*eta;
639 
640  default:
641  libmesh_error_msg("Invalid subtriangle lookup = " <<
642  subtriangle_lookup(p));
643  }
644  case 3:
645  switch (subtriangle_lookup(p))
646  {
647  case 0:
648  return -4 + 6*xi;
649  case 1:
650  return -9 + 13*xi + 8*eta;
651  case 2:
652  return -1 - 7*xi + 4*eta;
653 
654  default:
655  libmesh_error_msg("Invalid subtriangle lookup = " <<
656  subtriangle_lookup(p));
657  }
658  case 4:
659  switch (subtriangle_lookup(p))
660  {
661  case 0:
662  return 4*eta;
663  case 1:
664  return 1 - 2*xi + 3*eta;
665  case 2:
666  return -3 + 14*xi - eta;
667 
668  default:
669  libmesh_error_msg("Invalid subtriangle lookup = " <<
670  subtriangle_lookup(p));
671  }
672  case 5:
673  switch (subtriangle_lookup(p))
674  {
675  case 0:
676  return -2 + 6*xi;
677  case 1:
678  return -4 + 17./2.*xi + 7./2.*eta;
679  case 2:
680  return -2 + 13./2.*xi - 1./2.*eta;
681 
682  default:
683  libmesh_error_msg("Invalid subtriangle lookup = " <<
684  subtriangle_lookup(p));
685  }
686  case 6:
687  switch (subtriangle_lookup(p))
688  {
689  case 0:
690  return 4*eta;
691  case 1:
692  return 9 - 23./2.*xi - 23./2.*eta;
693  case 2:
694  return -1 + 5./2.*xi + 9./2.*eta;
695 
696  default:
697  libmesh_error_msg("Invalid subtriangle lookup = " <<
698  subtriangle_lookup(p));
699  }
700  case 7:
701  switch (subtriangle_lookup(p))
702  {
703  case 0:
704  return 0;
705  case 1:
706  return 7 - 17./2.*xi - 25./2.*eta;
707  case 2:
708  return 1 - 13./2.*xi + 7./2.*eta;
709 
710  default:
711  libmesh_error_msg("Invalid subtriangle lookup = " <<
712  subtriangle_lookup(p));
713  }
714  case 8:
715  switch (subtriangle_lookup(p))
716  {
717  case 0:
718  return 0;
719  case 1:
720  return -2 + 5./2.*xi + 7./2.*eta;
721  case 2:
722  return 1./2.*xi - 1./2.*eta;
723 
724  default:
725  libmesh_error_msg("Invalid subtriangle lookup = " <<
726  subtriangle_lookup(p));
727  }
728  case 9:
729  switch (subtriangle_lookup(p))
730  {
731  case 0:
732  return 0;
733  case 1:
734  return std::sqrt(2.) * (8 - 10*xi - 14*eta);
735  case 2:
736  return std::sqrt(2.) * (-2*xi + 2*eta);
737 
738  default:
739  libmesh_error_msg("Invalid subtriangle lookup = " <<
740  subtriangle_lookup(p));
741  }
742  case 10:
743  switch (subtriangle_lookup(p))
744  {
745  case 0:
746  return 0;
747  case 1:
748  return -4 + 4*xi + 8*eta;
749  case 2:
750  return -4 + 20*xi - 8*eta;
751 
752  default:
753  libmesh_error_msg("Invalid subtriangle lookup = " <<
754  subtriangle_lookup(p));
755  }
756  case 11:
757  switch (subtriangle_lookup(p))
758  {
759  case 0:
760  return -8*eta;
761  case 1:
762  return -12 + 16*xi + 12*eta;
763  case 2:
764  return 4 - 16*xi - 4*eta;
765 
766  default:
767  libmesh_error_msg("Invalid subtriangle lookup = " <<
768  subtriangle_lookup(p));
769  }
770 
771  default:
772  libmesh_error_msg("Invalid shape function index i = " <<
773  basis_num);
774  }
775 
776  // second derivative in xi-eta direction
777  case 1:
778  switch (basis_num)
779  {
780  case 0:
781  switch (subtriangle_lookup(p))
782  {
783  case 0:
784  return -6*eta;
785  case 1:
786  return -30 +42*xi
787  + 42*eta;
788  case 2:
789  return -6*xi;
790 
791  default:
792  libmesh_error_msg("Invalid subtriangle lookup = " <<
793  subtriangle_lookup(p));
794  }
795  case 1:
796  switch (subtriangle_lookup(p))
797  {
798  case 0:
799  return + 3*eta;
800  case 1:
801  return 15 - 21*xi - 21*eta;
802  case 2:
803  return 3*xi;
804 
805  default:
806  libmesh_error_msg("Invalid subtriangle lookup = " <<
807  subtriangle_lookup(p));
808  }
809  case 2:
810  switch (subtriangle_lookup(p))
811  {
812  case 0:
813  return 3*eta;
814  case 1:
815  return 15 - 21*xi - 21*eta;
816  case 2:
817  return 3*xi;
818 
819  default:
820  libmesh_error_msg("Invalid subtriangle lookup = " <<
821  subtriangle_lookup(p));
822  }
823  case 3:
824  switch (subtriangle_lookup(p))
825  {
826  case 0:
827  return -eta;
828  case 1:
829  return -4 + 8*xi + 3*eta;
830  case 2:
831  return -3 + 4*xi + 4*eta;
832 
833  default:
834  libmesh_error_msg("Invalid subtriangle lookup = " <<
835  subtriangle_lookup(p));
836  }
837  case 4:
838  switch (subtriangle_lookup(p))
839  {
840  case 0:
841  return -3 + 4*xi + 4*eta;
842  case 1:
843  return - 4 + 3*xi + 8*eta;
844  case 2:
845  return -xi;
846 
847  default:
848  libmesh_error_msg("Invalid subtriangle lookup = " <<
849  subtriangle_lookup(p));
850  }
851  case 5:
852  switch (subtriangle_lookup(p))
853  {
854  case 0:
855  return - 1./2.*eta;
856  case 1:
857  return -5./2. + 7./2.*xi + 7./2.*eta;
858  case 2:
859  return - 1./2.*xi;
860 
861  default:
862  libmesh_error_msg("Invalid subtriangle lookup = " <<
863  subtriangle_lookup(p));
864  }
865  case 6:
866  switch (subtriangle_lookup(p))
867  {
868  case 0:
869  return -1 + 4*xi + 7./2.*eta;
870  case 1:
871  return 19./2. - 23./2.*xi - 25./2.*eta;
872  case 2:
873  return 9./2.*xi;
874 
875  default:
876  libmesh_error_msg("Invalid subtriangle lookup = " <<
877  subtriangle_lookup(p));
878  }
879  case 7:
880  switch (subtriangle_lookup(p))
881  {
882  case 0:
883  return 9./2.*eta;
884  case 1:
885  return 19./2. - 25./2.*xi - 23./2.*eta;
886  case 2:
887  return -1 + 7./2.*xi + 4*eta;
888 
889  default:
890  libmesh_error_msg("Invalid subtriangle lookup = " <<
891  subtriangle_lookup(p));
892  }
893  case 8:
894  switch (subtriangle_lookup(p))
895  {
896  case 0:
897  return -1./2.*eta;
898  case 1:
899  return -5./2. + 7./2.*xi + 7./2.*eta;
900  case 2:
901  return -1./2.*xi;
902 
903  default:
904  libmesh_error_msg("Invalid subtriangle lookup = " <<
905  subtriangle_lookup(p));
906  }
907  case 9:
908  switch (subtriangle_lookup(p))
909  {
910  case 0:
911  return std::sqrt(2.) * (2*eta);
912  case 1:
913  return std::sqrt(2.) * (10 - 14*xi - 14*eta);
914  case 2:
915  return std::sqrt(2.) * (2*xi);
916 
917  default:
918  libmesh_error_msg("Invalid subtriangle lookup = " <<
919  subtriangle_lookup(p));
920  }
921  case 10:
922  switch (subtriangle_lookup(p))
923  {
924  case 0:
925  return -4*eta;
926  case 1:
927  return - 8 + 8*xi + 12*eta;
928  case 2:
929  return 4 - 8*xi - 8*eta;
930 
931  default:
932  libmesh_error_msg("Invalid subtriangle lookup = " <<
933  subtriangle_lookup(p));
934  }
935  case 11:
936  switch (subtriangle_lookup(p))
937  {
938  case 0:
939  return 4 - 8*xi - 8*eta;
940  case 1:
941  return -8 + 12*xi + 8*eta;
942  case 2:
943  return -4*xi;
944 
945  default:
946  libmesh_error_msg("Invalid subtriangle lookup = " <<
947  subtriangle_lookup(p));
948  }
949 
950  default:
951  libmesh_error_msg("Invalid shape function index i = " <<
952  basis_num);
953  }
954 
955  // second derivative in eta-eta direction
956  case 2:
957  switch (basis_num)
958  {
959  case 0:
960  switch (subtriangle_lookup(p))
961  {
962  case 0:
963  return -6 - 6*xi + 18*eta;
964  case 1:
965  return -30 + 42*xi + 42*eta;
966  case 2:
967  return -6 + 12*eta;
968 
969  default:
970  libmesh_error_msg("Invalid subtriangle lookup = " <<
971  subtriangle_lookup(p));
972  }
973  case 1:
974  switch (subtriangle_lookup(p))
975  {
976  case 0:
977  return 3*xi - 3*eta;
978  case 1:
979  return 12 - 21*xi - 15*eta;
980  case 2:
981  return 0;
982 
983  default:
984  libmesh_error_msg("Invalid subtriangle lookup = " <<
985  subtriangle_lookup(p));
986  }
987  case 2:
988  switch (subtriangle_lookup(p))
989  {
990  case 0:
991  return 6 + 3*xi - 15*eta;
992  case 1:
993  return 18 - 21.*xi - 27*eta;
994  case 2:
995  return 6 - 12*eta;
996 
997  default:
998  libmesh_error_msg("Invalid subtriangle lookup = " <<
999  subtriangle_lookup(p));
1000  }
1001  case 3:
1002  switch (subtriangle_lookup(p))
1003  {
1004  case 0:
1005  return -3 - xi + 14*eta;
1006  case 1:
1007  return 1 + 3*xi - 2*eta;
1008  case 2:
1009  return 4*xi;
1010 
1011  default:
1012  libmesh_error_msg("Invalid subtriangle lookup = " <<
1013  subtriangle_lookup(p));
1014  }
1015  case 4:
1016  switch (subtriangle_lookup(p))
1017  {
1018  case 0:
1019  return -1 + 4*xi - 7*eta;
1020  case 1:
1021  return -9 + 8*xi + 13*eta;
1022  case 2:
1023  return -4 + 6*eta;
1024 
1025  default:
1026  libmesh_error_msg("Invalid subtriangle lookup = " <<
1027  subtriangle_lookup(p));
1028  }
1029  case 5:
1030  switch (subtriangle_lookup(p))
1031  {
1032  case 0:
1033  return - 1./2.*xi + 1./2.*eta;
1034  case 1:
1035  return -2 + 7./2.*xi + 5./2.*eta;
1036  case 2:
1037  return 0;
1038 
1039  default:
1040  libmesh_error_msg("Invalid subtriangle lookup = " <<
1041  subtriangle_lookup(p));
1042  }
1043  case 6:
1044  switch (subtriangle_lookup(p))
1045  {
1046  case 0:
1047  return 1 + 7./2.*xi - 13./2.*eta;
1048  case 1:
1049  return 7 - 25./2.*xi - 17./2.*eta;
1050  case 2:
1051  return 0;
1052 
1053  default:
1054  libmesh_error_msg("Invalid subtriangle lookup = " <<
1055  subtriangle_lookup(p));
1056  }
1057  case 7:
1058  switch (subtriangle_lookup(p))
1059  {
1060  case 0:
1061  return -1 + 9./2.*xi + 5./2.*eta;
1062  case 1:
1063  return 9 - 23./2.*xi - 23./2.*eta;
1064  case 2:
1065  return 4*xi;
1066 
1067  default:
1068  libmesh_error_msg("Invalid subtriangle lookup = " <<
1069  subtriangle_lookup(p));
1070  }
1071  case 8:
1072  switch (subtriangle_lookup(p))
1073  {
1074  case 0:
1075  return -2 - 1./2.*xi + 13./2.*eta;
1076  case 1:
1077  return -4 + 7./2.*xi + 17./2.*eta;
1078  case 2:
1079  return -2 + 6*eta;
1080 
1081  default:
1082  libmesh_error_msg("Invalid subtriangle lookup = " <<
1083  subtriangle_lookup(p));
1084  }
1085  case 9:
1086  switch (subtriangle_lookup(p))
1087  {
1088  case 0:
1089  return std::sqrt(2.) * (2*xi - 2*eta);
1090  case 1:
1091  return std::sqrt(2.) * (8 - 14*xi - 10*eta);
1092  case 2:
1093  return 0;
1094 
1095  default:
1096  libmesh_error_msg("Invalid subtriangle lookup = " <<
1097  subtriangle_lookup(p));
1098  }
1099  case 10:
1100  switch (subtriangle_lookup(p))
1101  {
1102  case 0:
1103  return 4 - 4*xi - 16*eta;
1104  case 1:
1105  return -12 + 12*xi + 16*eta;
1106  case 2:
1107  return -8*xi;
1108 
1109  default:
1110  libmesh_error_msg("Invalid subtriangle lookup = " <<
1111  subtriangle_lookup(p));
1112  }
1113  case 11:
1114  switch (subtriangle_lookup(p))
1115  {
1116  case 0:
1117  return -4 - 8*xi + 20*eta;
1118  case 1:
1119  return -4 + 8*xi + 4*eta;
1120  case 2:
1121  return 0;
1122 
1123  default:
1124  libmesh_error_msg("Invalid subtriangle lookup = " <<
1125  subtriangle_lookup(p));
1126  }
1127 
1128  default:
1129  libmesh_error_msg("Invalid shape function index i = " <<
1130  basis_num);
1131  }
1132 
1133  default:
1134  libmesh_error_msg("Invalid shape function derivative j = " <<
1135  deriv_type);
1136  }
1137 }
1138 
1139 #endif
1140 
1141 
1142 
1143 Real clough_raw_shape_deriv(const unsigned int basis_num,
1144  const unsigned int deriv_type,
1145  const Point & p)
1146 {
1147  Real xi = p(0), eta = p(1);
1148 
1149  switch (deriv_type)
1150  {
1151  // derivative in xi direction
1152  case 0:
1153  switch (basis_num)
1154  {
1155  case 0:
1156  switch (subtriangle_lookup(p))
1157  {
1158  case 0:
1159  return -6*xi + 6*xi*xi
1160  - 3*eta*eta;
1161  case 1:
1162  return 9 - 30*xi + 21*xi*xi
1163  - 30*eta + 42*xi*eta
1164  + 21*eta*eta;
1165  case 2:
1166  return -6*xi + 9*xi*xi
1167  - 6*xi*eta;
1168 
1169  default:
1170  libmesh_error_msg("Invalid subtriangle lookup = " <<
1171  subtriangle_lookup(p));
1172  }
1173  case 1:
1174  switch (subtriangle_lookup(p))
1175  {
1176  case 0:
1177  return 6*xi - 6*xi*xi
1178  + 3./2.*eta*eta;
1179  case 1:
1180  return - 9./2. + 18*xi - 27./2.*xi*xi
1181  + 15*eta - 21*xi*eta
1182  - 21./2.*eta*eta;
1183  case 2:
1184  return 6*xi - 15./2.*xi*xi
1185  + 3*xi*eta;
1186 
1187  default:
1188  libmesh_error_msg("Invalid subtriangle lookup = " <<
1189  subtriangle_lookup(p));
1190  }
1191  case 2:
1192  switch (subtriangle_lookup(p))
1193  {
1194  case 0:
1195  return 3./2.*eta*eta;
1196  case 1:
1197  return - 9./2. + 12*xi - 15./2.*xi*xi
1198  + 15*eta - 21*xi*eta
1199  - 21./2.*eta*eta;
1200  case 2:
1201  return -3./2.*xi*xi
1202  + 3*xi*eta;
1203 
1204  default:
1205  libmesh_error_msg("Invalid subtriangle lookup = " <<
1206  subtriangle_lookup(p));
1207  }
1208  case 3:
1209  switch (subtriangle_lookup(p))
1210  {
1211  case 0:
1212  return 1 - 4*xi + 3*xi*xi
1213  - 1./2.*eta*eta;
1214  case 1:
1215  return 5./2. - 9*xi + 13./2.*xi*xi
1216  - 4*eta + 8*xi*eta
1217  + 3./2.*eta*eta;
1218  case 2:
1219  return 1 - xi - 7./2.*xi*xi
1220  - 3*eta + 4*xi*eta
1221  + 2*eta*eta;
1222 
1223  default:
1224  libmesh_error_msg("Invalid subtriangle lookup = " <<
1225  subtriangle_lookup(p));
1226  }
1227  case 4:
1228  switch (subtriangle_lookup(p))
1229  {
1230  case 0:
1231  return - 3*eta + 4*xi*eta
1232  + 2*eta*eta;
1233  case 1:
1234  return xi - xi*xi
1235  - 4*eta + 3*xi*eta
1236  + 4*eta*eta;
1237  case 2:
1238  return -3*xi + 7*xi*xi
1239  - xi*eta;
1240 
1241  default:
1242  libmesh_error_msg("Invalid subtriangle lookup = " <<
1243  subtriangle_lookup(p));
1244  }
1245  case 5:
1246  switch (subtriangle_lookup(p))
1247  {
1248  case 0:
1249  return -2*xi + 3*xi*xi
1250  - 1./4.*eta*eta;
1251  case 1:
1252  return 3./4. - 4*xi + 17./4.*xi*xi
1253  - 5./2.*eta + 7./2.*xi*eta
1254  + 7./4.*eta*eta;
1255  case 2:
1256  return -2*xi + 13./4.*xi*xi
1257  - 1./2.*xi*eta;
1258 
1259  default:
1260  libmesh_error_msg("Invalid subtriangle lookup = " <<
1261  subtriangle_lookup(p));
1262  }
1263  case 6:
1264  switch (subtriangle_lookup(p))
1265  {
1266  case 0:
1267  return -eta + 4*xi*eta
1268  + 7./4.*eta*eta;
1269  case 1:
1270  return -13./4. + 9*xi - 23./4.*xi*xi
1271  + 19./2.*eta - 23./2.*xi*eta
1272  - 25./4.*eta*eta;
1273  case 2:
1274  return -xi + 5./4.*xi*xi
1275  + 9./2.*xi*eta;
1276 
1277  default:
1278  libmesh_error_msg("Invalid subtriangle lookup = " <<
1279  subtriangle_lookup(p));
1280  }
1281  case 7:
1282  switch (subtriangle_lookup(p))
1283  {
1284  case 0:
1285  return 9./4.*eta*eta;
1286  case 1:
1287  return - 11./4. + 7*xi - 17./4.*xi*xi
1288  + 19./2.*eta - 25./2.*xi*eta
1289  - 23./4.*eta*eta;
1290  case 2:
1291  return xi - 13./4.*xi*xi
1292  - eta + 7./2.*xi*eta + 2*eta*eta;
1293 
1294  default:
1295  libmesh_error_msg("Invalid subtriangle lookup = " <<
1296  subtriangle_lookup(p));
1297  }
1298  case 8:
1299  switch (subtriangle_lookup(p))
1300  {
1301  case 0:
1302  return - 1./4.*eta*eta;
1303  case 1:
1304  return 3./4. - 2*xi + 5./4.*xi*xi
1305  - 5./2.*eta + 7./2.*xi*eta
1306  + 7./4.*eta*eta;
1307  case 2:
1308  return 1./4.*xi*xi
1309  - 1./2.*xi*eta;
1310 
1311  default:
1312  libmesh_error_msg("Invalid subtriangle lookup = " <<
1313  subtriangle_lookup(p));
1314  }
1315  case 9:
1316  switch (subtriangle_lookup(p))
1317  {
1318  case 0:
1319  return std::sqrt(2.) * eta*eta;
1320  case 1:
1321  return std::sqrt(2.) * (-3 + 8*xi - 5*xi*xi
1322  + 10*eta - 14*xi*eta
1323  - 7*eta*eta);
1324  case 2:
1325  return std::sqrt(2.) * (-xi*xi + 2*xi*eta);
1326 
1327  default:
1328  libmesh_error_msg("Invalid subtriangle lookup = " <<
1329  subtriangle_lookup(p));
1330  }
1331  case 10:
1332  switch (subtriangle_lookup(p))
1333  {
1334  case 0:
1335  return -2*eta*eta;
1336  case 1:
1337  return 2 - 4*xi + 2*xi*xi
1338  - 8*eta + 8*xi*eta
1339  + 6*eta*eta;
1340  case 2:
1341  return -4*xi + 10*xi*xi
1342  + 4*eta - 8*xi*eta
1343  - 4*eta*eta;
1344 
1345  default:
1346  libmesh_error_msg("Invalid subtriangle lookup = " <<
1347  subtriangle_lookup(p));
1348  }
1349  case 11:
1350  switch (subtriangle_lookup(p))
1351  {
1352  case 0:
1353  return 4*eta - 8*xi*eta
1354  - 4*eta*eta;
1355  case 1:
1356  return 4 - 12*xi + 8*xi*xi
1357  - 8*eta + 12*xi*eta
1358  + 4*eta*eta;
1359  case 2:
1360  return 4*xi - 8*xi*xi
1361  - 4*xi*eta;
1362 
1363  default:
1364  libmesh_error_msg("Invalid subtriangle lookup = " <<
1365  subtriangle_lookup(p));
1366  }
1367 
1368  default:
1369  libmesh_error_msg("Invalid shape function index i = " <<
1370  basis_num);
1371  }
1372 
1373  // derivative in eta direction
1374  case 1:
1375  switch (basis_num)
1376  {
1377  case 0:
1378  switch (subtriangle_lookup(p))
1379  {
1380  case 0:
1381  return - 6*eta - 6*xi*eta + 9*eta*eta;
1382  case 1:
1383  return 9 - 30*xi + 21*xi*xi
1384  - 30*eta + 42*xi*eta + 21*eta*eta;
1385  case 2:
1386  return - 3*xi*xi
1387  - 6*eta + 6*eta*eta;
1388 
1389  default:
1390  libmesh_error_msg("Invalid subtriangle lookup = " <<
1391  subtriangle_lookup(p));
1392  }
1393  case 1:
1394  switch (subtriangle_lookup(p))
1395  {
1396  case 0:
1397  return + 3*xi*eta
1398  - 3./2.*eta*eta;
1399  case 1:
1400  return - 9./2. + 15*xi - 21./2.*xi*xi
1401  + 12*eta - 21*xi*eta - 15./2.*eta*eta;
1402  case 2:
1403  return + 3./2.*xi*xi;
1404 
1405  default:
1406  libmesh_error_msg("Invalid subtriangle lookup = " <<
1407  subtriangle_lookup(p));
1408  }
1409  case 2:
1410  switch (subtriangle_lookup(p))
1411  {
1412  case 0:
1413  return 6*eta + 3*xi*eta - 15./2.*eta*eta;
1414  case 1:
1415  return - 9./2. + 15*xi - 21./2.*xi*xi
1416  + 18*eta - 21.*xi*eta - 27./2.*eta*eta;
1417  case 2:
1418  return 3./2.*xi*xi
1419  + 6*eta - 6*eta*eta;
1420 
1421  default:
1422  libmesh_error_msg("Invalid subtriangle lookup = " <<
1423  subtriangle_lookup(p));
1424  }
1425  case 3:
1426  switch (subtriangle_lookup(p))
1427  {
1428  case 0:
1429  return - 3*eta - xi*eta + 7*eta*eta;
1430  case 1:
1431  return - 4*xi + 4*xi*xi
1432  + eta + 3*xi*eta - eta*eta;
1433  case 2:
1434  return - 3*xi + 2*xi*xi
1435  + 4*xi*eta;
1436 
1437  default:
1438  libmesh_error_msg("Invalid subtriangle lookup = " <<
1439  subtriangle_lookup(p));
1440  }
1441  case 4:
1442  switch (subtriangle_lookup(p))
1443  {
1444  case 0:
1445  return 1 - 3*xi + 2*xi*xi
1446  - eta + 4*xi*eta - 7./2.*eta*eta;
1447  case 1:
1448  return 5./2. - 4*xi + 3./2.*xi*xi
1449  - 9.*eta + 8*xi*eta + 13./2.*eta*eta;
1450  case 2:
1451  return 1 - 1./2.*xi*xi - 4*eta + 3*eta*eta;
1452 
1453  default:
1454  libmesh_error_msg("Invalid subtriangle lookup = " <<
1455  subtriangle_lookup(p));
1456  }
1457  case 5:
1458  switch (subtriangle_lookup(p))
1459  {
1460  case 0:
1461  return - 1./2.*xi*eta + 1./4.*eta*eta;
1462  case 1:
1463  return 3./4. - 5./2.*xi + 7./4.*xi*xi
1464  - 2*eta + 7./2.*xi*eta + 5./4.*eta*eta;
1465  case 2:
1466  return - 1./4.*xi*xi;
1467 
1468  default:
1469  libmesh_error_msg("Invalid subtriangle lookup = " <<
1470  subtriangle_lookup(p));
1471  }
1472  case 6:
1473  switch (subtriangle_lookup(p))
1474  {
1475  case 0:
1476  return -xi + 2*xi*xi
1477  + eta + 7./2.*xi*eta - 13./4.*eta*eta;
1478  case 1:
1479  return - 11./4. + 19./2.*xi - 23./4.*xi*xi
1480  + 7*eta - 25./2.*xi*eta - 17./4.*eta*eta;
1481  case 2:
1482  return 9./4.*xi*xi;
1483 
1484  default:
1485  libmesh_error_msg("Invalid subtriangle lookup = " <<
1486  subtriangle_lookup(p));
1487  }
1488  case 7:
1489  switch (subtriangle_lookup(p))
1490  {
1491  case 0:
1492  return -eta + 9./2.*xi*eta + 5./4.*eta*eta;
1493  case 1:
1494  return - 13./4. + 19./2.*xi - 25./4.*xi*xi
1495  + 9*eta - 23./2.*xi*eta - 23./4.*eta*eta;
1496  case 2:
1497  return - xi + 7./4.*xi*xi + 4*xi*eta;
1498 
1499  default:
1500  libmesh_error_msg("Invalid subtriangle lookup = " <<
1501  subtriangle_lookup(p));
1502  }
1503  case 8:
1504  switch (subtriangle_lookup(p))
1505  {
1506  case 0:
1507  return -2*eta - 1./2.*xi*eta + 13./4.*eta*eta;
1508  case 1:
1509  return 3./4. - 5./2.*xi + 7./4.*xi*xi
1510  - 4*eta + 7./2.*xi*eta + 17./4.*eta*eta;
1511  case 2:
1512  return - 1./4.*xi*xi
1513  - 2*eta + 3*eta*eta;
1514 
1515  default:
1516  libmesh_error_msg("Invalid subtriangle lookup = " <<
1517  subtriangle_lookup(p));
1518  }
1519  case 9:
1520  switch (subtriangle_lookup(p))
1521  {
1522  case 0:
1523  return std::sqrt(2.) * (2*xi*eta - eta*eta);
1524  case 1:
1525  return std::sqrt(2.) * (- 3 + 10*xi - 7*xi*xi
1526  + 8*eta - 14*xi*eta - 5*eta*eta);
1527  case 2:
1528  return std::sqrt(2.) * (xi*xi);
1529 
1530  default:
1531  libmesh_error_msg("Invalid subtriangle lookup = " <<
1532  subtriangle_lookup(p));
1533  }
1534  case 10:
1535  switch (subtriangle_lookup(p))
1536  {
1537  case 0:
1538  return 4*eta - 4*xi*eta - 8*eta*eta;
1539  case 1:
1540  return 4 - 8*xi + 4*xi*xi
1541  - 12*eta + 12*xi*eta + 8*eta*eta;
1542  case 2:
1543  return 4*xi - 4*xi*xi
1544  - 8*xi*eta;
1545 
1546  default:
1547  libmesh_error_msg("Invalid subtriangle lookup = " <<
1548  subtriangle_lookup(p));
1549  }
1550  case 11:
1551  switch (subtriangle_lookup(p))
1552  {
1553  case 0:
1554  return 4*xi - 4*xi*xi
1555  - 4*eta - 8*xi*eta + 10.*eta*eta;
1556  case 1:
1557  return 2 - 8*xi + 6*xi*xi
1558  - 4*eta + 8*xi*eta + 2*eta*eta;
1559  case 2:
1560  return - 2*xi*xi;
1561 
1562  default:
1563  libmesh_error_msg("Invalid subtriangle lookup = " <<
1564  subtriangle_lookup(p));
1565  }
1566 
1567  default:
1568  libmesh_error_msg("Invalid shape function index i = " <<
1569  basis_num);
1570  }
1571 
1572  default:
1573  libmesh_error_msg("Invalid shape function derivative j = " <<
1574  deriv_type);
1575  }
1576 }
1577 
1578 Real clough_raw_shape(const unsigned int basis_num,
1579  const Point & p)
1580 {
1581  Real xi = p(0), eta = p(1);
1582 
1583  switch (basis_num)
1584  {
1585  case 0:
1586  switch (subtriangle_lookup(p))
1587  {
1588  case 0:
1589  return 1 - 3*xi*xi + 2*xi*xi*xi
1590  - 3*eta*eta - 3*xi*eta*eta + 3*eta*eta*eta;
1591  case 1:
1592  return -1 + 9*xi - 15*xi*xi + 7*xi*xi*xi
1593  + 9*eta - 30*xi*eta +21*xi*xi*eta
1594  - 15*eta*eta + 21*xi*eta*eta + 7*eta*eta*eta;
1595  case 2:
1596  return 1 - 3*xi*xi + 3*xi*xi*xi
1597  - 3*xi*xi*eta
1598  - 3*eta*eta + 2*eta*eta*eta;
1599 
1600  default:
1601  libmesh_error_msg("Invalid subtriangle lookup = " <<
1602  subtriangle_lookup(p));
1603  }
1604  case 1:
1605  switch (subtriangle_lookup(p))
1606  {
1607  case 0:
1608  return 3*xi*xi - 2*xi*xi*xi
1609  + 3./2.*xi*eta*eta
1610  - 1./2.*eta*eta*eta;
1611  case 1:
1612  return 1 - 9./2.*xi + 9*xi*xi - 9./2.*xi*xi*xi
1613  - 9./2.*eta + 15*xi*eta - 21./2.*xi*xi*eta
1614  + 6*eta*eta - 21./2.*xi*eta*eta - 5./2.*eta*eta*eta;
1615  case 2:
1616  return 3*xi*xi - 5./2.*xi*xi*xi
1617  + 3./2.*xi*xi*eta;
1618 
1619  default:
1620  libmesh_error_msg("Invalid subtriangle lookup = " <<
1621  subtriangle_lookup(p));
1622  }
1623  case 2:
1624  switch (subtriangle_lookup(p))
1625  {
1626  case 0:
1627  return 3*eta*eta + 3./2.*xi*eta*eta - 5./2.*eta*eta*eta;
1628  case 1:
1629  return 1 - 9./2.*xi + 6*xi*xi - 5./2.*xi*xi*xi
1630  - 9./2.*eta + 15*xi*eta - 21./2.*xi*xi*eta
1631  + 9*eta*eta - 21./2.*xi*eta*eta - 9./2.*eta*eta*eta;
1632  case 2:
1633  return -1./2.*xi*xi*xi
1634  + 3./2.*xi*xi*eta
1635  + 3*eta*eta - 2*eta*eta*eta;
1636 
1637  default:
1638  libmesh_error_msg("Invalid subtriangle lookup = " <<
1639  subtriangle_lookup(p));
1640  }
1641  case 3:
1642  switch (subtriangle_lookup(p))
1643  {
1644  case 0:
1645  return xi - 2*xi*xi + xi*xi*xi
1646  - 3./2.*eta*eta - 1./2.*xi*eta*eta + 7./3.*eta*eta*eta;
1647  case 1:
1648  return -1./6. + 5./2.*xi - 9./2.*xi*xi + 13./6.*xi*xi*xi
1649  - 4*xi*eta + 4*xi*xi*eta
1650  + 1./2.*eta*eta + 3./2.*xi*eta*eta - 1./3.*eta*eta*eta;
1651  case 2:
1652  return xi - 1./2.*xi*xi - 7./6.*xi*xi*xi
1653  - 3*xi*eta + 2*xi*xi*eta
1654  + 2*xi*eta*eta;
1655 
1656  default:
1657  libmesh_error_msg("Invalid subtriangle lookup = " <<
1658  subtriangle_lookup(p));
1659  }
1660  case 4:
1661  switch (subtriangle_lookup(p))
1662  {
1663  case 0:
1664  return eta - 3*xi*eta + 2*xi*xi*eta
1665  - 1./2.*eta*eta + 2*xi*eta*eta - 7./6.*eta*eta*eta;
1666  case 1:
1667  return -1./6. + 1./2.*xi*xi - 1./3.*xi*xi*xi
1668  + 5./2.*eta - 4*xi*eta + 3./2.*xi*xi*eta
1669  - 9./2.*eta*eta + 4*xi*eta*eta + 13./6.*eta*eta*eta;
1670  case 2:
1671  return -3./2.*xi*xi + 7./3.*xi*xi*xi
1672  + eta - 1./2.*xi*xi*eta - 2*eta*eta + eta*eta*eta;
1673 
1674  default:
1675  libmesh_error_msg("Invalid subtriangle lookup = " <<
1676  subtriangle_lookup(p));
1677  }
1678  case 5:
1679  switch (subtriangle_lookup(p))
1680  {
1681  case 0:
1682  return -xi*xi + xi*xi*xi
1683  - 1./4.*xi*eta*eta + 1./12.*eta*eta*eta;
1684  case 1:
1685  return -1./6. + 3./4.*xi - 2*xi*xi + 17./12.*xi*xi*xi
1686  + 3./4.*eta - 5./2.*xi*eta + 7./4.*xi*xi*eta
1687  - eta*eta + 7./4.*xi*eta*eta + 5./12.*eta*eta*eta;
1688  case 2:
1689  return -xi*xi + 13./12.*xi*xi*xi
1690  - 1./4.*xi*xi*eta;
1691 
1692  default:
1693  libmesh_error_msg("Invalid subtriangle lookup = " <<
1694  subtriangle_lookup(p));
1695  }
1696  case 6:
1697  switch (subtriangle_lookup(p))
1698  {
1699  case 0:
1700  return -xi*eta + 2*xi*xi*eta
1701  + 1./2.*eta*eta + 7./4.*xi*eta*eta - 13./12.*eta*eta*eta;
1702  case 1:
1703  return 2./3. - 13./4.*xi + 9./2.*xi*xi - 23./12.*xi*xi*xi
1704  - 11./4.*eta + 19./2.*xi*eta - 23./4.*xi*xi*eta
1705  + 7./2.*eta*eta - 25./4.*xi*eta*eta - 17./12.*eta*eta*eta;
1706  case 2:
1707  return -1./2.*xi*xi + 5./12.*xi*xi*xi
1708  + 9./4.*xi*xi*eta;
1709 
1710  default:
1711  libmesh_error_msg("Invalid subtriangle lookup = " <<
1712  subtriangle_lookup(p));
1713  }
1714  case 7:
1715  switch (subtriangle_lookup(p))
1716  {
1717  case 0:
1718  return -1./2.*eta*eta + 9./4.*xi*eta*eta + 5./12.*eta*eta*eta;
1719  case 1:
1720  return 2./3. - 11./4.*xi + 7./2.*xi*xi - 17./12.*xi*xi*xi
1721  - 13./4.*eta + 19./2.*xi*eta - 25./4.*xi*xi*eta
1722  + 9./2.*eta*eta - 23./4.*xi*eta*eta - 23./12.*eta*eta*eta;
1723  case 2:
1724  return 1./2.*xi*xi - 13./12.*xi*xi*xi
1725  - xi*eta + 7./4.*xi*xi*eta + 2*xi*eta*eta;
1726 
1727  default:
1728  libmesh_error_msg("Invalid subtriangle lookup = " <<
1729  subtriangle_lookup(p));
1730  }
1731  case 8:
1732  switch (subtriangle_lookup(p))
1733  {
1734  case 0:
1735  return -eta*eta - 1./4.*xi*eta*eta + 13./12.*eta*eta*eta;
1736  case 1:
1737  return -1./6. + 3./4.*xi - xi*xi + 5./12.*xi*xi*xi
1738  + 3./4.*eta - 5./2.*xi*eta + 7./4.*xi*xi*eta
1739  - 2*eta*eta + 7./4.*xi*eta*eta + 17./12.*eta*eta*eta;
1740  case 2:
1741  return 1./12.*xi*xi*xi
1742  - 1./4.*xi*xi*eta
1743  - eta*eta + eta*eta*eta;
1744 
1745  default:
1746  libmesh_error_msg("Invalid subtriangle lookup = " <<
1747  subtriangle_lookup(p));
1748  }
1749  case 9:
1750  switch (subtriangle_lookup(p))
1751  {
1752  case 0:
1753  return std::sqrt(2.) * (xi*eta*eta - 1./3.*eta*eta*eta);
1754  case 1:
1755  return std::sqrt(2.) * (2./3. - 3*xi + 4*xi*xi - 5./3.*xi*xi*xi
1756  - 3*eta + 10*xi*eta - 7*xi*xi*eta
1757  + 4*eta*eta - 7*xi*eta*eta - 5./3.*eta*eta*eta);
1758  case 2:
1759  return std::sqrt(2.) * (-1./3.*xi*xi*xi + xi*xi*eta);
1760 
1761  default:
1762  libmesh_error_msg("Invalid subtriangle lookup = " <<
1763  subtriangle_lookup(p));
1764  }
1765  case 10:
1766  switch (subtriangle_lookup(p))
1767  {
1768  case 0:
1769  return 2*eta*eta - 2*xi*eta*eta - 8./3.*eta*eta*eta;
1770  case 1:
1771  return -2./3. + 2*xi - 2*xi*xi + 2./3.*xi*xi*xi
1772  + 4*eta - 8*xi*eta + 4*xi*xi*eta
1773  - 6*eta*eta + 6*xi*eta*eta + 8./3.*eta*eta*eta;
1774  case 2:
1775  return -2*xi*xi + 10./3.*xi*xi*xi
1776  + 4*xi*eta - 4*xi*xi*eta
1777  - 4*xi*eta*eta;
1778 
1779  default:
1780  libmesh_error_msg("Invalid subtriangle lookup = " <<
1781  subtriangle_lookup(p));
1782  }
1783  case 11:
1784  switch (subtriangle_lookup(p))
1785  {
1786  case 0:
1787  return 4*xi*eta - 4*xi*xi*eta
1788  - 2*eta*eta - 4*xi*eta*eta + 10./3.*eta*eta*eta;
1789  case 1:
1790  return -2./3. + 4*xi - 6*xi*xi + 8./3.*xi*xi*xi
1791  + 2*eta - 8*xi*eta + 6*xi*xi*eta
1792  - 2*eta*eta + 4*xi*eta*eta + 2./3.*eta*eta*eta;
1793  case 2:
1794  return 2*xi*xi - 8./3.*xi*xi*xi
1795  - 2*xi*xi*eta;
1796 
1797  default:
1798  libmesh_error_msg("Invalid subtriangle lookup = " <<
1799  subtriangle_lookup(p));
1800  }
1801 
1802  default:
1803  libmesh_error_msg("Invalid shape function index i = " <<
1804  basis_num);
1805  }
1806 }
1807 
1808 
1809 } // end anonymous namespace
1810 
1811 
1812 namespace libMesh
1813 {
1814 
1815 
1816 template <>
1817 Real FE<2,CLOUGH>::shape(const ElemType,
1818  const Order,
1819  const unsigned int,
1820  const Point &)
1821 {
1822  libmesh_error_msg("Clough-Tocher elements require the real element \nto construct gradient-based degrees of freedom.");
1823  return 0.;
1824 }
1825 
1826 
1827 
1828 template <>
1829 Real FE<2,CLOUGH>::shape(const Elem * elem,
1830  const Order order,
1831  const unsigned int i,
1832  const Point & p,
1833  const bool add_p_level)
1834 {
1835  libmesh_assert(elem);
1836 
1837  clough_compute_coefs(elem);
1838 
1839  const ElemType type = elem->type();
1840 
1841  const Order totalorder =
1842  static_cast<Order>(order + add_p_level * elem->p_level());
1843 
1844  switch (totalorder)
1845  {
1846  // 2nd-order restricted Clough-Tocher element
1847  case SECOND:
1848  {
1849  // There may be a bug in the 2nd order case; the 3rd order
1850  // Clough-Tocher elements are pretty uniformly better anyways
1851  // so use those instead.
1852  libmesh_experimental();
1853  switch (type)
1854  {
1855  // C1 functions on the Clough-Tocher triangle.
1856  case TRI6:
1857  {
1858  libmesh_assert_less (i, 9);
1859  // FIXME: it would be nice to calculate (and cache)
1860  // clough_raw_shape(j,p) only once per triangle, not 1-7
1861  // times
1862  switch (i)
1863  {
1864  // Note: these DoF numbers are "scrambled" because my
1865  // initial numbering conventions didn't match libMesh
1866  case 0:
1867  return clough_raw_shape(0, p)
1868  + d1d2n * clough_raw_shape(10, p)
1869  + d1d3n * clough_raw_shape(11, p);
1870  case 3:
1871  return clough_raw_shape(1, p)
1872  + d2d3n * clough_raw_shape(11, p)
1873  + d2d1n * clough_raw_shape(9, p);
1874  case 6:
1875  return clough_raw_shape(2, p)
1876  + d3d1n * clough_raw_shape(9, p)
1877  + d3d2n * clough_raw_shape(10, p);
1878  case 1:
1879  return d1xd1x * clough_raw_shape(3, p)
1880  + d1xd1y * clough_raw_shape(4, p)
1881  + d1xd2n * clough_raw_shape(10, p)
1882  + d1xd3n * clough_raw_shape(11, p)
1883  + 0.5 * N01x * d3nd3n * clough_raw_shape(11, p)
1884  + 0.5 * N02x * d2nd2n * clough_raw_shape(10, p);
1885  case 2:
1886  return d1yd1y * clough_raw_shape(4, p)
1887  + d1yd1x * clough_raw_shape(3, p)
1888  + d1yd2n * clough_raw_shape(10, p)
1889  + d1yd3n * clough_raw_shape(11, p)
1890  + 0.5 * N01y * d3nd3n * clough_raw_shape(11, p)
1891  + 0.5 * N02y * d2nd2n * clough_raw_shape(10, p);
1892  case 4:
1893  return d2xd2x * clough_raw_shape(5, p)
1894  + d2xd2y * clough_raw_shape(6, p)
1895  + d2xd3n * clough_raw_shape(11, p)
1896  + d2xd1n * clough_raw_shape(9, p)
1897  + 0.5 * N10x * d3nd3n * clough_raw_shape(11, p)
1898  + 0.5 * N12x * d1nd1n * clough_raw_shape(9, p);
1899  case 5:
1900  return d2yd2y * clough_raw_shape(6, p)
1901  + d2yd2x * clough_raw_shape(5, p)
1902  + d2yd3n * clough_raw_shape(11, p)
1903  + d2yd1n * clough_raw_shape(9, p)
1904  + 0.5 * N10y * d3nd3n * clough_raw_shape(11, p)
1905  + 0.5 * N12y * d1nd1n * clough_raw_shape(9, p);
1906  case 7:
1907  return d3xd3x * clough_raw_shape(7, p)
1908  + d3xd3y * clough_raw_shape(8, p)
1909  + d3xd1n * clough_raw_shape(9, p)
1910  + d3xd2n * clough_raw_shape(10, p)
1911  + 0.5 * N20x * d2nd2n * clough_raw_shape(10, p)
1912  + 0.5 * N21x * d1nd1n * clough_raw_shape(9, p);
1913  case 8:
1914  return d3yd3y * clough_raw_shape(8, p)
1915  + d3yd3x * clough_raw_shape(7, p)
1916  + d3yd1n * clough_raw_shape(9, p)
1917  + d3yd2n * clough_raw_shape(10, p)
1918  + 0.5 * N20y * d2nd2n * clough_raw_shape(10, p)
1919  + 0.5 * N21y * d1nd1n * clough_raw_shape(9, p);
1920  default:
1921  libmesh_error_msg("Invalid shape function index i = " << i);
1922  }
1923  }
1924  default:
1925  libmesh_error_msg("ERROR: Unsupported element type = " << type);
1926  }
1927  }
1928  // 3rd-order Clough-Tocher element
1929  case THIRD:
1930  {
1931  switch (type)
1932  {
1933  // C1 functions on the Clough-Tocher triangle.
1934  case TRI6:
1935  {
1936  libmesh_assert_less (i, 12);
1937 
1938  // FIXME: it would be nice to calculate (and cache)
1939  // clough_raw_shape(j,p) only once per triangle, not 1-7
1940  // times
1941  switch (i)
1942  {
1943  // Note: these DoF numbers are "scrambled" because my
1944  // initial numbering conventions didn't match libMesh
1945  case 0:
1946  return clough_raw_shape(0, p)
1947  + d1d2n * clough_raw_shape(10, p)
1948  + d1d3n * clough_raw_shape(11, p);
1949  case 3:
1950  return clough_raw_shape(1, p)
1951  + d2d3n * clough_raw_shape(11, p)
1952  + d2d1n * clough_raw_shape(9, p);
1953  case 6:
1954  return clough_raw_shape(2, p)
1955  + d3d1n * clough_raw_shape(9, p)
1956  + d3d2n * clough_raw_shape(10, p);
1957  case 1:
1958  return d1xd1x * clough_raw_shape(3, p)
1959  + d1xd1y * clough_raw_shape(4, p)
1960  + d1xd2n * clough_raw_shape(10, p)
1961  + d1xd3n * clough_raw_shape(11, p);
1962  case 2:
1963  return d1yd1y * clough_raw_shape(4, p)
1964  + d1yd1x * clough_raw_shape(3, p)
1965  + d1yd2n * clough_raw_shape(10, p)
1966  + d1yd3n * clough_raw_shape(11, p);
1967  case 4:
1968  return d2xd2x * clough_raw_shape(5, p)
1969  + d2xd2y * clough_raw_shape(6, p)
1970  + d2xd3n * clough_raw_shape(11, p)
1971  + d2xd1n * clough_raw_shape(9, p);
1972  case 5:
1973  return d2yd2y * clough_raw_shape(6, p)
1974  + d2yd2x * clough_raw_shape(5, p)
1975  + d2yd3n * clough_raw_shape(11, p)
1976  + d2yd1n * clough_raw_shape(9, p);
1977  case 7:
1978  return d3xd3x * clough_raw_shape(7, p)
1979  + d3xd3y * clough_raw_shape(8, p)
1980  + d3xd1n * clough_raw_shape(9, p)
1981  + d3xd2n * clough_raw_shape(10, p);
1982  case 8:
1983  return d3yd3y * clough_raw_shape(8, p)
1984  + d3yd3x * clough_raw_shape(7, p)
1985  + d3yd1n * clough_raw_shape(9, p)
1986  + d3yd2n * clough_raw_shape(10, p);
1987  case 10:
1988  return d1nd1n * clough_raw_shape(9, p);
1989  case 11:
1990  return d2nd2n * clough_raw_shape(10, p);
1991  case 9:
1992  return d3nd3n * clough_raw_shape(11, p);
1993 
1994  default:
1995  libmesh_error_msg("Invalid shape function index i = " << i);
1996  }
1997  }
1998  default:
1999  libmesh_error_msg("ERROR: Unsupported element type = " << type);
2000  }
2001  }
2002  // by default throw an error
2003  default:
2004  libmesh_error_msg("ERROR: Unsupported polynomial order = " << order);
2005  }
2006 }
2007 
2008 
2009 
2010 template <>
2011 Real FE<2,CLOUGH>::shape_deriv(const ElemType,
2012  const Order,
2013  const unsigned int,
2014  const unsigned int,
2015  const Point &)
2016 {
2017  libmesh_error_msg("Clough-Tocher elements require the real element \nto construct gradient-based degrees of freedom.");
2018  return 0.;
2019 }
2020 
2021 
2022 
2023 template <>
2024 Real FE<2,CLOUGH>::shape_deriv(const Elem * elem,
2025  const Order order,
2026  const unsigned int i,
2027  const unsigned int j,
2028  const Point & p,
2029  const bool add_p_level)
2030 {
2031  libmesh_assert(elem);
2032 
2033  clough_compute_coefs(elem);
2034 
2035  const ElemType type = elem->type();
2036 
2037  const Order totalorder =
2038  static_cast<Order>(order + add_p_level * elem->p_level());
2039 
2040  switch (totalorder)
2041  {
2042  // 2nd-order restricted Clough-Tocher element
2043  case SECOND:
2044  {
2045  // There may be a bug in the 2nd order case; the 3rd order
2046  // Clough-Tocher elements are pretty uniformly better anyways
2047  // so use those instead.
2048  libmesh_experimental();
2049  switch (type)
2050  {
2051  // C1 functions on the Clough-Tocher triangle.
2052  case TRI6:
2053  {
2054  libmesh_assert_less (i, 9);
2055  // FIXME: it would be nice to calculate (and cache)
2056  // clough_raw_shape(j,p) only once per triangle, not 1-7
2057  // times
2058  switch (i)
2059  {
2060  // Note: these DoF numbers are "scrambled" because my
2061  // initial numbering conventions didn't match libMesh
2062  case 0:
2063  return clough_raw_shape_deriv(0, j, p)
2064  + d1d2n * clough_raw_shape_deriv(10, j, p)
2065  + d1d3n * clough_raw_shape_deriv(11, j, p);
2066  case 3:
2067  return clough_raw_shape_deriv(1, j, p)
2068  + d2d3n * clough_raw_shape_deriv(11, j, p)
2069  + d2d1n * clough_raw_shape_deriv(9, j, p);
2070  case 6:
2071  return clough_raw_shape_deriv(2, j, p)
2072  + d3d1n * clough_raw_shape_deriv(9, j, p)
2073  + d3d2n * clough_raw_shape_deriv(10, j, p);
2074  case 1:
2075  return d1xd1x * clough_raw_shape_deriv(3, j, p)
2076  + d1xd1y * clough_raw_shape_deriv(4, j, p)
2077  + d1xd2n * clough_raw_shape_deriv(10, j, p)
2078  + d1xd3n * clough_raw_shape_deriv(11, j, p)
2079  + 0.5 * N01x * d3nd3n * clough_raw_shape_deriv(11, j, p)
2080  + 0.5 * N02x * d2nd2n * clough_raw_shape_deriv(10, j, p);
2081  case 2:
2082  return d1yd1y * clough_raw_shape_deriv(4, j, p)
2083  + d1yd1x * clough_raw_shape_deriv(3, j, p)
2084  + d1yd2n * clough_raw_shape_deriv(10, j, p)
2085  + d1yd3n * clough_raw_shape_deriv(11, j, p)
2086  + 0.5 * N01y * d3nd3n * clough_raw_shape_deriv(11, j, p)
2087  + 0.5 * N02y * d2nd2n * clough_raw_shape_deriv(10, j, p);
2088  case 4:
2089  return d2xd2x * clough_raw_shape_deriv(5, j, p)
2090  + d2xd2y * clough_raw_shape_deriv(6, j, p)
2091  + d2xd3n * clough_raw_shape_deriv(11, j, p)
2092  + d2xd1n * clough_raw_shape_deriv(9, j, p)
2093  + 0.5 * N10x * d3nd3n * clough_raw_shape_deriv(11, j, p)
2094  + 0.5 * N12x * d1nd1n * clough_raw_shape_deriv(9, j, p);
2095  case 5:
2096  return d2yd2y * clough_raw_shape_deriv(6, j, p)
2097  + d2yd2x * clough_raw_shape_deriv(5, j, p)
2098  + d2yd3n * clough_raw_shape_deriv(11, j, p)
2099  + d2yd1n * clough_raw_shape_deriv(9, j, p)
2100  + 0.5 * N10y * d3nd3n * clough_raw_shape_deriv(11, j, p)
2101  + 0.5 * N12y * d1nd1n * clough_raw_shape_deriv(9, j, p);
2102  case 7:
2103  return d3xd3x * clough_raw_shape_deriv(7, j, p)
2104  + d3xd3y * clough_raw_shape_deriv(8, j, p)
2105  + d3xd1n * clough_raw_shape_deriv(9, j, p)
2106  + d3xd2n * clough_raw_shape_deriv(10, j, p)
2107  + 0.5 * N20x * d2nd2n * clough_raw_shape_deriv(10, j, p)
2108  + 0.5 * N21x * d1nd1n * clough_raw_shape_deriv(9, j, p);
2109  case 8:
2110  return d3yd3y * clough_raw_shape_deriv(8, j, p)
2111  + d3yd3x * clough_raw_shape_deriv(7, j, p)
2112  + d3yd1n * clough_raw_shape_deriv(9, j, p)
2113  + d3yd2n * clough_raw_shape_deriv(10, j, p)
2114  + 0.5 * N20y * d2nd2n * clough_raw_shape_deriv(10, j, p)
2115  + 0.5 * N21y * d1nd1n * clough_raw_shape_deriv(9, j, p);
2116  default:
2117  libmesh_error_msg("Invalid shape function index i = " << i);
2118  }
2119  }
2120  default:
2121  libmesh_error_msg("ERROR: Unsupported element type = " << type);
2122  }
2123  }
2124  // 3rd-order Clough-Tocher element
2125  case THIRD:
2126  {
2127  switch (type)
2128  {
2129  // C1 functions on the Clough-Tocher triangle.
2130  case TRI6:
2131  {
2132  libmesh_assert_less (i, 12);
2133 
2134  // FIXME: it would be nice to calculate (and cache)
2135  // clough_raw_shape(j,p) only once per triangle, not 1-7
2136  // times
2137  switch (i)
2138  {
2139  // Note: these DoF numbers are "scrambled" because my
2140  // initial numbering conventions didn't match libMesh
2141  case 0:
2142  return clough_raw_shape_deriv(0, j, p)
2143  + d1d2n * clough_raw_shape_deriv(10, j, p)
2144  + d1d3n * clough_raw_shape_deriv(11, j, p);
2145  case 3:
2146  return clough_raw_shape_deriv(1, j, p)
2147  + d2d3n * clough_raw_shape_deriv(11, j, p)
2148  + d2d1n * clough_raw_shape_deriv(9, j, p);
2149  case 6:
2150  return clough_raw_shape_deriv(2, j, p)
2151  + d3d1n * clough_raw_shape_deriv(9, j, p)
2152  + d3d2n * clough_raw_shape_deriv(10, j, p);
2153  case 1:
2154  return d1xd1x * clough_raw_shape_deriv(3, j, p)
2155  + d1xd1y * clough_raw_shape_deriv(4, j, p)
2156  + d1xd2n * clough_raw_shape_deriv(10, j, p)
2157  + d1xd3n * clough_raw_shape_deriv(11, j, p);
2158  case 2:
2159  return d1yd1y * clough_raw_shape_deriv(4, j, p)
2160  + d1yd1x * clough_raw_shape_deriv(3, j, p)
2161  + d1yd2n * clough_raw_shape_deriv(10, j, p)
2162  + d1yd3n * clough_raw_shape_deriv(11, j, p);
2163  case 4:
2164  return d2xd2x * clough_raw_shape_deriv(5, j, p)
2165  + d2xd2y * clough_raw_shape_deriv(6, j, p)
2166  + d2xd3n * clough_raw_shape_deriv(11, j, p)
2167  + d2xd1n * clough_raw_shape_deriv(9, j, p);
2168  case 5:
2169  return d2yd2y * clough_raw_shape_deriv(6, j, p)
2170  + d2yd2x * clough_raw_shape_deriv(5, j, p)
2171  + d2yd3n * clough_raw_shape_deriv(11, j, p)
2172  + d2yd1n * clough_raw_shape_deriv(9, j, p);
2173  case 7:
2174  return d3xd3x * clough_raw_shape_deriv(7, j, p)
2175  + d3xd3y * clough_raw_shape_deriv(8, j, p)
2176  + d3xd1n * clough_raw_shape_deriv(9, j, p)
2177  + d3xd2n * clough_raw_shape_deriv(10, j, p);
2178  case 8:
2179  return d3yd3y * clough_raw_shape_deriv(8, j, p)
2180  + d3yd3x * clough_raw_shape_deriv(7, j, p)
2181  + d3yd1n * clough_raw_shape_deriv(9, j, p)
2182  + d3yd2n * clough_raw_shape_deriv(10, j, p);
2183  case 10:
2184  return d1nd1n * clough_raw_shape_deriv(9, j, p);
2185  case 11:
2186  return d2nd2n * clough_raw_shape_deriv(10, j, p);
2187  case 9:
2188  return d3nd3n * clough_raw_shape_deriv(11, j, p);
2189 
2190  default:
2191  libmesh_error_msg("Invalid shape function index i = " << i);
2192  }
2193  }
2194  default:
2195  libmesh_error_msg("ERROR: Unsupported element type = " << type);
2196  }
2197  }
2198  // by default throw an error
2199  default:
2200  libmesh_error_msg("ERROR: Unsupported polynomial order = " << order);
2201  }
2202 }
2203 
2204 
2205 
2206 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
2207 
2208 template <>
2209 Real FE<2,CLOUGH>::shape_second_deriv(const ElemType,
2210  const Order,
2211  const unsigned int,
2212  const unsigned int,
2213  const Point &)
2214 {
2215  libmesh_error_msg("Clough-Tocher elements require the real element \nto construct gradient-based degrees of freedom.");
2216  return 0.;
2217 }
2218 
2219 
2220 template <>
2221 Real FE<2,CLOUGH>::shape_second_deriv(const Elem * elem,
2222  const Order order,
2223  const unsigned int i,
2224  const unsigned int j,
2225  const Point & p,
2226  const bool add_p_level)
2227 {
2228  libmesh_assert(elem);
2229 
2230  clough_compute_coefs(elem);
2231 
2232  const ElemType type = elem->type();
2233 
2234  const Order totalorder =
2235  static_cast<Order>(order + add_p_level * elem->p_level());
2236 
2237  switch (totalorder)
2238  {
2239  // 2nd-order restricted Clough-Tocher element
2240  case SECOND:
2241  {
2242  switch (type)
2243  {
2244  // C1 functions on the Clough-Tocher triangle.
2245  case TRI6:
2246  {
2247  libmesh_assert_less (i, 9);
2248  // FIXME: it would be nice to calculate (and cache)
2249  // clough_raw_shape(j,p) only once per triangle, not 1-7
2250  // times
2251  switch (i)
2252  {
2253  // Note: these DoF numbers are "scrambled" because my
2254  // initial numbering conventions didn't match libMesh
2255  case 0:
2256  return clough_raw_shape_second_deriv(0, j, p)
2257  + d1d2n * clough_raw_shape_second_deriv(10, j, p)
2258  + d1d3n * clough_raw_shape_second_deriv(11, j, p);
2259  case 3:
2260  return clough_raw_shape_second_deriv(1, j, p)
2261  + d2d3n * clough_raw_shape_second_deriv(11, j, p)
2262  + d2d1n * clough_raw_shape_second_deriv(9, j, p);
2263  case 6:
2264  return clough_raw_shape_second_deriv(2, j, p)
2265  + d3d1n * clough_raw_shape_second_deriv(9, j, p)
2266  + d3d2n * clough_raw_shape_second_deriv(10, j, p);
2267  case 1:
2268  return d1xd1x * clough_raw_shape_second_deriv(3, j, p)
2269  + d1xd1y * clough_raw_shape_second_deriv(4, j, p)
2270  + d1xd2n * clough_raw_shape_second_deriv(10, j, p)
2271  + d1xd3n * clough_raw_shape_second_deriv(11, j, p)
2272  + 0.5 * N01x * d3nd3n * clough_raw_shape_second_deriv(11, j, p)
2273  + 0.5 * N02x * d2nd2n * clough_raw_shape_second_deriv(10, j, p);
2274  case 2:
2275  return d1yd1y * clough_raw_shape_second_deriv(4, j, p)
2276  + d1yd1x * clough_raw_shape_second_deriv(3, j, p)
2277  + d1yd2n * clough_raw_shape_second_deriv(10, j, p)
2278  + d1yd3n * clough_raw_shape_second_deriv(11, j, p)
2279  + 0.5 * N01y * d3nd3n * clough_raw_shape_second_deriv(11, j, p)
2280  + 0.5 * N02y * d2nd2n * clough_raw_shape_second_deriv(10, j, p);
2281  case 4:
2282  return d2xd2x * clough_raw_shape_second_deriv(5, j, p)
2283  + d2xd2y * clough_raw_shape_second_deriv(6, j, p)
2284  + d2xd3n * clough_raw_shape_second_deriv(11, j, p)
2285  + d2xd1n * clough_raw_shape_second_deriv(9, j, p)
2286  + 0.5 * N10x * d3nd3n * clough_raw_shape_second_deriv(11, j, p)
2287  + 0.5 * N12x * d1nd1n * clough_raw_shape_second_deriv(9, j, p);
2288  case 5:
2289  return d2yd2y * clough_raw_shape_second_deriv(6, j, p)
2290  + d2yd2x * clough_raw_shape_second_deriv(5, j, p)
2291  + d2yd3n * clough_raw_shape_second_deriv(11, j, p)
2292  + d2yd1n * clough_raw_shape_second_deriv(9, j, p)
2293  + 0.5 * N10y * d3nd3n * clough_raw_shape_second_deriv(11, j, p)
2294  + 0.5 * N12y * d1nd1n * clough_raw_shape_second_deriv(9, j, p);
2295  case 7:
2296  return d3xd3x * clough_raw_shape_second_deriv(7, j, p)
2297  + d3xd3y * clough_raw_shape_second_deriv(8, j, p)
2298  + d3xd1n * clough_raw_shape_second_deriv(9, j, p)
2299  + d3xd2n * clough_raw_shape_second_deriv(10, j, p)
2300  + 0.5 * N20x * d2nd2n * clough_raw_shape_second_deriv(10, j, p)
2301  + 0.5 * N21x * d1nd1n * clough_raw_shape_second_deriv(9, j, p);
2302  case 8:
2303  return d3yd3y * clough_raw_shape_second_deriv(8, j, p)
2304  + d3yd3x * clough_raw_shape_second_deriv(7, j, p)
2305  + d3yd1n * clough_raw_shape_second_deriv(9, j, p)
2306  + d3yd2n * clough_raw_shape_second_deriv(10, j, p)
2307  + 0.5 * N20y * d2nd2n * clough_raw_shape_second_deriv(10, j, p)
2308  + 0.5 * N21y * d1nd1n * clough_raw_shape_second_deriv(9, j, p);
2309  default:
2310  libmesh_error_msg("Invalid shape function index i = " << i);
2311  }
2312  }
2313  default:
2314  libmesh_error_msg("ERROR: Unsupported element type = " << type);
2315  }
2316  }
2317  // 3rd-order Clough-Tocher element
2318  case THIRD:
2319  {
2320  switch (type)
2321  {
2322  // C1 functions on the Clough-Tocher triangle.
2323  case TRI6:
2324  {
2325  libmesh_assert_less (i, 12);
2326 
2327  // FIXME: it would be nice to calculate (and cache)
2328  // clough_raw_shape(j,p) only once per triangle, not 1-7
2329  // times
2330  switch (i)
2331  {
2332  // Note: these DoF numbers are "scrambled" because my
2333  // initial numbering conventions didn't match libMesh
2334  case 0:
2335  return clough_raw_shape_second_deriv(0, j, p)
2336  + d1d2n * clough_raw_shape_second_deriv(10, j, p)
2337  + d1d3n * clough_raw_shape_second_deriv(11, j, p);
2338  case 3:
2339  return clough_raw_shape_second_deriv(1, j, p)
2340  + d2d3n * clough_raw_shape_second_deriv(11, j, p)
2341  + d2d1n * clough_raw_shape_second_deriv(9, j, p);
2342  case 6:
2343  return clough_raw_shape_second_deriv(2, j, p)
2344  + d3d1n * clough_raw_shape_second_deriv(9, j, p)
2345  + d3d2n * clough_raw_shape_second_deriv(10, j, p);
2346  case 1:
2347  return d1xd1x * clough_raw_shape_second_deriv(3, j, p)
2348  + d1xd1y * clough_raw_shape_second_deriv(4, j, p)
2349  + d1xd2n * clough_raw_shape_second_deriv(10, j, p)
2350  + d1xd3n * clough_raw_shape_second_deriv(11, j, p);
2351  case 2:
2352  return d1yd1y * clough_raw_shape_second_deriv(4, j, p)
2353  + d1yd1x * clough_raw_shape_second_deriv(3, j, p)
2354  + d1yd2n * clough_raw_shape_second_deriv(10, j, p)
2355  + d1yd3n * clough_raw_shape_second_deriv(11, j, p);
2356  case 4:
2357  return d2xd2x * clough_raw_shape_second_deriv(5, j, p)
2358  + d2xd2y * clough_raw_shape_second_deriv(6, j, p)
2359  + d2xd3n * clough_raw_shape_second_deriv(11, j, p)
2360  + d2xd1n * clough_raw_shape_second_deriv(9, j, p);
2361  case 5:
2362  return d2yd2y * clough_raw_shape_second_deriv(6, j, p)
2363  + d2yd2x * clough_raw_shape_second_deriv(5, j, p)
2364  + d2yd3n * clough_raw_shape_second_deriv(11, j, p)
2365  + d2yd1n * clough_raw_shape_second_deriv(9, j, p);
2366  case 7:
2367  return d3xd3x * clough_raw_shape_second_deriv(7, j, p)
2368  + d3xd3y * clough_raw_shape_second_deriv(8, j, p)
2369  + d3xd1n * clough_raw_shape_second_deriv(9, j, p)
2370  + d3xd2n * clough_raw_shape_second_deriv(10, j, p);
2371  case 8:
2372  return d3yd3y * clough_raw_shape_second_deriv(8, j, p)
2373  + d3yd3x * clough_raw_shape_second_deriv(7, j, p)
2374  + d3yd1n * clough_raw_shape_second_deriv(9, j, p)
2375  + d3yd2n * clough_raw_shape_second_deriv(10, j, p);
2376  case 10:
2377  return d1nd1n * clough_raw_shape_second_deriv(9, j, p);
2378  case 11:
2379  return d2nd2n * clough_raw_shape_second_deriv(10, j, p);
2380  case 9:
2381  return d3nd3n * clough_raw_shape_second_deriv(11, j, p);
2382 
2383  default:
2384  libmesh_error_msg("Invalid shape function index i = " << i);
2385  }
2386  }
2387  default:
2388  libmesh_error_msg("ERROR: Unsupported element type = " << type);
2389  }
2390  }
2391  // by default throw an error
2392  default:
2393  libmesh_error_msg("ERROR: Unsupported polynomial order = " << order);
2394  }
2395 }
2396 
2397 #endif
2398 
2399 } // namespace libMesh
libMesh::FE::shape_second_deriv
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:67
libMesh::FEInterface::n_shape_functions
static unsigned int n_shape_functions(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:446
libMesh::n_threads
unsigned int n_threads()
Definition: libmesh_base.h:96
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
Definition: factoryfunction.C:55
libMesh::Order
Order
Definition: enum_order.h:40
libMesh::FEInterface::shape_deriv_function
static shape_deriv_ptr shape_deriv_function(const unsigned int dim, const FEType &fe_t)
Definition: fe_interface.C:916
std::sqrt
MetaPhysicL::DualNumber< T, D > sqrt(const MetaPhysicL::DualNumber< T, D > &in)
libMesh::Elem::p_level
unsigned int p_level() const
Definition: elem.h:2512
libMesh::SECOND
Definition: enum_order.h:43
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::FE::shape_deriv
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::DofObject::invalid_id
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
Definition: dof_object.h:421
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
libMesh::FEMap::map_fe_type
static FEFamily map_fe_type(const Elem &elem)
Definition: fe_map.C:47
libMesh::FEInterface::shape_deriv_ptr
Real(* shape_deriv_ptr)(const Elem *elem, const Order o, const unsigned int i, const unsigned int j, const Point &p, const bool add_p_level)
Definition: fe_interface.h:345
libMesh::TRI6
Definition: enum_elem_type.h:40
libMesh::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
libMesh::FEFamily
FEFamily
Definition: enum_fe_family.h:34
libMesh::DofObject::id
dof_id_type id() const
Definition: dof_object.h:767
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:100
libMesh::THIRD
Definition: enum_order.h:44
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:121
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
Generated on Sat Jan 25 2020 12:06:57 for libMesh by   doxygen 1.8.16