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