Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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 48347136 : void clough_compute_coefs(const Elem * elem, CloughCoefs & coefs)
69 : {
70 48347136 : const FEFamily mapping_family = FEMap::map_fe_type(*elem);
71 48347136 : 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 48347136 : FEInterface::n_shape_functions(map_fe_type, /*extra_order=*/0, elem);
77 :
78 : // Degrees of freedom are at vertices and edge midpoints
79 24317280 : std::vector<Point> dofpt;
80 48347136 : dofpt.push_back(Point(0,0));
81 48347136 : dofpt.push_back(Point(1,0));
82 48347136 : dofpt.push_back(Point(0,1));
83 48347136 : dofpt.push_back(Point(1/2.,1/2.));
84 48347136 : dofpt.push_back(Point(0,1/2.));
85 48347136 : dofpt.push_back(Point(1/2.,0));
86 :
87 : // Mapping functions - first derivatives at each dofpt
88 96981696 : std::vector<Real> dxdxi(6), dxdeta(6), dydxi(6), dydeta(6);
89 96981696 : std::vector<Real> dxidx(6), detadx(6), dxidy(6), detady(6);
90 :
91 : FEInterface::shape_deriv_ptr shape_deriv_ptr =
92 48347136 : FEInterface::shape_deriv_function(map_fe_type, elem);
93 :
94 338429952 : for (int p = 0; p != 6; ++p)
95 : {
96 : // libMesh::err << p << ' ' << dofpt[p];
97 2031510672 : for (int i = 0; i != n_mapping_shape_functions; ++i)
98 : {
99 : const Real ddxi = shape_deriv_ptr
100 2179397664 : (map_fe_type, elem, i, 0, dofpt[p], /*add_p_level=*/false);
101 : const Real ddeta = shape_deriv_ptr
102 2179397664 : (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 2179397664 : dxdxi[p] += elem->point(i)(0) * ddxi;
108 1741427856 : dydxi[p] += elem->point(i)(1) * ddxi;
109 1741427856 : dxdeta[p] += elem->point(i)(0) * ddeta;
110 2179397664 : 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 363034656 : const Real inv_jac = 1. / (dxdxi[p]*dydeta[p] -
124 363034656 : dxdeta[p]*dydxi[p]);
125 290082816 : dxidx[p] = dydeta[p] * inv_jac;
126 290082816 : dxidy[p] = - dxdeta[p] * inv_jac;
127 290082816 : detadx[p] = - dydxi[p] * inv_jac;
128 363034656 : detady[p] = dxdxi[p] * inv_jac;
129 : }
130 :
131 : // Calculate midpoint normal vectors
132 48347136 : Real N1x = dydeta[3] - dydxi[3];
133 48347136 : Real N1y = dxdxi[3] - dxdeta[3];
134 48347136 : Real Nlength = std::sqrt(static_cast<Real>(N1x*N1x + N1y*N1y));
135 48347136 : N1x /= Nlength; N1y /= Nlength;
136 :
137 48347136 : Real N2x = - dydeta[4];
138 48347136 : Real N2y = dxdeta[4];
139 48347136 : Nlength = std::sqrt(static_cast<Real>(N2x*N2x + N2y*N2y));
140 48347136 : N2x /= Nlength; N2y /= Nlength;
141 :
142 48347136 : Real N3x = dydxi[5];
143 48347136 : Real N3y = - dxdxi[5];
144 48347136 : Nlength = std::sqrt(static_cast<Real>(N3x*N3x + N3y*N3y));
145 48347136 : N3x /= Nlength; N3y /= Nlength;
146 :
147 : // Calculate corner normal vectors (used for reduced element)
148 48347136 : coefs.N01x = dydxi[0];
149 48347136 : coefs.N01y = - dxdxi[0];
150 48347136 : Nlength = std::sqrt(static_cast<Real>(coefs.N01x*coefs.N01x + coefs.N01y*coefs.N01y));
151 48347136 : coefs.N01x /= Nlength; coefs.N01y /= Nlength;
152 :
153 48347136 : coefs.N10x = dydxi[1];
154 48347136 : coefs.N10y = - dxdxi[1];
155 48347136 : Nlength = std::sqrt(static_cast<Real>(coefs.N10x*coefs.N10x + coefs.N10y*coefs.N10y));
156 48347136 : coefs.N10x /= Nlength; coefs.N10y /= Nlength;
157 :
158 48347136 : coefs.N02x = - dydeta[0];
159 48347136 : coefs.N02y = dxdeta[0];
160 48347136 : Nlength = std::sqrt(static_cast<Real>(coefs.N02x*coefs.N02x + coefs.N02y*coefs.N02y));
161 48347136 : coefs.N02x /= Nlength; coefs.N02y /= Nlength;
162 :
163 48347136 : coefs.N20x = - dydeta[2];
164 48347136 : coefs.N20y = dxdeta[2];
165 48347136 : Nlength = std::sqrt(static_cast<Real>(coefs.N20x*coefs.N20x + coefs.N20y*coefs.N20y));
166 48347136 : coefs.N20x /= Nlength; coefs.N20y /= Nlength;
167 :
168 48347136 : coefs.N12x = dydeta[1] - dydxi[1];
169 48347136 : coefs.N12y = dxdxi[1] - dxdeta[1];
170 48347136 : Nlength = std::sqrt(static_cast<Real>(coefs.N12x*coefs.N12x + coefs.N12y*coefs.N12y));
171 48347136 : coefs.N12x /= Nlength; coefs.N12y /= Nlength;
172 :
173 48347136 : coefs.N21x = dydeta[1] - dydxi[1];
174 48347136 : coefs.N21y = dxdxi[1] - dxdeta[1];
175 48347136 : Nlength = std::sqrt(static_cast<Real>(coefs.N21x*coefs.N21x + coefs.N21y*coefs.N21y));
176 48347136 : 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 72664416 : 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 42431664 : N1x = -N1x; N1y = -N1y;
199 42431664 : coefs.N12x = -coefs.N12x; coefs.N12y = -coefs.N12y;
200 42431664 : 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 72664416 : 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 10409448 : N2x = -N2x; N2y = -N2y;
217 10409448 : coefs.N02x = -coefs.N02x; coefs.N02y = -coefs.N02y;
218 10409448 : 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 72664416 : 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 19814748 : N3x = -N3x;
235 19814748 : N3y = -N3y;
236 19814748 : coefs.N01x = -coefs.N01x; coefs.N01y = -coefs.N01y;
237 19814748 : 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 60505776 : Real d1d2ndxi = clough_raw_shape_deriv(0, 0, dofpt[4]);
258 60505776 : Real d1d2ndeta = clough_raw_shape_deriv(0, 1, dofpt[4]);
259 48347136 : Real d1d2ndx = d1d2ndxi * dxidx[4] + d1d2ndeta * detadx[4];
260 48347136 : Real d1d2ndy = d1d2ndxi * dxidy[4] + d1d2ndeta * detady[4];
261 60505776 : Real d1d3ndxi = clough_raw_shape_deriv(0, 0, dofpt[5]);
262 60505776 : Real d1d3ndeta = clough_raw_shape_deriv(0, 1, dofpt[5]);
263 48347136 : Real d1d3ndx = d1d3ndxi * dxidx[5] + d1d3ndeta * detadx[5];
264 48347136 : Real d1d3ndy = d1d3ndxi * dxidy[5] + d1d3ndeta * detady[5];
265 60505776 : Real d2d3ndxi = clough_raw_shape_deriv(1, 0, dofpt[5]);
266 60505776 : Real d2d3ndeta = clough_raw_shape_deriv(1, 1, dofpt[5]);
267 48347136 : Real d2d3ndx = d2d3ndxi * dxidx[5] + d2d3ndeta * detadx[5];
268 48347136 : Real d2d3ndy = d2d3ndxi * dxidy[5] + d2d3ndeta * detady[5];
269 60505776 : Real d2d1ndxi = clough_raw_shape_deriv(1, 0, dofpt[3]);
270 60505776 : Real d2d1ndeta = clough_raw_shape_deriv(1, 1, dofpt[3]);
271 48347136 : Real d2d1ndx = d2d1ndxi * dxidx[3] + d2d1ndeta * detadx[3];
272 48347136 : Real d2d1ndy = d2d1ndxi * dxidy[3] + d2d1ndeta * detady[3];
273 60505776 : Real d3d1ndxi = clough_raw_shape_deriv(2, 0, dofpt[3]);
274 60505776 : Real d3d1ndeta = clough_raw_shape_deriv(2, 1, dofpt[3]);
275 48347136 : Real d3d1ndx = d3d1ndxi * dxidx[3] + d3d1ndeta * detadx[3];
276 48347136 : Real d3d1ndy = d3d1ndxi * dxidy[3] + d3d1ndeta * detady[3];
277 60505776 : Real d3d2ndxi = clough_raw_shape_deriv(2, 0, dofpt[4]);
278 60505776 : Real d3d2ndeta = clough_raw_shape_deriv(2, 1, dofpt[4]);
279 48347136 : Real d3d2ndx = d3d2ndxi * dxidx[4] + d3d2ndeta * detadx[4];
280 48347136 : Real d3d2ndy = d3d2ndxi * dxidy[4] + d3d2ndeta * detady[4];
281 60505776 : Real d1xd2ndxi = clough_raw_shape_deriv(3, 0, dofpt[4]);
282 60505776 : Real d1xd2ndeta = clough_raw_shape_deriv(3, 1, dofpt[4]);
283 48347136 : Real d1xd2ndx = d1xd2ndxi * dxidx[4] + d1xd2ndeta * detadx[4];
284 48347136 : Real d1xd2ndy = d1xd2ndxi * dxidy[4] + d1xd2ndeta * detady[4];
285 60505776 : Real d1xd3ndxi = clough_raw_shape_deriv(3, 0, dofpt[5]);
286 60505776 : Real d1xd3ndeta = clough_raw_shape_deriv(3, 1, dofpt[5]);
287 48347136 : Real d1xd3ndx = d1xd3ndxi * dxidx[5] + d1xd3ndeta * detadx[5];
288 48347136 : Real d1xd3ndy = d1xd3ndxi * dxidy[5] + d1xd3ndeta * detady[5];
289 60505776 : Real d1yd2ndxi = clough_raw_shape_deriv(4, 0, dofpt[4]);
290 60505776 : Real d1yd2ndeta = clough_raw_shape_deriv(4, 1, dofpt[4]);
291 48347136 : Real d1yd2ndx = d1yd2ndxi * dxidx[4] + d1yd2ndeta * detadx[4];
292 48347136 : Real d1yd2ndy = d1yd2ndxi * dxidy[4] + d1yd2ndeta * detady[4];
293 60505776 : Real d1yd3ndxi = clough_raw_shape_deriv(4, 0, dofpt[5]);
294 60505776 : Real d1yd3ndeta = clough_raw_shape_deriv(4, 1, dofpt[5]);
295 48347136 : Real d1yd3ndx = d1yd3ndxi * dxidx[5] + d1yd3ndeta * detadx[5];
296 48347136 : Real d1yd3ndy = d1yd3ndxi * dxidy[5] + d1yd3ndeta * detady[5];
297 60505776 : Real d2xd3ndxi = clough_raw_shape_deriv(5, 0, dofpt[5]);
298 60505776 : Real d2xd3ndeta = clough_raw_shape_deriv(5, 1, dofpt[5]);
299 48347136 : Real d2xd3ndx = d2xd3ndxi * dxidx[5] + d2xd3ndeta * detadx[5];
300 48347136 : Real d2xd3ndy = d2xd3ndxi * dxidy[5] + d2xd3ndeta * detady[5];
301 60505776 : Real d2xd1ndxi = clough_raw_shape_deriv(5, 0, dofpt[3]);
302 60505776 : Real d2xd1ndeta = clough_raw_shape_deriv(5, 1, dofpt[3]);
303 48347136 : Real d2xd1ndx = d2xd1ndxi * dxidx[3] + d2xd1ndeta * detadx[3];
304 48347136 : Real d2xd1ndy = d2xd1ndxi * dxidy[3] + d2xd1ndeta * detady[3];
305 60505776 : Real d2yd3ndxi = clough_raw_shape_deriv(6, 0, dofpt[5]);
306 60505776 : Real d2yd3ndeta = clough_raw_shape_deriv(6, 1, dofpt[5]);
307 48347136 : Real d2yd3ndx = d2yd3ndxi * dxidx[5] + d2yd3ndeta * detadx[5];
308 48347136 : Real d2yd3ndy = d2yd3ndxi * dxidy[5] + d2yd3ndeta * detady[5];
309 60505776 : Real d2yd1ndxi = clough_raw_shape_deriv(6, 0, dofpt[3]);
310 60505776 : Real d2yd1ndeta = clough_raw_shape_deriv(6, 1, dofpt[3]);
311 48347136 : Real d2yd1ndx = d2yd1ndxi * dxidx[3] + d2yd1ndeta * detadx[3];
312 48347136 : Real d2yd1ndy = d2yd1ndxi * dxidy[3] + d2yd1ndeta * detady[3];
313 60505776 : Real d3xd1ndxi = clough_raw_shape_deriv(7, 0, dofpt[3]);
314 60505776 : Real d3xd1ndeta = clough_raw_shape_deriv(7, 1, dofpt[3]);
315 48347136 : Real d3xd1ndx = d3xd1ndxi * dxidx[3] + d3xd1ndeta * detadx[3];
316 48347136 : Real d3xd1ndy = d3xd1ndxi * dxidy[3] + d3xd1ndeta * detady[3];
317 60505776 : Real d3xd2ndxi = clough_raw_shape_deriv(7, 0, dofpt[4]);
318 60505776 : Real d3xd2ndeta = clough_raw_shape_deriv(7, 1, dofpt[4]);
319 48347136 : Real d3xd2ndx = d3xd2ndxi * dxidx[4] + d3xd2ndeta * detadx[4];
320 48347136 : Real d3xd2ndy = d3xd2ndxi * dxidy[4] + d3xd2ndeta * detady[4];
321 60505776 : Real d3yd1ndxi = clough_raw_shape_deriv(8, 0, dofpt[3]);
322 60505776 : Real d3yd1ndeta = clough_raw_shape_deriv(8, 1, dofpt[3]);
323 48347136 : Real d3yd1ndx = d3yd1ndxi * dxidx[3] + d3yd1ndeta * detadx[3];
324 48347136 : Real d3yd1ndy = d3yd1ndxi * dxidy[3] + d3yd1ndeta * detady[3];
325 60505776 : Real d3yd2ndxi = clough_raw_shape_deriv(8, 0, dofpt[4]);
326 60505776 : Real d3yd2ndeta = clough_raw_shape_deriv(8, 1, dofpt[4]);
327 48347136 : Real d3yd2ndx = d3yd2ndxi * dxidx[4] + d3yd2ndeta * detadx[4];
328 48347136 : Real d3yd2ndy = d3yd2ndxi * dxidy[4] + d3yd2ndeta * detady[4];
329 60505776 : Real d1nd1ndxi = clough_raw_shape_deriv(9, 0, dofpt[3]);
330 60505776 : Real d1nd1ndeta = clough_raw_shape_deriv(9, 1, dofpt[3]);
331 48347136 : Real d1nd1ndx = d1nd1ndxi * dxidx[3] + d1nd1ndeta * detadx[3];
332 48347136 : Real d1nd1ndy = d1nd1ndxi * dxidy[3] + d1nd1ndeta * detady[3];
333 60505776 : Real d2nd2ndxi = clough_raw_shape_deriv(10, 0, dofpt[4]);
334 60505776 : Real d2nd2ndeta = clough_raw_shape_deriv(10, 1, dofpt[4]);
335 48347136 : Real d2nd2ndx = d2nd2ndxi * dxidx[4] + d2nd2ndeta * detadx[4];
336 48347136 : Real d2nd2ndy = d2nd2ndxi * dxidy[4] + d2nd2ndeta * detady[4];
337 60505776 : Real d3nd3ndxi = clough_raw_shape_deriv(11, 0, dofpt[5]);
338 60505776 : Real d3nd3ndeta = clough_raw_shape_deriv(11, 1, dofpt[5]);
339 48347136 : Real d3nd3ndx = d3nd3ndxi * dxidx[3] + d3nd3ndeta * detadx[3];
340 48347136 : Real d3nd3ndy = d3nd3ndxi * dxidy[3] + d3nd3ndeta * detady[3];
341 :
342 48347136 : Real d1xd1dxi = clough_raw_shape_deriv(3, 0, dofpt[0]);
343 48347136 : Real d1xd1deta = clough_raw_shape_deriv(3, 1, dofpt[0]);
344 48347136 : Real d1xd1dx = d1xd1dxi * dxidx[0] + d1xd1deta * detadx[0];
345 48347136 : Real d1xd1dy = d1xd1dxi * dxidy[0] + d1xd1deta * detady[0];
346 48347136 : Real d1yd1dxi = clough_raw_shape_deriv(4, 0, dofpt[0]);
347 48347136 : Real d1yd1deta = clough_raw_shape_deriv(4, 1, dofpt[0]);
348 48347136 : Real d1yd1dx = d1yd1dxi * dxidx[0] + d1yd1deta * detadx[0];
349 48347136 : Real d1yd1dy = d1yd1dxi * dxidy[0] + d1yd1deta * detady[0];
350 60505776 : Real d2xd2dxi = clough_raw_shape_deriv(5, 0, dofpt[1]);
351 60505776 : Real d2xd2deta = clough_raw_shape_deriv(5, 1, dofpt[1]);
352 48347136 : Real d2xd2dx = d2xd2dxi * dxidx[1] + d2xd2deta * detadx[1];
353 48347136 : 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 60505776 : Real d2yd2dxi = clough_raw_shape_deriv(6, 0, dofpt[1]);
374 60505776 : Real d2yd2deta = clough_raw_shape_deriv(6, 1, dofpt[1]);
375 48347136 : Real d2yd2dx = d2yd2dxi * dxidx[1] + d2yd2deta * detadx[1];
376 48347136 : Real d2yd2dy = d2yd2dxi * dxidy[1] + d2yd2deta * detady[1];
377 60505776 : Real d3xd3dxi = clough_raw_shape_deriv(7, 0, dofpt[2]);
378 60505776 : Real d3xd3deta = clough_raw_shape_deriv(7, 1, dofpt[2]);
379 48347136 : Real d3xd3dx = d3xd3dxi * dxidx[1] + d3xd3deta * detadx[1];
380 48347136 : Real d3xd3dy = d3xd3dxi * dxidy[1] + d3xd3deta * detady[1];
381 60505776 : Real d3yd3dxi = clough_raw_shape_deriv(8, 0, dofpt[2]);
382 60505776 : Real d3yd3deta = clough_raw_shape_deriv(8, 1, dofpt[2]);
383 48347136 : Real d3yd3dx = d3yd3dxi * dxidx[1] + d3yd3deta * detadx[1];
384 48347136 : Real d3yd3dy = d3yd3dxi * dxidy[1] + d3yd3deta * detady[1];
385 :
386 : // Calculate normal derivatives
387 :
388 48347136 : Real d1nd1ndn = d1nd1ndx * N1x + d1nd1ndy * N1y;
389 48347136 : Real d1xd2ndn = d1xd2ndx * N2x + d1xd2ndy * N2y;
390 48347136 : Real d1xd3ndn = d1xd3ndx * N3x + d1xd3ndy * N3y;
391 48347136 : Real d1yd2ndn = d1yd2ndx * N2x + d1yd2ndy * N2y;
392 48347136 : Real d1yd3ndn = d1yd3ndx * N3x + d1yd3ndy * N3y;
393 :
394 48347136 : Real d2nd2ndn = d2nd2ndx * N2x + d2nd2ndy * N2y;
395 48347136 : Real d2xd3ndn = d2xd3ndx * N3x + d2xd3ndy * N3y;
396 48347136 : Real d2xd1ndn = d2xd1ndx * N1x + d2xd1ndy * N1y;
397 48347136 : Real d2yd3ndn = d2yd3ndx * N3x + d2yd3ndy * N3y;
398 48347136 : Real d2yd1ndn = d2yd1ndx * N1x + d2yd1ndy * N1y;
399 :
400 48347136 : Real d3nd3ndn = d3nd3ndx * N3x + d3nd3ndy * N3y;
401 48347136 : Real d3xd1ndn = d3xd1ndx * N1x + d3xd1ndy * N1y;
402 48347136 : Real d3xd2ndn = d3xd2ndx * N2x + d3xd2ndy * N2y;
403 48347136 : Real d3yd1ndn = d3yd1ndx * N1x + d3yd1ndy * N1y;
404 48347136 : Real d3yd2ndn = d3yd2ndx * N2x + d3yd2ndy * N2y;
405 :
406 : // Calculate midpoint scaling factors
407 :
408 48347136 : coefs.d1nd1n = 1. / d1nd1ndn;
409 48347136 : coefs.d2nd2n = 1. / d2nd2ndn;
410 48347136 : coefs.d3nd3n = 1. / d3nd3ndn;
411 :
412 : // Calculate midpoint derivative adjustments to nodal value
413 : // interpolant functions
414 :
415 48347136 : coefs.d1d2n = -(d1d2ndx * N2x + d1d2ndy * N2y) / d2nd2ndn;
416 48347136 : coefs.d1d3n = -(d1d3ndx * N3x + d1d3ndy * N3y) / d3nd3ndn;
417 48347136 : coefs.d2d3n = -(d2d3ndx * N3x + d2d3ndy * N3y) / d3nd3ndn;
418 48347136 : coefs.d2d1n = -(d2d1ndx * N1x + d2d1ndy * N1y) / d1nd1ndn;
419 48347136 : coefs.d3d1n = -(d3d1ndx * N1x + d3d1ndy * N1y) / d1nd1ndn;
420 48347136 : coefs.d3d2n = -(d3d2ndx * N2x + d3d2ndy * N2y) / d2nd2ndn;
421 :
422 : // Calculate nodal derivative scaling factors
423 :
424 48347136 : coefs.d1xd1x = d1yd1dy / (d1yd1dy * d1xd1dx - d1xd1dy * d1yd1dx);
425 48347136 : coefs.d1xd1y = d1xd1dy / (d1xd1dy * d1yd1dx - d1xd1dx * d1yd1dy);
426 : // coefs.d1xd1y = - coefs.d1xd1x * (d1xd1dy / d1yd1dy);
427 48347136 : coefs.d1yd1y = d1xd1dx / (d1xd1dx * d1yd1dy - d1yd1dx * d1xd1dy);
428 48347136 : coefs.d1yd1x = d1yd1dx / (d1yd1dx * d1xd1dy - d1yd1dy * d1xd1dx);
429 : // coefs.d1yd1x = - coefs.d1yd1y * (d1yd1dx / d1xd1dx);
430 48347136 : coefs.d2xd2x = d2yd2dy / (d2yd2dy * d2xd2dx - d2xd2dy * d2yd2dx);
431 48347136 : coefs.d2xd2y = d2xd2dy / (d2xd2dy * d2yd2dx - d2xd2dx * d2yd2dy);
432 : // coefs.d2xd2y = - coefs.d2xd2x * (d2xd2dy / d2yd2dy);
433 48347136 : coefs.d2yd2y = d2xd2dx / (d2xd2dx * d2yd2dy - d2yd2dx * d2xd2dy);
434 48347136 : coefs.d2yd2x = d2yd2dx / (d2yd2dx * d2xd2dy - d2yd2dy * d2xd2dx);
435 : // coefs.d2yd2x = - coefs.d2yd2y * (d2yd2dx / d2xd2dx);
436 48347136 : coefs.d3xd3x = d3yd3dy / (d3yd3dy * d3xd3dx - d3xd3dy * d3yd3dx);
437 48347136 : coefs.d3xd3y = d3xd3dy / (d3xd3dy * d3yd3dx - d3xd3dx * d3yd3dy);
438 : // coefs.d3xd3y = - coefs.d3xd3x * (d3xd3dy / d3yd3dy);
439 48347136 : coefs.d3yd3y = d3xd3dx / (d3xd3dx * d3yd3dy - d3yd3dx * d3xd3dy);
440 48347136 : 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 48347136 : coefs.d1xd2n = -(coefs.d1xd1x * d1xd2ndn + coefs.d1xd1y * d1yd2ndn) / d2nd2ndn;
460 48347136 : coefs.d1yd2n = -(coefs.d1yd1y * d1yd2ndn + coefs.d1yd1x * d1xd2ndn) / d2nd2ndn;
461 48347136 : coefs.d1xd3n = -(coefs.d1xd1x * d1xd3ndn + coefs.d1xd1y * d1yd3ndn) / d3nd3ndn;
462 48347136 : coefs.d1yd3n = -(coefs.d1yd1y * d1yd3ndn + coefs.d1yd1x * d1xd3ndn) / d3nd3ndn;
463 48347136 : coefs.d2xd3n = -(coefs.d2xd2x * d2xd3ndn + coefs.d2xd2y * d2yd3ndn) / d3nd3ndn;
464 48347136 : coefs.d2yd3n = -(coefs.d2yd2y * d2yd3ndn + coefs.d2yd2x * d2xd3ndn) / d3nd3ndn;
465 48347136 : coefs.d2xd1n = -(coefs.d2xd2x * d2xd1ndn + coefs.d2xd2y * d2yd1ndn) / d1nd1ndn;
466 48347136 : coefs.d2yd1n = -(coefs.d2yd2y * d2yd1ndn + coefs.d2yd2x * d2xd1ndn) / d1nd1ndn;
467 48347136 : coefs.d3xd1n = -(coefs.d3xd3x * d3xd1ndn + coefs.d3xd3y * d3yd1ndn) / d1nd1ndn;
468 48347136 : coefs.d3yd1n = -(coefs.d3yd3y * d3yd1ndn + coefs.d3yd3x * d3xd1ndn) / d1nd1ndn;
469 48347136 : coefs.d3xd2n = -(coefs.d3xd3x * d3xd2ndn + coefs.d3xd3y * d3yd2ndn) / d2nd2ndn;
470 48347136 : 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 48347136 : }
515 :
516 :
517 2062744272 : unsigned char subtriangle_lookup(const Point & p)
518 : {
519 2755786752 : if ((p(0) >= p(1)) && (p(0) + 2 * p(1) <= 1))
520 279727164 : return 0;
521 1643488452 : if ((p(0) <= p(1)) && (p(1) + 2 * p(0) <= 1))
522 687415068 : return 2;
523 182349792 : return 1;
524 : }
525 :
526 :
527 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
528 : // Return shape function second derivatives on the unit right
529 : // triangle
530 65478456 : Real clough_raw_shape_second_deriv(const unsigned int basis_num,
531 : const unsigned int deriv_type,
532 : const Point & p)
533 : {
534 65478456 : Real xi = p(0), eta = p(1);
535 :
536 65478456 : switch (deriv_type)
537 : {
538 :
539 : // second derivative in xi-xi direction
540 21826152 : case 0:
541 21826152 : switch (basis_num)
542 : {
543 453776 : case 0:
544 453776 : switch (subtriangle_lookup(p))
545 : {
546 150810 : case 0:
547 201498 : return -6 + 12*xi;
548 202402 : case 1:
549 202402 : return -30 + 42*xi + 42*eta;
550 151468 : case 2:
551 202382 : return -6 + 18*xi - 6*eta;
552 :
553 0 : default:
554 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
555 : subtriangle_lookup(p));
556 : }
557 453776 : case 1:
558 453776 : switch (subtriangle_lookup(p))
559 : {
560 150810 : case 0:
561 201498 : return 6 - 12*xi;
562 202402 : case 1:
563 202402 : return 18 - 27*xi - 21*eta;
564 151468 : case 2:
565 202382 : return 6 - 15*xi + 3*eta;
566 :
567 0 : default:
568 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
569 : subtriangle_lookup(p));
570 : }
571 453776 : case 2:
572 453776 : switch (subtriangle_lookup(p))
573 : {
574 50688 : case 0:
575 50688 : return 0;
576 202402 : case 1:
577 202402 : return 12 - 15*xi - 21*eta;
578 151468 : case 2:
579 202382 : return -3*xi + 3*eta;
580 :
581 0 : default:
582 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
583 : subtriangle_lookup(p));
584 : }
585 907552 : case 3:
586 907552 : switch (subtriangle_lookup(p))
587 : {
588 301620 : case 0:
589 402996 : return -4 + 6*xi;
590 404804 : case 1:
591 404804 : return -9 + 13*xi + 8*eta;
592 302936 : case 2:
593 404764 : return -1 - 7*xi + 4*eta;
594 :
595 0 : default:
596 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
597 : subtriangle_lookup(p));
598 : }
599 907552 : case 4:
600 907552 : switch (subtriangle_lookup(p))
601 : {
602 301620 : case 0:
603 402996 : return 4*eta;
604 404804 : case 1:
605 404804 : return 1 - 2*xi + 3*eta;
606 302936 : case 2:
607 404764 : return -3 + 14*xi - eta;
608 :
609 0 : default:
610 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
611 : subtriangle_lookup(p));
612 : }
613 907552 : case 5:
614 907552 : switch (subtriangle_lookup(p))
615 : {
616 301620 : case 0:
617 402996 : return -2 + 6*xi;
618 404804 : case 1:
619 404804 : return -4 + 17./2.*xi + 7./2.*eta;
620 302936 : case 2:
621 404764 : return -2 + 13./2.*xi - 1./2.*eta;
622 :
623 0 : default:
624 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
625 : subtriangle_lookup(p));
626 : }
627 907552 : case 6:
628 907552 : switch (subtriangle_lookup(p))
629 : {
630 301620 : case 0:
631 402996 : return 4*eta;
632 404804 : case 1:
633 404804 : return 9 - 23./2.*xi - 23./2.*eta;
634 302936 : case 2:
635 404764 : return -1 + 5./2.*xi + 9./2.*eta;
636 :
637 0 : default:
638 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
639 : subtriangle_lookup(p));
640 : }
641 907552 : case 7:
642 907552 : switch (subtriangle_lookup(p))
643 : {
644 101376 : case 0:
645 101376 : return 0;
646 404804 : case 1:
647 404804 : return 7 - 17./2.*xi - 25./2.*eta;
648 302936 : case 2:
649 404764 : return 1 - 13./2.*xi + 7./2.*eta;
650 :
651 0 : default:
652 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
653 : subtriangle_lookup(p));
654 : }
655 907552 : case 8:
656 907552 : switch (subtriangle_lookup(p))
657 : {
658 101376 : case 0:
659 101376 : return 0;
660 404804 : case 1:
661 404804 : return -2 + 5./2.*xi + 7./2.*eta;
662 302936 : case 2:
663 404764 : return 1./2.*xi - 1./2.*eta;
664 :
665 0 : default:
666 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
667 : subtriangle_lookup(p));
668 : }
669 3176432 : case 9:
670 3176432 : switch (subtriangle_lookup(p))
671 : {
672 354816 : case 0:
673 354816 : return 0;
674 1416814 : case 1:
675 1416814 : return std::sqrt(Real(2)) * (8 - 10*xi - 14*eta);
676 1060276 : case 2:
677 1416674 : return std::sqrt(Real(2)) * (-2*xi + 2*eta);
678 :
679 0 : default:
680 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
681 : subtriangle_lookup(p));
682 : }
683 3176432 : case 10:
684 3176432 : switch (subtriangle_lookup(p))
685 : {
686 354816 : case 0:
687 354816 : return 0;
688 1416814 : case 1:
689 1416814 : return -4 + 4*xi + 8*eta;
690 1060276 : case 2:
691 1416674 : return -4 + 20*xi - 8*eta;
692 :
693 0 : default:
694 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
695 : subtriangle_lookup(p));
696 : }
697 3176432 : case 11:
698 3176432 : switch (subtriangle_lookup(p))
699 : {
700 1055670 : case 0:
701 1410486 : return -8*eta;
702 1416814 : case 1:
703 1416814 : return -12 + 16*xi + 12*eta;
704 1060276 : case 2:
705 1416674 : return 4 - 16*xi - 4*eta;
706 :
707 0 : default:
708 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
709 : subtriangle_lookup(p));
710 : }
711 :
712 0 : default:
713 0 : libmesh_error_msg("Invalid shape function index i = " <<
714 : basis_num);
715 : }
716 :
717 : // second derivative in xi-eta direction
718 21826152 : case 1:
719 21826152 : switch (basis_num)
720 : {
721 453776 : case 0:
722 453776 : switch (subtriangle_lookup(p))
723 : {
724 150810 : case 0:
725 201498 : return -6*eta;
726 202402 : case 1:
727 202402 : return -30 +42*xi
728 202402 : + 42*eta;
729 151468 : case 2:
730 202382 : return -6*xi;
731 :
732 0 : default:
733 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
734 : subtriangle_lookup(p));
735 : }
736 453776 : case 1:
737 453776 : switch (subtriangle_lookup(p))
738 : {
739 150810 : case 0:
740 201498 : return + 3*eta;
741 202402 : case 1:
742 202402 : return 15 - 21*xi - 21*eta;
743 151468 : case 2:
744 202382 : return 3*xi;
745 :
746 0 : default:
747 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
748 : subtriangle_lookup(p));
749 : }
750 453776 : case 2:
751 453776 : switch (subtriangle_lookup(p))
752 : {
753 150810 : case 0:
754 201498 : return 3*eta;
755 202402 : case 1:
756 202402 : return 15 - 21*xi - 21*eta;
757 151468 : case 2:
758 202382 : return 3*xi;
759 :
760 0 : default:
761 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
762 : subtriangle_lookup(p));
763 : }
764 907552 : case 3:
765 907552 : switch (subtriangle_lookup(p))
766 : {
767 301620 : case 0:
768 402996 : return -eta;
769 404804 : case 1:
770 404804 : return -4 + 8*xi + 3*eta;
771 302936 : case 2:
772 404764 : return -3 + 4*xi + 4*eta;
773 :
774 0 : default:
775 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
776 : subtriangle_lookup(p));
777 : }
778 907552 : case 4:
779 907552 : switch (subtriangle_lookup(p))
780 : {
781 301620 : case 0:
782 402996 : return -3 + 4*xi + 4*eta;
783 404804 : case 1:
784 404804 : return - 4 + 3*xi + 8*eta;
785 302936 : case 2:
786 404764 : return -xi;
787 :
788 0 : default:
789 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
790 : subtriangle_lookup(p));
791 : }
792 907552 : case 5:
793 907552 : switch (subtriangle_lookup(p))
794 : {
795 301620 : case 0:
796 402996 : return - 1./2.*eta;
797 404804 : case 1:
798 404804 : return -5./2. + 7./2.*xi + 7./2.*eta;
799 302936 : case 2:
800 404764 : return - 1./2.*xi;
801 :
802 0 : default:
803 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
804 : subtriangle_lookup(p));
805 : }
806 907552 : case 6:
807 907552 : switch (subtriangle_lookup(p))
808 : {
809 301620 : case 0:
810 402996 : return -1 + 4*xi + 7./2.*eta;
811 404804 : case 1:
812 404804 : return 19./2. - 23./2.*xi - 25./2.*eta;
813 302936 : case 2:
814 404764 : return 9./2.*xi;
815 :
816 0 : default:
817 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
818 : subtriangle_lookup(p));
819 : }
820 907552 : case 7:
821 907552 : switch (subtriangle_lookup(p))
822 : {
823 301620 : case 0:
824 402996 : return 9./2.*eta;
825 404804 : case 1:
826 404804 : return 19./2. - 25./2.*xi - 23./2.*eta;
827 302936 : case 2:
828 404764 : return -1 + 7./2.*xi + 4*eta;
829 :
830 0 : default:
831 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
832 : subtriangle_lookup(p));
833 : }
834 907552 : case 8:
835 907552 : switch (subtriangle_lookup(p))
836 : {
837 301620 : case 0:
838 402996 : return -1./2.*eta;
839 404804 : case 1:
840 404804 : return -5./2. + 7./2.*xi + 7./2.*eta;
841 302936 : case 2:
842 404764 : return -1./2.*xi;
843 :
844 0 : default:
845 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
846 : subtriangle_lookup(p));
847 : }
848 3176432 : case 9:
849 3176432 : switch (subtriangle_lookup(p))
850 : {
851 1055670 : case 0:
852 1410486 : return std::sqrt(Real(2)) * (2*eta);
853 1416814 : case 1:
854 1416814 : return std::sqrt(Real(2)) * (10 - 14*xi - 14*eta);
855 1060276 : case 2:
856 1416674 : return std::sqrt(Real(2)) * (2*xi);
857 :
858 0 : default:
859 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
860 : subtriangle_lookup(p));
861 : }
862 3176432 : case 10:
863 3176432 : switch (subtriangle_lookup(p))
864 : {
865 1055670 : case 0:
866 1410486 : return -4*eta;
867 1416814 : case 1:
868 1416814 : return - 8 + 8*xi + 12*eta;
869 1060276 : case 2:
870 1416674 : return 4 - 8*xi - 8*eta;
871 :
872 0 : default:
873 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
874 : subtriangle_lookup(p));
875 : }
876 3176432 : case 11:
877 3176432 : switch (subtriangle_lookup(p))
878 : {
879 1055670 : case 0:
880 1410486 : return 4 - 8*xi - 8*eta;
881 1416814 : case 1:
882 1416814 : return -8 + 12*xi + 8*eta;
883 1060276 : case 2:
884 1416674 : return -4*xi;
885 :
886 0 : default:
887 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
888 : subtriangle_lookup(p));
889 : }
890 :
891 0 : default:
892 0 : libmesh_error_msg("Invalid shape function index i = " <<
893 : basis_num);
894 : }
895 :
896 : // second derivative in eta-eta direction
897 21826152 : case 2:
898 21826152 : switch (basis_num)
899 : {
900 453776 : case 0:
901 453776 : switch (subtriangle_lookup(p))
902 : {
903 150810 : case 0:
904 201498 : return -6 - 6*xi + 18*eta;
905 202402 : case 1:
906 202402 : return -30 + 42*xi + 42*eta;
907 151468 : case 2:
908 202382 : return -6 + 12*eta;
909 :
910 0 : default:
911 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
912 : subtriangle_lookup(p));
913 : }
914 453776 : case 1:
915 453776 : switch (subtriangle_lookup(p))
916 : {
917 150810 : case 0:
918 201498 : return 3*xi - 3*eta;
919 202402 : case 1:
920 202402 : return 12 - 21*xi - 15*eta;
921 50914 : case 2:
922 50914 : return 0;
923 :
924 0 : default:
925 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
926 : subtriangle_lookup(p));
927 : }
928 453776 : case 2:
929 453776 : switch (subtriangle_lookup(p))
930 : {
931 150810 : case 0:
932 201498 : return 6 + 3*xi - 15*eta;
933 202402 : case 1:
934 202402 : return 18 - 21.*xi - 27*eta;
935 151468 : case 2:
936 202382 : return 6 - 12*eta;
937 :
938 0 : default:
939 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
940 : subtriangle_lookup(p));
941 : }
942 907552 : case 3:
943 907552 : switch (subtriangle_lookup(p))
944 : {
945 301620 : case 0:
946 402996 : return -3 - xi + 14*eta;
947 404804 : case 1:
948 404804 : return 1 + 3*xi - 2*eta;
949 302936 : case 2:
950 404764 : return 4*xi;
951 :
952 0 : default:
953 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
954 : subtriangle_lookup(p));
955 : }
956 907552 : case 4:
957 907552 : switch (subtriangle_lookup(p))
958 : {
959 301620 : case 0:
960 402996 : return -1 + 4*xi - 7*eta;
961 404804 : case 1:
962 404804 : return -9 + 8*xi + 13*eta;
963 302936 : case 2:
964 404764 : return -4 + 6*eta;
965 :
966 0 : default:
967 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
968 : subtriangle_lookup(p));
969 : }
970 907552 : case 5:
971 907552 : switch (subtriangle_lookup(p))
972 : {
973 301620 : case 0:
974 402996 : return - 1./2.*xi + 1./2.*eta;
975 404804 : case 1:
976 404804 : return -2 + 7./2.*xi + 5./2.*eta;
977 101828 : case 2:
978 101828 : return 0;
979 :
980 0 : default:
981 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
982 : subtriangle_lookup(p));
983 : }
984 907552 : case 6:
985 907552 : switch (subtriangle_lookup(p))
986 : {
987 301620 : case 0:
988 402996 : return 1 + 7./2.*xi - 13./2.*eta;
989 404804 : case 1:
990 404804 : return 7 - 25./2.*xi - 17./2.*eta;
991 101828 : case 2:
992 101828 : return 0;
993 :
994 0 : default:
995 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
996 : subtriangle_lookup(p));
997 : }
998 907552 : case 7:
999 907552 : switch (subtriangle_lookup(p))
1000 : {
1001 301620 : case 0:
1002 402996 : return -1 + 9./2.*xi + 5./2.*eta;
1003 404804 : case 1:
1004 404804 : return 9 - 23./2.*xi - 23./2.*eta;
1005 302936 : case 2:
1006 404764 : return 4*xi;
1007 :
1008 0 : default:
1009 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1010 : subtriangle_lookup(p));
1011 : }
1012 907552 : case 8:
1013 907552 : switch (subtriangle_lookup(p))
1014 : {
1015 301620 : case 0:
1016 402996 : return -2 - 1./2.*xi + 13./2.*eta;
1017 404804 : case 1:
1018 404804 : return -4 + 7./2.*xi + 17./2.*eta;
1019 302936 : case 2:
1020 404764 : return -2 + 6*eta;
1021 :
1022 0 : default:
1023 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1024 : subtriangle_lookup(p));
1025 : }
1026 3176432 : case 9:
1027 3176432 : switch (subtriangle_lookup(p))
1028 : {
1029 1055670 : case 0:
1030 1410486 : return std::sqrt(Real(2)) * (2*xi - 2*eta);
1031 1416814 : case 1:
1032 1416814 : return std::sqrt(Real(2)) * (8 - 14*xi - 10*eta);
1033 356398 : case 2:
1034 356398 : return 0;
1035 :
1036 0 : default:
1037 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1038 : subtriangle_lookup(p));
1039 : }
1040 3176432 : case 10:
1041 3176432 : switch (subtriangle_lookup(p))
1042 : {
1043 1055670 : case 0:
1044 1410486 : return 4 - 4*xi - 16*eta;
1045 1416814 : case 1:
1046 1416814 : return -12 + 12*xi + 16*eta;
1047 1060276 : case 2:
1048 1416674 : return -8*xi;
1049 :
1050 0 : default:
1051 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1052 : subtriangle_lookup(p));
1053 : }
1054 3176432 : case 11:
1055 3176432 : switch (subtriangle_lookup(p))
1056 : {
1057 1055670 : case 0:
1058 1410486 : return -4 - 8*xi + 20*eta;
1059 1416814 : case 1:
1060 1416814 : return -4 + 8*xi + 4*eta;
1061 356398 : case 2:
1062 356398 : return 0;
1063 :
1064 0 : default:
1065 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1066 : subtriangle_lookup(p));
1067 : }
1068 :
1069 0 : default:
1070 0 : libmesh_error_msg("Invalid shape function index i = " <<
1071 : basis_num);
1072 : }
1073 :
1074 0 : default:
1075 0 : libmesh_error_msg("Invalid shape function derivative j = " <<
1076 : deriv_type);
1077 : }
1078 : }
1079 :
1080 : #endif
1081 :
1082 :
1083 :
1084 2663376624 : Real clough_raw_shape_deriv(const unsigned int basis_num,
1085 : const unsigned int deriv_type,
1086 : const Point & p)
1087 : {
1088 2663376624 : Real xi = p(0), eta = p(1);
1089 :
1090 2663376624 : switch (deriv_type)
1091 : {
1092 : // derivative in xi direction
1093 1331688312 : case 0:
1094 1331688312 : switch (basis_num)
1095 : {
1096 72924176 : case 0:
1097 72924176 : switch (subtriangle_lookup(p))
1098 : {
1099 36372375 : case 0:
1100 48592775 : return -6*xi + 6*xi*xi
1101 48592775 : - 3*eta*eta;
1102 243857 : case 1:
1103 243857 : return 9 - 30*xi + 21*xi*xi
1104 243857 : - 30*eta + 42*xi*eta
1105 243857 : + 21*eta*eta;
1106 36369223 : case 2:
1107 48588630 : return -6*xi + 9*xi*xi
1108 48588630 : - 6*xi*eta;
1109 :
1110 0 : default:
1111 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1112 : subtriangle_lookup(p));
1113 : }
1114 72924176 : case 1:
1115 72924176 : switch (subtriangle_lookup(p))
1116 : {
1117 36372375 : case 0:
1118 48592775 : return 6*xi - 6*xi*xi
1119 48592775 : + 3./2.*eta*eta;
1120 48590993 : case 1:
1121 48590993 : return - 9./2. + 18*xi - 27./2.*xi*xi
1122 48590993 : + 15*eta - 21*xi*eta
1123 48590993 : - 21./2.*eta*eta;
1124 180727 : case 2:
1125 241494 : return 6*xi - 15./2.*xi*xi
1126 241494 : + 3*xi*eta;
1127 :
1128 0 : default:
1129 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1130 : subtriangle_lookup(p));
1131 : }
1132 72924176 : case 2:
1133 72924176 : switch (subtriangle_lookup(p))
1134 : {
1135 183879 : case 0:
1136 245639 : return 3./2.*eta*eta;
1137 48590993 : case 1:
1138 48590993 : return - 9./2. + 12*xi - 15./2.*xi*xi
1139 48590993 : + 15*eta - 21*xi*eta
1140 48590993 : - 21./2.*eta*eta;
1141 36369223 : case 2:
1142 48588630 : return -3./2.*xi*xi
1143 48588630 : + 3*xi*eta;
1144 :
1145 0 : default:
1146 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1147 : subtriangle_lookup(p));
1148 : }
1149 109659856 : case 3:
1150 109659856 : switch (subtriangle_lookup(p))
1151 : {
1152 72744750 : case 0:
1153 97185550 : return 1 - 4*xi + 3*xi*xi
1154 97185550 : - 1./2.*eta*eta;
1155 487714 : case 1:
1156 487714 : return 5./2. - 9*xi + 13./2.*xi*xi
1157 487714 : - 4*eta + 8*xi*eta
1158 487714 : + 3./2.*eta*eta;
1159 36549950 : case 2:
1160 48830124 : return 1 - xi - 7./2.*xi*xi
1161 48830124 : - 3*eta + 4*xi*eta
1162 48830124 : + 2*eta*eta;
1163 :
1164 0 : default:
1165 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1166 : subtriangle_lookup(p));
1167 : }
1168 109659856 : case 4:
1169 109659856 : switch (subtriangle_lookup(p))
1170 : {
1171 72744750 : case 0:
1172 97185550 : return - 3*eta + 4*xi*eta
1173 97185550 : + 2*eta*eta;
1174 487714 : case 1:
1175 487714 : return xi - xi*xi
1176 487714 : - 4*eta + 3*xi*eta
1177 487714 : + 4*eta*eta;
1178 36549950 : case 2:
1179 48830124 : return -3*xi + 7*xi*xi
1180 48830124 : - xi*eta;
1181 :
1182 0 : default:
1183 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1184 : subtriangle_lookup(p));
1185 : }
1186 109659856 : case 5:
1187 109659856 : switch (subtriangle_lookup(p))
1188 : {
1189 72744750 : case 0:
1190 97185550 : return -2*xi + 3*xi*xi
1191 97185550 : - 1./4.*eta*eta;
1192 48834850 : case 1:
1193 48834850 : return 3./4. - 4*xi + 17./4.*xi*xi
1194 48834850 : - 5./2.*eta + 7./2.*xi*eta
1195 48834850 : + 7./4.*eta*eta;
1196 361454 : case 2:
1197 482988 : return -2*xi + 13./4.*xi*xi
1198 482988 : - 1./2.*xi*eta;
1199 :
1200 0 : default:
1201 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1202 : subtriangle_lookup(p));
1203 : }
1204 109659856 : case 6:
1205 109659856 : switch (subtriangle_lookup(p))
1206 : {
1207 72744750 : case 0:
1208 97185550 : return -eta + 4*xi*eta
1209 97185550 : + 7./4.*eta*eta;
1210 48834850 : case 1:
1211 48834850 : return -13./4. + 9*xi - 23./4.*xi*xi
1212 48834850 : + 19./2.*eta - 23./2.*xi*eta
1213 48834850 : - 25./4.*eta*eta;
1214 361454 : case 2:
1215 482988 : return -xi + 5./4.*xi*xi
1216 482988 : + 9./2.*xi*eta;
1217 :
1218 0 : default:
1219 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1220 : subtriangle_lookup(p));
1221 : }
1222 109659856 : case 7:
1223 109659856 : switch (subtriangle_lookup(p))
1224 : {
1225 367758 : case 0:
1226 491278 : return 9./4.*eta*eta;
1227 48834850 : case 1:
1228 48834850 : return - 11./4. + 7*xi - 17./4.*xi*xi
1229 48834850 : + 19./2.*eta - 25./2.*xi*eta
1230 48834850 : - 23./4.*eta*eta;
1231 72738446 : case 2:
1232 97177260 : return xi - 13./4.*xi*xi
1233 97177260 : - eta + 7./2.*xi*eta + 2*eta*eta;
1234 :
1235 0 : default:
1236 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1237 : subtriangle_lookup(p));
1238 : }
1239 109659856 : case 8:
1240 109659856 : switch (subtriangle_lookup(p))
1241 : {
1242 367758 : case 0:
1243 491278 : return - 1./4.*eta*eta;
1244 48834850 : case 1:
1245 48834850 : return 3./4. - 2*xi + 5./4.*xi*xi
1246 48834850 : - 5./2.*eta + 7./2.*xi*eta
1247 48834850 : + 7./4.*eta*eta;
1248 72738446 : case 2:
1249 97177260 : return 1./4.*xi*xi
1250 97177260 : - 1./2.*xi*eta;
1251 :
1252 0 : default:
1253 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1254 : subtriangle_lookup(p));
1255 : }
1256 40018784 : case 9:
1257 40018784 : switch (subtriangle_lookup(p))
1258 : {
1259 1287153 : case 0:
1260 1719473 : return std::sqrt(Real(2)) * eta*eta;
1261 50054135 : case 1:
1262 50054135 : return std::sqrt(Real(2)) * (-3 + 8*xi - 5*xi*xi
1263 50054135 : + 10*eta - 14*xi*eta
1264 50054135 : - 7*eta*eta);
1265 1265089 : case 2:
1266 1690458 : return std::sqrt(Real(2)) * (-xi*xi + 2*xi*eta);
1267 :
1268 0 : default:
1269 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1270 : subtriangle_lookup(p));
1271 : }
1272 40018784 : case 10:
1273 40018784 : switch (subtriangle_lookup(p))
1274 : {
1275 1287153 : case 0:
1276 1719473 : return -2*eta*eta;
1277 1706999 : case 1:
1278 1706999 : return 2 - 4*xi + 2*xi*xi
1279 1706999 : - 8*eta + 8*xi*eta
1280 1706999 : + 6*eta*eta;
1281 37453585 : case 2:
1282 50037594 : return -4*xi + 10*xi*xi
1283 50037594 : + 4*eta - 8*xi*eta
1284 50037594 : - 4*eta*eta;
1285 :
1286 0 : default:
1287 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1288 : subtriangle_lookup(p));
1289 : }
1290 40018784 : case 11:
1291 40018784 : switch (subtriangle_lookup(p))
1292 : {
1293 37475649 : case 0:
1294 50066609 : return 4*eta - 8*xi*eta
1295 50066609 : - 4*eta*eta;
1296 1706999 : case 1:
1297 1706999 : return 4 - 12*xi + 8*xi*xi
1298 1706999 : - 8*eta + 12*xi*eta
1299 1706999 : + 4*eta*eta;
1300 1265089 : case 2:
1301 1690458 : return 4*xi - 8*xi*xi
1302 1690458 : - 4*xi*eta;
1303 :
1304 0 : default:
1305 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1306 : subtriangle_lookup(p));
1307 : }
1308 :
1309 0 : default:
1310 0 : libmesh_error_msg("Invalid shape function index i = " <<
1311 : basis_num);
1312 : }
1313 :
1314 : // derivative in eta direction
1315 1331688312 : case 1:
1316 1331688312 : switch (basis_num)
1317 : {
1318 72924176 : case 0:
1319 72924176 : switch (subtriangle_lookup(p))
1320 : {
1321 36372375 : case 0:
1322 48592775 : return - 6*eta - 6*xi*eta + 9*eta*eta;
1323 243857 : case 1:
1324 243857 : return 9 - 30*xi + 21*xi*xi
1325 243857 : - 30*eta + 42*xi*eta + 21*eta*eta;
1326 36369223 : case 2:
1327 48588630 : return - 3*xi*xi
1328 48588630 : - 6*eta + 6*eta*eta;
1329 :
1330 0 : default:
1331 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1332 : subtriangle_lookup(p));
1333 : }
1334 72924176 : case 1:
1335 72924176 : switch (subtriangle_lookup(p))
1336 : {
1337 36372375 : case 0:
1338 48592775 : return + 3*xi*eta
1339 48592775 : - 3./2.*eta*eta;
1340 48590993 : case 1:
1341 48590993 : return - 9./2. + 15*xi - 21./2.*xi*xi
1342 48590993 : + 12*eta - 21*xi*eta - 15./2.*eta*eta;
1343 180727 : case 2:
1344 241494 : return + 3./2.*xi*xi;
1345 :
1346 0 : default:
1347 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1348 : subtriangle_lookup(p));
1349 : }
1350 72924176 : case 2:
1351 72924176 : switch (subtriangle_lookup(p))
1352 : {
1353 183879 : case 0:
1354 245639 : return 6*eta + 3*xi*eta - 15./2.*eta*eta;
1355 48590993 : case 1:
1356 48590993 : return - 9./2. + 15*xi - 21./2.*xi*xi
1357 48590993 : + 18*eta - 21.*xi*eta - 27./2.*eta*eta;
1358 36369223 : case 2:
1359 48588630 : return 3./2.*xi*xi
1360 48588630 : + 6*eta - 6*eta*eta;
1361 :
1362 0 : default:
1363 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1364 : subtriangle_lookup(p));
1365 : }
1366 109659856 : case 3:
1367 109659856 : switch (subtriangle_lookup(p))
1368 : {
1369 72744750 : case 0:
1370 97185550 : return - 3*eta - xi*eta + 7*eta*eta;
1371 487714 : case 1:
1372 487714 : return - 4*xi + 4*xi*xi
1373 487714 : + eta + 3*xi*eta - eta*eta;
1374 36549950 : case 2:
1375 48830124 : return - 3*xi + 2*xi*xi
1376 48830124 : + 4*xi*eta;
1377 :
1378 0 : default:
1379 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1380 : subtriangle_lookup(p));
1381 : }
1382 109659856 : case 4:
1383 109659856 : switch (subtriangle_lookup(p))
1384 : {
1385 72744750 : case 0:
1386 97185550 : return 1 - 3*xi + 2*xi*xi
1387 97185550 : - eta + 4*xi*eta - 7./2.*eta*eta;
1388 487714 : case 1:
1389 487714 : return 5./2. - 4*xi + 3./2.*xi*xi
1390 487714 : - 9.*eta + 8*xi*eta + 13./2.*eta*eta;
1391 36549950 : case 2:
1392 48830124 : return 1 - 1./2.*xi*xi - 4*eta + 3*eta*eta;
1393 :
1394 0 : default:
1395 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1396 : subtriangle_lookup(p));
1397 : }
1398 109659856 : case 5:
1399 109659856 : switch (subtriangle_lookup(p))
1400 : {
1401 72744750 : case 0:
1402 97185550 : return - 1./2.*xi*eta + 1./4.*eta*eta;
1403 48834850 : case 1:
1404 48834850 : return 3./4. - 5./2.*xi + 7./4.*xi*xi
1405 48834850 : - 2*eta + 7./2.*xi*eta + 5./4.*eta*eta;
1406 361454 : case 2:
1407 482988 : return - 1./4.*xi*xi;
1408 :
1409 0 : default:
1410 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1411 : subtriangle_lookup(p));
1412 : }
1413 109659856 : case 6:
1414 109659856 : switch (subtriangle_lookup(p))
1415 : {
1416 72744750 : case 0:
1417 97185550 : return -xi + 2*xi*xi
1418 97185550 : + eta + 7./2.*xi*eta - 13./4.*eta*eta;
1419 48834850 : case 1:
1420 48834850 : return - 11./4. + 19./2.*xi - 23./4.*xi*xi
1421 48834850 : + 7*eta - 25./2.*xi*eta - 17./4.*eta*eta;
1422 361454 : case 2:
1423 482988 : return 9./4.*xi*xi;
1424 :
1425 0 : default:
1426 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1427 : subtriangle_lookup(p));
1428 : }
1429 109659856 : case 7:
1430 109659856 : switch (subtriangle_lookup(p))
1431 : {
1432 367758 : case 0:
1433 491278 : return -eta + 9./2.*xi*eta + 5./4.*eta*eta;
1434 48834850 : case 1:
1435 48834850 : return - 13./4. + 19./2.*xi - 25./4.*xi*xi
1436 48834850 : + 9*eta - 23./2.*xi*eta - 23./4.*eta*eta;
1437 72738446 : case 2:
1438 97177260 : return - xi + 7./4.*xi*xi + 4*xi*eta;
1439 :
1440 0 : default:
1441 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1442 : subtriangle_lookup(p));
1443 : }
1444 109659856 : case 8:
1445 109659856 : switch (subtriangle_lookup(p))
1446 : {
1447 367758 : case 0:
1448 491278 : return -2*eta - 1./2.*xi*eta + 13./4.*eta*eta;
1449 48834850 : case 1:
1450 48834850 : return 3./4. - 5./2.*xi + 7./4.*xi*xi
1451 48834850 : - 4*eta + 7./2.*xi*eta + 17./4.*eta*eta;
1452 72738446 : case 2:
1453 97177260 : return - 1./4.*xi*xi
1454 97177260 : - 2*eta + 3*eta*eta;
1455 :
1456 0 : default:
1457 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1458 : subtriangle_lookup(p));
1459 : }
1460 40018784 : case 9:
1461 40018784 : switch (subtriangle_lookup(p))
1462 : {
1463 1287153 : case 0:
1464 1719473 : return std::sqrt(Real(2)) * (2*xi*eta - eta*eta);
1465 50054135 : case 1:
1466 50054135 : return std::sqrt(Real(2)) * (- 3 + 10*xi - 7*xi*xi
1467 50054135 : + 8*eta - 14*xi*eta - 5*eta*eta);
1468 1265089 : case 2:
1469 1690458 : return std::sqrt(Real(2)) * (xi*xi);
1470 :
1471 0 : default:
1472 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1473 : subtriangle_lookup(p));
1474 : }
1475 40018784 : case 10:
1476 40018784 : switch (subtriangle_lookup(p))
1477 : {
1478 1287153 : case 0:
1479 1719473 : return 4*eta - 4*xi*eta - 8*eta*eta;
1480 1706999 : case 1:
1481 1706999 : return 4 - 8*xi + 4*xi*xi
1482 1706999 : - 12*eta + 12*xi*eta + 8*eta*eta;
1483 37453585 : case 2:
1484 50037594 : return 4*xi - 4*xi*xi
1485 50037594 : - 8*xi*eta;
1486 :
1487 0 : default:
1488 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1489 : subtriangle_lookup(p));
1490 : }
1491 40018784 : case 11:
1492 40018784 : switch (subtriangle_lookup(p))
1493 : {
1494 37475649 : case 0:
1495 50066609 : return 4*xi - 4*xi*xi
1496 50066609 : - 4*eta - 8*xi*eta + 10.*eta*eta;
1497 1706999 : case 1:
1498 1706999 : return 2 - 8*xi + 6*xi*xi
1499 1706999 : - 4*eta + 8*xi*eta + 2*eta*eta;
1500 1265089 : case 2:
1501 1690458 : return - 2*xi*xi;
1502 :
1503 0 : default:
1504 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1505 : subtriangle_lookup(p));
1506 : }
1507 :
1508 0 : default:
1509 0 : libmesh_error_msg("Invalid shape function index i = " <<
1510 : basis_num);
1511 : }
1512 :
1513 0 : default:
1514 0 : libmesh_error_msg("Invalid shape function derivative j = " <<
1515 : deriv_type);
1516 : }
1517 : }
1518 :
1519 26931672 : Real clough_raw_shape(const unsigned int basis_num,
1520 : const Point & p)
1521 : {
1522 26931672 : Real xi = p(0), eta = p(1);
1523 :
1524 26931672 : switch (basis_num)
1525 : {
1526 560012 : case 0:
1527 560012 : switch (subtriangle_lookup(p))
1528 : {
1529 191596 : case 0:
1530 255931 : return 1 - 3*xi*xi + 2*xi*xi*xi
1531 255931 : - 3*eta*eta - 3*xi*eta*eta + 3*eta*eta*eta;
1532 245301 : case 1:
1533 245301 : return -1 + 9*xi - 15*xi*xi + 7*xi*xi*xi
1534 245301 : + 9*eta - 30*xi*eta +21*xi*xi*eta
1535 245301 : - 15*eta*eta + 21*xi*eta*eta + 7*eta*eta*eta;
1536 184757 : case 2:
1537 246870 : return 1 - 3*xi*xi + 3*xi*xi*xi
1538 246870 : - 3*xi*xi*eta
1539 246870 : - 3*eta*eta + 2*eta*eta*eta;
1540 :
1541 0 : default:
1542 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1543 : subtriangle_lookup(p));
1544 : }
1545 560012 : case 1:
1546 560012 : switch (subtriangle_lookup(p))
1547 : {
1548 191596 : case 0:
1549 255931 : return 3*xi*xi - 2*xi*xi*xi
1550 255931 : + 3./2.*xi*eta*eta
1551 255931 : - 1./2.*eta*eta*eta;
1552 245301 : case 1:
1553 245301 : return 1 - 9./2.*xi + 9*xi*xi - 9./2.*xi*xi*xi
1554 245301 : - 9./2.*eta + 15*xi*eta - 21./2.*xi*xi*eta
1555 245301 : + 6*eta*eta - 21./2.*xi*eta*eta - 5./2.*eta*eta*eta;
1556 184757 : case 2:
1557 246870 : return 3*xi*xi - 5./2.*xi*xi*xi
1558 246870 : + 3./2.*xi*xi*eta;
1559 :
1560 0 : default:
1561 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1562 : subtriangle_lookup(p));
1563 : }
1564 560012 : case 2:
1565 560012 : switch (subtriangle_lookup(p))
1566 : {
1567 191596 : case 0:
1568 255931 : return 3*eta*eta + 3./2.*xi*eta*eta - 5./2.*eta*eta*eta;
1569 245301 : case 1:
1570 245301 : return 1 - 9./2.*xi + 6*xi*xi - 5./2.*xi*xi*xi
1571 245301 : - 9./2.*eta + 15*xi*eta - 21./2.*xi*xi*eta
1572 245301 : + 9*eta*eta - 21./2.*xi*eta*eta - 9./2.*eta*eta*eta;
1573 184757 : case 2:
1574 246870 : return -1./2.*xi*xi*xi
1575 246870 : + 3./2.*xi*xi*eta
1576 246870 : + 3*eta*eta - 2*eta*eta*eta;
1577 :
1578 0 : default:
1579 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1580 : subtriangle_lookup(p));
1581 : }
1582 1120024 : case 3:
1583 1120024 : switch (subtriangle_lookup(p))
1584 : {
1585 383192 : case 0:
1586 511862 : return xi - 2*xi*xi + xi*xi*xi
1587 511862 : - 3./2.*eta*eta - 1./2.*xi*eta*eta + 7./Real(3)*eta*eta*eta;
1588 490602 : case 1:
1589 490602 : return -1./Real(6) + 5./2.*xi - 9./2.*xi*xi + 13./Real(6)*xi*xi*xi
1590 490602 : - 4*xi*eta + 4*xi*xi*eta
1591 490602 : + 1./2.*eta*eta + 3./2.*xi*eta*eta - 1./Real(3)*eta*eta*eta;
1592 369514 : case 2:
1593 493740 : return xi - 1./2.*xi*xi - 7./Real(6)*xi*xi*xi
1594 493740 : - 3*xi*eta + 2*xi*xi*eta
1595 493740 : + 2*xi*eta*eta;
1596 :
1597 0 : default:
1598 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1599 : subtriangle_lookup(p));
1600 : }
1601 1120024 : case 4:
1602 1120024 : switch (subtriangle_lookup(p))
1603 : {
1604 383192 : case 0:
1605 511862 : return eta - 3*xi*eta + 2*xi*xi*eta
1606 511862 : - 1./2.*eta*eta + 2*xi*eta*eta - 7./Real(6)*eta*eta*eta;
1607 490602 : case 1:
1608 490602 : return -1./Real(6) + 1./2.*xi*xi - 1./Real(3)*xi*xi*xi
1609 490602 : + 5./2.*eta - 4*xi*eta + 3./2.*xi*xi*eta
1610 490602 : - 9./2.*eta*eta + 4*xi*eta*eta + 13./Real(6)*eta*eta*eta;
1611 369514 : case 2:
1612 493740 : return -3./2.*xi*xi + 7./Real(3)*xi*xi*xi
1613 493740 : + eta - 1./2.*xi*xi*eta - 2*eta*eta + eta*eta*eta;
1614 :
1615 0 : default:
1616 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1617 : subtriangle_lookup(p));
1618 : }
1619 1120024 : case 5:
1620 1120024 : switch (subtriangle_lookup(p))
1621 : {
1622 383192 : case 0:
1623 511862 : return -xi*xi + xi*xi*xi
1624 511862 : - 1./4.*xi*eta*eta + 1./Real(12)*eta*eta*eta;
1625 490602 : case 1:
1626 490602 : return -1./Real(6) + 3./4.*xi - 2*xi*xi + 17./Real(12)*xi*xi*xi
1627 490602 : + 3./4.*eta - 5./2.*xi*eta + 7./4.*xi*xi*eta
1628 490602 : - eta*eta + 7./4.*xi*eta*eta + 5./Real(12)*eta*eta*eta;
1629 369514 : case 2:
1630 493740 : return -xi*xi + 13./Real(12)*xi*xi*xi
1631 493740 : - 1./4.*xi*xi*eta;
1632 :
1633 0 : default:
1634 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1635 : subtriangle_lookup(p));
1636 : }
1637 1120024 : case 6:
1638 1120024 : switch (subtriangle_lookup(p))
1639 : {
1640 383192 : case 0:
1641 511862 : return -xi*eta + 2*xi*xi*eta
1642 511862 : + 1./2.*eta*eta + 7./4.*xi*eta*eta - 13./Real(12)*eta*eta*eta;
1643 490602 : case 1:
1644 490602 : return 2./Real(3) - 13./4.*xi + 9./2.*xi*xi - 23./Real(12)*xi*xi*xi
1645 490602 : - 11./4.*eta + 19./2.*xi*eta - 23./4.*xi*xi*eta
1646 490602 : + 7./2.*eta*eta - 25./4.*xi*eta*eta - 17./Real(12)*eta*eta*eta;
1647 369514 : case 2:
1648 493740 : return -1./2.*xi*xi + 5./Real(12)*xi*xi*xi
1649 493740 : + 9./4.*xi*xi*eta;
1650 :
1651 0 : default:
1652 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1653 : subtriangle_lookup(p));
1654 : }
1655 1120024 : case 7:
1656 1120024 : switch (subtriangle_lookup(p))
1657 : {
1658 383192 : case 0:
1659 511862 : return -1./2.*eta*eta + 9./4.*xi*eta*eta + 5./Real(12)*eta*eta*eta;
1660 490602 : case 1:
1661 490602 : return 2./Real(3) - 11./4.*xi + 7./2.*xi*xi - 17./Real(12)*xi*xi*xi
1662 490602 : - 13./4.*eta + 19./2.*xi*eta - 25./4.*xi*xi*eta
1663 490602 : + 9./2.*eta*eta - 23./4.*xi*eta*eta - 23./Real(12)*eta*eta*eta;
1664 369514 : case 2:
1665 493740 : return 1./2.*xi*xi - 13./Real(12)*xi*xi*xi
1666 493740 : - xi*eta + 7./4.*xi*xi*eta + 2*xi*eta*eta;
1667 :
1668 0 : default:
1669 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1670 : subtriangle_lookup(p));
1671 : }
1672 1120024 : case 8:
1673 1120024 : switch (subtriangle_lookup(p))
1674 : {
1675 383192 : case 0:
1676 511862 : return -eta*eta - 1./4.*xi*eta*eta + 13./Real(12)*eta*eta*eta;
1677 490602 : case 1:
1678 490602 : return -1./Real(6) + 3./4.*xi - xi*xi + 5./Real(12)*xi*xi*xi
1679 490602 : + 3./4.*eta - 5./2.*xi*eta + 7./4.*xi*xi*eta
1680 490602 : - 2*eta*eta + 7./4.*xi*eta*eta + 17./Real(12)*eta*eta*eta;
1681 369514 : case 2:
1682 493740 : return 1./Real(12)*xi*xi*xi
1683 493740 : - 1./4.*xi*xi*eta
1684 493740 : - eta*eta + eta*eta*eta;
1685 :
1686 0 : default:
1687 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1688 : subtriangle_lookup(p));
1689 : }
1690 3920084 : case 9:
1691 3920084 : switch (subtriangle_lookup(p))
1692 : {
1693 1341172 : case 0:
1694 1791517 : return std::sqrt(Real(2)) * (xi*eta*eta - 1./Real(3)*eta*eta*eta);
1695 1717107 : case 1:
1696 1717107 : return std::sqrt(Real(2)) * (2./Real(3) - 3*xi + 4*xi*xi - 5./Real(3)*xi*xi*xi
1697 1717107 : - 3*eta + 10*xi*eta - 7*xi*xi*eta
1698 1717107 : + 4*eta*eta - 7*xi*eta*eta - 5./Real(3)*eta*eta*eta);
1699 1293299 : case 2:
1700 1728090 : return std::sqrt(Real(2)) * (-1./Real(3)*xi*xi*xi + xi*xi*eta);
1701 :
1702 0 : default:
1703 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1704 : subtriangle_lookup(p));
1705 : }
1706 3920084 : case 10:
1707 3920084 : switch (subtriangle_lookup(p))
1708 : {
1709 1341172 : case 0:
1710 1791517 : return 2*eta*eta - 2*xi*eta*eta - 8./Real(3)*eta*eta*eta;
1711 1717107 : case 1:
1712 1717107 : return -2./Real(3) + 2*xi - 2*xi*xi + 2./Real(3)*xi*xi*xi
1713 1717107 : + 4*eta - 8*xi*eta + 4*xi*xi*eta
1714 1717107 : - 6*eta*eta + 6*xi*eta*eta + 8./Real(3)*eta*eta*eta;
1715 1293299 : case 2:
1716 1728090 : return -2*xi*xi + 10./Real(3)*xi*xi*xi
1717 1728090 : + 4*xi*eta - 4*xi*xi*eta
1718 1728090 : - 4*xi*eta*eta;
1719 :
1720 0 : default:
1721 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1722 : subtriangle_lookup(p));
1723 : }
1724 3920084 : case 11:
1725 3920084 : switch (subtriangle_lookup(p))
1726 : {
1727 1341172 : case 0:
1728 1791517 : return 4*xi*eta - 4*xi*xi*eta
1729 1791517 : - 2*eta*eta - 4*xi*eta*eta + 10./Real(3)*eta*eta*eta;
1730 1717107 : case 1:
1731 1717107 : return -2./Real(3) + 4*xi - 6*xi*xi + 8./Real(3)*xi*xi*xi
1732 1717107 : + 2*eta - 8*xi*eta + 6*xi*xi*eta
1733 1717107 : - 2*eta*eta + 4*xi*eta*eta + 2./Real(3)*eta*eta*eta;
1734 1293299 : case 2:
1735 1728090 : return 2*xi*xi - 8./Real(3)*xi*xi*xi
1736 1728090 : - 2*xi*xi*eta;
1737 :
1738 0 : default:
1739 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1740 : subtriangle_lookup(p));
1741 : }
1742 :
1743 0 : default:
1744 0 : 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 3993164 : LIBMESH_DEFAULT_VECTORIZED_FE(2,CLOUGH)
1758 :
1759 :
1760 : template <>
1761 8977224 : 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 2257080 : libmesh_assert(elem);
1768 :
1769 : CloughCoefs coefs;
1770 8977224 : clough_compute_coefs(elem, coefs);
1771 :
1772 8977224 : const ElemType type = elem->type();
1773 :
1774 : const Order totalorder =
1775 11234304 : order + add_p_level*elem->p_level();
1776 :
1777 8977224 : switch (totalorder)
1778 : {
1779 : // 2nd-order restricted Clough-Tocher element
1780 0 : 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 0 : switch (type)
1787 : {
1788 : // C1 functions on the Clough-Tocher triangle.
1789 0 : case TRI6:
1790 : case TRI7:
1791 : {
1792 0 : 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 0 : switch (i)
1797 : {
1798 : // Note: these DoF numbers are "scrambled" because my
1799 : // initial numbering conventions didn't match libMesh
1800 0 : case 0:
1801 0 : return clough_raw_shape(0, p)
1802 0 : + coefs.d1d2n * clough_raw_shape(10, p)
1803 0 : + coefs.d1d3n * clough_raw_shape(11, p);
1804 0 : case 3:
1805 0 : return clough_raw_shape(1, p)
1806 0 : + coefs.d2d3n * clough_raw_shape(11, p)
1807 0 : + coefs.d2d1n * clough_raw_shape(9, p);
1808 0 : case 6:
1809 0 : return clough_raw_shape(2, p)
1810 0 : + coefs.d3d1n * clough_raw_shape(9, p)
1811 0 : + coefs.d3d2n * clough_raw_shape(10, p);
1812 0 : case 1:
1813 0 : return coefs.d1xd1x * clough_raw_shape(3, p)
1814 0 : + coefs.d1xd1y * clough_raw_shape(4, p)
1815 0 : + coefs.d1xd2n * clough_raw_shape(10, p)
1816 0 : + coefs.d1xd3n * clough_raw_shape(11, p)
1817 0 : + 0.5 * coefs.N01x * coefs.d3nd3n * clough_raw_shape(11, p)
1818 0 : + 0.5 * coefs.N02x * coefs.d2nd2n * clough_raw_shape(10, p);
1819 0 : case 2:
1820 0 : return coefs.d1yd1y * clough_raw_shape(4, p)
1821 0 : + coefs.d1yd1x * clough_raw_shape(3, p)
1822 0 : + coefs.d1yd2n * clough_raw_shape(10, p)
1823 0 : + coefs.d1yd3n * clough_raw_shape(11, p)
1824 0 : + 0.5 * coefs.N01y * coefs.d3nd3n * clough_raw_shape(11, p)
1825 0 : + 0.5 * coefs.N02y * coefs.d2nd2n * clough_raw_shape(10, p);
1826 0 : case 4:
1827 0 : return coefs.d2xd2x * clough_raw_shape(5, p)
1828 0 : + coefs.d2xd2y * clough_raw_shape(6, p)
1829 0 : + coefs.d2xd3n * clough_raw_shape(11, p)
1830 0 : + coefs.d2xd1n * clough_raw_shape(9, p)
1831 0 : + 0.5 * coefs.N10x * coefs.d3nd3n * clough_raw_shape(11, p)
1832 0 : + 0.5 * coefs.N12x * coefs.d1nd1n * clough_raw_shape(9, p);
1833 0 : case 5:
1834 0 : return coefs.d2yd2y * clough_raw_shape(6, p)
1835 0 : + coefs.d2yd2x * clough_raw_shape(5, p)
1836 0 : + coefs.d2yd3n * clough_raw_shape(11, p)
1837 0 : + coefs.d2yd1n * clough_raw_shape(9, p)
1838 0 : + 0.5 * coefs.N10y * coefs.d3nd3n * clough_raw_shape(11, p)
1839 0 : + 0.5 * coefs.N12y * coefs.d1nd1n * clough_raw_shape(9, p);
1840 0 : case 7:
1841 0 : return coefs.d3xd3x * clough_raw_shape(7, p)
1842 0 : + coefs.d3xd3y * clough_raw_shape(8, p)
1843 0 : + coefs.d3xd1n * clough_raw_shape(9, p)
1844 0 : + coefs.d3xd2n * clough_raw_shape(10, p)
1845 0 : + 0.5 * coefs.N20x * coefs.d2nd2n * clough_raw_shape(10, p)
1846 0 : + 0.5 * coefs.N21x * coefs.d1nd1n * clough_raw_shape(9, p);
1847 0 : case 8:
1848 0 : return coefs.d3yd3y * clough_raw_shape(8, p)
1849 0 : + coefs.d3yd3x * clough_raw_shape(7, p)
1850 0 : + coefs.d3yd1n * clough_raw_shape(9, p)
1851 0 : + coefs.d3yd2n * clough_raw_shape(10, p)
1852 0 : + 0.5 * coefs.N20y * coefs.d2nd2n * clough_raw_shape(10, p)
1853 0 : + 0.5 * coefs.N21y * coefs.d1nd1n * clough_raw_shape(9, p);
1854 0 : default:
1855 0 : libmesh_error_msg("Invalid shape function index i = " << i);
1856 : }
1857 : }
1858 0 : default:
1859 0 : libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
1860 : }
1861 : }
1862 : // 3rd-order Clough-Tocher element
1863 8977224 : case THIRD:
1864 : {
1865 8977224 : switch (type)
1866 : {
1867 : // C1 functions on the Clough-Tocher triangle.
1868 8977224 : case TRI6:
1869 : case TRI7:
1870 : {
1871 2257080 : 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 8977224 : switch (i)
1877 : {
1878 : // Note: these DoF numbers are "scrambled" because my
1879 : // initial numbering conventions didn't match libMesh
1880 748102 : case 0:
1881 748102 : return clough_raw_shape(0, p)
1882 748102 : + coefs.d1d2n * clough_raw_shape(10, p)
1883 748102 : + coefs.d1d3n * clough_raw_shape(11, p);
1884 748102 : case 3:
1885 748102 : return clough_raw_shape(1, p)
1886 748102 : + coefs.d2d3n * clough_raw_shape(11, p)
1887 748102 : + coefs.d2d1n * clough_raw_shape(9, p);
1888 748102 : case 6:
1889 748102 : return clough_raw_shape(2, p)
1890 748102 : + coefs.d3d1n * clough_raw_shape(9, p)
1891 748102 : + coefs.d3d2n * clough_raw_shape(10, p);
1892 748102 : case 1:
1893 748102 : return coefs.d1xd1x * clough_raw_shape(3, p)
1894 748102 : + coefs.d1xd1y * clough_raw_shape(4, p)
1895 748102 : + coefs.d1xd2n * clough_raw_shape(10, p)
1896 748102 : + coefs.d1xd3n * clough_raw_shape(11, p);
1897 748102 : case 2:
1898 748102 : return coefs.d1yd1y * clough_raw_shape(4, p)
1899 748102 : + coefs.d1yd1x * clough_raw_shape(3, p)
1900 748102 : + coefs.d1yd2n * clough_raw_shape(10, p)
1901 748102 : + coefs.d1yd3n * clough_raw_shape(11, p);
1902 748102 : case 4:
1903 748102 : return coefs.d2xd2x * clough_raw_shape(5, p)
1904 748102 : + coefs.d2xd2y * clough_raw_shape(6, p)
1905 748102 : + coefs.d2xd3n * clough_raw_shape(11, p)
1906 748102 : + coefs.d2xd1n * clough_raw_shape(9, p);
1907 748102 : case 5:
1908 748102 : return coefs.d2yd2y * clough_raw_shape(6, p)
1909 748102 : + coefs.d2yd2x * clough_raw_shape(5, p)
1910 748102 : + coefs.d2yd3n * clough_raw_shape(11, p)
1911 748102 : + coefs.d2yd1n * clough_raw_shape(9, p);
1912 748102 : case 7:
1913 748102 : return coefs.d3xd3x * clough_raw_shape(7, p)
1914 748102 : + coefs.d3xd3y * clough_raw_shape(8, p)
1915 748102 : + coefs.d3xd1n * clough_raw_shape(9, p)
1916 748102 : + coefs.d3xd2n * clough_raw_shape(10, p);
1917 748102 : case 8:
1918 748102 : return coefs.d3yd3y * clough_raw_shape(8, p)
1919 748102 : + coefs.d3yd3x * clough_raw_shape(7, p)
1920 748102 : + coefs.d3yd1n * clough_raw_shape(9, p)
1921 748102 : + coefs.d3yd2n * clough_raw_shape(10, p);
1922 748102 : case 10:
1923 748102 : return coefs.d1nd1n * clough_raw_shape(9, p);
1924 748102 : case 11:
1925 748102 : return coefs.d2nd2n * clough_raw_shape(10, p);
1926 748102 : case 9:
1927 748102 : return coefs.d3nd3n * clough_raw_shape(11, p);
1928 :
1929 0 : default:
1930 0 : libmesh_error_msg("Invalid shape function index i = " << i);
1931 : }
1932 : }
1933 0 : default:
1934 0 : libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
1935 : }
1936 : }
1937 : // by default throw an error
1938 0 : default:
1939 0 : libmesh_error_msg("ERROR: Unsupported polynomial order = " << order);
1940 : }
1941 : }
1942 :
1943 :
1944 :
1945 : template <>
1946 0 : Real FE<2,CLOUGH>::shape(const ElemType,
1947 : const Order,
1948 : const unsigned int,
1949 : const Point &)
1950 : {
1951 0 : 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 0 : 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 0 : return FE<2,CLOUGH>::shape(elem, fet.order, i, p, add_p_level);
1964 : }
1965 :
1966 :
1967 :
1968 :
1969 : template <>
1970 17543760 : 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 4411344 : libmesh_assert(elem);
1978 :
1979 : CloughCoefs coefs;
1980 17543760 : clough_compute_coefs(elem, coefs);
1981 :
1982 17543760 : const ElemType type = elem->type();
1983 :
1984 : const Order totalorder =
1985 21955104 : order + add_p_level*elem->p_level();
1986 :
1987 17543760 : switch (totalorder)
1988 : {
1989 : // 2nd-order restricted Clough-Tocher element
1990 0 : 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 0 : switch (type)
1997 : {
1998 : // C1 functions on the Clough-Tocher triangle.
1999 0 : case TRI6:
2000 : case TRI7:
2001 : {
2002 0 : 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 0 : switch (i)
2007 : {
2008 : // Note: these DoF numbers are "scrambled" because my
2009 : // initial numbering conventions didn't match libMesh
2010 0 : case 0:
2011 0 : return clough_raw_shape_deriv(0, j, p)
2012 0 : + coefs.d1d2n * clough_raw_shape_deriv(10, j, p)
2013 0 : + coefs.d1d3n * clough_raw_shape_deriv(11, j, p);
2014 0 : case 3:
2015 0 : return clough_raw_shape_deriv(1, j, p)
2016 0 : + coefs.d2d3n * clough_raw_shape_deriv(11, j, p)
2017 0 : + coefs.d2d1n * clough_raw_shape_deriv(9, j, p);
2018 0 : case 6:
2019 0 : return clough_raw_shape_deriv(2, j, p)
2020 0 : + coefs.d3d1n * clough_raw_shape_deriv(9, j, p)
2021 0 : + coefs.d3d2n * clough_raw_shape_deriv(10, j, p);
2022 0 : case 1:
2023 0 : return coefs.d1xd1x * clough_raw_shape_deriv(3, j, p)
2024 0 : + coefs.d1xd1y * clough_raw_shape_deriv(4, j, p)
2025 0 : + coefs.d1xd2n * clough_raw_shape_deriv(10, j, p)
2026 0 : + coefs.d1xd3n * clough_raw_shape_deriv(11, j, p)
2027 0 : + 0.5 * coefs.N01x * coefs.d3nd3n * clough_raw_shape_deriv(11, j, p)
2028 0 : + 0.5 * coefs.N02x * coefs.d2nd2n * clough_raw_shape_deriv(10, j, p);
2029 0 : case 2:
2030 0 : return coefs.d1yd1y * clough_raw_shape_deriv(4, j, p)
2031 0 : + coefs.d1yd1x * clough_raw_shape_deriv(3, j, p)
2032 0 : + coefs.d1yd2n * clough_raw_shape_deriv(10, j, p)
2033 0 : + coefs.d1yd3n * clough_raw_shape_deriv(11, j, p)
2034 0 : + 0.5 * coefs.N01y * coefs.d3nd3n * clough_raw_shape_deriv(11, j, p)
2035 0 : + 0.5 * coefs.N02y * coefs.d2nd2n * clough_raw_shape_deriv(10, j, p);
2036 0 : case 4:
2037 0 : return coefs.d2xd2x * clough_raw_shape_deriv(5, j, p)
2038 0 : + coefs.d2xd2y * clough_raw_shape_deriv(6, j, p)
2039 0 : + coefs.d2xd3n * clough_raw_shape_deriv(11, j, p)
2040 0 : + coefs.d2xd1n * clough_raw_shape_deriv(9, j, p)
2041 0 : + 0.5 * coefs.N10x * coefs.d3nd3n * clough_raw_shape_deriv(11, j, p)
2042 0 : + 0.5 * coefs.N12x * coefs.d1nd1n * clough_raw_shape_deriv(9, j, p);
2043 0 : case 5:
2044 0 : return coefs.d2yd2y * clough_raw_shape_deriv(6, j, p)
2045 0 : + coefs.d2yd2x * clough_raw_shape_deriv(5, j, p)
2046 0 : + coefs.d2yd3n * clough_raw_shape_deriv(11, j, p)
2047 0 : + coefs.d2yd1n * clough_raw_shape_deriv(9, j, p)
2048 0 : + 0.5 * coefs.N10y * coefs.d3nd3n * clough_raw_shape_deriv(11, j, p)
2049 0 : + 0.5 * coefs.N12y * coefs.d1nd1n * clough_raw_shape_deriv(9, j, p);
2050 0 : case 7:
2051 0 : return coefs.d3xd3x * clough_raw_shape_deriv(7, j, p)
2052 0 : + coefs.d3xd3y * clough_raw_shape_deriv(8, j, p)
2053 0 : + coefs.d3xd1n * clough_raw_shape_deriv(9, j, p)
2054 0 : + coefs.d3xd2n * clough_raw_shape_deriv(10, j, p)
2055 0 : + 0.5 * coefs.N20x * coefs.d2nd2n * clough_raw_shape_deriv(10, j, p)
2056 0 : + 0.5 * coefs.N21x * coefs.d1nd1n * clough_raw_shape_deriv(9, j, p);
2057 0 : case 8:
2058 0 : return coefs.d3yd3y * clough_raw_shape_deriv(8, j, p)
2059 0 : + coefs.d3yd3x * clough_raw_shape_deriv(7, j, p)
2060 0 : + coefs.d3yd1n * clough_raw_shape_deriv(9, j, p)
2061 0 : + coefs.d3yd2n * clough_raw_shape_deriv(10, j, p)
2062 0 : + 0.5 * coefs.N20y * coefs.d2nd2n * clough_raw_shape_deriv(10, j, p)
2063 0 : + 0.5 * coefs.N21y * coefs.d1nd1n * clough_raw_shape_deriv(9, j, p);
2064 0 : default:
2065 0 : libmesh_error_msg("Invalid shape function index i = " << i);
2066 : }
2067 : }
2068 0 : default:
2069 0 : libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
2070 : }
2071 : }
2072 : // 3rd-order Clough-Tocher element
2073 17543760 : case THIRD:
2074 : {
2075 17543760 : switch (type)
2076 : {
2077 : // C1 functions on the Clough-Tocher triangle.
2078 17543760 : case TRI6:
2079 : case TRI7:
2080 : {
2081 4411344 : 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 17543760 : switch (i)
2087 : {
2088 : // Note: these DoF numbers are "scrambled" because my
2089 : // initial numbering conventions didn't match libMesh
2090 1461980 : case 0:
2091 1461980 : return clough_raw_shape_deriv(0, j, p)
2092 1461980 : + coefs.d1d2n * clough_raw_shape_deriv(10, j, p)
2093 1461980 : + coefs.d1d3n * clough_raw_shape_deriv(11, j, p);
2094 1461980 : case 3:
2095 1461980 : return clough_raw_shape_deriv(1, j, p)
2096 1461980 : + coefs.d2d3n * clough_raw_shape_deriv(11, j, p)
2097 1461980 : + coefs.d2d1n * clough_raw_shape_deriv(9, j, p);
2098 1461980 : case 6:
2099 1461980 : return clough_raw_shape_deriv(2, j, p)
2100 1461980 : + coefs.d3d1n * clough_raw_shape_deriv(9, j, p)
2101 1461980 : + coefs.d3d2n * clough_raw_shape_deriv(10, j, p);
2102 1461980 : case 1:
2103 1461980 : return coefs.d1xd1x * clough_raw_shape_deriv(3, j, p)
2104 1461980 : + coefs.d1xd1y * clough_raw_shape_deriv(4, j, p)
2105 1461980 : + coefs.d1xd2n * clough_raw_shape_deriv(10, j, p)
2106 1461980 : + coefs.d1xd3n * clough_raw_shape_deriv(11, j, p);
2107 1461980 : case 2:
2108 1461980 : return coefs.d1yd1y * clough_raw_shape_deriv(4, j, p)
2109 1461980 : + coefs.d1yd1x * clough_raw_shape_deriv(3, j, p)
2110 1461980 : + coefs.d1yd2n * clough_raw_shape_deriv(10, j, p)
2111 1461980 : + coefs.d1yd3n * clough_raw_shape_deriv(11, j, p);
2112 1461980 : case 4:
2113 1461980 : return coefs.d2xd2x * clough_raw_shape_deriv(5, j, p)
2114 1461980 : + coefs.d2xd2y * clough_raw_shape_deriv(6, j, p)
2115 1461980 : + coefs.d2xd3n * clough_raw_shape_deriv(11, j, p)
2116 1461980 : + coefs.d2xd1n * clough_raw_shape_deriv(9, j, p);
2117 1461980 : case 5:
2118 1461980 : return coefs.d2yd2y * clough_raw_shape_deriv(6, j, p)
2119 1461980 : + coefs.d2yd2x * clough_raw_shape_deriv(5, j, p)
2120 1461980 : + coefs.d2yd3n * clough_raw_shape_deriv(11, j, p)
2121 1461980 : + coefs.d2yd1n * clough_raw_shape_deriv(9, j, p);
2122 1461980 : case 7:
2123 1461980 : return coefs.d3xd3x * clough_raw_shape_deriv(7, j, p)
2124 1461980 : + coefs.d3xd3y * clough_raw_shape_deriv(8, j, p)
2125 1461980 : + coefs.d3xd1n * clough_raw_shape_deriv(9, j, p)
2126 1461980 : + coefs.d3xd2n * clough_raw_shape_deriv(10, j, p);
2127 1461980 : case 8:
2128 1461980 : return coefs.d3yd3y * clough_raw_shape_deriv(8, j, p)
2129 1461980 : + coefs.d3yd3x * clough_raw_shape_deriv(7, j, p)
2130 1461980 : + coefs.d3yd1n * clough_raw_shape_deriv(9, j, p)
2131 1461980 : + coefs.d3yd2n * clough_raw_shape_deriv(10, j, p);
2132 1461980 : case 10:
2133 1461980 : return coefs.d1nd1n * clough_raw_shape_deriv(9, j, p);
2134 1461980 : case 11:
2135 1461980 : return coefs.d2nd2n * clough_raw_shape_deriv(10, j, p);
2136 1461980 : case 9:
2137 1461980 : return coefs.d3nd3n * clough_raw_shape_deriv(11, j, p);
2138 :
2139 0 : default:
2140 0 : libmesh_error_msg("Invalid shape function index i = " << i);
2141 : }
2142 : }
2143 0 : default:
2144 0 : libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
2145 : }
2146 : }
2147 : // by default throw an error
2148 0 : default:
2149 0 : libmesh_error_msg("ERROR: Unsupported polynomial order = " << order);
2150 : }
2151 : }
2152 :
2153 :
2154 : template <>
2155 0 : Real FE<2,CLOUGH>::shape_deriv(const ElemType,
2156 : const Order,
2157 : const unsigned int,
2158 : const unsigned int,
2159 : const Point &)
2160 : {
2161 0 : 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 0 : 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 0 : 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 21826152 : 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 5490216 : libmesh_assert(elem);
2191 :
2192 : CloughCoefs coefs;
2193 21826152 : clough_compute_coefs(elem, coefs);
2194 :
2195 21826152 : const ElemType type = elem->type();
2196 :
2197 : const Order totalorder =
2198 27316368 : order + add_p_level*elem->p_level();
2199 :
2200 21826152 : switch (totalorder)
2201 : {
2202 : // 2nd-order restricted Clough-Tocher element
2203 0 : case SECOND:
2204 : {
2205 0 : switch (type)
2206 : {
2207 : // C1 functions on the Clough-Tocher triangle.
2208 0 : case TRI6:
2209 : case TRI7:
2210 : {
2211 0 : 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 0 : switch (i)
2216 : {
2217 : // Note: these DoF numbers are "scrambled" because my
2218 : // initial numbering conventions didn't match libMesh
2219 0 : case 0:
2220 0 : return clough_raw_shape_second_deriv(0, j, p)
2221 0 : + coefs.d1d2n * clough_raw_shape_second_deriv(10, j, p)
2222 0 : + coefs.d1d3n * clough_raw_shape_second_deriv(11, j, p);
2223 0 : case 3:
2224 0 : return clough_raw_shape_second_deriv(1, j, p)
2225 0 : + coefs.d2d3n * clough_raw_shape_second_deriv(11, j, p)
2226 0 : + coefs.d2d1n * clough_raw_shape_second_deriv(9, j, p);
2227 0 : case 6:
2228 0 : return clough_raw_shape_second_deriv(2, j, p)
2229 0 : + coefs.d3d1n * clough_raw_shape_second_deriv(9, j, p)
2230 0 : + coefs.d3d2n * clough_raw_shape_second_deriv(10, j, p);
2231 0 : case 1:
2232 0 : return coefs.d1xd1x * clough_raw_shape_second_deriv(3, j, p)
2233 0 : + coefs.d1xd1y * clough_raw_shape_second_deriv(4, j, p)
2234 0 : + coefs.d1xd2n * clough_raw_shape_second_deriv(10, j, p)
2235 0 : + coefs.d1xd3n * clough_raw_shape_second_deriv(11, j, p)
2236 0 : + 0.5 * coefs.N01x * coefs.d3nd3n * clough_raw_shape_second_deriv(11, j, p)
2237 0 : + 0.5 * coefs.N02x * coefs.d2nd2n * clough_raw_shape_second_deriv(10, j, p);
2238 0 : case 2:
2239 0 : return coefs.d1yd1y * clough_raw_shape_second_deriv(4, j, p)
2240 0 : + coefs.d1yd1x * clough_raw_shape_second_deriv(3, j, p)
2241 0 : + coefs.d1yd2n * clough_raw_shape_second_deriv(10, j, p)
2242 0 : + coefs.d1yd3n * clough_raw_shape_second_deriv(11, j, p)
2243 0 : + 0.5 * coefs.N01y * coefs.d3nd3n * clough_raw_shape_second_deriv(11, j, p)
2244 0 : + 0.5 * coefs.N02y * coefs.d2nd2n * clough_raw_shape_second_deriv(10, j, p);
2245 0 : case 4:
2246 0 : return coefs.d2xd2x * clough_raw_shape_second_deriv(5, j, p)
2247 0 : + coefs.d2xd2y * clough_raw_shape_second_deriv(6, j, p)
2248 0 : + coefs.d2xd3n * clough_raw_shape_second_deriv(11, j, p)
2249 0 : + coefs.d2xd1n * clough_raw_shape_second_deriv(9, j, p)
2250 0 : + 0.5 * coefs.N10x * coefs.d3nd3n * clough_raw_shape_second_deriv(11, j, p)
2251 0 : + 0.5 * coefs.N12x * coefs.d1nd1n * clough_raw_shape_second_deriv(9, j, p);
2252 0 : case 5:
2253 0 : return coefs.d2yd2y * clough_raw_shape_second_deriv(6, j, p)
2254 0 : + coefs.d2yd2x * clough_raw_shape_second_deriv(5, j, p)
2255 0 : + coefs.d2yd3n * clough_raw_shape_second_deriv(11, j, p)
2256 0 : + coefs.d2yd1n * clough_raw_shape_second_deriv(9, j, p)
2257 0 : + 0.5 * coefs.N10y * coefs.d3nd3n * clough_raw_shape_second_deriv(11, j, p)
2258 0 : + 0.5 * coefs.N12y * coefs.d1nd1n * clough_raw_shape_second_deriv(9, j, p);
2259 0 : case 7:
2260 0 : return coefs.d3xd3x * clough_raw_shape_second_deriv(7, j, p)
2261 0 : + coefs.d3xd3y * clough_raw_shape_second_deriv(8, j, p)
2262 0 : + coefs.d3xd1n * clough_raw_shape_second_deriv(9, j, p)
2263 0 : + coefs.d3xd2n * clough_raw_shape_second_deriv(10, j, p)
2264 0 : + 0.5 * coefs.N20x * coefs.d2nd2n * clough_raw_shape_second_deriv(10, j, p)
2265 0 : + 0.5 * coefs.N21x * coefs.d1nd1n * clough_raw_shape_second_deriv(9, j, p);
2266 0 : case 8:
2267 0 : return coefs.d3yd3y * clough_raw_shape_second_deriv(8, j, p)
2268 0 : + coefs.d3yd3x * clough_raw_shape_second_deriv(7, j, p)
2269 0 : + coefs.d3yd1n * clough_raw_shape_second_deriv(9, j, p)
2270 0 : + coefs.d3yd2n * clough_raw_shape_second_deriv(10, j, p)
2271 0 : + 0.5 * coefs.N20y * coefs.d2nd2n * clough_raw_shape_second_deriv(10, j, p)
2272 0 : + 0.5 * coefs.N21y * coefs.d1nd1n * clough_raw_shape_second_deriv(9, j, p);
2273 0 : default:
2274 0 : libmesh_error_msg("Invalid shape function index i = " << i);
2275 : }
2276 : }
2277 0 : default:
2278 0 : libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
2279 : }
2280 : }
2281 : // 3rd-order Clough-Tocher element
2282 21826152 : case THIRD:
2283 : {
2284 21826152 : switch (type)
2285 : {
2286 : // C1 functions on the Clough-Tocher triangle.
2287 21826152 : case TRI6:
2288 : case TRI7:
2289 : {
2290 5490216 : 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 21826152 : switch (i)
2296 : {
2297 : // Note: these DoF numbers are "scrambled" because my
2298 : // initial numbering conventions didn't match libMesh
2299 1818846 : case 0:
2300 1818846 : return clough_raw_shape_second_deriv(0, j, p)
2301 1818846 : + coefs.d1d2n * clough_raw_shape_second_deriv(10, j, p)
2302 1818846 : + coefs.d1d3n * clough_raw_shape_second_deriv(11, j, p);
2303 1818846 : case 3:
2304 1818846 : return clough_raw_shape_second_deriv(1, j, p)
2305 1818846 : + coefs.d2d3n * clough_raw_shape_second_deriv(11, j, p)
2306 1818846 : + coefs.d2d1n * clough_raw_shape_second_deriv(9, j, p);
2307 1818846 : case 6:
2308 1818846 : return clough_raw_shape_second_deriv(2, j, p)
2309 1818846 : + coefs.d3d1n * clough_raw_shape_second_deriv(9, j, p)
2310 1818846 : + coefs.d3d2n * clough_raw_shape_second_deriv(10, j, p);
2311 1818846 : case 1:
2312 1818846 : return coefs.d1xd1x * clough_raw_shape_second_deriv(3, j, p)
2313 1818846 : + coefs.d1xd1y * clough_raw_shape_second_deriv(4, j, p)
2314 1818846 : + coefs.d1xd2n * clough_raw_shape_second_deriv(10, j, p)
2315 1818846 : + coefs.d1xd3n * clough_raw_shape_second_deriv(11, j, p);
2316 1818846 : case 2:
2317 1818846 : return coefs.d1yd1y * clough_raw_shape_second_deriv(4, j, p)
2318 1818846 : + coefs.d1yd1x * clough_raw_shape_second_deriv(3, j, p)
2319 1818846 : + coefs.d1yd2n * clough_raw_shape_second_deriv(10, j, p)
2320 1818846 : + coefs.d1yd3n * clough_raw_shape_second_deriv(11, j, p);
2321 1818846 : case 4:
2322 1818846 : return coefs.d2xd2x * clough_raw_shape_second_deriv(5, j, p)
2323 1818846 : + coefs.d2xd2y * clough_raw_shape_second_deriv(6, j, p)
2324 1818846 : + coefs.d2xd3n * clough_raw_shape_second_deriv(11, j, p)
2325 1818846 : + coefs.d2xd1n * clough_raw_shape_second_deriv(9, j, p);
2326 1818846 : case 5:
2327 1818846 : return coefs.d2yd2y * clough_raw_shape_second_deriv(6, j, p)
2328 1818846 : + coefs.d2yd2x * clough_raw_shape_second_deriv(5, j, p)
2329 1818846 : + coefs.d2yd3n * clough_raw_shape_second_deriv(11, j, p)
2330 1818846 : + coefs.d2yd1n * clough_raw_shape_second_deriv(9, j, p);
2331 1818846 : case 7:
2332 1818846 : return coefs.d3xd3x * clough_raw_shape_second_deriv(7, j, p)
2333 1818846 : + coefs.d3xd3y * clough_raw_shape_second_deriv(8, j, p)
2334 1818846 : + coefs.d3xd1n * clough_raw_shape_second_deriv(9, j, p)
2335 1818846 : + coefs.d3xd2n * clough_raw_shape_second_deriv(10, j, p);
2336 1818846 : case 8:
2337 1818846 : return coefs.d3yd3y * clough_raw_shape_second_deriv(8, j, p)
2338 1818846 : + coefs.d3yd3x * clough_raw_shape_second_deriv(7, j, p)
2339 1818846 : + coefs.d3yd1n * clough_raw_shape_second_deriv(9, j, p)
2340 1818846 : + coefs.d3yd2n * clough_raw_shape_second_deriv(10, j, p);
2341 1818846 : case 10:
2342 1818846 : return coefs.d1nd1n * clough_raw_shape_second_deriv(9, j, p);
2343 1818846 : case 11:
2344 1818846 : return coefs.d2nd2n * clough_raw_shape_second_deriv(10, j, p);
2345 1818846 : case 9:
2346 1818846 : return coefs.d3nd3n * clough_raw_shape_second_deriv(11, j, p);
2347 :
2348 0 : default:
2349 0 : libmesh_error_msg("Invalid shape function index i = " << i);
2350 : }
2351 : }
2352 0 : default:
2353 0 : libmesh_error_msg("ERROR: Unsupported element type = " << Utility::enum_to_string(type));
2354 : }
2355 : }
2356 : // by default throw an error
2357 0 : default:
2358 0 : libmesh_error_msg("ERROR: Unsupported polynomial order = " << order);
2359 : }
2360 : }
2361 :
2362 :
2363 : template <>
2364 0 : 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 0 : 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 0 : 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 0 : 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
|
