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