fe_clough_shape_2D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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 = 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 
1126 
1127 
1128 Real clough_raw_shape_deriv(const unsigned int basis_num,
1129  const unsigned int deriv_type,
1130  const Point & p)
1131 {
1132  Real xi = p(0), eta = p(1);
1133 
1134  switch (deriv_type)
1135  {
1136  // derivative in xi direction
1137  case 0:
1138  switch (basis_num)
1139  {
1140  case 0:
1141  switch (subtriangle_lookup(p))
1142  {
1143  case 0:
1144  return -6*xi + 6*xi*xi
1145  - 3*eta*eta;
1146  case 1:
1147  return 9 - 30*xi + 21*xi*xi
1148  - 30*eta + 42*xi*eta
1149  + 21*eta*eta;
1150  case 2:
1151  return -6*xi + 9*xi*xi
1152  - 6*xi*eta;
1153 
1154  default:
1155  libmesh_error_msg("Invalid subtriangle lookup = " <<
1156  subtriangle_lookup(p));
1157  }
1158  case 1:
1159  switch (subtriangle_lookup(p))
1160  {
1161  case 0:
1162  return 6*xi - 6*xi*xi
1163  + 3./2.*eta*eta;
1164  case 1:
1165  return - 9./2. + 18*xi - 27./2.*xi*xi
1166  + 15*eta - 21*xi*eta
1167  - 21./2.*eta*eta;
1168  case 2:
1169  return 6*xi - 15./2.*xi*xi
1170  + 3*xi*eta;
1171 
1172  default:
1173  libmesh_error_msg("Invalid subtriangle lookup = " <<
1174  subtriangle_lookup(p));
1175  }
1176  case 2:
1177  switch (subtriangle_lookup(p))
1178  {
1179  case 0:
1180  return 3./2.*eta*eta;
1181  case 1:
1182  return - 9./2. + 12*xi - 15./2.*xi*xi
1183  + 15*eta - 21*xi*eta
1184  - 21./2.*eta*eta;
1185  case 2:
1186  return -3./2.*xi*xi
1187  + 3*xi*eta;
1188 
1189  default:
1190  libmesh_error_msg("Invalid subtriangle lookup = " <<
1191  subtriangle_lookup(p));
1192  }
1193  case 3:
1194  switch (subtriangle_lookup(p))
1195  {
1196  case 0:
1197  return 1 - 4*xi + 3*xi*xi
1198  - 1./2.*eta*eta;
1199  case 1:
1200  return 5./2. - 9*xi + 13./2.*xi*xi
1201  - 4*eta + 8*xi*eta
1202  + 3./2.*eta*eta;
1203  case 2:
1204  return 1 - xi - 7./2.*xi*xi
1205  - 3*eta + 4*xi*eta
1206  + 2*eta*eta;
1207 
1208  default:
1209  libmesh_error_msg("Invalid subtriangle lookup = " <<
1210  subtriangle_lookup(p));
1211  }
1212  case 4:
1213  switch (subtriangle_lookup(p))
1214  {
1215  case 0:
1216  return - 3*eta + 4*xi*eta
1217  + 2*eta*eta;
1218  case 1:
1219  return xi - xi*xi
1220  - 4*eta + 3*xi*eta
1221  + 4*eta*eta;
1222  case 2:
1223  return -3*xi + 7*xi*xi
1224  - xi*eta;
1225 
1226  default:
1227  libmesh_error_msg("Invalid subtriangle lookup = " <<
1228  subtriangle_lookup(p));
1229  }
1230  case 5:
1231  switch (subtriangle_lookup(p))
1232  {
1233  case 0:
1234  return -2*xi + 3*xi*xi
1235  - 1./4.*eta*eta;
1236  case 1:
1237  return 3./4. - 4*xi + 17./4.*xi*xi
1238  - 5./2.*eta + 7./2.*xi*eta
1239  + 7./4.*eta*eta;
1240  case 2:
1241  return -2*xi + 13./4.*xi*xi
1242  - 1./2.*xi*eta;
1243 
1244  default:
1245  libmesh_error_msg("Invalid subtriangle lookup = " <<
1246  subtriangle_lookup(p));
1247  }
1248  case 6:
1249  switch (subtriangle_lookup(p))
1250  {
1251  case 0:
1252  return -eta + 4*xi*eta
1253  + 7./4.*eta*eta;
1254  case 1:
1255  return -13./4. + 9*xi - 23./4.*xi*xi
1256  + 19./2.*eta - 23./2.*xi*eta
1257  - 25./4.*eta*eta;
1258  case 2:
1259  return -xi + 5./4.*xi*xi
1260  + 9./2.*xi*eta;
1261 
1262  default:
1263  libmesh_error_msg("Invalid subtriangle lookup = " <<
1264  subtriangle_lookup(p));
1265  }
1266  case 7:
1267  switch (subtriangle_lookup(p))
1268  {
1269  case 0:
1270  return 9./4.*eta*eta;
1271  case 1:
1272  return - 11./4. + 7*xi - 17./4.*xi*xi
1273  + 19./2.*eta - 25./2.*xi*eta
1274  - 23./4.*eta*eta;
1275  case 2:
1276  return xi - 13./4.*xi*xi
1277  - eta + 7./2.*xi*eta + 2*eta*eta;
1278 
1279  default:
1280  libmesh_error_msg("Invalid subtriangle lookup = " <<
1281  subtriangle_lookup(p));
1282  }
1283  case 8:
1284  switch (subtriangle_lookup(p))
1285  {
1286  case 0:
1287  return - 1./4.*eta*eta;
1288  case 1:
1289  return 3./4. - 2*xi + 5./4.*xi*xi
1290  - 5./2.*eta + 7./2.*xi*eta
1291  + 7./4.*eta*eta;
1292  case 2:
1293  return 1./4.*xi*xi
1294  - 1./2.*xi*eta;
1295 
1296  default:
1297  libmesh_error_msg("Invalid subtriangle lookup = " <<
1298  subtriangle_lookup(p));
1299  }
1300  case 9:
1301  switch (subtriangle_lookup(p))
1302  {
1303  case 0:
1304  return std::sqrt(2.) * eta*eta;
1305  case 1:
1306  return std::sqrt(2.) * (-3 + 8*xi - 5*xi*xi
1307  + 10*eta - 14*xi*eta
1308  - 7*eta*eta);
1309  case 2:
1310  return std::sqrt(2.) * (-xi*xi + 2*xi*eta);
1311 
1312  default:
1313  libmesh_error_msg("Invalid subtriangle lookup = " <<
1314  subtriangle_lookup(p));
1315  }
1316  case 10:
1317  switch (subtriangle_lookup(p))
1318  {
1319  case 0:
1320  return -2*eta*eta;
1321  case 1:
1322  return 2 - 4*xi + 2*xi*xi
1323  - 8*eta + 8*xi*eta
1324  + 6*eta*eta;
1325  case 2:
1326  return -4*xi + 10*xi*xi
1327  + 4*eta - 8*xi*eta
1328  - 4*eta*eta;
1329 
1330  default:
1331  libmesh_error_msg("Invalid subtriangle lookup = " <<
1332  subtriangle_lookup(p));
1333  }
1334  case 11:
1335  switch (subtriangle_lookup(p))
1336  {
1337  case 0:
1338  return 4*eta - 8*xi*eta
1339  - 4*eta*eta;
1340  case 1:
1341  return 4 - 12*xi + 8*xi*xi
1342  - 8*eta + 12*xi*eta
1343  + 4*eta*eta;
1344  case 2:
1345  return 4*xi - 8*xi*xi
1346  - 4*xi*eta;
1347 
1348  default:
1349  libmesh_error_msg("Invalid subtriangle lookup = " <<
1350  subtriangle_lookup(p));
1351  }
1352 
1353  default:
1354  libmesh_error_msg("Invalid shape function index i = " <<
1355  basis_num);
1356  }
1357 
1358  // derivative in eta direction
1359  case 1:
1360  switch (basis_num)
1361  {
1362  case 0:
1363  switch (subtriangle_lookup(p))
1364  {
1365  case 0:
1366  return - 6*eta - 6*xi*eta + 9*eta*eta;
1367  case 1:
1368  return 9 - 30*xi + 21*xi*xi
1369  - 30*eta + 42*xi*eta + 21*eta*eta;
1370  case 2:
1371  return - 3*xi*xi
1372  - 6*eta + 6*eta*eta;
1373 
1374  default:
1375  libmesh_error_msg("Invalid subtriangle lookup = " <<
1376  subtriangle_lookup(p));
1377  }
1378  case 1:
1379  switch (subtriangle_lookup(p))
1380  {
1381  case 0:
1382  return + 3*xi*eta
1383  - 3./2.*eta*eta;
1384  case 1:
1385  return - 9./2. + 15*xi - 21./2.*xi*xi
1386  + 12*eta - 21*xi*eta - 15./2.*eta*eta;
1387  case 2:
1388  return + 3./2.*xi*xi;
1389 
1390  default:
1391  libmesh_error_msg("Invalid subtriangle lookup = " <<
1392  subtriangle_lookup(p));
1393  }
1394  case 2:
1395  switch (subtriangle_lookup(p))
1396  {
1397  case 0:
1398  return 6*eta + 3*xi*eta - 15./2.*eta*eta;
1399  case 1:
1400  return - 9./2. + 15*xi - 21./2.*xi*xi
1401  + 18*eta - 21.*xi*eta - 27./2.*eta*eta;
1402  case 2:
1403  return 3./2.*xi*xi
1404  + 6*eta - 6*eta*eta;
1405 
1406  default:
1407  libmesh_error_msg("Invalid subtriangle lookup = " <<
1408  subtriangle_lookup(p));
1409  }
1410  case 3:
1411  switch (subtriangle_lookup(p))
1412  {
1413  case 0:
1414  return - 3*eta - xi*eta + 7*eta*eta;
1415  case 1:
1416  return - 4*xi + 4*xi*xi
1417  + eta + 3*xi*eta - eta*eta;
1418  case 2:
1419  return - 3*xi + 2*xi*xi
1420  + 4*xi*eta;
1421 
1422  default:
1423  libmesh_error_msg("Invalid subtriangle lookup = " <<
1424  subtriangle_lookup(p));
1425  }
1426  case 4:
1427  switch (subtriangle_lookup(p))
1428  {
1429  case 0:
1430  return 1 - 3*xi + 2*xi*xi
1431  - eta + 4*xi*eta - 7./2.*eta*eta;
1432  case 1:
1433  return 5./2. - 4*xi + 3./2.*xi*xi
1434  - 9.*eta + 8*xi*eta + 13./2.*eta*eta;
1435  case 2:
1436  return 1 - 1./2.*xi*xi - 4*eta + 3*eta*eta;
1437 
1438  default:
1439  libmesh_error_msg("Invalid subtriangle lookup = " <<
1440  subtriangle_lookup(p));
1441  }
1442  case 5:
1443  switch (subtriangle_lookup(p))
1444  {
1445  case 0:
1446  return - 1./2.*xi*eta + 1./4.*eta*eta;
1447  case 1:
1448  return 3./4. - 5./2.*xi + 7./4.*xi*xi
1449  - 2*eta + 7./2.*xi*eta + 5./4.*eta*eta;
1450  case 2:
1451  return - 1./4.*xi*xi;
1452 
1453  default:
1454  libmesh_error_msg("Invalid subtriangle lookup = " <<
1455  subtriangle_lookup(p));
1456  }
1457  case 6:
1458  switch (subtriangle_lookup(p))
1459  {
1460  case 0:
1461  return -xi + 2*xi*xi
1462  + eta + 7./2.*xi*eta - 13./4.*eta*eta;
1463  case 1:
1464  return - 11./4. + 19./2.*xi - 23./4.*xi*xi
1465  + 7*eta - 25./2.*xi*eta - 17./4.*eta*eta;
1466  case 2:
1467  return 9./4.*xi*xi;
1468 
1469  default:
1470  libmesh_error_msg("Invalid subtriangle lookup = " <<
1471  subtriangle_lookup(p));
1472  }
1473  case 7:
1474  switch (subtriangle_lookup(p))
1475  {
1476  case 0:
1477  return -eta + 9./2.*xi*eta + 5./4.*eta*eta;
1478  case 1:
1479  return - 13./4. + 19./2.*xi - 25./4.*xi*xi
1480  + 9*eta - 23./2.*xi*eta - 23./4.*eta*eta;
1481  case 2:
1482  return - xi + 7./4.*xi*xi + 4*xi*eta;
1483 
1484  default:
1485  libmesh_error_msg("Invalid subtriangle lookup = " <<
1486  subtriangle_lookup(p));
1487  }
1488  case 8:
1489  switch (subtriangle_lookup(p))
1490  {
1491  case 0:
1492  return -2*eta - 1./2.*xi*eta + 13./4.*eta*eta;
1493  case 1:
1494  return 3./4. - 5./2.*xi + 7./4.*xi*xi
1495  - 4*eta + 7./2.*xi*eta + 17./4.*eta*eta;
1496  case 2:
1497  return - 1./4.*xi*xi
1498  - 2*eta + 3*eta*eta;
1499 
1500  default:
1501  libmesh_error_msg("Invalid subtriangle lookup = " <<
1502  subtriangle_lookup(p));
1503  }
1504  case 9:
1505  switch (subtriangle_lookup(p))
1506  {
1507  case 0:
1508  return std::sqrt(2.) * (2*xi*eta - eta*eta);
1509  case 1:
1510  return std::sqrt(2.) * (- 3 + 10*xi - 7*xi*xi
1511  + 8*eta - 14*xi*eta - 5*eta*eta);
1512  case 2:
1513  return std::sqrt(2.) * (xi*xi);
1514 
1515  default:
1516  libmesh_error_msg("Invalid subtriangle lookup = " <<
1517  subtriangle_lookup(p));
1518  }
1519  case 10:
1520  switch (subtriangle_lookup(p))
1521  {
1522  case 0:
1523  return 4*eta - 4*xi*eta - 8*eta*eta;
1524  case 1:
1525  return 4 - 8*xi + 4*xi*xi
1526  - 12*eta + 12*xi*eta + 8*eta*eta;
1527  case 2:
1528  return 4*xi - 4*xi*xi
1529  - 8*xi*eta;
1530 
1531  default:
1532  libmesh_error_msg("Invalid subtriangle lookup = " <<
1533  subtriangle_lookup(p));
1534  }
1535  case 11:
1536  switch (subtriangle_lookup(p))
1537  {
1538  case 0:
1539  return 4*xi - 4*xi*xi
1540  - 4*eta - 8*xi*eta + 10.*eta*eta;
1541  case 1:
1542  return 2 - 8*xi + 6*xi*xi
1543  - 4*eta + 8*xi*eta + 2*eta*eta;
1544  case 2:
1545  return - 2*xi*xi;
1546 
1547  default:
1548  libmesh_error_msg("Invalid subtriangle lookup = " <<
1549  subtriangle_lookup(p));
1550  }
1551 
1552  default:
1553  libmesh_error_msg("Invalid shape function index i = " <<
1554  basis_num);
1555  }
1556 
1557  default:
1558  libmesh_error_msg("Invalid shape function derivative j = " <<
1559  deriv_type);
1560  }
1561 }
1562 
1563 Real clough_raw_shape(const unsigned int basis_num,
1564  const Point & p)
1565 {
1566  Real xi = p(0), eta = p(1);
1567 
1568  switch (basis_num)
1569  {
1570  case 0:
1571  switch (subtriangle_lookup(p))
1572  {
1573  case 0:
1574  return 1 - 3*xi*xi + 2*xi*xi*xi
1575  - 3*eta*eta - 3*xi*eta*eta + 3*eta*eta*eta;
1576  case 1:
1577  return -1 + 9*xi - 15*xi*xi + 7*xi*xi*xi
1578  + 9*eta - 30*xi*eta +21*xi*xi*eta
1579  - 15*eta*eta + 21*xi*eta*eta + 7*eta*eta*eta;
1580  case 2:
1581  return 1 - 3*xi*xi + 3*xi*xi*xi
1582  - 3*xi*xi*eta
1583  - 3*eta*eta + 2*eta*eta*eta;
1584 
1585  default:
1586  libmesh_error_msg("Invalid subtriangle lookup = " <<
1587  subtriangle_lookup(p));
1588  }
1589  case 1:
1590  switch (subtriangle_lookup(p))
1591  {
1592  case 0:
1593  return 3*xi*xi - 2*xi*xi*xi
1594  + 3./2.*xi*eta*eta
1595  - 1./2.*eta*eta*eta;
1596  case 1:
1597  return 1 - 9./2.*xi + 9*xi*xi - 9./2.*xi*xi*xi
1598  - 9./2.*eta + 15*xi*eta - 21./2.*xi*xi*eta
1599  + 6*eta*eta - 21./2.*xi*eta*eta - 5./2.*eta*eta*eta;
1600  case 2:
1601  return 3*xi*xi - 5./2.*xi*xi*xi
1602  + 3./2.*xi*xi*eta;
1603 
1604  default:
1605  libmesh_error_msg("Invalid subtriangle lookup = " <<
1606  subtriangle_lookup(p));
1607  }
1608  case 2:
1609  switch (subtriangle_lookup(p))
1610  {
1611  case 0:
1612  return 3*eta*eta + 3./2.*xi*eta*eta - 5./2.*eta*eta*eta;
1613  case 1:
1614  return 1 - 9./2.*xi + 6*xi*xi - 5./2.*xi*xi*xi
1615  - 9./2.*eta + 15*xi*eta - 21./2.*xi*xi*eta
1616  + 9*eta*eta - 21./2.*xi*eta*eta - 9./2.*eta*eta*eta;
1617  case 2:
1618  return -1./2.*xi*xi*xi
1619  + 3./2.*xi*xi*eta
1620  + 3*eta*eta - 2*eta*eta*eta;
1621 
1622  default:
1623  libmesh_error_msg("Invalid subtriangle lookup = " <<
1624  subtriangle_lookup(p));
1625  }
1626  case 3:
1627  switch (subtriangle_lookup(p))
1628  {
1629  case 0:
1630  return xi - 2*xi*xi + xi*xi*xi
1631  - 3./2.*eta*eta - 1./2.*xi*eta*eta + 7./3.*eta*eta*eta;
1632  case 1:
1633  return -1./6. + 5./2.*xi - 9./2.*xi*xi + 13./6.*xi*xi*xi
1634  - 4*xi*eta + 4*xi*xi*eta
1635  + 1./2.*eta*eta + 3./2.*xi*eta*eta - 1./3.*eta*eta*eta;
1636  case 2:
1637  return xi - 1./2.*xi*xi - 7./6.*xi*xi*xi
1638  - 3*xi*eta + 2*xi*xi*eta
1639  + 2*xi*eta*eta;
1640 
1641  default:
1642  libmesh_error_msg("Invalid subtriangle lookup = " <<
1643  subtriangle_lookup(p));
1644  }
1645  case 4:
1646  switch (subtriangle_lookup(p))
1647  {
1648  case 0:
1649  return eta - 3*xi*eta + 2*xi*xi*eta
1650  - 1./2.*eta*eta + 2*xi*eta*eta - 7./6.*eta*eta*eta;
1651  case 1:
1652  return -1./6. + 1./2.*xi*xi - 1./3.*xi*xi*xi
1653  + 5./2.*eta - 4*xi*eta + 3./2.*xi*xi*eta
1654  - 9./2.*eta*eta + 4*xi*eta*eta + 13./6.*eta*eta*eta;
1655  case 2:
1656  return -3./2.*xi*xi + 7./3.*xi*xi*xi
1657  + eta - 1./2.*xi*xi*eta - 2*eta*eta + eta*eta*eta;
1658 
1659  default:
1660  libmesh_error_msg("Invalid subtriangle lookup = " <<
1661  subtriangle_lookup(p));
1662  }
1663  case 5:
1664  switch (subtriangle_lookup(p))
1665  {
1666  case 0:
1667  return -xi*xi + xi*xi*xi
1668  - 1./4.*xi*eta*eta + 1./12.*eta*eta*eta;
1669  case 1:
1670  return -1./6. + 3./4.*xi - 2*xi*xi + 17./12.*xi*xi*xi
1671  + 3./4.*eta - 5./2.*xi*eta + 7./4.*xi*xi*eta
1672  - eta*eta + 7./4.*xi*eta*eta + 5./12.*eta*eta*eta;
1673  case 2:
1674  return -xi*xi + 13./12.*xi*xi*xi
1675  - 1./4.*xi*xi*eta;
1676 
1677  default:
1678  libmesh_error_msg("Invalid subtriangle lookup = " <<
1679  subtriangle_lookup(p));
1680  }
1681  case 6:
1682  switch (subtriangle_lookup(p))
1683  {
1684  case 0:
1685  return -xi*eta + 2*xi*xi*eta
1686  + 1./2.*eta*eta + 7./4.*xi*eta*eta - 13./12.*eta*eta*eta;
1687  case 1:
1688  return 2./3. - 13./4.*xi + 9./2.*xi*xi - 23./12.*xi*xi*xi
1689  - 11./4.*eta + 19./2.*xi*eta - 23./4.*xi*xi*eta
1690  + 7./2.*eta*eta - 25./4.*xi*eta*eta - 17./12.*eta*eta*eta;
1691  case 2:
1692  return -1./2.*xi*xi + 5./12.*xi*xi*xi
1693  + 9./4.*xi*xi*eta;
1694 
1695  default:
1696  libmesh_error_msg("Invalid subtriangle lookup = " <<
1697  subtriangle_lookup(p));
1698  }
1699  case 7:
1700  switch (subtriangle_lookup(p))
1701  {
1702  case 0:
1703  return -1./2.*eta*eta + 9./4.*xi*eta*eta + 5./12.*eta*eta*eta;
1704  case 1:
1705  return 2./3. - 11./4.*xi + 7./2.*xi*xi - 17./12.*xi*xi*xi
1706  - 13./4.*eta + 19./2.*xi*eta - 25./4.*xi*xi*eta
1707  + 9./2.*eta*eta - 23./4.*xi*eta*eta - 23./12.*eta*eta*eta;
1708  case 2:
1709  return 1./2.*xi*xi - 13./12.*xi*xi*xi
1710  - xi*eta + 7./4.*xi*xi*eta + 2*xi*eta*eta;
1711 
1712  default:
1713  libmesh_error_msg("Invalid subtriangle lookup = " <<
1714  subtriangle_lookup(p));
1715  }
1716  case 8:
1717  switch (subtriangle_lookup(p))
1718  {
1719  case 0:
1720  return -eta*eta - 1./4.*xi*eta*eta + 13./12.*eta*eta*eta;
1721  case 1:
1722  return -1./6. + 3./4.*xi - xi*xi + 5./12.*xi*xi*xi
1723  + 3./4.*eta - 5./2.*xi*eta + 7./4.*xi*xi*eta
1724  - 2*eta*eta + 7./4.*xi*eta*eta + 17./12.*eta*eta*eta;
1725  case 2:
1726  return 1./12.*xi*xi*xi
1727  - 1./4.*xi*xi*eta
1728  - eta*eta + eta*eta*eta;
1729 
1730  default:
1731  libmesh_error_msg("Invalid subtriangle lookup = " <<
1732  subtriangle_lookup(p));
1733  }
1734  case 9:
1735  switch (subtriangle_lookup(p))
1736  {
1737  case 0:
1738  return std::sqrt(2.) * (xi*eta*eta - 1./3.*eta*eta*eta);
1739  case 1:
1740  return std::sqrt(2.) * (2./3. - 3*xi + 4*xi*xi - 5./3.*xi*xi*xi
1741  - 3*eta + 10*xi*eta - 7*xi*xi*eta
1742  + 4*eta*eta - 7*xi*eta*eta - 5./3.*eta*eta*eta);
1743  case 2:
1744  return std::sqrt(2.) * (-1./3.*xi*xi*xi + xi*xi*eta);
1745 
1746  default:
1747  libmesh_error_msg("Invalid subtriangle lookup = " <<
1748  subtriangle_lookup(p));
1749  }
1750  case 10:
1751  switch (subtriangle_lookup(p))
1752  {
1753  case 0:
1754  return 2*eta*eta - 2*xi*eta*eta - 8./3.*eta*eta*eta;
1755  case 1:
1756  return -2./3. + 2*xi - 2*xi*xi + 2./3.*xi*xi*xi
1757  + 4*eta - 8*xi*eta + 4*xi*xi*eta
1758  - 6*eta*eta + 6*xi*eta*eta + 8./3.*eta*eta*eta;
1759  case 2:
1760  return -2*xi*xi + 10./3.*xi*xi*xi
1761  + 4*xi*eta - 4*xi*xi*eta
1762  - 4*xi*eta*eta;
1763 
1764  default:
1765  libmesh_error_msg("Invalid subtriangle lookup = " <<
1766  subtriangle_lookup(p));
1767  }
1768  case 11:
1769  switch (subtriangle_lookup(p))
1770  {
1771  case 0:
1772  return 4*xi*eta - 4*xi*xi*eta
1773  - 2*eta*eta - 4*xi*eta*eta + 10./3.*eta*eta*eta;
1774  case 1:
1775  return -2./3. + 4*xi - 6*xi*xi + 8./3.*xi*xi*xi
1776  + 2*eta - 8*xi*eta + 6*xi*xi*eta
1777  - 2*eta*eta + 4*xi*eta*eta + 2./3.*eta*eta*eta;
1778  case 2:
1779  return 2*xi*xi - 8./3.*xi*xi*xi
1780  - 2*xi*xi*eta;
1781 
1782  default:
1783  libmesh_error_msg("Invalid subtriangle lookup = " <<
1784  subtriangle_lookup(p));
1785  }
1786 
1787  default:
1788  libmesh_error_msg("Invalid shape function index i = " <<
1789  basis_num);
1790  }
1791 }
1792 
1793 
1794 } // end anonymous namespace
1795 
1796 
1797 namespace libMesh
1798 {
1799 
1800 
1801 template <>
1802 Real FE<2,CLOUGH>::shape(const ElemType,
1803  const Order,
1804  const unsigned int,
1805  const Point &)
1806 {
1807  libmesh_error_msg("Clough-Tocher elements require the real element \nto construct gradient-based degrees of freedom.");
1808  return 0.;
1809 }
1810 
1811 
1812 
1813 template <>
1814 Real FE<2,CLOUGH>::shape(const Elem * elem,
1815  const Order order,
1816  const unsigned int i,
1817  const Point & p)
1818 {
1819  libmesh_assert(elem);
1820 
1821  clough_compute_coefs(elem);
1822 
1823  const ElemType type = elem->type();
1824 
1825  const Order totalorder = static_cast<Order>(order + elem->p_level());
1826 
1827  switch (totalorder)
1828  {
1829  // 2nd-order restricted Clough-Tocher element
1830  case SECOND:
1831  {
1832  // There may be a bug in the 2nd order case; the 3rd order
1833  // Clough-Tocher elements are pretty uniformly better anyways
1834  // so use those instead.
1835  libmesh_experimental();
1836  switch (type)
1837  {
1838  // C1 functions on the Clough-Tocher triangle.
1839  case TRI6:
1840  {
1841  libmesh_assert_less (i, 9);
1842  // FIXME: it would be nice to calculate (and cache)
1843  // clough_raw_shape(j,p) only once per triangle, not 1-7
1844  // times
1845  switch (i)
1846  {
1847  // Note: these DoF numbers are "scrambled" because my
1848  // initial numbering conventions didn't match libMesh
1849  case 0:
1850  return clough_raw_shape(0, p)
1851  + d1d2n * clough_raw_shape(10, p)
1852  + d1d3n * clough_raw_shape(11, p);
1853  case 3:
1854  return clough_raw_shape(1, p)
1855  + d2d3n * clough_raw_shape(11, p)
1856  + d2d1n * clough_raw_shape(9, p);
1857  case 6:
1858  return clough_raw_shape(2, p)
1859  + d3d1n * clough_raw_shape(9, p)
1860  + d3d2n * clough_raw_shape(10, p);
1861  case 1:
1862  return d1xd1x * clough_raw_shape(3, p)
1863  + d1xd1y * clough_raw_shape(4, p)
1864  + d1xd2n * clough_raw_shape(10, p)
1865  + d1xd3n * clough_raw_shape(11, p)
1866  + 0.5 * N01x * d3nd3n * clough_raw_shape(11, p)
1867  + 0.5 * N02x * d2nd2n * clough_raw_shape(10, p);
1868  case 2:
1869  return d1yd1y * clough_raw_shape(4, p)
1870  + d1yd1x * clough_raw_shape(3, p)
1871  + d1yd2n * clough_raw_shape(10, p)
1872  + d1yd3n * clough_raw_shape(11, p)
1873  + 0.5 * N01y * d3nd3n * clough_raw_shape(11, p)
1874  + 0.5 * N02y * d2nd2n * clough_raw_shape(10, p);
1875  case 4:
1876  return d2xd2x * clough_raw_shape(5, p)
1877  + d2xd2y * clough_raw_shape(6, p)
1878  + d2xd3n * clough_raw_shape(11, p)
1879  + d2xd1n * clough_raw_shape(9, p)
1880  + 0.5 * N10x * d3nd3n * clough_raw_shape(11, p)
1881  + 0.5 * N12x * d1nd1n * clough_raw_shape(9, p);
1882  case 5:
1883  return d2yd2y * clough_raw_shape(6, p)
1884  + d2yd2x * clough_raw_shape(5, p)
1885  + d2yd3n * clough_raw_shape(11, p)
1886  + d2yd1n * clough_raw_shape(9, p)
1887  + 0.5 * N10y * d3nd3n * clough_raw_shape(11, p)
1888  + 0.5 * N12y * d1nd1n * clough_raw_shape(9, p);
1889  case 7:
1890  return d3xd3x * clough_raw_shape(7, p)
1891  + d3xd3y * clough_raw_shape(8, p)
1892  + d3xd1n * clough_raw_shape(9, p)
1893  + d3xd2n * clough_raw_shape(10, p)
1894  + 0.5 * N20x * d2nd2n * clough_raw_shape(10, p)
1895  + 0.5 * N21x * d1nd1n * clough_raw_shape(9, p);
1896  case 8:
1897  return d3yd3y * clough_raw_shape(8, p)
1898  + d3yd3x * clough_raw_shape(7, p)
1899  + d3yd1n * clough_raw_shape(9, p)
1900  + d3yd2n * clough_raw_shape(10, p)
1901  + 0.5 * N20y * d2nd2n * clough_raw_shape(10, p)
1902  + 0.5 * N21y * d1nd1n * clough_raw_shape(9, p);
1903  default:
1904  libmesh_error_msg("Invalid shape function index i = " << i);
1905  }
1906  }
1907  default:
1908  libmesh_error_msg("ERROR: Unsupported element type = " << type);
1909  }
1910  }
1911  // 3rd-order Clough-Tocher element
1912  case THIRD:
1913  {
1914  switch (type)
1915  {
1916  // C1 functions on the Clough-Tocher triangle.
1917  case TRI6:
1918  {
1919  libmesh_assert_less (i, 12);
1920 
1921  // FIXME: it would be nice to calculate (and cache)
1922  // clough_raw_shape(j,p) only once per triangle, not 1-7
1923  // times
1924  switch (i)
1925  {
1926  // Note: these DoF numbers are "scrambled" because my
1927  // initial numbering conventions didn't match libMesh
1928  case 0:
1929  return clough_raw_shape(0, p)
1930  + d1d2n * clough_raw_shape(10, p)
1931  + d1d3n * clough_raw_shape(11, p);
1932  case 3:
1933  return clough_raw_shape(1, p)
1934  + d2d3n * clough_raw_shape(11, p)
1935  + d2d1n * clough_raw_shape(9, p);
1936  case 6:
1937  return clough_raw_shape(2, p)
1938  + d3d1n * clough_raw_shape(9, p)
1939  + d3d2n * clough_raw_shape(10, p);
1940  case 1:
1941  return d1xd1x * clough_raw_shape(3, p)
1942  + d1xd1y * clough_raw_shape(4, p)
1943  + d1xd2n * clough_raw_shape(10, p)
1944  + d1xd3n * clough_raw_shape(11, p);
1945  case 2:
1946  return d1yd1y * clough_raw_shape(4, p)
1947  + d1yd1x * clough_raw_shape(3, p)
1948  + d1yd2n * clough_raw_shape(10, p)
1949  + d1yd3n * clough_raw_shape(11, p);
1950  case 4:
1951  return d2xd2x * clough_raw_shape(5, p)
1952  + d2xd2y * clough_raw_shape(6, p)
1953  + d2xd3n * clough_raw_shape(11, p)
1954  + d2xd1n * clough_raw_shape(9, p);
1955  case 5:
1956  return d2yd2y * clough_raw_shape(6, p)
1957  + d2yd2x * clough_raw_shape(5, p)
1958  + d2yd3n * clough_raw_shape(11, p)
1959  + d2yd1n * clough_raw_shape(9, p);
1960  case 7:
1961  return d3xd3x * clough_raw_shape(7, p)
1962  + d3xd3y * clough_raw_shape(8, p)
1963  + d3xd1n * clough_raw_shape(9, p)
1964  + d3xd2n * clough_raw_shape(10, p);
1965  case 8:
1966  return d3yd3y * clough_raw_shape(8, p)
1967  + d3yd3x * clough_raw_shape(7, p)
1968  + d3yd1n * clough_raw_shape(9, p)
1969  + d3yd2n * clough_raw_shape(10, p);
1970  case 10:
1971  return d1nd1n * clough_raw_shape(9, p);
1972  case 11:
1973  return d2nd2n * clough_raw_shape(10, p);
1974  case 9:
1975  return d3nd3n * clough_raw_shape(11, p);
1976 
1977  default:
1978  libmesh_error_msg("Invalid shape function index i = " << i);
1979  }
1980  }
1981  default:
1982  libmesh_error_msg("ERROR: Unsupported element type = " << type);
1983  }
1984  }
1985  // by default throw an error
1986  default:
1987  libmesh_error_msg("ERROR: Unsupported polynomial order = " << order);
1988  }
1989 }
1990 
1991 
1992 
1993 template <>
1994 Real FE<2,CLOUGH>::shape_deriv(const ElemType,
1995  const Order,
1996  const unsigned int,
1997  const unsigned int,
1998  const Point &)
1999 {
2000  libmesh_error_msg("Clough-Tocher elements require the real element \nto construct gradient-based degrees of freedom.");
2001  return 0.;
2002 }
2003 
2004 
2005 
2006 template <>
2007 Real FE<2,CLOUGH>::shape_deriv(const Elem * elem,
2008  const Order order,
2009  const unsigned int i,
2010  const unsigned int j,
2011  const Point & p)
2012 {
2013  libmesh_assert(elem);
2014 
2015  clough_compute_coefs(elem);
2016 
2017  const ElemType type = elem->type();
2018 
2019  const Order totalorder = static_cast<Order>(order + elem->p_level());
2020 
2021  switch (totalorder)
2022  {
2023  // 2nd-order restricted Clough-Tocher element
2024  case SECOND:
2025  {
2026  // There may be a bug in the 2nd order case; the 3rd order
2027  // Clough-Tocher elements are pretty uniformly better anyways
2028  // so use those instead.
2029  libmesh_experimental();
2030  switch (type)
2031  {
2032  // C1 functions on the Clough-Tocher triangle.
2033  case TRI6:
2034  {
2035  libmesh_assert_less (i, 9);
2036  // FIXME: it would be nice to calculate (and cache)
2037  // clough_raw_shape(j,p) only once per triangle, not 1-7
2038  // times
2039  switch (i)
2040  {
2041  // Note: these DoF numbers are "scrambled" because my
2042  // initial numbering conventions didn't match libMesh
2043  case 0:
2044  return clough_raw_shape_deriv(0, j, p)
2045  + d1d2n * clough_raw_shape_deriv(10, j, p)
2046  + d1d3n * clough_raw_shape_deriv(11, j, p);
2047  case 3:
2048  return clough_raw_shape_deriv(1, j, p)
2049  + d2d3n * clough_raw_shape_deriv(11, j, p)
2050  + d2d1n * clough_raw_shape_deriv(9, j, p);
2051  case 6:
2052  return clough_raw_shape_deriv(2, j, p)
2053  + d3d1n * clough_raw_shape_deriv(9, j, p)
2054  + d3d2n * clough_raw_shape_deriv(10, j, p);
2055  case 1:
2056  return d1xd1x * clough_raw_shape_deriv(3, j, p)
2057  + d1xd1y * clough_raw_shape_deriv(4, j, p)
2058  + d1xd2n * clough_raw_shape_deriv(10, j, p)
2059  + d1xd3n * clough_raw_shape_deriv(11, j, p)
2060  + 0.5 * N01x * d3nd3n * clough_raw_shape_deriv(11, j, p)
2061  + 0.5 * N02x * d2nd2n * clough_raw_shape_deriv(10, j, p);
2062  case 2:
2063  return d1yd1y * clough_raw_shape_deriv(4, j, p)
2064  + d1yd1x * clough_raw_shape_deriv(3, j, p)
2065  + d1yd2n * clough_raw_shape_deriv(10, j, p)
2066  + d1yd3n * clough_raw_shape_deriv(11, j, p)
2067  + 0.5 * N01y * d3nd3n * clough_raw_shape_deriv(11, j, p)
2068  + 0.5 * N02y * d2nd2n * clough_raw_shape_deriv(10, j, p);
2069  case 4:
2070  return d2xd2x * clough_raw_shape_deriv(5, j, p)
2071  + d2xd2y * clough_raw_shape_deriv(6, j, p)
2072  + d2xd3n * clough_raw_shape_deriv(11, j, p)
2073  + d2xd1n * clough_raw_shape_deriv(9, j, p)
2074  + 0.5 * N10x * d3nd3n * clough_raw_shape_deriv(11, j, p)
2075  + 0.5 * N12x * d1nd1n * clough_raw_shape_deriv(9, j, p);
2076  case 5:
2077  return d2yd2y * clough_raw_shape_deriv(6, j, p)
2078  + d2yd2x * clough_raw_shape_deriv(5, j, p)
2079  + d2yd3n * clough_raw_shape_deriv(11, j, p)
2080  + d2yd1n * clough_raw_shape_deriv(9, j, p)
2081  + 0.5 * N10y * d3nd3n * clough_raw_shape_deriv(11, j, p)
2082  + 0.5 * N12y * d1nd1n * clough_raw_shape_deriv(9, j, p);
2083  case 7:
2084  return d3xd3x * clough_raw_shape_deriv(7, j, p)
2085  + d3xd3y * clough_raw_shape_deriv(8, j, p)
2086  + d3xd1n * clough_raw_shape_deriv(9, j, p)
2087  + d3xd2n * clough_raw_shape_deriv(10, j, p)
2088  + 0.5 * N20x * d2nd2n * clough_raw_shape_deriv(10, j, p)
2089  + 0.5 * N21x * d1nd1n * clough_raw_shape_deriv(9, j, p);
2090  case 8:
2091  return d3yd3y * clough_raw_shape_deriv(8, j, p)
2092  + d3yd3x * clough_raw_shape_deriv(7, j, p)
2093  + d3yd1n * clough_raw_shape_deriv(9, j, p)
2094  + d3yd2n * clough_raw_shape_deriv(10, j, p)
2095  + 0.5 * N20y * d2nd2n * clough_raw_shape_deriv(10, j, p)
2096  + 0.5 * N21y * d1nd1n * clough_raw_shape_deriv(9, j, p);
2097  default:
2098  libmesh_error_msg("Invalid shape function index i = " << i);
2099  }
2100  }
2101  default:
2102  libmesh_error_msg("ERROR: Unsupported element type = " << type);
2103  }
2104  }
2105  // 3rd-order Clough-Tocher element
2106  case THIRD:
2107  {
2108  switch (type)
2109  {
2110  // C1 functions on the Clough-Tocher triangle.
2111  case TRI6:
2112  {
2113  libmesh_assert_less (i, 12);
2114 
2115  // FIXME: it would be nice to calculate (and cache)
2116  // clough_raw_shape(j,p) only once per triangle, not 1-7
2117  // times
2118  switch (i)
2119  {
2120  // Note: these DoF numbers are "scrambled" because my
2121  // initial numbering conventions didn't match libMesh
2122  case 0:
2123  return clough_raw_shape_deriv(0, j, p)
2124  + d1d2n * clough_raw_shape_deriv(10, j, p)
2125  + d1d3n * clough_raw_shape_deriv(11, j, p);
2126  case 3:
2127  return clough_raw_shape_deriv(1, j, p)
2128  + d2d3n * clough_raw_shape_deriv(11, j, p)
2129  + d2d1n * clough_raw_shape_deriv(9, j, p);
2130  case 6:
2131  return clough_raw_shape_deriv(2, j, p)
2132  + d3d1n * clough_raw_shape_deriv(9, j, p)
2133  + d3d2n * clough_raw_shape_deriv(10, j, p);
2134  case 1:
2135  return d1xd1x * clough_raw_shape_deriv(3, j, p)
2136  + d1xd1y * clough_raw_shape_deriv(4, j, p)
2137  + d1xd2n * clough_raw_shape_deriv(10, j, p)
2138  + d1xd3n * clough_raw_shape_deriv(11, j, p);
2139  case 2:
2140  return d1yd1y * clough_raw_shape_deriv(4, j, p)
2141  + d1yd1x * clough_raw_shape_deriv(3, j, p)
2142  + d1yd2n * clough_raw_shape_deriv(10, j, p)
2143  + d1yd3n * clough_raw_shape_deriv(11, j, p);
2144  case 4:
2145  return d2xd2x * clough_raw_shape_deriv(5, j, p)
2146  + d2xd2y * clough_raw_shape_deriv(6, j, p)
2147  + d2xd3n * clough_raw_shape_deriv(11, j, p)
2148  + d2xd1n * clough_raw_shape_deriv(9, j, p);
2149  case 5:
2150  return d2yd2y * clough_raw_shape_deriv(6, j, p)
2151  + d2yd2x * clough_raw_shape_deriv(5, j, p)
2152  + d2yd3n * clough_raw_shape_deriv(11, j, p)
2153  + d2yd1n * clough_raw_shape_deriv(9, j, p);
2154  case 7:
2155  return d3xd3x * clough_raw_shape_deriv(7, j, p)
2156  + d3xd3y * clough_raw_shape_deriv(8, j, p)
2157  + d3xd1n * clough_raw_shape_deriv(9, j, p)
2158  + d3xd2n * clough_raw_shape_deriv(10, j, p);
2159  case 8:
2160  return d3yd3y * clough_raw_shape_deriv(8, j, p)
2161  + d3yd3x * clough_raw_shape_deriv(7, j, p)
2162  + d3yd1n * clough_raw_shape_deriv(9, j, p)
2163  + d3yd2n * clough_raw_shape_deriv(10, j, p);
2164  case 10:
2165  return d1nd1n * clough_raw_shape_deriv(9, j, p);
2166  case 11:
2167  return d2nd2n * clough_raw_shape_deriv(10, j, p);
2168  case 9:
2169  return d3nd3n * clough_raw_shape_deriv(11, j, p);
2170 
2171  default:
2172  libmesh_error_msg("Invalid shape function index i = " << i);
2173  }
2174  }
2175  default:
2176  libmesh_error_msg("ERROR: Unsupported element type = " << type);
2177  }
2178  }
2179  // by default throw an error
2180  default:
2181  libmesh_error_msg("ERROR: Unsupported polynomial order = " << order);
2182  }
2183 }
2184 
2185 
2186 
2187 template <>
2188 Real FE<2,CLOUGH>::shape_second_deriv(const Elem * elem,
2189  const Order order,
2190  const unsigned int i,
2191  const unsigned int j,
2192  const Point & p)
2193 {
2194  libmesh_assert(elem);
2195 
2196  clough_compute_coefs(elem);
2197 
2198  const ElemType type = elem->type();
2199 
2200  const Order totalorder = static_cast<Order>(order + elem->p_level());
2201 
2202  switch (totalorder)
2203  {
2204  // 2nd-order restricted Clough-Tocher element
2205  case SECOND:
2206  {
2207  switch (type)
2208  {
2209  // C1 functions on the Clough-Tocher triangle.
2210  case TRI6:
2211  {
2212  libmesh_assert_less (i, 9);
2213  // FIXME: it would be nice to calculate (and cache)
2214  // clough_raw_shape(j,p) only once per triangle, not 1-7
2215  // times
2216  switch (i)
2217  {
2218  // Note: these DoF numbers are "scrambled" because my
2219  // initial numbering conventions didn't match libMesh
2220  case 0:
2221  return clough_raw_shape_second_deriv(0, j, p)
2222  + d1d2n * clough_raw_shape_second_deriv(10, j, p)
2223  + d1d3n * clough_raw_shape_second_deriv(11, j, p);
2224  case 3:
2225  return clough_raw_shape_second_deriv(1, j, p)
2226  + d2d3n * clough_raw_shape_second_deriv(11, j, p)
2227  + d2d1n * clough_raw_shape_second_deriv(9, j, p);
2228  case 6:
2229  return clough_raw_shape_second_deriv(2, j, p)
2230  + d3d1n * clough_raw_shape_second_deriv(9, j, p)
2231  + d3d2n * clough_raw_shape_second_deriv(10, j, p);
2232  case 1:
2233  return d1xd1x * clough_raw_shape_second_deriv(3, j, p)
2234  + d1xd1y * clough_raw_shape_second_deriv(4, j, p)
2235  + d1xd2n * clough_raw_shape_second_deriv(10, j, p)
2236  + d1xd3n * clough_raw_shape_second_deriv(11, j, p)
2237  + 0.5 * N01x * d3nd3n * clough_raw_shape_second_deriv(11, j, p)
2238  + 0.5 * N02x * d2nd2n * clough_raw_shape_second_deriv(10, j, p);
2239  case 2:
2240  return d1yd1y * clough_raw_shape_second_deriv(4, j, p)
2241  + d1yd1x * clough_raw_shape_second_deriv(3, j, p)
2242  + d1yd2n * clough_raw_shape_second_deriv(10, j, p)
2243  + d1yd3n * clough_raw_shape_second_deriv(11, j, p)
2244  + 0.5 * N01y * d3nd3n * clough_raw_shape_second_deriv(11, j, p)
2245  + 0.5 * N02y * d2nd2n * clough_raw_shape_second_deriv(10, j, p);
2246  case 4:
2247  return d2xd2x * clough_raw_shape_second_deriv(5, j, p)
2248  + d2xd2y * clough_raw_shape_second_deriv(6, j, p)
2249  + d2xd3n * clough_raw_shape_second_deriv(11, j, p)
2250  + d2xd1n * clough_raw_shape_second_deriv(9, j, p)
2251  + 0.5 * N10x * d3nd3n * clough_raw_shape_second_deriv(11, j, p)
2252  + 0.5 * N12x * d1nd1n * clough_raw_shape_second_deriv(9, j, p);
2253  case 5:
2254  return d2yd2y * clough_raw_shape_second_deriv(6, j, p)
2255  + d2yd2x * clough_raw_shape_second_deriv(5, j, p)
2256  + d2yd3n * clough_raw_shape_second_deriv(11, j, p)
2257  + d2yd1n * clough_raw_shape_second_deriv(9, j, p)
2258  + 0.5 * N10y * d3nd3n * clough_raw_shape_second_deriv(11, j, p)
2259  + 0.5 * N12y * d1nd1n * clough_raw_shape_second_deriv(9, j, p);
2260  case 7:
2261  return d3xd3x * clough_raw_shape_second_deriv(7, j, p)
2262  + d3xd3y * clough_raw_shape_second_deriv(8, j, p)
2263  + d3xd1n * clough_raw_shape_second_deriv(9, j, p)
2264  + d3xd2n * clough_raw_shape_second_deriv(10, j, p)
2265  + 0.5 * N20x * d2nd2n * clough_raw_shape_second_deriv(10, j, p)
2266  + 0.5 * N21x * d1nd1n * clough_raw_shape_second_deriv(9, j, p);
2267  case 8:
2268  return d3yd3y * clough_raw_shape_second_deriv(8, j, p)
2269  + d3yd3x * clough_raw_shape_second_deriv(7, j, p)
2270  + d3yd1n * clough_raw_shape_second_deriv(9, j, p)
2271  + d3yd2n * clough_raw_shape_second_deriv(10, j, p)
2272  + 0.5 * N20y * d2nd2n * clough_raw_shape_second_deriv(10, j, p)
2273  + 0.5 * N21y * d1nd1n * clough_raw_shape_second_deriv(9, j, p);
2274  default:
2275  libmesh_error_msg("Invalid shape function index i = " << i);
2276  }
2277  }
2278  default:
2279  libmesh_error_msg("ERROR: Unsupported element type = " << type);
2280  }
2281  }
2282  // 3rd-order Clough-Tocher element
2283  case THIRD:
2284  {
2285  switch (type)
2286  {
2287  // C1 functions on the Clough-Tocher triangle.
2288  case TRI6:
2289  {
2290  libmesh_assert_less (i, 12);
2291 
2292  // FIXME: it would be nice to calculate (and cache)
2293  // clough_raw_shape(j,p) only once per triangle, not 1-7
2294  // times
2295  switch (i)
2296  {
2297  // Note: these DoF numbers are "scrambled" because my
2298  // initial numbering conventions didn't match libMesh
2299  case 0:
2300  return clough_raw_shape_second_deriv(0, j, p)
2301  + d1d2n * clough_raw_shape_second_deriv(10, j, p)
2302  + d1d3n * clough_raw_shape_second_deriv(11, j, p);
2303  case 3:
2304  return clough_raw_shape_second_deriv(1, j, p)
2305  + d2d3n * clough_raw_shape_second_deriv(11, j, p)
2306  + d2d1n * clough_raw_shape_second_deriv(9, j, p);
2307  case 6:
2308  return clough_raw_shape_second_deriv(2, j, p)
2309  + d3d1n * clough_raw_shape_second_deriv(9, j, p)
2310  + d3d2n * clough_raw_shape_second_deriv(10, j, p);
2311  case 1:
2312  return d1xd1x * clough_raw_shape_second_deriv(3, j, p)
2313  + d1xd1y * clough_raw_shape_second_deriv(4, j, p)
2314  + d1xd2n * clough_raw_shape_second_deriv(10, j, p)
2315  + d1xd3n * clough_raw_shape_second_deriv(11, j, p);
2316  case 2:
2317  return d1yd1y * clough_raw_shape_second_deriv(4, j, p)
2318  + d1yd1x * clough_raw_shape_second_deriv(3, j, p)
2319  + d1yd2n * clough_raw_shape_second_deriv(10, j, p)
2320  + d1yd3n * clough_raw_shape_second_deriv(11, j, p);
2321  case 4:
2322  return d2xd2x * clough_raw_shape_second_deriv(5, j, p)
2323  + d2xd2y * clough_raw_shape_second_deriv(6, j, p)
2324  + d2xd3n * clough_raw_shape_second_deriv(11, j, p)
2325  + d2xd1n * clough_raw_shape_second_deriv(9, j, p);
2326  case 5:
2327  return d2yd2y * clough_raw_shape_second_deriv(6, j, p)
2328  + d2yd2x * clough_raw_shape_second_deriv(5, j, p)
2329  + d2yd3n * clough_raw_shape_second_deriv(11, j, p)
2330  + d2yd1n * clough_raw_shape_second_deriv(9, j, p);
2331  case 7:
2332  return d3xd3x * clough_raw_shape_second_deriv(7, j, p)
2333  + d3xd3y * clough_raw_shape_second_deriv(8, j, p)
2334  + d3xd1n * clough_raw_shape_second_deriv(9, j, p)
2335  + d3xd2n * clough_raw_shape_second_deriv(10, j, p);
2336  case 8:
2337  return d3yd3y * clough_raw_shape_second_deriv(8, j, p)
2338  + d3yd3x * clough_raw_shape_second_deriv(7, j, p)
2339  + d3yd1n * clough_raw_shape_second_deriv(9, j, p)
2340  + d3yd2n * clough_raw_shape_second_deriv(10, j, p);
2341  case 10:
2342  return d1nd1n * clough_raw_shape_second_deriv(9, j, p);
2343  case 11:
2344  return d2nd2n * clough_raw_shape_second_deriv(10, j, p);
2345  case 9:
2346  return d3nd3n * clough_raw_shape_second_deriv(11, j, p);
2347 
2348  default:
2349  libmesh_error_msg("Invalid shape function index i = " << i);
2350  }
2351  }
2352  default:
2353  libmesh_error_msg("ERROR: Unsupported element type = " << type);
2354  }
2355  }
2356  // by default throw an error
2357  default:
2358  libmesh_error_msg("ERROR: Unsupported polynomial order = " << order);
2359  }
2360 }
2361 
2362 } // namespace libMesh
libMesh::n_threads
unsigned int n_threads()
Definition: libmesh_base.h:96
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::Elem
The base class for all geometric element types.
Definition: elem.h:100
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::Elem::p_level
unsigned int p_level() const
Definition: elem.h:2555
libMesh::FE::n_shape_functions
virtual unsigned int n_shape_functions() const override
Definition: fe.C:50
libMesh::DofObject::id
dof_id_type id() const
Definition: dof_object.h:655
libMesh::DofObject::invalid_id
static const dof_id_type invalid_id
Definition: dof_object.h:347
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Point
A geometric point in (x,y,z) space.
Definition: point.h:38
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:64
generated by doxygen