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