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