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