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 162308676 : void clough_compute_coefs(const Elem * elem, CloughCoefs & coefs)
69 : {
70 162308676 : const FEFamily mapping_family = FEMap::map_fe_type(*elem);
71 162308676 : 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 162308676 : FEInterface::n_shape_functions(map_fe_type, /*extra_order=*/0, elem);
77 :
78 : // Degrees of freedom are at vertices and edge midpoints
79 24340320 : std::vector<Point> dofpt;
80 162308676 : dofpt.push_back(Point(0,0));
81 162308676 : dofpt.push_back(Point(1,0));
82 162308676 : dofpt.push_back(Point(0,1));
83 162308676 : dofpt.push_back(Point(1/2.,1/2.));
84 162308676 : dofpt.push_back(Point(0,1/2.));
85 162308676 : dofpt.push_back(Point(1/2.,0));
86 :
87 : // Mapping functions - first derivatives at each dofpt
88 210989316 : std::vector<Real> dxdxi(6), dxdeta(6), dydxi(6), dydeta(6);
89 210989316 : std::vector<Real> dxidx(6), detadx(6), dxidy(6), detady(6);
90 :
91 : FEInterface::shape_deriv_ptr shape_deriv_ptr =
92 162308676 : FEInterface::shape_deriv_function(map_fe_type, elem);
93 :
94 1136160732 : for (int p = 0; p != 6; ++p)
95 : {
96 : // libMesh::err << p << ' ' << dofpt[p];
97 6820445448 : for (int i = 0; i != n_mapping_shape_functions; ++i)
98 : {
99 : const Real ddxi = shape_deriv_ptr
100 6285029760 : (map_fe_type, elem, i, 0, dofpt[p], /*add_p_level=*/false);
101 : const Real ddeta = shape_deriv_ptr
102 6285029760 : (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 6285029760 : dxdxi[p] += elem->point(i)(0) * ddxi;
108 5846593392 : dydxi[p] += elem->point(i)(1) * ddxi;
109 5846593392 : dxdeta[p] += elem->point(i)(0) * ddeta;
110 6285029760 : 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 1046873016 : const Real inv_jac = 1. / (dxdxi[p]*dydeta[p] -
124 1046873016 : dxdeta[p]*dydxi[p]);
125 973852056 : dxidx[p] = dydeta[p] * inv_jac;
126 973852056 : dxidy[p] = - dxdeta[p] * inv_jac;
127 973852056 : detadx[p] = - dydxi[p] * inv_jac;
128 1046873016 : detady[p] = dxdxi[p] * inv_jac;
129 : }
130 :
131 : // Calculate midpoint normal vectors
132 162308676 : Real N1x = dydeta[3] - dydxi[3];
133 162308676 : Real N1y = dxdxi[3] - dxdeta[3];
134 162308676 : Real Nlength = std::sqrt(static_cast<Real>(N1x*N1x + N1y*N1y));
135 162308676 : N1x /= Nlength; N1y /= Nlength;
136 :
137 162308676 : Real N2x = - dydeta[4];
138 162308676 : Real N2y = dxdeta[4];
139 162308676 : Nlength = std::sqrt(static_cast<Real>(N2x*N2x + N2y*N2y));
140 162308676 : N2x /= Nlength; N2y /= Nlength;
141 :
142 162308676 : Real N3x = dydxi[5];
143 162308676 : Real N3y = - dxdxi[5];
144 162308676 : Nlength = std::sqrt(static_cast<Real>(N3x*N3x + N3y*N3y));
145 162308676 : N3x /= Nlength; N3y /= Nlength;
146 :
147 : // Calculate corner normal vectors (used for reduced element)
148 162308676 : coefs.N01x = dydxi[0];
149 162308676 : coefs.N01y = - dxdxi[0];
150 162308676 : Nlength = std::sqrt(static_cast<Real>(coefs.N01x*coefs.N01x + coefs.N01y*coefs.N01y));
151 162308676 : coefs.N01x /= Nlength; coefs.N01y /= Nlength;
152 :
153 162308676 : coefs.N10x = dydxi[1];
154 162308676 : coefs.N10y = - dxdxi[1];
155 162308676 : Nlength = std::sqrt(static_cast<Real>(coefs.N10x*coefs.N10x + coefs.N10y*coefs.N10y));
156 162308676 : coefs.N10x /= Nlength; coefs.N10y /= Nlength;
157 :
158 162308676 : coefs.N02x = - dydeta[0];
159 162308676 : coefs.N02y = dxdeta[0];
160 162308676 : Nlength = std::sqrt(static_cast<Real>(coefs.N02x*coefs.N02x + coefs.N02y*coefs.N02y));
161 162308676 : coefs.N02x /= Nlength; coefs.N02y /= Nlength;
162 :
163 162308676 : coefs.N20x = - dydeta[2];
164 162308676 : coefs.N20y = dxdeta[2];
165 162308676 : Nlength = std::sqrt(static_cast<Real>(coefs.N20x*coefs.N20x + coefs.N20y*coefs.N20y));
166 162308676 : coefs.N20x /= Nlength; coefs.N20y /= Nlength;
167 :
168 162308676 : coefs.N12x = dydeta[1] - dydxi[1];
169 162308676 : coefs.N12y = dxdxi[1] - dxdeta[1];
170 162308676 : Nlength = std::sqrt(static_cast<Real>(coefs.N12x*coefs.N12x + coefs.N12y*coefs.N12y));
171 162308676 : coefs.N12x /= Nlength; coefs.N12y /= Nlength;
172 :
173 162308676 : coefs.N21x = dydeta[1] - dydxi[1];
174 162308676 : coefs.N21y = dxdxi[1] - dxdeta[1];
175 162308676 : Nlength = std::sqrt(static_cast<Real>(coefs.N21x*coefs.N21x + coefs.N21y*coefs.N21y));
176 162308676 : 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 186648996 : 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 143942748 : N1x = -N1x; N1y = -N1y;
199 143942748 : coefs.N12x = -coefs.N12x; coefs.N12y = -coefs.N12y;
200 143942748 : 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 186648996 : 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 36539808 : N2x = -N2x; N2y = -N2y;
217 36539808 : coefs.N02x = -coefs.N02x; coefs.N02y = -coefs.N02y;
218 36539808 : 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 186648996 : 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 63860412 : N3x = -N3x;
235 63860412 : N3y = -N3y;
236 63860412 : coefs.N01x = -coefs.N01x; coefs.N01y = -coefs.N01y;
237 63860412 : 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 174478836 : Real d1d2ndxi = clough_raw_shape_deriv(0, 0, dofpt[4]);
258 174478836 : Real d1d2ndeta = clough_raw_shape_deriv(0, 1, dofpt[4]);
259 162308676 : Real d1d2ndx = d1d2ndxi * dxidx[4] + d1d2ndeta * detadx[4];
260 162308676 : Real d1d2ndy = d1d2ndxi * dxidy[4] + d1d2ndeta * detady[4];
261 174478836 : Real d1d3ndxi = clough_raw_shape_deriv(0, 0, dofpt[5]);
262 174478836 : Real d1d3ndeta = clough_raw_shape_deriv(0, 1, dofpt[5]);
263 162308676 : Real d1d3ndx = d1d3ndxi * dxidx[5] + d1d3ndeta * detadx[5];
264 162308676 : Real d1d3ndy = d1d3ndxi * dxidy[5] + d1d3ndeta * detady[5];
265 174478836 : Real d2d3ndxi = clough_raw_shape_deriv(1, 0, dofpt[5]);
266 174478836 : Real d2d3ndeta = clough_raw_shape_deriv(1, 1, dofpt[5]);
267 162308676 : Real d2d3ndx = d2d3ndxi * dxidx[5] + d2d3ndeta * detadx[5];
268 162308676 : Real d2d3ndy = d2d3ndxi * dxidy[5] + d2d3ndeta * detady[5];
269 174478836 : Real d2d1ndxi = clough_raw_shape_deriv(1, 0, dofpt[3]);
270 174478836 : Real d2d1ndeta = clough_raw_shape_deriv(1, 1, dofpt[3]);
271 162308676 : Real d2d1ndx = d2d1ndxi * dxidx[3] + d2d1ndeta * detadx[3];
272 162308676 : Real d2d1ndy = d2d1ndxi * dxidy[3] + d2d1ndeta * detady[3];
273 174478836 : Real d3d1ndxi = clough_raw_shape_deriv(2, 0, dofpt[3]);
274 174478836 : Real d3d1ndeta = clough_raw_shape_deriv(2, 1, dofpt[3]);
275 162308676 : Real d3d1ndx = d3d1ndxi * dxidx[3] + d3d1ndeta * detadx[3];
276 162308676 : Real d3d1ndy = d3d1ndxi * dxidy[3] + d3d1ndeta * detady[3];
277 174478836 : Real d3d2ndxi = clough_raw_shape_deriv(2, 0, dofpt[4]);
278 174478836 : Real d3d2ndeta = clough_raw_shape_deriv(2, 1, dofpt[4]);
279 162308676 : Real d3d2ndx = d3d2ndxi * dxidx[4] + d3d2ndeta * detadx[4];
280 162308676 : Real d3d2ndy = d3d2ndxi * dxidy[4] + d3d2ndeta * detady[4];
281 174478836 : Real d1xd2ndxi = clough_raw_shape_deriv(3, 0, dofpt[4]);
282 174478836 : Real d1xd2ndeta = clough_raw_shape_deriv(3, 1, dofpt[4]);
283 162308676 : Real d1xd2ndx = d1xd2ndxi * dxidx[4] + d1xd2ndeta * detadx[4];
284 162308676 : Real d1xd2ndy = d1xd2ndxi * dxidy[4] + d1xd2ndeta * detady[4];
285 174478836 : Real d1xd3ndxi = clough_raw_shape_deriv(3, 0, dofpt[5]);
286 174478836 : Real d1xd3ndeta = clough_raw_shape_deriv(3, 1, dofpt[5]);
287 162308676 : Real d1xd3ndx = d1xd3ndxi * dxidx[5] + d1xd3ndeta * detadx[5];
288 162308676 : Real d1xd3ndy = d1xd3ndxi * dxidy[5] + d1xd3ndeta * detady[5];
289 174478836 : Real d1yd2ndxi = clough_raw_shape_deriv(4, 0, dofpt[4]);
290 174478836 : Real d1yd2ndeta = clough_raw_shape_deriv(4, 1, dofpt[4]);
291 162308676 : Real d1yd2ndx = d1yd2ndxi * dxidx[4] + d1yd2ndeta * detadx[4];
292 162308676 : Real d1yd2ndy = d1yd2ndxi * dxidy[4] + d1yd2ndeta * detady[4];
293 174478836 : Real d1yd3ndxi = clough_raw_shape_deriv(4, 0, dofpt[5]);
294 174478836 : Real d1yd3ndeta = clough_raw_shape_deriv(4, 1, dofpt[5]);
295 162308676 : Real d1yd3ndx = d1yd3ndxi * dxidx[5] + d1yd3ndeta * detadx[5];
296 162308676 : Real d1yd3ndy = d1yd3ndxi * dxidy[5] + d1yd3ndeta * detady[5];
297 174478836 : Real d2xd3ndxi = clough_raw_shape_deriv(5, 0, dofpt[5]);
298 174478836 : Real d2xd3ndeta = clough_raw_shape_deriv(5, 1, dofpt[5]);
299 162308676 : Real d2xd3ndx = d2xd3ndxi * dxidx[5] + d2xd3ndeta * detadx[5];
300 162308676 : Real d2xd3ndy = d2xd3ndxi * dxidy[5] + d2xd3ndeta * detady[5];
301 174478836 : Real d2xd1ndxi = clough_raw_shape_deriv(5, 0, dofpt[3]);
302 174478836 : Real d2xd1ndeta = clough_raw_shape_deriv(5, 1, dofpt[3]);
303 162308676 : Real d2xd1ndx = d2xd1ndxi * dxidx[3] + d2xd1ndeta * detadx[3];
304 162308676 : Real d2xd1ndy = d2xd1ndxi * dxidy[3] + d2xd1ndeta * detady[3];
305 174478836 : Real d2yd3ndxi = clough_raw_shape_deriv(6, 0, dofpt[5]);
306 174478836 : Real d2yd3ndeta = clough_raw_shape_deriv(6, 1, dofpt[5]);
307 162308676 : Real d2yd3ndx = d2yd3ndxi * dxidx[5] + d2yd3ndeta * detadx[5];
308 162308676 : Real d2yd3ndy = d2yd3ndxi * dxidy[5] + d2yd3ndeta * detady[5];
309 174478836 : Real d2yd1ndxi = clough_raw_shape_deriv(6, 0, dofpt[3]);
310 174478836 : Real d2yd1ndeta = clough_raw_shape_deriv(6, 1, dofpt[3]);
311 162308676 : Real d2yd1ndx = d2yd1ndxi * dxidx[3] + d2yd1ndeta * detadx[3];
312 162308676 : Real d2yd1ndy = d2yd1ndxi * dxidy[3] + d2yd1ndeta * detady[3];
313 174478836 : Real d3xd1ndxi = clough_raw_shape_deriv(7, 0, dofpt[3]);
314 174478836 : Real d3xd1ndeta = clough_raw_shape_deriv(7, 1, dofpt[3]);
315 162308676 : Real d3xd1ndx = d3xd1ndxi * dxidx[3] + d3xd1ndeta * detadx[3];
316 162308676 : Real d3xd1ndy = d3xd1ndxi * dxidy[3] + d3xd1ndeta * detady[3];
317 174478836 : Real d3xd2ndxi = clough_raw_shape_deriv(7, 0, dofpt[4]);
318 174478836 : Real d3xd2ndeta = clough_raw_shape_deriv(7, 1, dofpt[4]);
319 162308676 : Real d3xd2ndx = d3xd2ndxi * dxidx[4] + d3xd2ndeta * detadx[4];
320 162308676 : Real d3xd2ndy = d3xd2ndxi * dxidy[4] + d3xd2ndeta * detady[4];
321 174478836 : Real d3yd1ndxi = clough_raw_shape_deriv(8, 0, dofpt[3]);
322 174478836 : Real d3yd1ndeta = clough_raw_shape_deriv(8, 1, dofpt[3]);
323 162308676 : Real d3yd1ndx = d3yd1ndxi * dxidx[3] + d3yd1ndeta * detadx[3];
324 162308676 : Real d3yd1ndy = d3yd1ndxi * dxidy[3] + d3yd1ndeta * detady[3];
325 174478836 : Real d3yd2ndxi = clough_raw_shape_deriv(8, 0, dofpt[4]);
326 174478836 : Real d3yd2ndeta = clough_raw_shape_deriv(8, 1, dofpt[4]);
327 162308676 : Real d3yd2ndx = d3yd2ndxi * dxidx[4] + d3yd2ndeta * detadx[4];
328 162308676 : Real d3yd2ndy = d3yd2ndxi * dxidy[4] + d3yd2ndeta * detady[4];
329 174478836 : Real d1nd1ndxi = clough_raw_shape_deriv(9, 0, dofpt[3]);
330 174478836 : Real d1nd1ndeta = clough_raw_shape_deriv(9, 1, dofpt[3]);
331 162308676 : Real d1nd1ndx = d1nd1ndxi * dxidx[3] + d1nd1ndeta * detadx[3];
332 162308676 : Real d1nd1ndy = d1nd1ndxi * dxidy[3] + d1nd1ndeta * detady[3];
333 174478836 : Real d2nd2ndxi = clough_raw_shape_deriv(10, 0, dofpt[4]);
334 174478836 : Real d2nd2ndeta = clough_raw_shape_deriv(10, 1, dofpt[4]);
335 162308676 : Real d2nd2ndx = d2nd2ndxi * dxidx[4] + d2nd2ndeta * detadx[4];
336 162308676 : Real d2nd2ndy = d2nd2ndxi * dxidy[4] + d2nd2ndeta * detady[4];
337 174478836 : Real d3nd3ndxi = clough_raw_shape_deriv(11, 0, dofpt[5]);
338 174478836 : Real d3nd3ndeta = clough_raw_shape_deriv(11, 1, dofpt[5]);
339 162308676 : Real d3nd3ndx = d3nd3ndxi * dxidx[3] + d3nd3ndeta * detadx[3];
340 162308676 : Real d3nd3ndy = d3nd3ndxi * dxidy[3] + d3nd3ndeta * detady[3];
341 :
342 162308676 : Real d1xd1dxi = clough_raw_shape_deriv(3, 0, dofpt[0]);
343 162308676 : Real d1xd1deta = clough_raw_shape_deriv(3, 1, dofpt[0]);
344 162308676 : Real d1xd1dx = d1xd1dxi * dxidx[0] + d1xd1deta * detadx[0];
345 162308676 : Real d1xd1dy = d1xd1dxi * dxidy[0] + d1xd1deta * detady[0];
346 162308676 : Real d1yd1dxi = clough_raw_shape_deriv(4, 0, dofpt[0]);
347 162308676 : Real d1yd1deta = clough_raw_shape_deriv(4, 1, dofpt[0]);
348 162308676 : Real d1yd1dx = d1yd1dxi * dxidx[0] + d1yd1deta * detadx[0];
349 162308676 : Real d1yd1dy = d1yd1dxi * dxidy[0] + d1yd1deta * detady[0];
350 174478836 : Real d2xd2dxi = clough_raw_shape_deriv(5, 0, dofpt[1]);
351 174478836 : Real d2xd2deta = clough_raw_shape_deriv(5, 1, dofpt[1]);
352 162308676 : Real d2xd2dx = d2xd2dxi * dxidx[1] + d2xd2deta * detadx[1];
353 162308676 : 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 174478836 : Real d2yd2dxi = clough_raw_shape_deriv(6, 0, dofpt[1]);
374 174478836 : Real d2yd2deta = clough_raw_shape_deriv(6, 1, dofpt[1]);
375 162308676 : Real d2yd2dx = d2yd2dxi * dxidx[1] + d2yd2deta * detadx[1];
376 162308676 : Real d2yd2dy = d2yd2dxi * dxidy[1] + d2yd2deta * detady[1];
377 174478836 : Real d3xd3dxi = clough_raw_shape_deriv(7, 0, dofpt[2]);
378 174478836 : Real d3xd3deta = clough_raw_shape_deriv(7, 1, dofpt[2]);
379 162308676 : Real d3xd3dx = d3xd3dxi * dxidx[1] + d3xd3deta * detadx[1];
380 162308676 : Real d3xd3dy = d3xd3dxi * dxidy[1] + d3xd3deta * detady[1];
381 174478836 : Real d3yd3dxi = clough_raw_shape_deriv(8, 0, dofpt[2]);
382 174478836 : Real d3yd3deta = clough_raw_shape_deriv(8, 1, dofpt[2]);
383 162308676 : Real d3yd3dx = d3yd3dxi * dxidx[1] + d3yd3deta * detadx[1];
384 162308676 : Real d3yd3dy = d3yd3dxi * dxidy[1] + d3yd3deta * detady[1];
385 :
386 : // Calculate normal derivatives
387 :
388 162308676 : Real d1nd1ndn = d1nd1ndx * N1x + d1nd1ndy * N1y;
389 162308676 : Real d1xd2ndn = d1xd2ndx * N2x + d1xd2ndy * N2y;
390 162308676 : Real d1xd3ndn = d1xd3ndx * N3x + d1xd3ndy * N3y;
391 162308676 : Real d1yd2ndn = d1yd2ndx * N2x + d1yd2ndy * N2y;
392 162308676 : Real d1yd3ndn = d1yd3ndx * N3x + d1yd3ndy * N3y;
393 :
394 162308676 : Real d2nd2ndn = d2nd2ndx * N2x + d2nd2ndy * N2y;
395 162308676 : Real d2xd3ndn = d2xd3ndx * N3x + d2xd3ndy * N3y;
396 162308676 : Real d2xd1ndn = d2xd1ndx * N1x + d2xd1ndy * N1y;
397 162308676 : Real d2yd3ndn = d2yd3ndx * N3x + d2yd3ndy * N3y;
398 162308676 : Real d2yd1ndn = d2yd1ndx * N1x + d2yd1ndy * N1y;
399 :
400 162308676 : Real d3nd3ndn = d3nd3ndx * N3x + d3nd3ndy * N3y;
401 162308676 : Real d3xd1ndn = d3xd1ndx * N1x + d3xd1ndy * N1y;
402 162308676 : Real d3xd2ndn = d3xd2ndx * N2x + d3xd2ndy * N2y;
403 162308676 : Real d3yd1ndn = d3yd1ndx * N1x + d3yd1ndy * N1y;
404 162308676 : Real d3yd2ndn = d3yd2ndx * N2x + d3yd2ndy * N2y;
405 :
406 : // Calculate midpoint scaling factors
407 :
408 162308676 : coefs.d1nd1n = 1. / d1nd1ndn;
409 162308676 : coefs.d2nd2n = 1. / d2nd2ndn;
410 162308676 : coefs.d3nd3n = 1. / d3nd3ndn;
411 :
412 : // Calculate midpoint derivative adjustments to nodal value
413 : // interpolant functions
414 :
415 162308676 : coefs.d1d2n = -(d1d2ndx * N2x + d1d2ndy * N2y) / d2nd2ndn;
416 162308676 : coefs.d1d3n = -(d1d3ndx * N3x + d1d3ndy * N3y) / d3nd3ndn;
417 162308676 : coefs.d2d3n = -(d2d3ndx * N3x + d2d3ndy * N3y) / d3nd3ndn;
418 162308676 : coefs.d2d1n = -(d2d1ndx * N1x + d2d1ndy * N1y) / d1nd1ndn;
419 162308676 : coefs.d3d1n = -(d3d1ndx * N1x + d3d1ndy * N1y) / d1nd1ndn;
420 162308676 : coefs.d3d2n = -(d3d2ndx * N2x + d3d2ndy * N2y) / d2nd2ndn;
421 :
422 : // Calculate nodal derivative scaling factors
423 :
424 162308676 : coefs.d1xd1x = d1yd1dy / (d1yd1dy * d1xd1dx - d1xd1dy * d1yd1dx);
425 162308676 : coefs.d1xd1y = d1xd1dy / (d1xd1dy * d1yd1dx - d1xd1dx * d1yd1dy);
426 : // coefs.d1xd1y = - coefs.d1xd1x * (d1xd1dy / d1yd1dy);
427 162308676 : coefs.d1yd1y = d1xd1dx / (d1xd1dx * d1yd1dy - d1yd1dx * d1xd1dy);
428 162308676 : coefs.d1yd1x = d1yd1dx / (d1yd1dx * d1xd1dy - d1yd1dy * d1xd1dx);
429 : // coefs.d1yd1x = - coefs.d1yd1y * (d1yd1dx / d1xd1dx);
430 162308676 : coefs.d2xd2x = d2yd2dy / (d2yd2dy * d2xd2dx - d2xd2dy * d2yd2dx);
431 162308676 : coefs.d2xd2y = d2xd2dy / (d2xd2dy * d2yd2dx - d2xd2dx * d2yd2dy);
432 : // coefs.d2xd2y = - coefs.d2xd2x * (d2xd2dy / d2yd2dy);
433 162308676 : coefs.d2yd2y = d2xd2dx / (d2xd2dx * d2yd2dy - d2yd2dx * d2xd2dy);
434 162308676 : coefs.d2yd2x = d2yd2dx / (d2yd2dx * d2xd2dy - d2yd2dy * d2xd2dx);
435 : // coefs.d2yd2x = - coefs.d2yd2y * (d2yd2dx / d2xd2dx);
436 162308676 : coefs.d3xd3x = d3yd3dy / (d3yd3dy * d3xd3dx - d3xd3dy * d3yd3dx);
437 162308676 : coefs.d3xd3y = d3xd3dy / (d3xd3dy * d3yd3dx - d3xd3dx * d3yd3dy);
438 : // coefs.d3xd3y = - coefs.d3xd3x * (d3xd3dy / d3yd3dy);
439 162308676 : coefs.d3yd3y = d3xd3dx / (d3xd3dx * d3yd3dy - d3yd3dx * d3xd3dy);
440 162308676 : 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 162308676 : coefs.d1xd2n = -(coefs.d1xd1x * d1xd2ndn + coefs.d1xd1y * d1yd2ndn) / d2nd2ndn;
460 162308676 : coefs.d1yd2n = -(coefs.d1yd1y * d1yd2ndn + coefs.d1yd1x * d1xd2ndn) / d2nd2ndn;
461 162308676 : coefs.d1xd3n = -(coefs.d1xd1x * d1xd3ndn + coefs.d1xd1y * d1yd3ndn) / d3nd3ndn;
462 162308676 : coefs.d1yd3n = -(coefs.d1yd1y * d1yd3ndn + coefs.d1yd1x * d1xd3ndn) / d3nd3ndn;
463 162308676 : coefs.d2xd3n = -(coefs.d2xd2x * d2xd3ndn + coefs.d2xd2y * d2yd3ndn) / d3nd3ndn;
464 162308676 : coefs.d2yd3n = -(coefs.d2yd2y * d2yd3ndn + coefs.d2yd2x * d2xd3ndn) / d3nd3ndn;
465 162308676 : coefs.d2xd1n = -(coefs.d2xd2x * d2xd1ndn + coefs.d2xd2y * d2yd1ndn) / d1nd1ndn;
466 162308676 : coefs.d2yd1n = -(coefs.d2yd2y * d2yd1ndn + coefs.d2yd2x * d2xd1ndn) / d1nd1ndn;
467 162308676 : coefs.d3xd1n = -(coefs.d3xd3x * d3xd1ndn + coefs.d3xd3y * d3yd1ndn) / d1nd1ndn;
468 162308676 : coefs.d3yd1n = -(coefs.d3yd3y * d3yd1ndn + coefs.d3yd3x * d3xd1ndn) / d1nd1ndn;
469 162308676 : coefs.d3xd2n = -(coefs.d3xd3x * d3xd2ndn + coefs.d3xd3y * d3yd2ndn) / d2nd2ndn;
470 162308676 : 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 162308676 : }
515 :
516 :
517 8557895412 : unsigned char subtriangle_lookup(const Point & p)
518 : {
519 9251594532 : if ((p(0) >= p(1)) && (p(0) + 2 * p(1) <= 1))
520 279991404 : return 0;
521 5517491532 : if ((p(0) <= p(1)) && (p(1) + 2 * p(0) <= 1))
522 2852232480 : return 2;
523 182528244 : return 1;
524 : }
525 :
526 :
527 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
528 : // Return shape function second derivatives on the unit right
529 : // triangle
530 221085072 : Real clough_raw_shape_second_deriv(const unsigned int basis_num,
531 : const unsigned int deriv_type,
532 : const Point & p)
533 : {
534 221085072 : Real xi = p(0), eta = p(1);
535 :
536 221085072 : switch (deriv_type)
537 : {
538 :
539 : // second derivative in xi-xi direction
540 73695024 : case 0:
541 73695024 : switch (basis_num)
542 : {
543 1894506 : case 0:
544 1894506 : switch (subtriangle_lookup(p))
545 : {
546 629868 : case 0:
547 680580 : return -6 + 12*xi;
548 683252 : case 1:
549 683252 : return -30 + 42*xi + 42*eta;
550 632314 : case 2:
551 683252 : return -6 + 18*xi - 6*eta;
552 :
553 0 : default:
554 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
555 : subtriangle_lookup(p));
556 : }
557 1894506 : case 1:
558 1894506 : switch (subtriangle_lookup(p))
559 : {
560 629868 : case 0:
561 680580 : return 6 - 12*xi;
562 683252 : case 1:
563 683252 : return 18 - 27*xi - 21*eta;
564 632314 : case 2:
565 683252 : return 6 - 15*xi + 3*eta;
566 :
567 0 : default:
568 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
569 : subtriangle_lookup(p));
570 : }
571 1894506 : case 2:
572 1894506 : switch (subtriangle_lookup(p))
573 : {
574 50712 : case 0:
575 50712 : return 0;
576 683252 : case 1:
577 683252 : return 12 - 15*xi - 21*eta;
578 632314 : case 2:
579 683252 : return -3*xi + 3*eta;
580 :
581 0 : default:
582 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
583 : subtriangle_lookup(p));
584 : }
585 3789012 : case 3:
586 3789012 : switch (subtriangle_lookup(p))
587 : {
588 1259736 : case 0:
589 1361160 : return -4 + 6*xi;
590 1366504 : case 1:
591 1366504 : return -9 + 13*xi + 8*eta;
592 1264628 : case 2:
593 1366504 : return -1 - 7*xi + 4*eta;
594 :
595 0 : default:
596 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
597 : subtriangle_lookup(p));
598 : }
599 3789012 : case 4:
600 3789012 : switch (subtriangle_lookup(p))
601 : {
602 1259736 : case 0:
603 1361160 : return 4*eta;
604 1366504 : case 1:
605 1366504 : return 1 - 2*xi + 3*eta;
606 1264628 : case 2:
607 1366504 : return -3 + 14*xi - eta;
608 :
609 0 : default:
610 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
611 : subtriangle_lookup(p));
612 : }
613 3789012 : case 5:
614 3789012 : switch (subtriangle_lookup(p))
615 : {
616 1259736 : case 0:
617 1361160 : return -2 + 6*xi;
618 1366504 : case 1:
619 1366504 : return -4 + 17./2.*xi + 7./2.*eta;
620 1264628 : case 2:
621 1366504 : 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 3789012 : case 6:
628 3789012 : switch (subtriangle_lookup(p))
629 : {
630 1259736 : case 0:
631 1361160 : return 4*eta;
632 1366504 : case 1:
633 1366504 : return 9 - 23./2.*xi - 23./2.*eta;
634 1264628 : case 2:
635 1366504 : 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 3789012 : case 7:
642 3789012 : switch (subtriangle_lookup(p))
643 : {
644 101424 : case 0:
645 101424 : return 0;
646 1366504 : case 1:
647 1366504 : return 7 - 17./2.*xi - 25./2.*eta;
648 1264628 : case 2:
649 1366504 : 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 3789012 : case 8:
656 3789012 : switch (subtriangle_lookup(p))
657 : {
658 101424 : case 0:
659 101424 : return 0;
660 1366504 : case 1:
661 1366504 : return -2 + 5./2.*xi + 7./2.*eta;
662 1264628 : case 2:
663 1366504 : 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 13261542 : case 9:
670 13261542 : switch (subtriangle_lookup(p))
671 : {
672 354984 : case 0:
673 354984 : return 0;
674 4782764 : case 1:
675 4782764 : return std::sqrt(Real(2)) * (8 - 10*xi - 14*eta);
676 4426198 : case 2:
677 4782764 : 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 13261542 : case 10:
684 13261542 : switch (subtriangle_lookup(p))
685 : {
686 354984 : case 0:
687 354984 : return 0;
688 4782764 : case 1:
689 4782764 : return -4 + 4*xi + 8*eta;
690 4426198 : case 2:
691 4782764 : return -4 + 20*xi - 8*eta;
692 :
693 0 : default:
694 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
695 : subtriangle_lookup(p));
696 : }
697 13261542 : case 11:
698 13261542 : switch (subtriangle_lookup(p))
699 : {
700 4409076 : case 0:
701 4764060 : return -8*eta;
702 4782764 : case 1:
703 4782764 : return -12 + 16*xi + 12*eta;
704 4426198 : case 2:
705 4782764 : 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 73695024 : case 1:
719 73695024 : switch (basis_num)
720 : {
721 1894506 : case 0:
722 1894506 : switch (subtriangle_lookup(p))
723 : {
724 629868 : case 0:
725 680580 : return -6*eta;
726 683252 : case 1:
727 683252 : return -30 +42*xi
728 683252 : + 42*eta;
729 632314 : case 2:
730 683252 : return -6*xi;
731 :
732 0 : default:
733 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
734 : subtriangle_lookup(p));
735 : }
736 1894506 : case 1:
737 1894506 : switch (subtriangle_lookup(p))
738 : {
739 629868 : case 0:
740 680580 : return + 3*eta;
741 683252 : case 1:
742 683252 : return 15 - 21*xi - 21*eta;
743 632314 : case 2:
744 683252 : return 3*xi;
745 :
746 0 : default:
747 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
748 : subtriangle_lookup(p));
749 : }
750 1894506 : case 2:
751 1894506 : switch (subtriangle_lookup(p))
752 : {
753 629868 : case 0:
754 680580 : return 3*eta;
755 683252 : case 1:
756 683252 : return 15 - 21*xi - 21*eta;
757 632314 : case 2:
758 683252 : return 3*xi;
759 :
760 0 : default:
761 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
762 : subtriangle_lookup(p));
763 : }
764 3789012 : case 3:
765 3789012 : switch (subtriangle_lookup(p))
766 : {
767 1259736 : case 0:
768 1361160 : return -eta;
769 1366504 : case 1:
770 1366504 : return -4 + 8*xi + 3*eta;
771 1264628 : case 2:
772 1366504 : return -3 + 4*xi + 4*eta;
773 :
774 0 : default:
775 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
776 : subtriangle_lookup(p));
777 : }
778 3789012 : case 4:
779 3789012 : switch (subtriangle_lookup(p))
780 : {
781 1259736 : case 0:
782 1361160 : return -3 + 4*xi + 4*eta;
783 1366504 : case 1:
784 1366504 : return - 4 + 3*xi + 8*eta;
785 1264628 : case 2:
786 1366504 : return -xi;
787 :
788 0 : default:
789 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
790 : subtriangle_lookup(p));
791 : }
792 3789012 : case 5:
793 3789012 : switch (subtriangle_lookup(p))
794 : {
795 1259736 : case 0:
796 1361160 : return - 1./2.*eta;
797 1366504 : case 1:
798 1366504 : return -5./2. + 7./2.*xi + 7./2.*eta;
799 1264628 : case 2:
800 1366504 : return - 1./2.*xi;
801 :
802 0 : default:
803 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
804 : subtriangle_lookup(p));
805 : }
806 3789012 : case 6:
807 3789012 : switch (subtriangle_lookup(p))
808 : {
809 1259736 : case 0:
810 1361160 : return -1 + 4*xi + 7./2.*eta;
811 1366504 : case 1:
812 1366504 : return 19./2. - 23./2.*xi - 25./2.*eta;
813 1264628 : case 2:
814 1366504 : return 9./2.*xi;
815 :
816 0 : default:
817 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
818 : subtriangle_lookup(p));
819 : }
820 3789012 : case 7:
821 3789012 : switch (subtriangle_lookup(p))
822 : {
823 1259736 : case 0:
824 1361160 : return 9./2.*eta;
825 1366504 : case 1:
826 1366504 : return 19./2. - 25./2.*xi - 23./2.*eta;
827 1264628 : case 2:
828 1366504 : 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 3789012 : case 8:
835 3789012 : switch (subtriangle_lookup(p))
836 : {
837 1259736 : case 0:
838 1361160 : return -1./2.*eta;
839 1366504 : case 1:
840 1366504 : return -5./2. + 7./2.*xi + 7./2.*eta;
841 1264628 : case 2:
842 1366504 : return -1./2.*xi;
843 :
844 0 : default:
845 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
846 : subtriangle_lookup(p));
847 : }
848 13261542 : case 9:
849 13261542 : switch (subtriangle_lookup(p))
850 : {
851 4409076 : case 0:
852 4764060 : return std::sqrt(Real(2)) * (2*eta);
853 4782764 : case 1:
854 4782764 : return std::sqrt(Real(2)) * (10 - 14*xi - 14*eta);
855 4426198 : case 2:
856 4782764 : 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 13261542 : case 10:
863 13261542 : switch (subtriangle_lookup(p))
864 : {
865 4409076 : case 0:
866 4764060 : return -4*eta;
867 4782764 : case 1:
868 4782764 : return - 8 + 8*xi + 12*eta;
869 4426198 : case 2:
870 4782764 : return 4 - 8*xi - 8*eta;
871 :
872 0 : default:
873 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
874 : subtriangle_lookup(p));
875 : }
876 13261542 : case 11:
877 13261542 : switch (subtriangle_lookup(p))
878 : {
879 4409076 : case 0:
880 4764060 : return 4 - 8*xi - 8*eta;
881 4782764 : case 1:
882 4782764 : return -8 + 12*xi + 8*eta;
883 4426198 : case 2:
884 4782764 : 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 73695024 : case 2:
898 73695024 : switch (basis_num)
899 : {
900 1894506 : case 0:
901 1894506 : switch (subtriangle_lookup(p))
902 : {
903 629868 : case 0:
904 680580 : return -6 - 6*xi + 18*eta;
905 683252 : case 1:
906 683252 : return -30 + 42*xi + 42*eta;
907 632314 : case 2:
908 683252 : return -6 + 12*eta;
909 :
910 0 : default:
911 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
912 : subtriangle_lookup(p));
913 : }
914 1894506 : case 1:
915 1894506 : switch (subtriangle_lookup(p))
916 : {
917 629868 : case 0:
918 680580 : return 3*xi - 3*eta;
919 683252 : case 1:
920 683252 : return 12 - 21*xi - 15*eta;
921 50938 : case 2:
922 50938 : return 0;
923 :
924 0 : default:
925 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
926 : subtriangle_lookup(p));
927 : }
928 1894506 : case 2:
929 1894506 : switch (subtriangle_lookup(p))
930 : {
931 629868 : case 0:
932 680580 : return 6 + 3*xi - 15*eta;
933 683252 : case 1:
934 683252 : return 18 - 21.*xi - 27*eta;
935 632314 : case 2:
936 683252 : return 6 - 12*eta;
937 :
938 0 : default:
939 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
940 : subtriangle_lookup(p));
941 : }
942 3789012 : case 3:
943 3789012 : switch (subtriangle_lookup(p))
944 : {
945 1259736 : case 0:
946 1361160 : return -3 - xi + 14*eta;
947 1366504 : case 1:
948 1366504 : return 1 + 3*xi - 2*eta;
949 1264628 : case 2:
950 1366504 : return 4*xi;
951 :
952 0 : default:
953 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
954 : subtriangle_lookup(p));
955 : }
956 3789012 : case 4:
957 3789012 : switch (subtriangle_lookup(p))
958 : {
959 1259736 : case 0:
960 1361160 : return -1 + 4*xi - 7*eta;
961 1366504 : case 1:
962 1366504 : return -9 + 8*xi + 13*eta;
963 1264628 : case 2:
964 1366504 : return -4 + 6*eta;
965 :
966 0 : default:
967 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
968 : subtriangle_lookup(p));
969 : }
970 3789012 : case 5:
971 3789012 : switch (subtriangle_lookup(p))
972 : {
973 1259736 : case 0:
974 1361160 : return - 1./2.*xi + 1./2.*eta;
975 1366504 : case 1:
976 1366504 : return -2 + 7./2.*xi + 5./2.*eta;
977 101876 : case 2:
978 101876 : return 0;
979 :
980 0 : default:
981 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
982 : subtriangle_lookup(p));
983 : }
984 3789012 : case 6:
985 3789012 : switch (subtriangle_lookup(p))
986 : {
987 1259736 : case 0:
988 1361160 : return 1 + 7./2.*xi - 13./2.*eta;
989 1366504 : case 1:
990 1366504 : return 7 - 25./2.*xi - 17./2.*eta;
991 101876 : case 2:
992 101876 : return 0;
993 :
994 0 : default:
995 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
996 : subtriangle_lookup(p));
997 : }
998 3789012 : case 7:
999 3789012 : switch (subtriangle_lookup(p))
1000 : {
1001 1259736 : case 0:
1002 1361160 : return -1 + 9./2.*xi + 5./2.*eta;
1003 1366504 : case 1:
1004 1366504 : return 9 - 23./2.*xi - 23./2.*eta;
1005 1264628 : case 2:
1006 1366504 : return 4*xi;
1007 :
1008 0 : default:
1009 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1010 : subtriangle_lookup(p));
1011 : }
1012 3789012 : case 8:
1013 3789012 : switch (subtriangle_lookup(p))
1014 : {
1015 1259736 : case 0:
1016 1361160 : return -2 - 1./2.*xi + 13./2.*eta;
1017 1366504 : case 1:
1018 1366504 : return -4 + 7./2.*xi + 17./2.*eta;
1019 1264628 : case 2:
1020 1366504 : return -2 + 6*eta;
1021 :
1022 0 : default:
1023 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1024 : subtriangle_lookup(p));
1025 : }
1026 13261542 : case 9:
1027 13261542 : switch (subtriangle_lookup(p))
1028 : {
1029 4409076 : case 0:
1030 4764060 : return std::sqrt(Real(2)) * (2*xi - 2*eta);
1031 4782764 : case 1:
1032 4782764 : return std::sqrt(Real(2)) * (8 - 14*xi - 10*eta);
1033 356566 : case 2:
1034 356566 : return 0;
1035 :
1036 0 : default:
1037 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1038 : subtriangle_lookup(p));
1039 : }
1040 13261542 : case 10:
1041 13261542 : switch (subtriangle_lookup(p))
1042 : {
1043 4409076 : case 0:
1044 4764060 : return 4 - 4*xi - 16*eta;
1045 4782764 : case 1:
1046 4782764 : return -12 + 12*xi + 16*eta;
1047 4426198 : case 2:
1048 4782764 : return -8*xi;
1049 :
1050 0 : default:
1051 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1052 : subtriangle_lookup(p));
1053 : }
1054 13261542 : case 11:
1055 13261542 : switch (subtriangle_lookup(p))
1056 : {
1057 4409076 : case 0:
1058 4764060 : return -4 - 8*xi + 20*eta;
1059 4782764 : case 1:
1060 4782764 : return -4 + 8*xi + 4*eta;
1061 356566 : case 2:
1062 356566 : 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 8940697344 : Real clough_raw_shape_deriv(const unsigned int basis_num,
1085 : const unsigned int deriv_type,
1086 : const Point & p)
1087 : {
1088 8940697344 : Real xi = p(0), eta = p(1);
1089 :
1090 8940697344 : switch (deriv_type)
1091 : {
1092 : // derivative in xi direction
1093 4470348672 : case 0:
1094 4470348672 : switch (basis_num)
1095 : {
1096 302537823 : case 0:
1097 302537823 : switch (subtriangle_lookup(p))
1098 : {
1099 150898182 : case 0:
1100 163130178 : return -6*xi + 6*xi*xi
1101 163130178 : - 3*eta*eta;
1102 813509 : case 1:
1103 813509 : return 9 - 30*xi + 21*xi*xi
1104 813509 : - 30*eta + 42*xi*eta
1105 813509 : + 21*eta*eta;
1106 150887546 : case 2:
1107 163118510 : return -6*xi + 9*xi*xi
1108 163118510 : - 6*xi*eta;
1109 :
1110 0 : default:
1111 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1112 : subtriangle_lookup(p));
1113 : }
1114 302537823 : case 1:
1115 302537823 : switch (subtriangle_lookup(p))
1116 : {
1117 150898182 : case 0:
1118 163130178 : return 6*xi - 6*xi*xi
1119 163130178 : + 3./2.*eta*eta;
1120 163122185 : case 1:
1121 163122185 : return - 9./2. + 18*xi - 27./2.*xi*xi
1122 163122185 : + 15*eta - 21*xi*eta
1123 163122185 : - 21./2.*eta*eta;
1124 749030 : case 2:
1125 809834 : return 6*xi - 15./2.*xi*xi
1126 809834 : + 3*xi*eta;
1127 :
1128 0 : default:
1129 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1130 : subtriangle_lookup(p));
1131 : }
1132 302537823 : case 2:
1133 302537823 : switch (subtriangle_lookup(p))
1134 : {
1135 759666 : case 0:
1136 821502 : return 3./2.*eta*eta;
1137 163122185 : case 1:
1138 163122185 : return - 9./2. + 12*xi - 15./2.*xi*xi
1139 163122185 : + 15*eta - 21*xi*eta
1140 163122185 : - 21./2.*eta*eta;
1141 150887546 : case 2:
1142 163118510 : return -3./2.*xi*xi
1143 163118510 : + 3*xi*eta;
1144 :
1145 0 : default:
1146 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1147 : subtriangle_lookup(p));
1148 : }
1149 454937130 : case 3:
1150 454937130 : switch (subtriangle_lookup(p))
1151 : {
1152 301796364 : case 0:
1153 326260356 : return 1 - 4*xi + 3*xi*xi
1154 326260356 : - 1./2.*eta*eta;
1155 1627018 : case 1:
1156 1627018 : return 5./2. - 9*xi + 13./2.*xi*xi
1157 1627018 : - 4*eta + 8*xi*eta
1158 1627018 : + 3./2.*eta*eta;
1159 151636576 : case 2:
1160 163928344 : return 1 - xi - 7./2.*xi*xi
1161 163928344 : - 3*eta + 4*xi*eta
1162 163928344 : + 2*eta*eta;
1163 :
1164 0 : default:
1165 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1166 : subtriangle_lookup(p));
1167 : }
1168 454937130 : case 4:
1169 454937130 : switch (subtriangle_lookup(p))
1170 : {
1171 301796364 : case 0:
1172 326260356 : return - 3*eta + 4*xi*eta
1173 326260356 : + 2*eta*eta;
1174 1627018 : case 1:
1175 1627018 : return xi - xi*xi
1176 1627018 : - 4*eta + 3*xi*eta
1177 1627018 : + 4*eta*eta;
1178 151636576 : case 2:
1179 163928344 : return -3*xi + 7*xi*xi
1180 163928344 : - xi*eta;
1181 :
1182 0 : default:
1183 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1184 : subtriangle_lookup(p));
1185 : }
1186 454937130 : case 5:
1187 454937130 : switch (subtriangle_lookup(p))
1188 : {
1189 301796364 : case 0:
1190 326260356 : return -2*xi + 3*xi*xi
1191 326260356 : - 1./4.*eta*eta;
1192 163935694 : case 1:
1193 163935694 : return 3./4. - 4*xi + 17./4.*xi*xi
1194 163935694 : - 5./2.*eta + 7./2.*xi*eta
1195 163935694 : + 7./4.*eta*eta;
1196 1498060 : case 2:
1197 1619668 : return -2*xi + 13./4.*xi*xi
1198 1619668 : - 1./2.*xi*eta;
1199 :
1200 0 : default:
1201 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1202 : subtriangle_lookup(p));
1203 : }
1204 454937130 : case 6:
1205 454937130 : switch (subtriangle_lookup(p))
1206 : {
1207 301796364 : case 0:
1208 326260356 : return -eta + 4*xi*eta
1209 326260356 : + 7./4.*eta*eta;
1210 163935694 : case 1:
1211 163935694 : return -13./4. + 9*xi - 23./4.*xi*xi
1212 163935694 : + 19./2.*eta - 23./2.*xi*eta
1213 163935694 : - 25./4.*eta*eta;
1214 1498060 : case 2:
1215 1619668 : return -xi + 5./4.*xi*xi
1216 1619668 : + 9./2.*xi*eta;
1217 :
1218 0 : default:
1219 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1220 : subtriangle_lookup(p));
1221 : }
1222 454937130 : case 7:
1223 454937130 : switch (subtriangle_lookup(p))
1224 : {
1225 1519332 : case 0:
1226 1643004 : return 9./4.*eta*eta;
1227 163935694 : case 1:
1228 163935694 : return - 11./4. + 7*xi - 17./4.*xi*xi
1229 163935694 : + 19./2.*eta - 25./2.*xi*eta
1230 163935694 : - 23./4.*eta*eta;
1231 301775092 : case 2:
1232 326237020 : return xi - 13./4.*xi*xi
1233 326237020 : - 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 454937130 : case 8:
1240 454937130 : switch (subtriangle_lookup(p))
1241 : {
1242 1519332 : case 0:
1243 1643004 : return - 1./4.*eta*eta;
1244 163935694 : case 1:
1245 163935694 : return 3./4. - 2*xi + 5./4.*xi*xi
1246 163935694 : - 5./2.*eta + 7./2.*xi*eta
1247 163935694 : + 7./4.*eta*eta;
1248 301775092 : case 2:
1249 326237020 : return 1./4.*xi*xi
1250 326237020 : - 1./2.*xi*eta;
1251 :
1252 0 : default:
1253 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1254 : subtriangle_lookup(p));
1255 : }
1256 165964053 : case 9:
1257 165964053 : switch (subtriangle_lookup(p))
1258 : {
1259 5317662 : case 0:
1260 5750514 : return std::sqrt(Real(2)) * eta*eta;
1261 168003239 : case 1:
1262 168003239 : return std::sqrt(Real(2)) * (-3 + 8*xi - 5*xi*xi
1263 168003239 : + 10*eta - 14*xi*eta
1264 168003239 : - 7*eta*eta);
1265 5243210 : case 2:
1266 5668838 : 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 165964053 : case 10:
1273 165964053 : switch (subtriangle_lookup(p))
1274 : {
1275 5317662 : case 0:
1276 5750514 : return -2*eta*eta;
1277 5694563 : case 1:
1278 5694563 : return 2 - 4*xi + 2*xi*xi
1279 5694563 : - 8*eta + 8*xi*eta
1280 5694563 : + 6*eta*eta;
1281 155381726 : case 2:
1282 167977514 : return -4*xi + 10*xi*xi
1283 167977514 : + 4*eta - 8*xi*eta
1284 167977514 : - 4*eta*eta;
1285 :
1286 0 : default:
1287 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1288 : subtriangle_lookup(p));
1289 : }
1290 165964053 : case 11:
1291 165964053 : switch (subtriangle_lookup(p))
1292 : {
1293 155456178 : case 0:
1294 168059190 : return 4*eta - 8*xi*eta
1295 168059190 : - 4*eta*eta;
1296 5694563 : case 1:
1297 5694563 : return 4 - 12*xi + 8*xi*xi
1298 5694563 : - 8*eta + 12*xi*eta
1299 5694563 : + 4*eta*eta;
1300 5243210 : case 2:
1301 5668838 : return 4*xi - 8*xi*xi
1302 5668838 : - 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 4470348672 : case 1:
1316 4470348672 : switch (basis_num)
1317 : {
1318 302537823 : case 0:
1319 302537823 : switch (subtriangle_lookup(p))
1320 : {
1321 150898182 : case 0:
1322 163130178 : return - 6*eta - 6*xi*eta + 9*eta*eta;
1323 813509 : case 1:
1324 813509 : return 9 - 30*xi + 21*xi*xi
1325 813509 : - 30*eta + 42*xi*eta + 21*eta*eta;
1326 150887546 : case 2:
1327 163118510 : return - 3*xi*xi
1328 163118510 : - 6*eta + 6*eta*eta;
1329 :
1330 0 : default:
1331 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1332 : subtriangle_lookup(p));
1333 : }
1334 302537823 : case 1:
1335 302537823 : switch (subtriangle_lookup(p))
1336 : {
1337 150898182 : case 0:
1338 163130178 : return + 3*xi*eta
1339 163130178 : - 3./2.*eta*eta;
1340 163122185 : case 1:
1341 163122185 : return - 9./2. + 15*xi - 21./2.*xi*xi
1342 163122185 : + 12*eta - 21*xi*eta - 15./2.*eta*eta;
1343 749030 : case 2:
1344 809834 : return + 3./2.*xi*xi;
1345 :
1346 0 : default:
1347 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1348 : subtriangle_lookup(p));
1349 : }
1350 302537823 : case 2:
1351 302537823 : switch (subtriangle_lookup(p))
1352 : {
1353 759666 : case 0:
1354 821502 : return 6*eta + 3*xi*eta - 15./2.*eta*eta;
1355 163122185 : case 1:
1356 163122185 : return - 9./2. + 15*xi - 21./2.*xi*xi
1357 163122185 : + 18*eta - 21.*xi*eta - 27./2.*eta*eta;
1358 150887546 : case 2:
1359 163118510 : return 3./2.*xi*xi
1360 163118510 : + 6*eta - 6*eta*eta;
1361 :
1362 0 : default:
1363 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1364 : subtriangle_lookup(p));
1365 : }
1366 454937130 : case 3:
1367 454937130 : switch (subtriangle_lookup(p))
1368 : {
1369 301796364 : case 0:
1370 326260356 : return - 3*eta - xi*eta + 7*eta*eta;
1371 1627018 : case 1:
1372 1627018 : return - 4*xi + 4*xi*xi
1373 1627018 : + eta + 3*xi*eta - eta*eta;
1374 151636576 : case 2:
1375 163928344 : return - 3*xi + 2*xi*xi
1376 163928344 : + 4*xi*eta;
1377 :
1378 0 : default:
1379 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1380 : subtriangle_lookup(p));
1381 : }
1382 454937130 : case 4:
1383 454937130 : switch (subtriangle_lookup(p))
1384 : {
1385 301796364 : case 0:
1386 326260356 : return 1 - 3*xi + 2*xi*xi
1387 326260356 : - eta + 4*xi*eta - 7./2.*eta*eta;
1388 1627018 : case 1:
1389 1627018 : return 5./2. - 4*xi + 3./2.*xi*xi
1390 1627018 : - 9.*eta + 8*xi*eta + 13./2.*eta*eta;
1391 151636576 : case 2:
1392 163928344 : 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 454937130 : case 5:
1399 454937130 : switch (subtriangle_lookup(p))
1400 : {
1401 301796364 : case 0:
1402 326260356 : return - 1./2.*xi*eta + 1./4.*eta*eta;
1403 163935694 : case 1:
1404 163935694 : return 3./4. - 5./2.*xi + 7./4.*xi*xi
1405 163935694 : - 2*eta + 7./2.*xi*eta + 5./4.*eta*eta;
1406 1498060 : case 2:
1407 1619668 : return - 1./4.*xi*xi;
1408 :
1409 0 : default:
1410 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1411 : subtriangle_lookup(p));
1412 : }
1413 454937130 : case 6:
1414 454937130 : switch (subtriangle_lookup(p))
1415 : {
1416 301796364 : case 0:
1417 326260356 : return -xi + 2*xi*xi
1418 326260356 : + eta + 7./2.*xi*eta - 13./4.*eta*eta;
1419 163935694 : case 1:
1420 163935694 : return - 11./4. + 19./2.*xi - 23./4.*xi*xi
1421 163935694 : + 7*eta - 25./2.*xi*eta - 17./4.*eta*eta;
1422 1498060 : case 2:
1423 1619668 : return 9./4.*xi*xi;
1424 :
1425 0 : default:
1426 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1427 : subtriangle_lookup(p));
1428 : }
1429 454937130 : case 7:
1430 454937130 : switch (subtriangle_lookup(p))
1431 : {
1432 1519332 : case 0:
1433 1643004 : return -eta + 9./2.*xi*eta + 5./4.*eta*eta;
1434 163935694 : case 1:
1435 163935694 : return - 13./4. + 19./2.*xi - 25./4.*xi*xi
1436 163935694 : + 9*eta - 23./2.*xi*eta - 23./4.*eta*eta;
1437 301775092 : case 2:
1438 326237020 : 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 454937130 : case 8:
1445 454937130 : switch (subtriangle_lookup(p))
1446 : {
1447 1519332 : case 0:
1448 1643004 : return -2*eta - 1./2.*xi*eta + 13./4.*eta*eta;
1449 163935694 : case 1:
1450 163935694 : return 3./4. - 5./2.*xi + 7./4.*xi*xi
1451 163935694 : - 4*eta + 7./2.*xi*eta + 17./4.*eta*eta;
1452 301775092 : case 2:
1453 326237020 : return - 1./4.*xi*xi
1454 326237020 : - 2*eta + 3*eta*eta;
1455 :
1456 0 : default:
1457 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1458 : subtriangle_lookup(p));
1459 : }
1460 165964053 : case 9:
1461 165964053 : switch (subtriangle_lookup(p))
1462 : {
1463 5317662 : case 0:
1464 5750514 : return std::sqrt(Real(2)) * (2*xi*eta - eta*eta);
1465 168003239 : case 1:
1466 168003239 : return std::sqrt(Real(2)) * (- 3 + 10*xi - 7*xi*xi
1467 168003239 : + 8*eta - 14*xi*eta - 5*eta*eta);
1468 5243210 : case 2:
1469 5668838 : 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 165964053 : case 10:
1476 165964053 : switch (subtriangle_lookup(p))
1477 : {
1478 5317662 : case 0:
1479 5750514 : return 4*eta - 4*xi*eta - 8*eta*eta;
1480 5694563 : case 1:
1481 5694563 : return 4 - 8*xi + 4*xi*xi
1482 5694563 : - 12*eta + 12*xi*eta + 8*eta*eta;
1483 155381726 : case 2:
1484 167977514 : return 4*xi - 4*xi*xi
1485 167977514 : - 8*xi*eta;
1486 :
1487 0 : default:
1488 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1489 : subtriangle_lookup(p));
1490 : }
1491 165964053 : case 11:
1492 165964053 : switch (subtriangle_lookup(p))
1493 : {
1494 155456178 : case 0:
1495 168059190 : return 4*xi - 4*xi*xi
1496 168059190 : - 4*eta - 8*xi*eta + 10.*eta*eta;
1497 5694563 : case 1:
1498 5694563 : return 2 - 8*xi + 6*xi*xi
1499 5694563 : - 4*eta + 8*xi*eta + 2*eta*eta;
1500 5243210 : case 2:
1501 5668838 : 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 89812116 : Real clough_raw_shape(const unsigned int basis_num,
1520 : const Point & p)
1521 : {
1522 89812116 : Real xi = p(0), eta = p(1);
1523 :
1524 89812116 : switch (basis_num)
1525 : {
1526 2306443 : case 0:
1527 2306443 : switch (subtriangle_lookup(p))
1528 : {
1529 787293 : case 0:
1530 851704 : return 1 - 3*xi*xi + 2*xi*xi*xi
1531 851704 : - 3*eta*eta - 3*xi*eta*eta + 3*eta*eta*eta;
1532 816507 : case 1:
1533 816507 : return -1 + 9*xi - 15*xi*xi + 7*xi*xi*xi
1534 816507 : + 9*eta - 30*xi*eta +21*xi*xi*eta
1535 816507 : - 15*eta*eta + 21*xi*eta*eta + 7*eta*eta*eta;
1536 764420 : case 2:
1537 826570 : return 1 - 3*xi*xi + 3*xi*xi*xi
1538 826570 : - 3*xi*xi*eta
1539 826570 : - 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 2306443 : case 1:
1546 2306443 : switch (subtriangle_lookup(p))
1547 : {
1548 787293 : case 0:
1549 851704 : return 3*xi*xi - 2*xi*xi*xi
1550 851704 : + 3./2.*xi*eta*eta
1551 851704 : - 1./2.*eta*eta*eta;
1552 816507 : case 1:
1553 816507 : return 1 - 9./2.*xi + 9*xi*xi - 9./2.*xi*xi*xi
1554 816507 : - 9./2.*eta + 15*xi*eta - 21./2.*xi*xi*eta
1555 816507 : + 6*eta*eta - 21./2.*xi*eta*eta - 5./2.*eta*eta*eta;
1556 764420 : case 2:
1557 826570 : return 3*xi*xi - 5./2.*xi*xi*xi
1558 826570 : + 3./2.*xi*xi*eta;
1559 :
1560 0 : default:
1561 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1562 : subtriangle_lookup(p));
1563 : }
1564 2306443 : case 2:
1565 2306443 : switch (subtriangle_lookup(p))
1566 : {
1567 787293 : case 0:
1568 851704 : return 3*eta*eta + 3./2.*xi*eta*eta - 5./2.*eta*eta*eta;
1569 816507 : case 1:
1570 816507 : return 1 - 9./2.*xi + 6*xi*xi - 5./2.*xi*xi*xi
1571 816507 : - 9./2.*eta + 15*xi*eta - 21./2.*xi*xi*eta
1572 816507 : + 9*eta*eta - 21./2.*xi*eta*eta - 9./2.*eta*eta*eta;
1573 764420 : case 2:
1574 826570 : return -1./2.*xi*xi*xi
1575 826570 : + 3./2.*xi*xi*eta
1576 826570 : + 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 4612886 : case 3:
1583 4612886 : switch (subtriangle_lookup(p))
1584 : {
1585 1574586 : case 0:
1586 1703408 : return xi - 2*xi*xi + xi*xi*xi
1587 1703408 : - 3./2.*eta*eta - 1./2.*xi*eta*eta + 7./Real(3)*eta*eta*eta;
1588 1633014 : case 1:
1589 1633014 : return -1./Real(6) + 5./2.*xi - 9./2.*xi*xi + 13./Real(6)*xi*xi*xi
1590 1633014 : - 4*xi*eta + 4*xi*xi*eta
1591 1633014 : + 1./2.*eta*eta + 3./2.*xi*eta*eta - 1./Real(3)*eta*eta*eta;
1592 1528840 : case 2:
1593 1653140 : return xi - 1./2.*xi*xi - 7./Real(6)*xi*xi*xi
1594 1653140 : - 3*xi*eta + 2*xi*xi*eta
1595 1653140 : + 2*xi*eta*eta;
1596 :
1597 0 : default:
1598 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1599 : subtriangle_lookup(p));
1600 : }
1601 4612886 : case 4:
1602 4612886 : switch (subtriangle_lookup(p))
1603 : {
1604 1574586 : case 0:
1605 1703408 : return eta - 3*xi*eta + 2*xi*xi*eta
1606 1703408 : - 1./2.*eta*eta + 2*xi*eta*eta - 7./Real(6)*eta*eta*eta;
1607 1633014 : case 1:
1608 1633014 : return -1./Real(6) + 1./2.*xi*xi - 1./Real(3)*xi*xi*xi
1609 1633014 : + 5./2.*eta - 4*xi*eta + 3./2.*xi*xi*eta
1610 1633014 : - 9./2.*eta*eta + 4*xi*eta*eta + 13./Real(6)*eta*eta*eta;
1611 1528840 : case 2:
1612 1653140 : return -3./2.*xi*xi + 7./Real(3)*xi*xi*xi
1613 1653140 : + 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 4612886 : case 5:
1620 4612886 : switch (subtriangle_lookup(p))
1621 : {
1622 1574586 : case 0:
1623 1703408 : return -xi*xi + xi*xi*xi
1624 1703408 : - 1./4.*xi*eta*eta + 1./Real(12)*eta*eta*eta;
1625 1633014 : case 1:
1626 1633014 : return -1./Real(6) + 3./4.*xi - 2*xi*xi + 17./Real(12)*xi*xi*xi
1627 1633014 : + 3./4.*eta - 5./2.*xi*eta + 7./4.*xi*xi*eta
1628 1633014 : - eta*eta + 7./4.*xi*eta*eta + 5./Real(12)*eta*eta*eta;
1629 1528840 : case 2:
1630 1653140 : return -xi*xi + 13./Real(12)*xi*xi*xi
1631 1653140 : - 1./4.*xi*xi*eta;
1632 :
1633 0 : default:
1634 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1635 : subtriangle_lookup(p));
1636 : }
1637 4612886 : case 6:
1638 4612886 : switch (subtriangle_lookup(p))
1639 : {
1640 1574586 : case 0:
1641 1703408 : return -xi*eta + 2*xi*xi*eta
1642 1703408 : + 1./2.*eta*eta + 7./4.*xi*eta*eta - 13./Real(12)*eta*eta*eta;
1643 1633014 : case 1:
1644 1633014 : return 2./Real(3) - 13./4.*xi + 9./2.*xi*xi - 23./Real(12)*xi*xi*xi
1645 1633014 : - 11./4.*eta + 19./2.*xi*eta - 23./4.*xi*xi*eta
1646 1633014 : + 7./2.*eta*eta - 25./4.*xi*eta*eta - 17./Real(12)*eta*eta*eta;
1647 1528840 : case 2:
1648 1653140 : return -1./2.*xi*xi + 5./Real(12)*xi*xi*xi
1649 1653140 : + 9./4.*xi*xi*eta;
1650 :
1651 0 : default:
1652 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1653 : subtriangle_lookup(p));
1654 : }
1655 4612886 : case 7:
1656 4612886 : switch (subtriangle_lookup(p))
1657 : {
1658 1574586 : case 0:
1659 1703408 : return -1./2.*eta*eta + 9./4.*xi*eta*eta + 5./Real(12)*eta*eta*eta;
1660 1633014 : case 1:
1661 1633014 : return 2./Real(3) - 11./4.*xi + 7./2.*xi*xi - 17./Real(12)*xi*xi*xi
1662 1633014 : - 13./4.*eta + 19./2.*xi*eta - 25./4.*xi*xi*eta
1663 1633014 : + 9./2.*eta*eta - 23./4.*xi*eta*eta - 23./Real(12)*eta*eta*eta;
1664 1528840 : case 2:
1665 1653140 : return 1./2.*xi*xi - 13./Real(12)*xi*xi*xi
1666 1653140 : - 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 4612886 : case 8:
1673 4612886 : switch (subtriangle_lookup(p))
1674 : {
1675 1574586 : case 0:
1676 1703408 : return -eta*eta - 1./4.*xi*eta*eta + 13./Real(12)*eta*eta*eta;
1677 1633014 : case 1:
1678 1633014 : return -1./Real(6) + 3./4.*xi - xi*xi + 5./Real(12)*xi*xi*xi
1679 1633014 : + 3./4.*eta - 5./2.*xi*eta + 7./4.*xi*xi*eta
1680 1633014 : - 2*eta*eta + 7./4.*xi*eta*eta + 17./Real(12)*eta*eta*eta;
1681 1528840 : case 2:
1682 1653140 : return 1./Real(12)*xi*xi*xi
1683 1653140 : - 1./4.*xi*xi*eta
1684 1653140 : - eta*eta + eta*eta*eta;
1685 :
1686 0 : default:
1687 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1688 : subtriangle_lookup(p));
1689 : }
1690 16145101 : case 9:
1691 16145101 : switch (subtriangle_lookup(p))
1692 : {
1693 5511051 : case 0:
1694 5961928 : return std::sqrt(Real(2)) * (xi*eta*eta - 1./Real(3)*eta*eta*eta);
1695 5715549 : case 1:
1696 5715549 : return std::sqrt(Real(2)) * (2./Real(3) - 3*xi + 4*xi*xi - 5./Real(3)*xi*xi*xi
1697 5715549 : - 3*eta + 10*xi*eta - 7*xi*xi*eta
1698 5715549 : + 4*eta*eta - 7*xi*eta*eta - 5./Real(3)*eta*eta*eta);
1699 5350940 : case 2:
1700 5785990 : 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 16145101 : case 10:
1707 16145101 : switch (subtriangle_lookup(p))
1708 : {
1709 5511051 : case 0:
1710 5961928 : return 2*eta*eta - 2*xi*eta*eta - 8./Real(3)*eta*eta*eta;
1711 5715549 : case 1:
1712 5715549 : return -2./Real(3) + 2*xi - 2*xi*xi + 2./Real(3)*xi*xi*xi
1713 5715549 : + 4*eta - 8*xi*eta + 4*xi*xi*eta
1714 5715549 : - 6*eta*eta + 6*xi*eta*eta + 8./Real(3)*eta*eta*eta;
1715 5350940 : case 2:
1716 5785990 : return -2*xi*xi + 10./Real(3)*xi*xi*xi
1717 5785990 : + 4*xi*eta - 4*xi*xi*eta
1718 5785990 : - 4*xi*eta*eta;
1719 :
1720 0 : default:
1721 0 : libmesh_error_msg("Invalid subtriangle lookup = " <<
1722 : subtriangle_lookup(p));
1723 : }
1724 16145101 : case 11:
1725 16145101 : switch (subtriangle_lookup(p))
1726 : {
1727 5511051 : case 0:
1728 5961928 : return 4*xi*eta - 4*xi*xi*eta
1729 5961928 : - 2*eta*eta - 4*xi*eta*eta + 10./Real(3)*eta*eta*eta;
1730 5715549 : case 1:
1731 5715549 : return -2./Real(3) + 4*xi - 6*xi*xi + 8./Real(3)*xi*xi*xi
1732 5715549 : + 2*eta - 8*xi*eta + 6*xi*xi*eta
1733 5715549 : - 2*eta*eta + 4*xi*eta*eta + 2./Real(3)*eta*eta*eta;
1734 5350940 : case 2:
1735 5785990 : return 2*xi*xi - 8./Real(3)*xi*xi*xi
1736 5785990 : - 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 12732740 : LIBMESH_DEFAULT_VECTORIZED_FE(2,CLOUGH)
1758 :
1759 :
1760 : template <>
1761 29937372 : 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 2260056 : libmesh_assert(elem);
1768 :
1769 : CloughCoefs coefs;
1770 29937372 : clough_compute_coefs(elem, coefs);
1771 :
1772 29937372 : const ElemType type = elem->type();
1773 :
1774 : const Order totalorder =
1775 32197428 : order + add_p_level*elem->p_level();
1776 :
1777 29937372 : 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 29937372 : case THIRD:
1864 : {
1865 29937372 : switch (type)
1866 : {
1867 : // C1 functions on the Clough-Tocher triangle.
1868 29937372 : case TRI6:
1869 : case TRI7:
1870 : {
1871 2260056 : 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 29937372 : switch (i)
1877 : {
1878 : // Note: these DoF numbers are "scrambled" because my
1879 : // initial numbering conventions didn't match libMesh
1880 2494781 : case 0:
1881 2494781 : return clough_raw_shape(0, p)
1882 2494781 : + coefs.d1d2n * clough_raw_shape(10, p)
1883 2494781 : + coefs.d1d3n * clough_raw_shape(11, p);
1884 2494781 : case 3:
1885 2494781 : return clough_raw_shape(1, p)
1886 2494781 : + coefs.d2d3n * clough_raw_shape(11, p)
1887 2494781 : + coefs.d2d1n * clough_raw_shape(9, p);
1888 2494781 : case 6:
1889 2494781 : return clough_raw_shape(2, p)
1890 2494781 : + coefs.d3d1n * clough_raw_shape(9, p)
1891 2494781 : + coefs.d3d2n * clough_raw_shape(10, p);
1892 2494781 : case 1:
1893 2494781 : return coefs.d1xd1x * clough_raw_shape(3, p)
1894 2494781 : + coefs.d1xd1y * clough_raw_shape(4, p)
1895 2494781 : + coefs.d1xd2n * clough_raw_shape(10, p)
1896 2494781 : + coefs.d1xd3n * clough_raw_shape(11, p);
1897 2494781 : case 2:
1898 2494781 : return coefs.d1yd1y * clough_raw_shape(4, p)
1899 2494781 : + coefs.d1yd1x * clough_raw_shape(3, p)
1900 2494781 : + coefs.d1yd2n * clough_raw_shape(10, p)
1901 2494781 : + coefs.d1yd3n * clough_raw_shape(11, p);
1902 2494781 : case 4:
1903 2494781 : return coefs.d2xd2x * clough_raw_shape(5, p)
1904 2494781 : + coefs.d2xd2y * clough_raw_shape(6, p)
1905 2494781 : + coefs.d2xd3n * clough_raw_shape(11, p)
1906 2494781 : + coefs.d2xd1n * clough_raw_shape(9, p);
1907 2494781 : case 5:
1908 2494781 : return coefs.d2yd2y * clough_raw_shape(6, p)
1909 2494781 : + coefs.d2yd2x * clough_raw_shape(5, p)
1910 2494781 : + coefs.d2yd3n * clough_raw_shape(11, p)
1911 2494781 : + coefs.d2yd1n * clough_raw_shape(9, p);
1912 2494781 : case 7:
1913 2494781 : return coefs.d3xd3x * clough_raw_shape(7, p)
1914 2494781 : + coefs.d3xd3y * clough_raw_shape(8, p)
1915 2494781 : + coefs.d3xd1n * clough_raw_shape(9, p)
1916 2494781 : + coefs.d3xd2n * clough_raw_shape(10, p);
1917 2494781 : case 8:
1918 2494781 : return coefs.d3yd3y * clough_raw_shape(8, p)
1919 2494781 : + coefs.d3yd3x * clough_raw_shape(7, p)
1920 2494781 : + coefs.d3yd1n * clough_raw_shape(9, p)
1921 2494781 : + coefs.d3yd2n * clough_raw_shape(10, p);
1922 2494781 : case 10:
1923 2494781 : return coefs.d1nd1n * clough_raw_shape(9, p);
1924 2494781 : case 11:
1925 2494781 : return coefs.d2nd2n * clough_raw_shape(10, p);
1926 2494781 : case 9:
1927 2494781 : 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 58676280 : 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 4417296 : libmesh_assert(elem);
1978 :
1979 : CloughCoefs coefs;
1980 58676280 : clough_compute_coefs(elem, coefs);
1981 :
1982 58676280 : const ElemType type = elem->type();
1983 :
1984 : const Order totalorder =
1985 63093576 : order + add_p_level*elem->p_level();
1986 :
1987 58676280 : 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 58676280 : case THIRD:
2074 : {
2075 58676280 : switch (type)
2076 : {
2077 : // C1 functions on the Clough-Tocher triangle.
2078 58676280 : case TRI6:
2079 : case TRI7:
2080 : {
2081 4417296 : 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 58676280 : switch (i)
2087 : {
2088 : // Note: these DoF numbers are "scrambled" because my
2089 : // initial numbering conventions didn't match libMesh
2090 4889690 : case 0:
2091 4889690 : return clough_raw_shape_deriv(0, j, p)
2092 4889690 : + coefs.d1d2n * clough_raw_shape_deriv(10, j, p)
2093 4889690 : + coefs.d1d3n * clough_raw_shape_deriv(11, j, p);
2094 4889690 : case 3:
2095 4889690 : return clough_raw_shape_deriv(1, j, p)
2096 4889690 : + coefs.d2d3n * clough_raw_shape_deriv(11, j, p)
2097 4889690 : + coefs.d2d1n * clough_raw_shape_deriv(9, j, p);
2098 4889690 : case 6:
2099 4889690 : return clough_raw_shape_deriv(2, j, p)
2100 4889690 : + coefs.d3d1n * clough_raw_shape_deriv(9, j, p)
2101 4889690 : + coefs.d3d2n * clough_raw_shape_deriv(10, j, p);
2102 4889690 : case 1:
2103 4889690 : return coefs.d1xd1x * clough_raw_shape_deriv(3, j, p)
2104 4889690 : + coefs.d1xd1y * clough_raw_shape_deriv(4, j, p)
2105 4889690 : + coefs.d1xd2n * clough_raw_shape_deriv(10, j, p)
2106 4889690 : + coefs.d1xd3n * clough_raw_shape_deriv(11, j, p);
2107 4889690 : case 2:
2108 4889690 : return coefs.d1yd1y * clough_raw_shape_deriv(4, j, p)
2109 4889690 : + coefs.d1yd1x * clough_raw_shape_deriv(3, j, p)
2110 4889690 : + coefs.d1yd2n * clough_raw_shape_deriv(10, j, p)
2111 4889690 : + coefs.d1yd3n * clough_raw_shape_deriv(11, j, p);
2112 4889690 : case 4:
2113 4889690 : return coefs.d2xd2x * clough_raw_shape_deriv(5, j, p)
2114 4889690 : + coefs.d2xd2y * clough_raw_shape_deriv(6, j, p)
2115 4889690 : + coefs.d2xd3n * clough_raw_shape_deriv(11, j, p)
2116 4889690 : + coefs.d2xd1n * clough_raw_shape_deriv(9, j, p);
2117 4889690 : case 5:
2118 4889690 : return coefs.d2yd2y * clough_raw_shape_deriv(6, j, p)
2119 4889690 : + coefs.d2yd2x * clough_raw_shape_deriv(5, j, p)
2120 4889690 : + coefs.d2yd3n * clough_raw_shape_deriv(11, j, p)
2121 4889690 : + coefs.d2yd1n * clough_raw_shape_deriv(9, j, p);
2122 4889690 : case 7:
2123 4889690 : return coefs.d3xd3x * clough_raw_shape_deriv(7, j, p)
2124 4889690 : + coefs.d3xd3y * clough_raw_shape_deriv(8, j, p)
2125 4889690 : + coefs.d3xd1n * clough_raw_shape_deriv(9, j, p)
2126 4889690 : + coefs.d3xd2n * clough_raw_shape_deriv(10, j, p);
2127 4889690 : case 8:
2128 4889690 : return coefs.d3yd3y * clough_raw_shape_deriv(8, j, p)
2129 4889690 : + coefs.d3yd3x * clough_raw_shape_deriv(7, j, p)
2130 4889690 : + coefs.d3yd1n * clough_raw_shape_deriv(9, j, p)
2131 4889690 : + coefs.d3yd2n * clough_raw_shape_deriv(10, j, p);
2132 4889690 : case 10:
2133 4889690 : return coefs.d1nd1n * clough_raw_shape_deriv(9, j, p);
2134 4889690 : case 11:
2135 4889690 : return coefs.d2nd2n * clough_raw_shape_deriv(10, j, p);
2136 4889690 : case 9:
2137 4889690 : 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 73695024 : 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 5492808 : libmesh_assert(elem);
2191 :
2192 : CloughCoefs coefs;
2193 73695024 : clough_compute_coefs(elem, coefs);
2194 :
2195 73695024 : const ElemType type = elem->type();
2196 :
2197 : const Order totalorder =
2198 79187832 : order + add_p_level*elem->p_level();
2199 :
2200 73695024 : 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 73695024 : case THIRD:
2283 : {
2284 73695024 : switch (type)
2285 : {
2286 : // C1 functions on the Clough-Tocher triangle.
2287 73695024 : case TRI6:
2288 : case TRI7:
2289 : {
2290 5492808 : 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 73695024 : switch (i)
2296 : {
2297 : // Note: these DoF numbers are "scrambled" because my
2298 : // initial numbering conventions didn't match libMesh
2299 6141252 : case 0:
2300 6141252 : return clough_raw_shape_second_deriv(0, j, p)
2301 6141252 : + coefs.d1d2n * clough_raw_shape_second_deriv(10, j, p)
2302 6141252 : + coefs.d1d3n * clough_raw_shape_second_deriv(11, j, p);
2303 6141252 : case 3:
2304 6141252 : return clough_raw_shape_second_deriv(1, j, p)
2305 6141252 : + coefs.d2d3n * clough_raw_shape_second_deriv(11, j, p)
2306 6141252 : + coefs.d2d1n * clough_raw_shape_second_deriv(9, j, p);
2307 6141252 : case 6:
2308 6141252 : return clough_raw_shape_second_deriv(2, j, p)
2309 6141252 : + coefs.d3d1n * clough_raw_shape_second_deriv(9, j, p)
2310 6141252 : + coefs.d3d2n * clough_raw_shape_second_deriv(10, j, p);
2311 6141252 : case 1:
2312 6141252 : return coefs.d1xd1x * clough_raw_shape_second_deriv(3, j, p)
2313 6141252 : + coefs.d1xd1y * clough_raw_shape_second_deriv(4, j, p)
2314 6141252 : + coefs.d1xd2n * clough_raw_shape_second_deriv(10, j, p)
2315 6141252 : + coefs.d1xd3n * clough_raw_shape_second_deriv(11, j, p);
2316 6141252 : case 2:
2317 6141252 : return coefs.d1yd1y * clough_raw_shape_second_deriv(4, j, p)
2318 6141252 : + coefs.d1yd1x * clough_raw_shape_second_deriv(3, j, p)
2319 6141252 : + coefs.d1yd2n * clough_raw_shape_second_deriv(10, j, p)
2320 6141252 : + coefs.d1yd3n * clough_raw_shape_second_deriv(11, j, p);
2321 6141252 : case 4:
2322 6141252 : return coefs.d2xd2x * clough_raw_shape_second_deriv(5, j, p)
2323 6141252 : + coefs.d2xd2y * clough_raw_shape_second_deriv(6, j, p)
2324 6141252 : + coefs.d2xd3n * clough_raw_shape_second_deriv(11, j, p)
2325 6141252 : + coefs.d2xd1n * clough_raw_shape_second_deriv(9, j, p);
2326 6141252 : case 5:
2327 6141252 : return coefs.d2yd2y * clough_raw_shape_second_deriv(6, j, p)
2328 6141252 : + coefs.d2yd2x * clough_raw_shape_second_deriv(5, j, p)
2329 6141252 : + coefs.d2yd3n * clough_raw_shape_second_deriv(11, j, p)
2330 6141252 : + coefs.d2yd1n * clough_raw_shape_second_deriv(9, j, p);
2331 6141252 : case 7:
2332 6141252 : return coefs.d3xd3x * clough_raw_shape_second_deriv(7, j, p)
2333 6141252 : + coefs.d3xd3y * clough_raw_shape_second_deriv(8, j, p)
2334 6141252 : + coefs.d3xd1n * clough_raw_shape_second_deriv(9, j, p)
2335 6141252 : + coefs.d3xd2n * clough_raw_shape_second_deriv(10, j, p);
2336 6141252 : case 8:
2337 6141252 : return coefs.d3yd3y * clough_raw_shape_second_deriv(8, j, p)
2338 6141252 : + coefs.d3yd3x * clough_raw_shape_second_deriv(7, j, p)
2339 6141252 : + coefs.d3yd1n * clough_raw_shape_second_deriv(9, j, p)
2340 6141252 : + coefs.d3yd2n * clough_raw_shape_second_deriv(10, j, p);
2341 6141252 : case 10:
2342 6141252 : return coefs.d1nd1n * clough_raw_shape_second_deriv(9, j, p);
2343 6141252 : case 11:
2344 6141252 : return coefs.d2nd2n * clough_raw_shape_second_deriv(10, j, p);
2345 6141252 : case 9:
2346 6141252 : 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
|
