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

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo