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