fe_clough_shape_2D.C
Go to the documentation of this file.
00001 // The libMesh Finite Element Library. 00002 // Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner 00003 00004 // This library is free software; you can redistribute it and/or 00005 // modify it under the terms of the GNU Lesser General Public 00006 // License as published by the Free Software Foundation; either 00007 // version 2.1 of the License, or (at your option) any later version. 00008 00009 // This library is distributed in the hope that it will be useful, 00010 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00011 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00012 // Lesser General Public License for more details. 00013 00014 // You should have received a copy of the GNU Lesser General Public 00015 // License along with this library; if not, write to the Free Software 00016 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00017 00018 00019 // C++ inlcludes 00020 00021 // Local includes 00022 #include "libmesh/fe.h" 00023 #include "libmesh/elem.h" 00024 00025 00026 // Anonymous namespace for persistant variables. 00027 // This allows us to cache the global-to-local mapping transformation 00028 // FIXME: This should also screw up multithreading royally 00029 namespace 00030 { 00031 using namespace libMesh; 00032 00033 static dof_id_type old_elem_id = DofObject::invalid_id; 00034 // Coefficient naming: d(1)d(2n) is the coefficient of the 00035 // global shape function corresponding to value 1 in terms of the 00036 // local shape function corresponding to normal derivative 2 00037 static Real d1d2n, d1d3n, d2d3n, d2d1n, d3d1n, d3d2n; 00038 static Real d1xd1x, d1xd1y, d1xd2n, d1xd3n; 00039 static Real d1yd1x, d1yd1y, d1yd2n, d1yd3n; 00040 static Real d2xd2x, d2xd2y, d2xd3n, d2xd1n; 00041 static Real d2yd2x, d2yd2y, d2yd3n, d2yd1n; 00042 static Real d3xd3x, d3xd3y, d3xd1n, d3xd2n; 00043 static Real d3yd3x, d3yd3y, d3yd1n, d3yd2n; 00044 static Real d1nd1n, d2nd2n, d3nd3n; 00045 // Normal vector naming: N01x is the x component of the 00046 // unit vector at point 0 normal to (possibly curved) side 01 00047 static Real N01x, N01y, N10x, N10y; 00048 static Real N02x, N02y, N20x, N20y; 00049 static Real N21x, N21y, N12x, N12y; 00050 00051 Real clough_raw_shape_second_deriv(const unsigned int basis_num, 00052 const unsigned int deriv_type, 00053 const Point& p); 00054 Real clough_raw_shape_deriv(const unsigned int basis_num, 00055 const unsigned int deriv_type, 00056 const Point& p); 00057 Real clough_raw_shape(const unsigned int basis_num, 00058 const Point& p); 00059 unsigned char subtriangle_lookup(const Point& p); 00060 00061 00062 // Compute the static coefficients for an element 00063 void clough_compute_coefs(const Elem* elem) 00064 { 00065 // Using static globals for old_elem_id, etc. will fail 00066 // horribly with more than one thread. 00067 libmesh_assert_equal_to (libMesh::n_threads(), 1); 00068 00069 // Coefficients are cached from old elements 00070 if (elem->id() == old_elem_id) 00071 return; 00072 00073 old_elem_id = elem->id(); 00074 00075 const Order mapping_order (elem->default_order()); 00076 const ElemType mapping_elem_type (elem->type()); 00077 const int n_mapping_shape_functions = 00078 FE<2,LAGRANGE>::n_shape_functions(mapping_elem_type, 00079 mapping_order); 00080 00081 // Degrees of freedom are at vertices and edge midpoints 00082 std::vector<Point> dofpt; 00083 dofpt.push_back(Point(0,0)); 00084 dofpt.push_back(Point(1,0)); 00085 dofpt.push_back(Point(0,1)); 00086 dofpt.push_back(Point(1/2.,1/2.)); 00087 dofpt.push_back(Point(0,1/2.)); 00088 dofpt.push_back(Point(1/2.,0)); 00089 00090 // Mapping functions - first derivatives at each dofpt 00091 std::vector<Real> dxdxi(6), dxdeta(6), dydxi(6), dydeta(6); 00092 std::vector<Real> dxidx(6), detadx(6), dxidy(6), detady(6); 00093 00094 for (int p = 0; p != 6; ++p) 00095 { 00096 // libMesh::err << p << ' ' << dofpt[p]; 00097 for (int i = 0; i != n_mapping_shape_functions; ++i) 00098 { 00099 const Real ddxi = FE<2,LAGRANGE>::shape_deriv 00100 (mapping_elem_type, mapping_order, i, 0, dofpt[p]); 00101 const Real ddeta = FE<2,LAGRANGE>::shape_deriv 00102 (mapping_elem_type, mapping_order, i, 1, dofpt[p]); 00103 00104 // libMesh::err << ddxi << ' '; 00105 // libMesh::err << ddeta << std::endl; 00106 00107 dxdxi[p] += elem->point(i)(0) * ddxi; 00108 dydxi[p] += elem->point(i)(1) * ddxi; 00109 dxdeta[p] += elem->point(i)(0) * ddeta; 00110 dydeta[p] += elem->point(i)(1) * ddeta; 00111 } 00112 00113 // for (int i = 0; i != 12; ++i) 00114 // libMesh::err << i << ' ' << clough_raw_shape(i, dofpt[p]) << std::endl; 00115 00116 // libMesh::err << elem->point(p)(0) << ' '; 00117 // libMesh::err << elem->point(p)(1) << ' '; 00118 // libMesh::err << dxdxi[p] << ' '; 00119 // libMesh::err << dydxi[p] << ' '; 00120 // libMesh::err << dxdeta[p] << ' '; 00121 // libMesh::err << dydeta[p] << std::endl << std::endl; 00122 00123 const Real inv_jac = 1. / (dxdxi[p]*dydeta[p] - 00124 dxdeta[p]*dydxi[p]); 00125 dxidx[p] = dydeta[p] * inv_jac; 00126 dxidy[p] = - dxdeta[p] * inv_jac; 00127 detadx[p] = - dydxi[p] * inv_jac; 00128 detady[p] = dxdxi[p] * inv_jac; 00129 } 00130 00131 // Calculate midpoint normal vectors 00132 Real N1x = dydeta[3] - dydxi[3]; 00133 Real N1y = dxdxi[3] - dxdeta[3]; 00134 Real Nlength = std::sqrt(static_cast<Real>(N1x*N1x + N1y*N1y)); 00135 N1x /= Nlength; N1y /= Nlength; 00136 00137 Real N2x = - dydeta[4]; 00138 Real N2y = dxdeta[4]; 00139 Nlength = std::sqrt(static_cast<Real>(N2x*N2x + N2y*N2y)); 00140 N2x /= Nlength; N2y /= Nlength; 00141 00142 Real N3x = dydxi[5]; 00143 Real N3y = - dxdxi[5]; 00144 Nlength = std::sqrt(static_cast<Real>(N3x*N3x + N3y*N3y)); 00145 N3x /= Nlength; N3y /= Nlength; 00146 00147 // Calculate corner normal vectors (used for reduced element) 00148 N01x = dydxi[0]; 00149 N01y = - dxdxi[0]; 00150 Nlength = std::sqrt(static_cast<Real>(N01x*N01x + N01y*N01y)); 00151 N01x /= Nlength; N01y /= Nlength; 00152 00153 N10x = dydxi[1]; 00154 N10y = - dxdxi[1]; 00155 Nlength = std::sqrt(static_cast<Real>(N10x*N10x + N10y*N10y)); 00156 N10x /= Nlength; N10y /= Nlength; 00157 00158 N02x = - dydeta[0]; 00159 N02y = dxdeta[0]; 00160 Nlength = std::sqrt(static_cast<Real>(N02x*N02x + N02y*N02y)); 00161 N02x /= Nlength; N02y /= Nlength; 00162 00163 N20x = - dydeta[2]; 00164 N20y = dxdeta[2]; 00165 Nlength = std::sqrt(static_cast<Real>(N20x*N20x + N20y*N20y)); 00166 N20x /= Nlength; N20y /= Nlength; 00167 00168 N12x = dydeta[1] - dydxi[1]; 00169 N12y = dxdxi[1] - dxdeta[1]; 00170 Nlength = std::sqrt(static_cast<Real>(N12x*N12x + N12y*N12y)); 00171 N12x /= Nlength; N12y /= Nlength; 00172 00173 N21x = dydeta[1] - dydxi[1]; 00174 N21y = dxdxi[1] - dxdeta[1]; 00175 Nlength = std::sqrt(static_cast<Real>(N21x*N21x + N21y*N21y)); 00176 N21x /= Nlength; N21y /= Nlength; 00177 00178 // for (int i=0; i != 6; ++i) { 00179 // libMesh::err << elem->node(i) << ' '; 00180 // } 00181 // libMesh::err << std::endl; 00182 00183 // for (int i=0; i != 6; ++i) { 00184 // libMesh::err << elem->point(i)(0) << ' '; 00185 // libMesh::err << elem->point(i)(1) << ' '; 00186 // } 00187 // libMesh::err << std::endl; 00188 00189 00190 // give normal vectors a globally unique orientation 00191 00192 if (elem->point(2) < elem->point(1)) 00193 { 00194 // libMesh::err << "Flipping nodes " << elem->node(2); 00195 // libMesh::err << " and " << elem->node(1); 00196 // libMesh::err << " around node " << elem->node(4); 00197 // libMesh::err << std::endl; 00198 N1x = -N1x; N1y = -N1y; 00199 N12x = -N12x; N12y = -N12y; 00200 N21x = -N21x; N21y = -N21y; 00201 } 00202 else 00203 { 00204 // libMesh::err << "Not flipping nodes " << elem->node(2); 00205 // libMesh::err << " and " << elem->node(1); 00206 // libMesh::err << " around node " << elem->node(4); 00207 // libMesh::err << std::endl; 00208 } 00209 if (elem->point(0) < elem->point(2)) 00210 { 00211 // libMesh::err << "Flipping nodes " << elem->node(0); 00212 // libMesh::err << " and " << elem->node(2); 00213 // libMesh::err << " around node " << elem->node(5); 00214 // libMesh::err << std::endl; 00215 // libMesh::err << N2x << ' ' << N2y << std::endl; 00216 N2x = -N2x; N2y = -N2y; 00217 N02x = -N02x; N02y = -N02y; 00218 N20x = -N20x; N20y = -N20y; 00219 // libMesh::err << N2x << ' ' << N2y << std::endl; 00220 } 00221 else 00222 { 00223 // libMesh::err << "Not flipping nodes " << elem->node(0); 00224 // libMesh::err << " and " << elem->node(2); 00225 // libMesh::err << " around node " << elem->node(5); 00226 // libMesh::err << std::endl; 00227 } 00228 if (elem->point(1) < elem->point(0)) 00229 { 00230 // libMesh::err << "Flipping nodes " << elem->node(1); 00231 // libMesh::err << " and " << elem->node(0); 00232 // libMesh::err << " around node " << elem->node(3); 00233 // libMesh::err << std::endl; 00234 N3x = -N3x; 00235 N3y = -N3y; 00236 N01x = -N01x; N01y = -N01y; 00237 N10x = -N10x; N10y = -N10y; 00238 } 00239 else 00240 { 00241 // libMesh::err << "Not flipping nodes " << elem->node(1); 00242 // libMesh::err << " and " << elem->node(0); 00243 // libMesh::err << " around node " << elem->node(3); 00244 // libMesh::err << std::endl; 00245 } 00246 00247 // libMesh::err << N2x << ' ' << N2y << std::endl; 00248 00249 // Cache basis function gradients 00250 // FIXME: the raw_shape calls shouldn't be done on every element! 00251 // FIXME: I should probably be looping, too... 00252 // Gradient naming: d(1)d(2n)d(xi) is the xi component of the 00253 // gradient of the 00254 // local basis function corresponding to value 1 at the node 00255 // corresponding to normal vector 2 00256 00257 Real d1d2ndxi = clough_raw_shape_deriv(0, 0, dofpt[4]); 00258 Real d1d2ndeta = clough_raw_shape_deriv(0, 1, dofpt[4]); 00259 Real d1d2ndx = d1d2ndxi * dxidx[4] + d1d2ndeta * detadx[4]; 00260 Real d1d2ndy = d1d2ndxi * dxidy[4] + d1d2ndeta * detady[4]; 00261 Real d1d3ndxi = clough_raw_shape_deriv(0, 0, dofpt[5]); 00262 Real d1d3ndeta = clough_raw_shape_deriv(0, 1, dofpt[5]); 00263 Real d1d3ndx = d1d3ndxi * dxidx[5] + d1d3ndeta * detadx[5]; 00264 Real d1d3ndy = d1d3ndxi * dxidy[5] + d1d3ndeta * detady[5]; 00265 Real d2d3ndxi = clough_raw_shape_deriv(1, 0, dofpt[5]); 00266 Real d2d3ndeta = clough_raw_shape_deriv(1, 1, dofpt[5]); 00267 Real d2d3ndx = d2d3ndxi * dxidx[5] + d2d3ndeta * detadx[5]; 00268 Real d2d3ndy = d2d3ndxi * dxidy[5] + d2d3ndeta * detady[5]; 00269 Real d2d1ndxi = clough_raw_shape_deriv(1, 0, dofpt[3]); 00270 Real d2d1ndeta = clough_raw_shape_deriv(1, 1, dofpt[3]); 00271 Real d2d1ndx = d2d1ndxi * dxidx[3] + d2d1ndeta * detadx[3]; 00272 Real d2d1ndy = d2d1ndxi * dxidy[3] + d2d1ndeta * detady[3]; 00273 Real d3d1ndxi = clough_raw_shape_deriv(2, 0, dofpt[3]); 00274 Real d3d1ndeta = clough_raw_shape_deriv(2, 1, dofpt[3]); 00275 Real d3d1ndx = d3d1ndxi * dxidx[3] + d3d1ndeta * detadx[3]; 00276 Real d3d1ndy = d3d1ndxi * dxidy[3] + d3d1ndeta * detady[3]; 00277 Real d3d2ndxi = clough_raw_shape_deriv(2, 0, dofpt[4]); 00278 Real d3d2ndeta = clough_raw_shape_deriv(2, 1, dofpt[4]); 00279 Real d3d2ndx = d3d2ndxi * dxidx[4] + d3d2ndeta * detadx[4]; 00280 Real d3d2ndy = d3d2ndxi * dxidy[4] + d3d2ndeta * detady[4]; 00281 Real d1xd2ndxi = clough_raw_shape_deriv(3, 0, dofpt[4]); 00282 Real d1xd2ndeta = clough_raw_shape_deriv(3, 1, dofpt[4]); 00283 Real d1xd2ndx = d1xd2ndxi * dxidx[4] + d1xd2ndeta * detadx[4]; 00284 Real d1xd2ndy = d1xd2ndxi * dxidy[4] + d1xd2ndeta * detady[4]; 00285 Real d1xd3ndxi = clough_raw_shape_deriv(3, 0, dofpt[5]); 00286 Real d1xd3ndeta = clough_raw_shape_deriv(3, 1, dofpt[5]); 00287 Real d1xd3ndx = d1xd3ndxi * dxidx[5] + d1xd3ndeta * detadx[5]; 00288 Real d1xd3ndy = d1xd3ndxi * dxidy[5] + d1xd3ndeta * detady[5]; 00289 Real d1yd2ndxi = clough_raw_shape_deriv(4, 0, dofpt[4]); 00290 Real d1yd2ndeta = clough_raw_shape_deriv(4, 1, dofpt[4]); 00291 Real d1yd2ndx = d1yd2ndxi * dxidx[4] + d1yd2ndeta * detadx[4]; 00292 Real d1yd2ndy = d1yd2ndxi * dxidy[4] + d1yd2ndeta * detady[4]; 00293 Real d1yd3ndxi = clough_raw_shape_deriv(4, 0, dofpt[5]); 00294 Real d1yd3ndeta = clough_raw_shape_deriv(4, 1, dofpt[5]); 00295 Real d1yd3ndx = d1yd3ndxi * dxidx[5] + d1yd3ndeta * detadx[5]; 00296 Real d1yd3ndy = d1yd3ndxi * dxidy[5] + d1yd3ndeta * detady[5]; 00297 Real d2xd3ndxi = clough_raw_shape_deriv(5, 0, dofpt[5]); 00298 Real d2xd3ndeta = clough_raw_shape_deriv(5, 1, dofpt[5]); 00299 Real d2xd3ndx = d2xd3ndxi * dxidx[5] + d2xd3ndeta * detadx[5]; 00300 Real d2xd3ndy = d2xd3ndxi * dxidy[5] + d2xd3ndeta * detady[5]; 00301 Real d2xd1ndxi = clough_raw_shape_deriv(5, 0, dofpt[3]); 00302 Real d2xd1ndeta = clough_raw_shape_deriv(5, 1, dofpt[3]); 00303 Real d2xd1ndx = d2xd1ndxi * dxidx[3] + d2xd1ndeta * detadx[3]; 00304 Real d2xd1ndy = d2xd1ndxi * dxidy[3] + d2xd1ndeta * detady[3]; 00305 Real d2yd3ndxi = clough_raw_shape_deriv(6, 0, dofpt[5]); 00306 Real d2yd3ndeta = clough_raw_shape_deriv(6, 1, dofpt[5]); 00307 Real d2yd3ndx = d2yd3ndxi * dxidx[5] + d2yd3ndeta * detadx[5]; 00308 Real d2yd3ndy = d2yd3ndxi * dxidy[5] + d2yd3ndeta * detady[5]; 00309 Real d2yd1ndxi = clough_raw_shape_deriv(6, 0, dofpt[3]); 00310 Real d2yd1ndeta = clough_raw_shape_deriv(6, 1, dofpt[3]); 00311 Real d2yd1ndx = d2yd1ndxi * dxidx[3] + d2yd1ndeta * detadx[3]; 00312 Real d2yd1ndy = d2yd1ndxi * dxidy[3] + d2yd1ndeta * detady[3]; 00313 Real d3xd1ndxi = clough_raw_shape_deriv(7, 0, dofpt[3]); 00314 Real d3xd1ndeta = clough_raw_shape_deriv(7, 1, dofpt[3]); 00315 Real d3xd1ndx = d3xd1ndxi * dxidx[3] + d3xd1ndeta * detadx[3]; 00316 Real d3xd1ndy = d3xd1ndxi * dxidy[3] + d3xd1ndeta * detady[3]; 00317 Real d3xd2ndxi = clough_raw_shape_deriv(7, 0, dofpt[4]); 00318 Real d3xd2ndeta = clough_raw_shape_deriv(7, 1, dofpt[4]); 00319 Real d3xd2ndx = d3xd2ndxi * dxidx[4] + d3xd2ndeta * detadx[4]; 00320 Real d3xd2ndy = d3xd2ndxi * dxidy[4] + d3xd2ndeta * detady[4]; 00321 Real d3yd1ndxi = clough_raw_shape_deriv(8, 0, dofpt[3]); 00322 Real d3yd1ndeta = clough_raw_shape_deriv(8, 1, dofpt[3]); 00323 Real d3yd1ndx = d3yd1ndxi * dxidx[3] + d3yd1ndeta * detadx[3]; 00324 Real d3yd1ndy = d3yd1ndxi * dxidy[3] + d3yd1ndeta * detady[3]; 00325 Real d3yd2ndxi = clough_raw_shape_deriv(8, 0, dofpt[4]); 00326 Real d3yd2ndeta = clough_raw_shape_deriv(8, 1, dofpt[4]); 00327 Real d3yd2ndx = d3yd2ndxi * dxidx[4] + d3yd2ndeta * detadx[4]; 00328 Real d3yd2ndy = d3yd2ndxi * dxidy[4] + d3yd2ndeta * detady[4]; 00329 Real d1nd1ndxi = clough_raw_shape_deriv(9, 0, dofpt[3]); 00330 Real d1nd1ndeta = clough_raw_shape_deriv(9, 1, dofpt[3]); 00331 Real d1nd1ndx = d1nd1ndxi * dxidx[3] + d1nd1ndeta * detadx[3]; 00332 Real d1nd1ndy = d1nd1ndxi * dxidy[3] + d1nd1ndeta * detady[3]; 00333 Real d2nd2ndxi = clough_raw_shape_deriv(10, 0, dofpt[4]); 00334 Real d2nd2ndeta = clough_raw_shape_deriv(10, 1, dofpt[4]); 00335 Real d2nd2ndx = d2nd2ndxi * dxidx[4] + d2nd2ndeta * detadx[4]; 00336 Real d2nd2ndy = d2nd2ndxi * dxidy[4] + d2nd2ndeta * detady[4]; 00337 Real d3nd3ndxi = clough_raw_shape_deriv(11, 0, dofpt[5]); 00338 Real d3nd3ndeta = clough_raw_shape_deriv(11, 1, dofpt[5]); 00339 Real d3nd3ndx = d3nd3ndxi * dxidx[3] + d3nd3ndeta * detadx[3]; 00340 Real d3nd3ndy = d3nd3ndxi * dxidy[3] + d3nd3ndeta * detady[3]; 00341 00342 Real d1xd1dxi = clough_raw_shape_deriv(3, 0, dofpt[0]); 00343 Real d1xd1deta = clough_raw_shape_deriv(3, 1, dofpt[0]); 00344 Real d1xd1dx = d1xd1dxi * dxidx[0] + d1xd1deta * detadx[0]; 00345 Real d1xd1dy = d1xd1dxi * dxidy[0] + d1xd1deta * detady[0]; 00346 Real d1yd1dxi = clough_raw_shape_deriv(4, 0, dofpt[0]); 00347 Real d1yd1deta = clough_raw_shape_deriv(4, 1, dofpt[0]); 00348 Real d1yd1dx = d1yd1dxi * dxidx[0] + d1yd1deta * detadx[0]; 00349 Real d1yd1dy = d1yd1dxi * dxidy[0] + d1yd1deta * detady[0]; 00350 Real d2xd2dxi = clough_raw_shape_deriv(5, 0, dofpt[1]); 00351 Real d2xd2deta = clough_raw_shape_deriv(5, 1, dofpt[1]); 00352 Real d2xd2dx = d2xd2dxi * dxidx[1] + d2xd2deta * detadx[1]; 00353 Real d2xd2dy = d2xd2dxi * dxidy[1] + d2xd2deta * detady[1]; 00354 00355 // libMesh::err << dofpt[4](0) << ' '; 00356 // libMesh::err << dofpt[4](1) << ' '; 00357 // libMesh::err << (int)subtriangle_lookup(dofpt[5]) << ' '; 00358 // libMesh::err << dxdxi[4] << ' '; 00359 // libMesh::err << dxdeta[4] << ' '; 00360 // libMesh::err << dydxi[4] << ' '; 00361 // libMesh::err << dydeta[4] << ' '; 00362 // libMesh::err << dxidx[4] << ' '; 00363 // libMesh::err << dxidy[4] << ' '; 00364 // libMesh::err << detadx[4] << ' '; 00365 // libMesh::err << detady[4] << ' '; 00366 // libMesh::err << N1x << ' '; 00367 // libMesh::err << N1y << ' '; 00368 // libMesh::err << d2yd1ndxi << ' '; 00369 // libMesh::err << d2yd1ndeta << ' '; 00370 // libMesh::err << d2yd1ndx << ' '; 00371 // libMesh::err << d2yd1ndy << std::endl; 00372 00373 Real d2yd2dxi = clough_raw_shape_deriv(6, 0, dofpt[1]); 00374 Real d2yd2deta = clough_raw_shape_deriv(6, 1, dofpt[1]); 00375 Real d2yd2dx = d2yd2dxi * dxidx[1] + d2yd2deta * detadx[1]; 00376 Real d2yd2dy = d2yd2dxi * dxidy[1] + d2yd2deta * detady[1]; 00377 Real d3xd3dxi = clough_raw_shape_deriv(7, 0, dofpt[2]); 00378 Real d3xd3deta = clough_raw_shape_deriv(7, 1, dofpt[2]); 00379 Real d3xd3dx = d3xd3dxi * dxidx[1] + d3xd3deta * detadx[1]; 00380 Real d3xd3dy = d3xd3dxi * dxidy[1] + d3xd3deta * detady[1]; 00381 Real d3yd3dxi = clough_raw_shape_deriv(8, 0, dofpt[2]); 00382 Real d3yd3deta = clough_raw_shape_deriv(8, 1, dofpt[2]); 00383 Real d3yd3dx = d3yd3dxi * dxidx[1] + d3yd3deta * detadx[1]; 00384 Real d3yd3dy = d3yd3dxi * dxidy[1] + d3yd3deta * detady[1]; 00385 00386 // Calculate normal derivatives 00387 00388 Real d1nd1ndn = d1nd1ndx * N1x + d1nd1ndy * N1y; 00389 Real d1xd2ndn = d1xd2ndx * N2x + d1xd2ndy * N2y; 00390 Real d1xd3ndn = d1xd3ndx * N3x + d1xd3ndy * N3y; 00391 Real d1yd2ndn = d1yd2ndx * N2x + d1yd2ndy * N2y; 00392 Real d1yd3ndn = d1yd3ndx * N3x + d1yd3ndy * N3y; 00393 00394 Real d2nd2ndn = d2nd2ndx * N2x + d2nd2ndy * N2y; 00395 Real d2xd3ndn = d2xd3ndx * N3x + d2xd3ndy * N3y; 00396 Real d2xd1ndn = d2xd1ndx * N1x + d2xd1ndy * N1y; 00397 Real d2yd3ndn = d2yd3ndx * N3x + d2yd3ndy * N3y; 00398 Real d2yd1ndn = d2yd1ndx * N1x + d2yd1ndy * N1y; 00399 00400 Real d3nd3ndn = d3nd3ndx * N3x + d3nd3ndy * N3y; 00401 Real d3xd1ndn = d3xd1ndx * N1x + d3xd1ndy * N1y; 00402 Real d3xd2ndn = d3xd2ndx * N2x + d3xd2ndy * N2y; 00403 Real d3yd1ndn = d3yd1ndx * N1x + d3yd1ndy * N1y; 00404 Real d3yd2ndn = d3yd2ndx * N2x + d3yd2ndy * N2y; 00405 00406 // Calculate midpoint scaling factors 00407 00408 d1nd1n = 1. / d1nd1ndn; 00409 d2nd2n = 1. / d2nd2ndn; 00410 d3nd3n = 1. / d3nd3ndn; 00411 00412 // Calculate midpoint derivative adjustments to nodal value 00413 // interpolant functions 00414 00415 d1d2n = -(d1d2ndx * N2x + d1d2ndy * N2y) / d2nd2ndn; 00416 d1d3n = -(d1d3ndx * N3x + d1d3ndy * N3y) / d3nd3ndn; 00417 d2d3n = -(d2d3ndx * N3x + d2d3ndy * N3y) / d3nd3ndn; 00418 d2d1n = -(d2d1ndx * N1x + d2d1ndy * N1y) / d1nd1ndn; 00419 d3d1n = -(d3d1ndx * N1x + d3d1ndy * N1y) / d1nd1ndn; 00420 d3d2n = -(d3d2ndx * N2x + d3d2ndy * N2y) / d2nd2ndn; 00421 00422 // Calculate nodal derivative scaling factors 00423 00424 d1xd1x = 1. / (d1xd1dx - d1xd1dy * d1yd1dx / d1yd1dy); 00425 d1xd1y = 1. / (d1yd1dx - d1xd1dx * d1yd1dy / d1xd1dy); 00426 // d1xd1y = - d1xd1x * (d1xd1dy / d1yd1dy); 00427 d1yd1y = 1. / (d1yd1dy - d1yd1dx * d1xd1dy / d1xd1dx); 00428 d1yd1x = 1. / (d1xd1dy - d1yd1dy * d1xd1dx / d1yd1dx); 00429 // d1yd1x = - d1yd1y * (d1yd1dx / d1xd1dx); 00430 d2xd2x = 1. / (d2xd2dx - d2xd2dy * d2yd2dx / d2yd2dy); 00431 d2xd2y = 1. / (d2yd2dx - d2xd2dx * d2yd2dy / d2xd2dy); 00432 // d2xd2y = - d2xd2x * (d2xd2dy / d2yd2dy); 00433 d2yd2y = 1. / (d2yd2dy - d2yd2dx * d2xd2dy / d2xd2dx); 00434 d2yd2x = 1. / (d2xd2dy - d2yd2dy * d2xd2dx / d2yd2dx); 00435 // d2yd2x = - d2yd2y * (d2yd2dx / d2xd2dx); 00436 d3xd3x = 1. / (d3xd3dx - d3xd3dy * d3yd3dx / d3yd3dy); 00437 d3xd3y = 1. / (d3yd3dx - d3xd3dx * d3yd3dy / d3xd3dy); 00438 // d3xd3y = - d3xd3x * (d3xd3dy / d3yd3dy); 00439 d3yd3y = 1. / (d3yd3dy - d3yd3dx * d3xd3dy / d3xd3dx); 00440 d3yd3x = 1. / (d3xd3dy - d3yd3dy * d3xd3dx / d3yd3dx); 00441 // d3yd3x = - d3yd3y * (d3yd3dx / d3xd3dx); 00442 00443 // libMesh::err << d1xd1dx << ' '; 00444 // libMesh::err << d1xd1dy << ' '; 00445 // libMesh::err << d1yd1dx << ' '; 00446 // libMesh::err << d1yd1dy << ' '; 00447 // libMesh::err << d2xd2dx << ' '; 00448 // libMesh::err << d2xd2dy << ' '; 00449 // libMesh::err << d2yd2dx << ' '; 00450 // libMesh::err << d2yd2dy << ' '; 00451 // libMesh::err << d3xd3dx << ' '; 00452 // libMesh::err << d3xd3dy << ' '; 00453 // libMesh::err << d3yd3dx << ' '; 00454 // libMesh::err << d3yd3dy << std::endl; 00455 00456 // Calculate midpoint derivative adjustments to nodal derivative 00457 // interpolant functions 00458 00459 d1xd2n = -(d1xd1x * d1xd2ndn + d1xd1y * d1yd2ndn) / d2nd2ndn; 00460 d1yd2n = -(d1yd1y * d1yd2ndn + d1yd1x * d1xd2ndn) / d2nd2ndn; 00461 d1xd3n = -(d1xd1x * d1xd3ndn + d1xd1y * d1yd3ndn) / d3nd3ndn; 00462 d1yd3n = -(d1yd1y * d1yd3ndn + d1yd1x * d1xd3ndn) / d3nd3ndn; 00463 d2xd3n = -(d2xd2x * d2xd3ndn + d2xd2y * d2yd3ndn) / d3nd3ndn; 00464 d2yd3n = -(d2yd2y * d2yd3ndn + d2yd2x * d2xd3ndn) / d3nd3ndn; 00465 d2xd1n = -(d2xd2x * d2xd1ndn + d2xd2y * d2yd1ndn) / d1nd1ndn; 00466 d2yd1n = -(d2yd2y * d2yd1ndn + d2yd2x * d2xd1ndn) / d1nd1ndn; 00467 d3xd1n = -(d3xd3x * d3xd1ndn + d3xd3y * d3yd1ndn) / d1nd1ndn; 00468 d3yd1n = -(d3yd3y * d3yd1ndn + d3yd3x * d3xd1ndn) / d1nd1ndn; 00469 d3xd2n = -(d3xd3x * d3xd2ndn + d3xd3y * d3yd2ndn) / d2nd2ndn; 00470 d3yd2n = -(d3yd3y * d3yd2ndn + d3yd3x * d3xd2ndn) / d2nd2ndn; 00471 00472 // Cross your fingers 00473 // libMesh::err << d1nd1ndn << ' '; 00474 // libMesh::err << d2xd1ndn << ' '; 00475 // libMesh::err << d2yd1ndn << ' '; 00476 // libMesh::err << std::endl; 00477 00478 // libMesh::err << "Transform variables: "; 00479 // libMesh::err << d1nd1n << ' '; 00480 // libMesh::err << d2nd2n << ' '; 00481 // libMesh::err << d3nd3n << ' '; 00482 // libMesh::err << d1d2n << ' '; 00483 // libMesh::err << d1d3n << ' '; 00484 // libMesh::err << d2d3n << ' '; 00485 // libMesh::err << d2d1n << ' '; 00486 // libMesh::err << d3d1n << ' '; 00487 // libMesh::err << d3d2n << std::endl; 00488 // libMesh::err << d1xd1x << ' '; 00489 // libMesh::err << d1xd1y << ' '; 00490 // libMesh::err << d1yd1x << ' '; 00491 // libMesh::err << d1yd1y << ' '; 00492 // libMesh::err << d2xd2x << ' '; 00493 // libMesh::err << d2xd2y << ' '; 00494 // libMesh::err << d2yd2x << ' '; 00495 // libMesh::err << d2yd2y << ' '; 00496 // libMesh::err << d3xd3x << ' '; 00497 // libMesh::err << d3xd3y << ' '; 00498 // libMesh::err << d3yd3x << ' '; 00499 // libMesh::err << d3yd3y << std::endl; 00500 // libMesh::err << d1xd2n << ' '; 00501 // libMesh::err << d1yd2n << ' '; 00502 // libMesh::err << d1xd3n << ' '; 00503 // libMesh::err << d1yd3n << ' '; 00504 // libMesh::err << d2xd3n << ' '; 00505 // libMesh::err << d2yd3n << ' '; 00506 // libMesh::err << d2xd1n << ' '; 00507 // libMesh::err << d2yd1n << ' '; 00508 // libMesh::err << d3xd1n << ' '; 00509 // libMesh::err << d3yd1n << ' '; 00510 // libMesh::err << d3xd2n << ' '; 00511 // libMesh::err << d3yd2n << ' '; 00512 // libMesh::err << std::endl; 00513 // libMesh::err << std::endl; 00514 } 00515 00516 00517 unsigned char subtriangle_lookup(const Point& p) 00518 { 00519 if ((p(0) >= p(1)) && (p(0) + 2 * p(1) <= 1)) 00520 return 0; 00521 if ((p(0) <= p(1)) && (p(1) + 2 * p(0) <= 1)) 00522 return 2; 00523 return 1; 00524 } 00525 00526 // Return shape function second derivatives on the unit right 00527 // triangle 00528 Real clough_raw_shape_second_deriv(const unsigned int basis_num, 00529 const unsigned int deriv_type, 00530 const Point& p) 00531 { 00532 Real xi = p(0), eta = p(1); 00533 00534 switch (deriv_type) 00535 { 00536 00537 // second derivative in xi-xi direction 00538 case 0: 00539 switch (basis_num) 00540 { 00541 case 0: 00542 switch (subtriangle_lookup(p)) 00543 { 00544 case 0: 00545 return -6 + 12*xi; 00546 case 1: 00547 return -30 + 42*xi + 42*eta; 00548 case 2: 00549 return -6 + 18*xi - 6*eta; 00550 00551 default: 00552 libmesh_error_msg("Invalid subtriangle lookup = " << 00553 subtriangle_lookup(p)); 00554 } 00555 case 1: 00556 switch (subtriangle_lookup(p)) 00557 { 00558 case 0: 00559 return 6 - 12*xi; 00560 case 1: 00561 return 18 - 27*xi - 21*eta; 00562 case 2: 00563 return 6 - 15*xi + 3*eta; 00564 00565 default: 00566 libmesh_error_msg("Invalid subtriangle lookup = " << 00567 subtriangle_lookup(p)); 00568 } 00569 case 2: 00570 switch (subtriangle_lookup(p)) 00571 { 00572 case 0: 00573 return 0; 00574 case 1: 00575 return 12 - 15*xi - 21*eta; 00576 case 2: 00577 return -3*xi + 3*eta; 00578 00579 default: 00580 libmesh_error_msg("Invalid subtriangle lookup = " << 00581 subtriangle_lookup(p)); 00582 } 00583 case 3: 00584 switch (subtriangle_lookup(p)) 00585 { 00586 case 0: 00587 return -4 + 6*xi; 00588 case 1: 00589 return -9 + 13*xi + 8*eta; 00590 case 2: 00591 return -1 - 7*xi + 4*eta; 00592 00593 default: 00594 libmesh_error_msg("Invalid subtriangle lookup = " << 00595 subtriangle_lookup(p)); 00596 } 00597 case 4: 00598 switch (subtriangle_lookup(p)) 00599 { 00600 case 0: 00601 return 4*eta; 00602 case 1: 00603 return 1 - 2*xi + 3*eta; 00604 case 2: 00605 return -3 + 14*xi - eta; 00606 00607 default: 00608 libmesh_error_msg("Invalid subtriangle lookup = " << 00609 subtriangle_lookup(p)); 00610 } 00611 case 5: 00612 switch (subtriangle_lookup(p)) 00613 { 00614 case 0: 00615 return -2 + 6*xi; 00616 case 1: 00617 return -4 + 17./2.*xi + 7./2.*eta; 00618 case 2: 00619 return -2 + 13./2.*xi - 1./2.*eta; 00620 00621 default: 00622 libmesh_error_msg("Invalid subtriangle lookup = " << 00623 subtriangle_lookup(p)); 00624 } 00625 case 6: 00626 switch (subtriangle_lookup(p)) 00627 { 00628 case 0: 00629 return 4*eta; 00630 case 1: 00631 return 9 - 23./2.*xi - 23./2.*eta; 00632 case 2: 00633 return -1 + 5./2.*xi + 9./2.*eta; 00634 00635 default: 00636 libmesh_error_msg("Invalid subtriangle lookup = " << 00637 subtriangle_lookup(p)); 00638 } 00639 case 7: 00640 switch (subtriangle_lookup(p)) 00641 { 00642 case 0: 00643 return 0; 00644 case 1: 00645 return 7 - 17./2.*xi - 25./2.*eta; 00646 case 2: 00647 return 1 - 13./2.*xi + 7./2.*eta; 00648 00649 default: 00650 libmesh_error_msg("Invalid subtriangle lookup = " << 00651 subtriangle_lookup(p)); 00652 } 00653 case 8: 00654 switch (subtriangle_lookup(p)) 00655 { 00656 case 0: 00657 return 0; 00658 case 1: 00659 return -2 + 5./2.*xi + 7./2.*eta; 00660 case 2: 00661 return 1./2.*xi - 1./2.*eta; 00662 00663 default: 00664 libmesh_error_msg("Invalid subtriangle lookup = " << 00665 subtriangle_lookup(p)); 00666 } 00667 case 9: 00668 switch (subtriangle_lookup(p)) 00669 { 00670 case 0: 00671 return 0; 00672 case 1: 00673 return std::sqrt(2.) * (8 - 10*xi - 14*eta); 00674 case 2: 00675 return std::sqrt(2.) * (-2*xi + 2*eta); 00676 00677 default: 00678 libmesh_error_msg("Invalid subtriangle lookup = " << 00679 subtriangle_lookup(p)); 00680 } 00681 case 10: 00682 switch (subtriangle_lookup(p)) 00683 { 00684 case 0: 00685 return 0; 00686 case 1: 00687 return -4 + 4*xi + 8*eta; 00688 case 2: 00689 return -4 + 20*xi - 8*eta; 00690 00691 default: 00692 libmesh_error_msg("Invalid subtriangle lookup = " << 00693 subtriangle_lookup(p)); 00694 } 00695 case 11: 00696 switch (subtriangle_lookup(p)) 00697 { 00698 case 0: 00699 return -8*eta; 00700 case 1: 00701 return -12 + 16*xi + 12*eta; 00702 case 2: 00703 return 4 - 16*xi - 4*eta; 00704 00705 default: 00706 libmesh_error_msg("Invalid subtriangle lookup = " << 00707 subtriangle_lookup(p)); 00708 } 00709 00710 default: 00711 libmesh_error_msg("Invalid shape function index i = " << 00712 basis_num); 00713 } 00714 00715 // second derivative in xi-eta direction 00716 case 1: 00717 switch (basis_num) 00718 { 00719 case 0: 00720 switch (subtriangle_lookup(p)) 00721 { 00722 case 0: 00723 return -6*eta; 00724 case 1: 00725 return -30 +42*xi 00726 + 42*eta; 00727 case 2: 00728 return -6*xi; 00729 00730 default: 00731 libmesh_error_msg("Invalid subtriangle lookup = " << 00732 subtriangle_lookup(p)); 00733 } 00734 case 1: 00735 switch (subtriangle_lookup(p)) 00736 { 00737 case 0: 00738 return + 3*eta; 00739 case 1: 00740 return 15 - 21*xi - 21*eta; 00741 case 2: 00742 return 3*xi; 00743 00744 default: 00745 libmesh_error_msg("Invalid subtriangle lookup = " << 00746 subtriangle_lookup(p)); 00747 } 00748 case 2: 00749 switch (subtriangle_lookup(p)) 00750 { 00751 case 0: 00752 return 3*eta; 00753 case 1: 00754 return 15 - 21*xi - 21*eta; 00755 case 2: 00756 return 3*xi; 00757 00758 default: 00759 libmesh_error_msg("Invalid subtriangle lookup = " << 00760 subtriangle_lookup(p)); 00761 } 00762 case 3: 00763 switch (subtriangle_lookup(p)) 00764 { 00765 case 0: 00766 return -eta; 00767 case 1: 00768 return -4 + 8*xi + 3*eta; 00769 case 2: 00770 return -3 + 4*xi + 4*eta; 00771 00772 default: 00773 libmesh_error_msg("Invalid subtriangle lookup = " << 00774 subtriangle_lookup(p)); 00775 } 00776 case 4: 00777 switch (subtriangle_lookup(p)) 00778 { 00779 case 0: 00780 return -3 + 4*xi + 4*eta; 00781 case 1: 00782 return - 4 + 3*xi + 8*eta; 00783 case 2: 00784 return -xi; 00785 00786 default: 00787 libmesh_error_msg("Invalid subtriangle lookup = " << 00788 subtriangle_lookup(p)); 00789 } 00790 case 5: 00791 switch (subtriangle_lookup(p)) 00792 { 00793 case 0: 00794 return - 1./2.*eta; 00795 case 1: 00796 return -5./2. + 7./2.*xi + 7./2.*eta; 00797 case 2: 00798 return - 1./2.*xi; 00799 00800 default: 00801 libmesh_error_msg("Invalid subtriangle lookup = " << 00802 subtriangle_lookup(p)); 00803 } 00804 case 6: 00805 switch (subtriangle_lookup(p)) 00806 { 00807 case 0: 00808 return -1 + 4*xi + 7./2.*eta; 00809 case 1: 00810 return 19./2. - 23./2.*xi - 25./2.*eta; 00811 case 2: 00812 return 9./2.*xi; 00813 00814 default: 00815 libmesh_error_msg("Invalid subtriangle lookup = " << 00816 subtriangle_lookup(p)); 00817 } 00818 case 7: 00819 switch (subtriangle_lookup(p)) 00820 { 00821 case 0: 00822 return 9./2.*eta; 00823 case 1: 00824 return 19./2. - 25./2.*xi - 23./2.*eta; 00825 case 2: 00826 return -1 + 7./2.*xi + 4*eta; 00827 00828 default: 00829 libmesh_error_msg("Invalid subtriangle lookup = " << 00830 subtriangle_lookup(p)); 00831 } 00832 case 8: 00833 switch (subtriangle_lookup(p)) 00834 { 00835 case 0: 00836 return -1./2.*eta; 00837 case 1: 00838 return -5./2. + 7./2.*xi + 7./2.*eta; 00839 case 2: 00840 return -1./2.*xi; 00841 00842 default: 00843 libmesh_error_msg("Invalid subtriangle lookup = " << 00844 subtriangle_lookup(p)); 00845 } 00846 case 9: 00847 switch (subtriangle_lookup(p)) 00848 { 00849 case 0: 00850 return std::sqrt(2.) * (2*eta); 00851 case 1: 00852 return std::sqrt(2.) * (10 - 14*xi - 14*eta); 00853 case 2: 00854 return std::sqrt(2.) * (2*xi); 00855 00856 default: 00857 libmesh_error_msg("Invalid subtriangle lookup = " << 00858 subtriangle_lookup(p)); 00859 } 00860 case 10: 00861 switch (subtriangle_lookup(p)) 00862 { 00863 case 0: 00864 return -4*eta; 00865 case 1: 00866 return - 8 + 8*xi + 12*eta; 00867 case 2: 00868 return 4 - 8*xi - 8*eta; 00869 00870 default: 00871 libmesh_error_msg("Invalid subtriangle lookup = " << 00872 subtriangle_lookup(p)); 00873 } 00874 case 11: 00875 switch (subtriangle_lookup(p)) 00876 { 00877 case 0: 00878 return 4 - 8*xi - 8*eta; 00879 case 1: 00880 return -8 + 12*xi + 8*eta; 00881 case 2: 00882 return -4*xi; 00883 00884 default: 00885 libmesh_error_msg("Invalid subtriangle lookup = " << 00886 subtriangle_lookup(p)); 00887 } 00888 00889 default: 00890 libmesh_error_msg("Invalid shape function index i = " << 00891 basis_num); 00892 } 00893 00894 // second derivative in eta-eta direction 00895 case 2: 00896 switch (basis_num) 00897 { 00898 case 0: 00899 switch (subtriangle_lookup(p)) 00900 { 00901 case 0: 00902 return -6 - 6*xi + 18*eta; 00903 case 1: 00904 return -30 + 42*xi + 42*eta; 00905 case 2: 00906 return -6 + 12*eta; 00907 00908 default: 00909 libmesh_error_msg("Invalid subtriangle lookup = " << 00910 subtriangle_lookup(p)); 00911 } 00912 case 1: 00913 switch (subtriangle_lookup(p)) 00914 { 00915 case 0: 00916 return 3*xi - 3*eta; 00917 case 1: 00918 return 12 - 21*xi - 15*eta; 00919 case 2: 00920 return 0; 00921 00922 default: 00923 libmesh_error_msg("Invalid subtriangle lookup = " << 00924 subtriangle_lookup(p)); 00925 } 00926 case 2: 00927 switch (subtriangle_lookup(p)) 00928 { 00929 case 0: 00930 return 6 + 3*xi - 15*eta; 00931 case 1: 00932 return 18 - 21.*xi - 27*eta; 00933 case 2: 00934 return 6 - 12*eta; 00935 00936 default: 00937 libmesh_error_msg("Invalid subtriangle lookup = " << 00938 subtriangle_lookup(p)); 00939 } 00940 case 3: 00941 switch (subtriangle_lookup(p)) 00942 { 00943 case 0: 00944 return -3 - xi + 14*eta; 00945 case 1: 00946 return 1 + 3*xi - 2*eta; 00947 case 2: 00948 return 4*xi; 00949 00950 default: 00951 libmesh_error_msg("Invalid subtriangle lookup = " << 00952 subtriangle_lookup(p)); 00953 } 00954 case 4: 00955 switch (subtriangle_lookup(p)) 00956 { 00957 case 0: 00958 return -1 + 4*xi - 7*eta; 00959 case 1: 00960 return -9 + 8*xi + 13*eta; 00961 case 2: 00962 return -4 + 6*eta; 00963 00964 default: 00965 libmesh_error_msg("Invalid subtriangle lookup = " << 00966 subtriangle_lookup(p)); 00967 } 00968 case 5: 00969 switch (subtriangle_lookup(p)) 00970 { 00971 case 0: 00972 return - 1./2.*xi + 1./2.*eta; 00973 case 1: 00974 return -2 + 7./2.*xi + 5./2.*eta; 00975 case 2: 00976 return 0; 00977 00978 default: 00979 libmesh_error_msg("Invalid subtriangle lookup = " << 00980 subtriangle_lookup(p)); 00981 } 00982 case 6: 00983 switch (subtriangle_lookup(p)) 00984 { 00985 case 0: 00986 return 1 + 7./2.*xi - 13./2.*eta; 00987 case 1: 00988 return 7 - 25./2.*xi - 17./2.*eta; 00989 case 2: 00990 return 0; 00991 00992 default: 00993 libmesh_error_msg("Invalid subtriangle lookup = " << 00994 subtriangle_lookup(p)); 00995 } 00996 case 7: 00997 switch (subtriangle_lookup(p)) 00998 { 00999 case 0: 01000 return -1 + 9./2.*xi + 5./2.*eta; 01001 case 1: 01002 return 9 - 23./2.*xi - 23./2.*eta; 01003 case 2: 01004 return 4*xi; 01005 01006 default: 01007 libmesh_error_msg("Invalid subtriangle lookup = " << 01008 subtriangle_lookup(p)); 01009 } 01010 case 8: 01011 switch (subtriangle_lookup(p)) 01012 { 01013 case 0: 01014 return -2 - 1./2.*xi + 13./2.*eta; 01015 case 1: 01016 return -4 + 7./2.*xi + 17./2.*eta; 01017 case 2: 01018 return -2 + 6*eta; 01019 01020 default: 01021 libmesh_error_msg("Invalid subtriangle lookup = " << 01022 subtriangle_lookup(p)); 01023 } 01024 case 9: 01025 switch (subtriangle_lookup(p)) 01026 { 01027 case 0: 01028 return std::sqrt(2.) * (2*xi - 2*eta); 01029 case 1: 01030 return std::sqrt(2.) * (8 - 14*xi - 10*eta); 01031 case 2: 01032 return 0; 01033 01034 default: 01035 libmesh_error_msg("Invalid subtriangle lookup = " << 01036 subtriangle_lookup(p)); 01037 } 01038 case 10: 01039 switch (subtriangle_lookup(p)) 01040 { 01041 case 0: 01042 return 4 - 4*xi - 16*eta; 01043 case 1: 01044 return -12 + 12*xi + 16*eta; 01045 case 2: 01046 return -8*xi; 01047 01048 default: 01049 libmesh_error_msg("Invalid subtriangle lookup = " << 01050 subtriangle_lookup(p)); 01051 } 01052 case 11: 01053 switch (subtriangle_lookup(p)) 01054 { 01055 case 0: 01056 return -4 - 8*xi + 20*eta; 01057 case 1: 01058 return -4 + 8*xi + 4*eta; 01059 case 2: 01060 return 0; 01061 01062 default: 01063 libmesh_error_msg("Invalid subtriangle lookup = " << 01064 subtriangle_lookup(p)); 01065 } 01066 01067 default: 01068 libmesh_error_msg("Invalid shape function index i = " << 01069 basis_num); 01070 } 01071 01072 default: 01073 libmesh_error_msg("Invalid shape function derivative j = " << 01074 deriv_type); 01075 } 01076 01077 libmesh_error_msg("We'll never get here!"); 01078 return 0.; 01079 } 01080 01081 01082 01083 Real clough_raw_shape_deriv(const unsigned int basis_num, 01084 const unsigned int deriv_type, 01085 const Point& p) 01086 { 01087 Real xi = p(0), eta = p(1); 01088 01089 switch (deriv_type) 01090 { 01091 // derivative in xi direction 01092 case 0: 01093 switch (basis_num) 01094 { 01095 case 0: 01096 switch (subtriangle_lookup(p)) 01097 { 01098 case 0: 01099 return -6*xi + 6*xi*xi 01100 - 3*eta*eta; 01101 case 1: 01102 return 9 - 30*xi + 21*xi*xi 01103 - 30*eta + 42*xi*eta 01104 + 21*eta*eta; 01105 case 2: 01106 return -6*xi + 9*xi*xi 01107 - 6*xi*eta; 01108 01109 default: 01110 libmesh_error_msg("Invalid subtriangle lookup = " << 01111 subtriangle_lookup(p)); 01112 } 01113 case 1: 01114 switch (subtriangle_lookup(p)) 01115 { 01116 case 0: 01117 return 6*xi - 6*xi*xi 01118 + 3./2.*eta*eta; 01119 case 1: 01120 return - 9./2. + 18*xi - 27./2.*xi*xi 01121 + 15*eta - 21*xi*eta 01122 - 21./2.*eta*eta; 01123 case 2: 01124 return 6*xi - 15./2.*xi*xi 01125 + 3*xi*eta; 01126 01127 default: 01128 libmesh_error_msg("Invalid subtriangle lookup = " << 01129 subtriangle_lookup(p)); 01130 } 01131 case 2: 01132 switch (subtriangle_lookup(p)) 01133 { 01134 case 0: 01135 return 3./2.*eta*eta; 01136 case 1: 01137 return - 9./2. + 12*xi - 15./2.*xi*xi 01138 + 15*eta - 21*xi*eta 01139 - 21./2.*eta*eta; 01140 case 2: 01141 return -3./2.*xi*xi 01142 + 3*xi*eta; 01143 01144 default: 01145 libmesh_error_msg("Invalid subtriangle lookup = " << 01146 subtriangle_lookup(p)); 01147 } 01148 case 3: 01149 switch (subtriangle_lookup(p)) 01150 { 01151 case 0: 01152 return 1 - 4*xi + 3*xi*xi 01153 - 1./2.*eta*eta; 01154 case 1: 01155 return 5./2. - 9*xi + 13./2.*xi*xi 01156 - 4*eta + 8*xi*eta 01157 + 3./2.*eta*eta; 01158 case 2: 01159 return 1 - xi - 7./2.*xi*xi 01160 - 3*eta + 4*xi*eta 01161 + 2*eta*eta; 01162 01163 default: 01164 libmesh_error_msg("Invalid subtriangle lookup = " << 01165 subtriangle_lookup(p)); 01166 } 01167 case 4: 01168 switch (subtriangle_lookup(p)) 01169 { 01170 case 0: 01171 return - 3*eta + 4*xi*eta 01172 + 2*eta*eta; 01173 case 1: 01174 return xi - xi*xi 01175 - 4*eta + 3*xi*eta 01176 + 4*eta*eta; 01177 case 2: 01178 return -3*xi + 7*xi*xi 01179 - xi*eta; 01180 01181 default: 01182 libmesh_error_msg("Invalid subtriangle lookup = " << 01183 subtriangle_lookup(p)); 01184 } 01185 case 5: 01186 switch (subtriangle_lookup(p)) 01187 { 01188 case 0: 01189 return -2*xi + 3*xi*xi 01190 - 1./4.*eta*eta; 01191 case 1: 01192 return 3./4. - 4*xi + 17./4.*xi*xi 01193 - 5./2.*eta + 7./2.*xi*eta 01194 + 7./4.*eta*eta; 01195 case 2: 01196 return -2*xi + 13./4.*xi*xi 01197 - 1./2.*xi*eta; 01198 01199 default: 01200 libmesh_error_msg("Invalid subtriangle lookup = " << 01201 subtriangle_lookup(p)); 01202 } 01203 case 6: 01204 switch (subtriangle_lookup(p)) 01205 { 01206 case 0: 01207 return -eta + 4*xi*eta 01208 + 7./4.*eta*eta; 01209 case 1: 01210 return -13./4. + 9*xi - 23./4.*xi*xi 01211 + 19./2.*eta - 23./2.*xi*eta 01212 - 25./4.*eta*eta; 01213 case 2: 01214 return -xi + 5./4.*xi*xi 01215 + 9./2.*xi*eta; 01216 01217 default: 01218 libmesh_error_msg("Invalid subtriangle lookup = " << 01219 subtriangle_lookup(p)); 01220 } 01221 case 7: 01222 switch (subtriangle_lookup(p)) 01223 { 01224 case 0: 01225 return 9./4.*eta*eta; 01226 case 1: 01227 return - 11./4. + 7*xi - 17./4.*xi*xi 01228 + 19./2.*eta - 25./2.*xi*eta 01229 - 23./4.*eta*eta; 01230 case 2: 01231 return xi - 13./4.*xi*xi 01232 - eta + 7./2.*xi*eta + 2*eta*eta; 01233 01234 default: 01235 libmesh_error_msg("Invalid subtriangle lookup = " << 01236 subtriangle_lookup(p)); 01237 } 01238 case 8: 01239 switch (subtriangle_lookup(p)) 01240 { 01241 case 0: 01242 return - 1./4.*eta*eta; 01243 case 1: 01244 return 3./4. - 2*xi + 5./4.*xi*xi 01245 - 5./2.*eta + 7./2.*xi*eta 01246 + 7./4.*eta*eta; 01247 case 2: 01248 return 1./4.*xi*xi 01249 - 1./2.*xi*eta; 01250 01251 default: 01252 libmesh_error_msg("Invalid subtriangle lookup = " << 01253 subtriangle_lookup(p)); 01254 } 01255 case 9: 01256 switch (subtriangle_lookup(p)) 01257 { 01258 case 0: 01259 return std::sqrt(2.) * eta*eta; 01260 case 1: 01261 return std::sqrt(2.) * (-3 + 8*xi - 5*xi*xi 01262 + 10*eta - 14*xi*eta 01263 - 7*eta*eta); 01264 case 2: 01265 return std::sqrt(2.) * (-xi*xi + 2*xi*eta); 01266 01267 default: 01268 libmesh_error_msg("Invalid subtriangle lookup = " << 01269 subtriangle_lookup(p)); 01270 } 01271 case 10: 01272 switch (subtriangle_lookup(p)) 01273 { 01274 case 0: 01275 return -2*eta*eta; 01276 case 1: 01277 return 2 - 4*xi + 2*xi*xi 01278 - 8*eta + 8*xi*eta 01279 + 6*eta*eta; 01280 case 2: 01281 return -4*xi + 10*xi*xi 01282 + 4*eta - 8*xi*eta 01283 - 4*eta*eta; 01284 01285 default: 01286 libmesh_error_msg("Invalid subtriangle lookup = " << 01287 subtriangle_lookup(p)); 01288 } 01289 case 11: 01290 switch (subtriangle_lookup(p)) 01291 { 01292 case 0: 01293 return 4*eta - 8*xi*eta 01294 - 4*eta*eta; 01295 case 1: 01296 return 4 - 12*xi + 8*xi*xi 01297 - 8*eta + 12*xi*eta 01298 + 4*eta*eta; 01299 case 2: 01300 return 4*xi - 8*xi*xi 01301 - 4*xi*eta; 01302 01303 default: 01304 libmesh_error_msg("Invalid subtriangle lookup = " << 01305 subtriangle_lookup(p)); 01306 } 01307 01308 default: 01309 libmesh_error_msg("Invalid shape function index i = " << 01310 basis_num); 01311 } 01312 01313 // derivative in eta direction 01314 case 1: 01315 switch (basis_num) 01316 { 01317 case 0: 01318 switch (subtriangle_lookup(p)) 01319 { 01320 case 0: 01321 return - 6*eta - 6*xi*eta + 9*eta*eta; 01322 case 1: 01323 return 9 - 30*xi + 21*xi*xi 01324 - 30*eta + 42*xi*eta + 21*eta*eta; 01325 case 2: 01326 return - 3*xi*xi 01327 - 6*eta + 6*eta*eta; 01328 01329 default: 01330 libmesh_error_msg("Invalid subtriangle lookup = " << 01331 subtriangle_lookup(p)); 01332 } 01333 case 1: 01334 switch (subtriangle_lookup(p)) 01335 { 01336 case 0: 01337 return + 3*xi*eta 01338 - 3./2.*eta*eta; 01339 case 1: 01340 return - 9./2. + 15*xi - 21./2.*xi*xi 01341 + 12*eta - 21*xi*eta - 15./2.*eta*eta; 01342 case 2: 01343 return + 3./2.*xi*xi; 01344 01345 default: 01346 libmesh_error_msg("Invalid subtriangle lookup = " << 01347 subtriangle_lookup(p)); 01348 } 01349 case 2: 01350 switch (subtriangle_lookup(p)) 01351 { 01352 case 0: 01353 return 6*eta + 3*xi*eta - 15./2.*eta*eta; 01354 case 1: 01355 return - 9./2. + 15*xi - 21./2.*xi*xi 01356 + 18*eta - 21.*xi*eta - 27./2.*eta*eta; 01357 case 2: 01358 return 3./2.*xi*xi 01359 + 6*eta - 6*eta*eta; 01360 01361 default: 01362 libmesh_error_msg("Invalid subtriangle lookup = " << 01363 subtriangle_lookup(p)); 01364 } 01365 case 3: 01366 switch (subtriangle_lookup(p)) 01367 { 01368 case 0: 01369 return - 3*eta - xi*eta + 7*eta*eta; 01370 case 1: 01371 return - 4*xi + 4*xi*xi 01372 + eta + 3*xi*eta - eta*eta; 01373 case 2: 01374 return - 3*xi + 2*xi*xi 01375 + 4*xi*eta; 01376 01377 default: 01378 libmesh_error_msg("Invalid subtriangle lookup = " << 01379 subtriangle_lookup(p)); 01380 } 01381 case 4: 01382 switch (subtriangle_lookup(p)) 01383 { 01384 case 0: 01385 return 1 - 3*xi + 2*xi*xi 01386 - eta + 4*xi*eta - 7./2.*eta*eta; 01387 case 1: 01388 return 5./2. - 4*xi + 3./2.*xi*xi 01389 - 9.*eta + 8*xi*eta + 13./2.*eta*eta; 01390 case 2: 01391 return 1 - 1./2.*xi*xi - 4*eta + 3*eta*eta; 01392 01393 default: 01394 libmesh_error_msg("Invalid subtriangle lookup = " << 01395 subtriangle_lookup(p)); 01396 } 01397 case 5: 01398 switch (subtriangle_lookup(p)) 01399 { 01400 case 0: 01401 return - 1./2.*xi*eta + 1./4.*eta*eta; 01402 case 1: 01403 return 3./4. - 5./2.*xi + 7./4.*xi*xi 01404 - 2*eta + 7./2.*xi*eta + 5./4.*eta*eta; 01405 case 2: 01406 return - 1./4.*xi*xi; 01407 01408 default: 01409 libmesh_error_msg("Invalid subtriangle lookup = " << 01410 subtriangle_lookup(p)); 01411 } 01412 case 6: 01413 switch (subtriangle_lookup(p)) 01414 { 01415 case 0: 01416 return -xi + 2*xi*xi 01417 + eta + 7./2.*xi*eta - 13./4.*eta*eta; 01418 case 1: 01419 return - 11./4. + 19./2.*xi - 23./4.*xi*xi 01420 + 7*eta - 25./2.*xi*eta - 17./4.*eta*eta; 01421 case 2: 01422 return 9./4.*xi*xi; 01423 01424 default: 01425 libmesh_error_msg("Invalid subtriangle lookup = " << 01426 subtriangle_lookup(p)); 01427 } 01428 case 7: 01429 switch (subtriangle_lookup(p)) 01430 { 01431 case 0: 01432 return -eta + 9./2.*xi*eta + 5./4.*eta*eta; 01433 case 1: 01434 return - 13./4. + 19./2.*xi - 25./4.*xi*xi 01435 + 9*eta - 23./2.*xi*eta - 23./4.*eta*eta; 01436 case 2: 01437 return - xi + 7./4.*xi*xi + 4*xi*eta; 01438 01439 default: 01440 libmesh_error_msg("Invalid subtriangle lookup = " << 01441 subtriangle_lookup(p)); 01442 } 01443 case 8: 01444 switch (subtriangle_lookup(p)) 01445 { 01446 case 0: 01447 return -2*eta - 1./2.*xi*eta + 13./4.*eta*eta; 01448 case 1: 01449 return 3./4. - 5./2.*xi + 7./4.*xi*xi 01450 - 4*eta + 7./2.*xi*eta + 17./4.*eta*eta; 01451 case 2: 01452 return - 1./4.*xi*xi 01453 - 2*eta + 3*eta*eta; 01454 01455 default: 01456 libmesh_error_msg("Invalid subtriangle lookup = " << 01457 subtriangle_lookup(p)); 01458 } 01459 case 9: 01460 switch (subtriangle_lookup(p)) 01461 { 01462 case 0: 01463 return std::sqrt(2.) * (2*xi*eta - eta*eta); 01464 case 1: 01465 return std::sqrt(2.) * (- 3 + 10*xi - 7*xi*xi 01466 + 8*eta - 14*xi*eta - 5*eta*eta); 01467 case 2: 01468 return std::sqrt(2.) * (xi*xi); 01469 01470 default: 01471 libmesh_error_msg("Invalid subtriangle lookup = " << 01472 subtriangle_lookup(p)); 01473 } 01474 case 10: 01475 switch (subtriangle_lookup(p)) 01476 { 01477 case 0: 01478 return 4*eta - 4*xi*eta - 8*eta*eta; 01479 case 1: 01480 return 4 - 8*xi + 4*xi*xi 01481 - 12*eta + 12*xi*eta + 8*eta*eta; 01482 case 2: 01483 return 4*xi - 4*xi*xi 01484 - 8*xi*eta; 01485 01486 default: 01487 libmesh_error_msg("Invalid subtriangle lookup = " << 01488 subtriangle_lookup(p)); 01489 } 01490 case 11: 01491 switch (subtriangle_lookup(p)) 01492 { 01493 case 0: 01494 return 4*xi - 4*xi*xi 01495 - 4*eta - 8*xi*eta + 10.*eta*eta; 01496 case 1: 01497 return 2 - 8*xi + 6*xi*xi 01498 - 4*eta + 8*xi*eta + 2*eta*eta; 01499 case 2: 01500 return - 2*xi*xi; 01501 01502 default: 01503 libmesh_error_msg("Invalid subtriangle lookup = " << 01504 subtriangle_lookup(p)); 01505 } 01506 01507 default: 01508 libmesh_error_msg("Invalid shape function index i = " << 01509 basis_num); 01510 } 01511 01512 default: 01513 libmesh_error_msg("Invalid shape function derivative j = " << 01514 deriv_type); 01515 } 01516 01517 libmesh_error_msg("We'll never get here!"); 01518 return 0.; 01519 } 01520 01521 Real clough_raw_shape(const unsigned int basis_num, 01522 const Point& p) 01523 { 01524 Real xi = p(0), eta = p(1); 01525 01526 switch (basis_num) 01527 { 01528 case 0: 01529 switch (subtriangle_lookup(p)) 01530 { 01531 case 0: 01532 return 1 - 3*xi*xi + 2*xi*xi*xi 01533 - 3*eta*eta - 3*xi*eta*eta + 3*eta*eta*eta; 01534 case 1: 01535 return -1 + 9*xi - 15*xi*xi + 7*xi*xi*xi 01536 + 9*eta - 30*xi*eta +21*xi*xi*eta 01537 - 15*eta*eta + 21*xi*eta*eta + 7*eta*eta*eta; 01538 case 2: 01539 return 1 - 3*xi*xi + 3*xi*xi*xi 01540 - 3*xi*xi*eta 01541 - 3*eta*eta + 2*eta*eta*eta; 01542 01543 default: 01544 libmesh_error_msg("Invalid subtriangle lookup = " << 01545 subtriangle_lookup(p)); 01546 } 01547 case 1: 01548 switch (subtriangle_lookup(p)) 01549 { 01550 case 0: 01551 return 3*xi*xi - 2*xi*xi*xi 01552 + 3./2.*xi*eta*eta 01553 - 1./2.*eta*eta*eta; 01554 case 1: 01555 return 1 - 9./2.*xi + 9*xi*xi - 9./2.*xi*xi*xi 01556 - 9./2.*eta + 15*xi*eta - 21./2.*xi*xi*eta 01557 + 6*eta*eta - 21./2.*xi*eta*eta - 5./2.*eta*eta*eta; 01558 case 2: 01559 return 3*xi*xi - 5./2.*xi*xi*xi 01560 + 3./2.*xi*xi*eta; 01561 01562 default: 01563 libmesh_error_msg("Invalid subtriangle lookup = " << 01564 subtriangle_lookup(p)); 01565 } 01566 case 2: 01567 switch (subtriangle_lookup(p)) 01568 { 01569 case 0: 01570 return 3*eta*eta + 3./2.*xi*eta*eta - 5./2.*eta*eta*eta; 01571 case 1: 01572 return 1 - 9./2.*xi + 6*xi*xi - 5./2.*xi*xi*xi 01573 - 9./2.*eta + 15*xi*eta - 21./2.*xi*xi*eta 01574 + 9*eta*eta - 21./2.*xi*eta*eta - 9./2.*eta*eta*eta; 01575 case 2: 01576 return -1./2.*xi*xi*xi 01577 + 3./2.*xi*xi*eta 01578 + 3*eta*eta - 2*eta*eta*eta; 01579 01580 default: 01581 libmesh_error_msg("Invalid subtriangle lookup = " << 01582 subtriangle_lookup(p)); 01583 } 01584 case 3: 01585 switch (subtriangle_lookup(p)) 01586 { 01587 case 0: 01588 return xi - 2*xi*xi + xi*xi*xi 01589 - 3./2.*eta*eta - 1./2.*xi*eta*eta + 7./3.*eta*eta*eta; 01590 case 1: 01591 return -1./6. + 5./2.*xi - 9./2.*xi*xi + 13./6.*xi*xi*xi 01592 - 4*xi*eta + 4*xi*xi*eta 01593 + 1./2.*eta*eta + 3./2.*xi*eta*eta - 1./3.*eta*eta*eta; 01594 case 2: 01595 return xi - 1./2.*xi*xi - 7./6.*xi*xi*xi 01596 - 3*xi*eta + 2*xi*xi*eta 01597 + 2*xi*eta*eta; 01598 01599 default: 01600 libmesh_error_msg("Invalid subtriangle lookup = " << 01601 subtriangle_lookup(p)); 01602 } 01603 case 4: 01604 switch (subtriangle_lookup(p)) 01605 { 01606 case 0: 01607 return eta - 3*xi*eta + 2*xi*xi*eta 01608 - 1./2.*eta*eta + 2*xi*eta*eta - 7./6.*eta*eta*eta; 01609 case 1: 01610 return -1./6. + 1./2.*xi*xi - 1./3.*xi*xi*xi 01611 + 5./2.*eta - 4*xi*eta + 3./2.*xi*xi*eta 01612 - 9./2.*eta*eta + 4*xi*eta*eta + 13./6.*eta*eta*eta; 01613 case 2: 01614 return -3./2.*xi*xi + 7./3.*xi*xi*xi 01615 + eta - 1./2.*xi*xi*eta - 2*eta*eta + eta*eta*eta; 01616 01617 default: 01618 libmesh_error_msg("Invalid subtriangle lookup = " << 01619 subtriangle_lookup(p)); 01620 } 01621 case 5: 01622 switch (subtriangle_lookup(p)) 01623 { 01624 case 0: 01625 return -xi*xi + xi*xi*xi 01626 - 1./4.*xi*eta*eta + 1./12.*eta*eta*eta; 01627 case 1: 01628 return -1./6. + 3./4.*xi - 2*xi*xi + 17./12.*xi*xi*xi 01629 + 3./4.*eta - 5./2.*xi*eta + 7./4.*xi*xi*eta 01630 - eta*eta + 7./4.*xi*eta*eta + 5./12.*eta*eta*eta; 01631 case 2: 01632 return -xi*xi + 13./12.*xi*xi*xi 01633 - 1./4.*xi*xi*eta; 01634 01635 default: 01636 libmesh_error_msg("Invalid subtriangle lookup = " << 01637 subtriangle_lookup(p)); 01638 } 01639 case 6: 01640 switch (subtriangle_lookup(p)) 01641 { 01642 case 0: 01643 return -xi*eta + 2*xi*xi*eta 01644 + 1./2.*eta*eta + 7./4.*xi*eta*eta - 13./12.*eta*eta*eta; 01645 case 1: 01646 return 2./3. - 13./4.*xi + 9./2.*xi*xi - 23./12.*xi*xi*xi 01647 - 11./4.*eta + 19./2.*xi*eta - 23./4.*xi*xi*eta 01648 + 7./2.*eta*eta - 25./4.*xi*eta*eta - 17./12.*eta*eta*eta; 01649 case 2: 01650 return -1./2.*xi*xi + 5./12.*xi*xi*xi 01651 + 9./4.*xi*xi*eta; 01652 01653 default: 01654 libmesh_error_msg("Invalid subtriangle lookup = " << 01655 subtriangle_lookup(p)); 01656 } 01657 case 7: 01658 switch (subtriangle_lookup(p)) 01659 { 01660 case 0: 01661 return -1./2.*eta*eta + 9./4.*xi*eta*eta + 5./12.*eta*eta*eta; 01662 case 1: 01663 return 2./3. - 11./4.*xi + 7./2.*xi*xi - 17./12.*xi*xi*xi 01664 - 13./4.*eta + 19./2.*xi*eta - 25./4.*xi*xi*eta 01665 + 9./2.*eta*eta - 23./4.*xi*eta*eta - 23./12.*eta*eta*eta; 01666 case 2: 01667 return 1./2.*xi*xi - 13./12.*xi*xi*xi 01668 - xi*eta + 7./4.*xi*xi*eta + 2*xi*eta*eta; 01669 01670 default: 01671 libmesh_error_msg("Invalid subtriangle lookup = " << 01672 subtriangle_lookup(p)); 01673 } 01674 case 8: 01675 switch (subtriangle_lookup(p)) 01676 { 01677 case 0: 01678 return -eta*eta - 1./4.*xi*eta*eta + 13./12.*eta*eta*eta; 01679 case 1: 01680 return -1./6. + 3./4.*xi - xi*xi + 5./12.*xi*xi*xi 01681 + 3./4.*eta - 5./2.*xi*eta + 7./4.*xi*xi*eta 01682 - 2*eta*eta + 7./4.*xi*eta*eta + 17./12.*eta*eta*eta; 01683 case 2: 01684 return 1./12.*xi*xi*xi 01685 - 1./4.*xi*xi*eta 01686 - eta*eta + eta*eta*eta; 01687 01688 default: 01689 libmesh_error_msg("Invalid subtriangle lookup = " << 01690 subtriangle_lookup(p)); 01691 } 01692 case 9: 01693 switch (subtriangle_lookup(p)) 01694 { 01695 case 0: 01696 return std::sqrt(2.) * (xi*eta*eta - 1./3.*eta*eta*eta); 01697 case 1: 01698 return std::sqrt(2.) * (2./3. - 3*xi + 4*xi*xi - 5./3.*xi*xi*xi 01699 - 3*eta + 10*xi*eta - 7*xi*xi*eta 01700 + 4*eta*eta - 7*xi*eta*eta - 5./3.*eta*eta*eta); 01701 case 2: 01702 return std::sqrt(2.) * (-1./3.*xi*xi*xi + xi*xi*eta); 01703 01704 default: 01705 libmesh_error_msg("Invalid subtriangle lookup = " << 01706 subtriangle_lookup(p)); 01707 } 01708 case 10: 01709 switch (subtriangle_lookup(p)) 01710 { 01711 case 0: 01712 return 2*eta*eta - 2*xi*eta*eta - 8./3.*eta*eta*eta; 01713 case 1: 01714 return -2./3. + 2*xi - 2*xi*xi + 2./3.*xi*xi*xi 01715 + 4*eta - 8*xi*eta + 4*xi*xi*eta 01716 - 6*eta*eta + 6*xi*eta*eta + 8./3.*eta*eta*eta; 01717 case 2: 01718 return -2*xi*xi + 10./3.*xi*xi*xi 01719 + 4*xi*eta - 4*xi*xi*eta 01720 - 4*xi*eta*eta; 01721 01722 default: 01723 libmesh_error_msg("Invalid subtriangle lookup = " << 01724 subtriangle_lookup(p)); 01725 } 01726 case 11: 01727 switch (subtriangle_lookup(p)) 01728 { 01729 case 0: 01730 return 4*xi*eta - 4*xi*xi*eta 01731 - 2*eta*eta - 4*xi*eta*eta + 10./3.*eta*eta*eta; 01732 case 1: 01733 return -2./3. + 4*xi - 6*xi*xi + 8./3.*xi*xi*xi 01734 + 2*eta - 8*xi*eta + 6*xi*xi*eta 01735 - 2*eta*eta + 4*xi*eta*eta + 2./3.*eta*eta*eta; 01736 case 2: 01737 return 2*xi*xi - 8./3.*xi*xi*xi 01738 - 2*xi*xi*eta; 01739 01740 default: 01741 libmesh_error_msg("Invalid subtriangle lookup = " << 01742 subtriangle_lookup(p)); 01743 } 01744 01745 default: 01746 libmesh_error_msg("Invalid shape function index i = " << 01747 basis_num); 01748 } 01749 01750 libmesh_error_msg("We'll never get here!"); 01751 return 0.; 01752 } 01753 01754 01755 } // end anonymous namespace 01756 01757 01758 namespace libMesh 01759 { 01760 01761 01762 template <> 01763 Real FE<2,CLOUGH>::shape(const ElemType, 01764 const Order, 01765 const unsigned int, 01766 const Point&) 01767 { 01768 libmesh_error_msg("Clough-Tocher elements require the real element \nto construct gradient-based degrees of freedom."); 01769 return 0.; 01770 } 01771 01772 01773 01774 template <> 01775 Real FE<2,CLOUGH>::shape(const Elem* elem, 01776 const Order order, 01777 const unsigned int i, 01778 const Point& p) 01779 { 01780 libmesh_assert(elem); 01781 01782 clough_compute_coefs(elem); 01783 01784 const ElemType type = elem->type(); 01785 01786 const Order totalorder = static_cast<Order>(order + elem->p_level()); 01787 01788 switch (totalorder) 01789 { 01790 // 2nd-order restricted Clough-Tocher element 01791 case SECOND: 01792 { 01793 // There may be a bug in the 2nd order case; the 3rd order 01794 // Clough-Tocher elements are pretty uniformly better anyways 01795 // so use those instead. 01796 libmesh_experimental(); 01797 switch (type) 01798 { 01799 // C1 functions on the Clough-Tocher triangle. 01800 case TRI6: 01801 { 01802 libmesh_assert_less (i, 9); 01803 // FIXME: it would be nice to calculate (and cache) 01804 // clough_raw_shape(j,p) only once per triangle, not 1-7 01805 // times 01806 switch (i) 01807 { 01808 // Note: these DoF numbers are "scrambled" because my 01809 // initial numbering conventions didn't match libMesh 01810 case 0: 01811 return clough_raw_shape(0, p) 01812 + d1d2n * clough_raw_shape(10, p) 01813 + d1d3n * clough_raw_shape(11, p); 01814 case 3: 01815 return clough_raw_shape(1, p) 01816 + d2d3n * clough_raw_shape(11, p) 01817 + d2d1n * clough_raw_shape(9, p); 01818 case 6: 01819 return clough_raw_shape(2, p) 01820 + d3d1n * clough_raw_shape(9, p) 01821 + d3d2n * clough_raw_shape(10, p); 01822 case 1: 01823 return d1xd1x * clough_raw_shape(3, p) 01824 + d1xd1y * clough_raw_shape(4, p) 01825 + d1xd2n * clough_raw_shape(10, p) 01826 + d1xd3n * clough_raw_shape(11, p) 01827 + 0.5 * N01x * d3nd3n * clough_raw_shape(11, p) 01828 + 0.5 * N02x * d2nd2n * clough_raw_shape(10, p); 01829 case 2: 01830 return d1yd1y * clough_raw_shape(4, p) 01831 + d1yd1x * clough_raw_shape(3, p) 01832 + d1yd2n * clough_raw_shape(10, p) 01833 + d1yd3n * clough_raw_shape(11, p) 01834 + 0.5 * N01y * d3nd3n * clough_raw_shape(11, p) 01835 + 0.5 * N02y * d2nd2n * clough_raw_shape(10, p); 01836 case 4: 01837 return d2xd2x * clough_raw_shape(5, p) 01838 + d2xd2y * clough_raw_shape(6, p) 01839 + d2xd3n * clough_raw_shape(11, p) 01840 + d2xd1n * clough_raw_shape(9, p) 01841 + 0.5 * N10x * d3nd3n * clough_raw_shape(11, p) 01842 + 0.5 * N12x * d1nd1n * clough_raw_shape(9, p); 01843 case 5: 01844 return d2yd2y * clough_raw_shape(6, p) 01845 + d2yd2x * clough_raw_shape(5, p) 01846 + d2yd3n * clough_raw_shape(11, p) 01847 + d2yd1n * clough_raw_shape(9, p) 01848 + 0.5 * N10y * d3nd3n * clough_raw_shape(11, p) 01849 + 0.5 * N12y * d1nd1n * clough_raw_shape(9, p); 01850 case 7: 01851 return d3xd3x * clough_raw_shape(7, p) 01852 + d3xd3y * clough_raw_shape(8, p) 01853 + d3xd1n * clough_raw_shape(9, p) 01854 + d3xd2n * clough_raw_shape(10, p) 01855 + 0.5 * N20x * d2nd2n * clough_raw_shape(10, p) 01856 + 0.5 * N21x * d1nd1n * clough_raw_shape(9, p); 01857 case 8: 01858 return d3yd3y * clough_raw_shape(8, p) 01859 + d3yd3x * clough_raw_shape(7, p) 01860 + d3yd1n * clough_raw_shape(9, p) 01861 + d3yd2n * clough_raw_shape(10, p) 01862 + 0.5 * N20y * d2nd2n * clough_raw_shape(10, p) 01863 + 0.5 * N21y * d1nd1n * clough_raw_shape(9, p); 01864 default: 01865 libmesh_error_msg("Invalid shape function index i = " << i); 01866 } 01867 } 01868 default: 01869 libmesh_error_msg("ERROR: Unsupported element type = " << type); 01870 } 01871 } 01872 // 3rd-order Clough-Tocher element 01873 case THIRD: 01874 { 01875 switch (type) 01876 { 01877 // C1 functions on the Clough-Tocher triangle. 01878 case TRI6: 01879 { 01880 libmesh_assert_less (i, 12); 01881 01882 // FIXME: it would be nice to calculate (and cache) 01883 // clough_raw_shape(j,p) only once per triangle, not 1-7 01884 // times 01885 switch (i) 01886 { 01887 // Note: these DoF numbers are "scrambled" because my 01888 // initial numbering conventions didn't match libMesh 01889 case 0: 01890 return clough_raw_shape(0, p) 01891 + d1d2n * clough_raw_shape(10, p) 01892 + d1d3n * clough_raw_shape(11, p); 01893 case 3: 01894 return clough_raw_shape(1, p) 01895 + d2d3n * clough_raw_shape(11, p) 01896 + d2d1n * clough_raw_shape(9, p); 01897 case 6: 01898 return clough_raw_shape(2, p) 01899 + d3d1n * clough_raw_shape(9, p) 01900 + d3d2n * clough_raw_shape(10, p); 01901 case 1: 01902 return d1xd1x * clough_raw_shape(3, p) 01903 + d1xd1y * clough_raw_shape(4, p) 01904 + d1xd2n * clough_raw_shape(10, p) 01905 + d1xd3n * clough_raw_shape(11, p); 01906 case 2: 01907 return d1yd1y * clough_raw_shape(4, p) 01908 + d1yd1x * clough_raw_shape(3, p) 01909 + d1yd2n * clough_raw_shape(10, p) 01910 + d1yd3n * clough_raw_shape(11, p); 01911 case 4: 01912 return d2xd2x * clough_raw_shape(5, p) 01913 + d2xd2y * clough_raw_shape(6, p) 01914 + d2xd3n * clough_raw_shape(11, p) 01915 + d2xd1n * clough_raw_shape(9, p); 01916 case 5: 01917 return d2yd2y * clough_raw_shape(6, p) 01918 + d2yd2x * clough_raw_shape(5, p) 01919 + d2yd3n * clough_raw_shape(11, p) 01920 + d2yd1n * clough_raw_shape(9, p); 01921 case 7: 01922 return d3xd3x * clough_raw_shape(7, p) 01923 + d3xd3y * clough_raw_shape(8, p) 01924 + d3xd1n * clough_raw_shape(9, p) 01925 + d3xd2n * clough_raw_shape(10, p); 01926 case 8: 01927 return d3yd3y * clough_raw_shape(8, p) 01928 + d3yd3x * clough_raw_shape(7, p) 01929 + d3yd1n * clough_raw_shape(9, p) 01930 + d3yd2n * clough_raw_shape(10, p); 01931 case 10: 01932 return d1nd1n * clough_raw_shape(9, p); 01933 case 11: 01934 return d2nd2n * clough_raw_shape(10, p); 01935 case 9: 01936 return d3nd3n * clough_raw_shape(11, p); 01937 01938 default: 01939 libmesh_error_msg("Invalid shape function index i = " << i); 01940 } 01941 } 01942 default: 01943 libmesh_error_msg("ERROR: Unsupported element type = " << type); 01944 } 01945 } 01946 // by default throw an error 01947 default: 01948 libmesh_error_msg("ERROR: Unsupported polynomial order = " << order); 01949 } 01950 01951 libmesh_error_msg("We'll never get here!"); 01952 return 0.; 01953 } 01954 01955 01956 01957 template <> 01958 Real FE<2,CLOUGH>::shape_deriv(const ElemType, 01959 const Order, 01960 const unsigned int, 01961 const unsigned int, 01962 const Point&) 01963 { 01964 libmesh_error_msg("Clough-Tocher elements require the real element \nto construct gradient-based degrees of freedom."); 01965 return 0.; 01966 } 01967 01968 01969 01970 template <> 01971 Real FE<2,CLOUGH>::shape_deriv(const Elem* elem, 01972 const Order order, 01973 const unsigned int i, 01974 const unsigned int j, 01975 const Point& p) 01976 { 01977 libmesh_assert(elem); 01978 01979 clough_compute_coefs(elem); 01980 01981 const ElemType type = elem->type(); 01982 01983 const Order totalorder = static_cast<Order>(order + elem->p_level()); 01984 01985 switch (totalorder) 01986 { 01987 // 2nd-order restricted Clough-Tocher element 01988 case SECOND: 01989 { 01990 // There may be a bug in the 2nd order case; the 3rd order 01991 // Clough-Tocher elements are pretty uniformly better anyways 01992 // so use those instead. 01993 libmesh_experimental(); 01994 switch (type) 01995 { 01996 // C1 functions on the Clough-Tocher triangle. 01997 case TRI6: 01998 { 01999 libmesh_assert_less (i, 9); 02000 // FIXME: it would be nice to calculate (and cache) 02001 // clough_raw_shape(j,p) only once per triangle, not 1-7 02002 // times 02003 switch (i) 02004 { 02005 // Note: these DoF numbers are "scrambled" because my 02006 // initial numbering conventions didn't match libMesh 02007 case 0: 02008 return clough_raw_shape_deriv(0, j, p) 02009 + d1d2n * clough_raw_shape_deriv(10, j, p) 02010 + d1d3n * clough_raw_shape_deriv(11, j, p); 02011 case 3: 02012 return clough_raw_shape_deriv(1, j, p) 02013 + d2d3n * clough_raw_shape_deriv(11, j, p) 02014 + d2d1n * clough_raw_shape_deriv(9, j, p); 02015 case 6: 02016 return clough_raw_shape_deriv(2, j, p) 02017 + d3d1n * clough_raw_shape_deriv(9, j, p) 02018 + d3d2n * clough_raw_shape_deriv(10, j, p); 02019 case 1: 02020 return d1xd1x * clough_raw_shape_deriv(3, j, p) 02021 + d1xd1y * clough_raw_shape_deriv(4, j, p) 02022 + d1xd2n * clough_raw_shape_deriv(10, j, p) 02023 + d1xd3n * clough_raw_shape_deriv(11, j, p) 02024 + 0.5 * N01x * d3nd3n * clough_raw_shape_deriv(11, j, p) 02025 + 0.5 * N02x * d2nd2n * clough_raw_shape_deriv(10, j, p); 02026 case 2: 02027 return d1yd1y * clough_raw_shape_deriv(4, j, p) 02028 + d1yd1x * clough_raw_shape_deriv(3, j, p) 02029 + d1yd2n * clough_raw_shape_deriv(10, j, p) 02030 + d1yd3n * clough_raw_shape_deriv(11, j, p) 02031 + 0.5 * N01y * d3nd3n * clough_raw_shape_deriv(11, j, p) 02032 + 0.5 * N02y * d2nd2n * clough_raw_shape_deriv(10, j, p); 02033 case 4: 02034 return d2xd2x * clough_raw_shape_deriv(5, j, p) 02035 + d2xd2y * clough_raw_shape_deriv(6, j, p) 02036 + d2xd3n * clough_raw_shape_deriv(11, j, p) 02037 + d2xd1n * clough_raw_shape_deriv(9, j, p) 02038 + 0.5 * N10x * d3nd3n * clough_raw_shape_deriv(11, j, p) 02039 + 0.5 * N12x * d1nd1n * clough_raw_shape_deriv(9, j, p); 02040 case 5: 02041 return d2yd2y * clough_raw_shape_deriv(6, j, p) 02042 + d2yd2x * clough_raw_shape_deriv(5, j, p) 02043 + d2yd3n * clough_raw_shape_deriv(11, j, p) 02044 + d2yd1n * clough_raw_shape_deriv(9, j, p) 02045 + 0.5 * N10y * d3nd3n * clough_raw_shape_deriv(11, j, p) 02046 + 0.5 * N12y * d1nd1n * clough_raw_shape_deriv(9, j, p); 02047 case 7: 02048 return d3xd3x * clough_raw_shape_deriv(7, j, p) 02049 + d3xd3y * clough_raw_shape_deriv(8, j, p) 02050 + d3xd1n * clough_raw_shape_deriv(9, j, p) 02051 + d3xd2n * clough_raw_shape_deriv(10, j, p) 02052 + 0.5 * N20x * d2nd2n * clough_raw_shape_deriv(10, j, p) 02053 + 0.5 * N21x * d1nd1n * clough_raw_shape_deriv(9, j, p); 02054 case 8: 02055 return d3yd3y * clough_raw_shape_deriv(8, j, p) 02056 + d3yd3x * clough_raw_shape_deriv(7, j, p) 02057 + d3yd1n * clough_raw_shape_deriv(9, j, p) 02058 + d3yd2n * clough_raw_shape_deriv(10, j, p) 02059 + 0.5 * N20y * d2nd2n * clough_raw_shape_deriv(10, j, p) 02060 + 0.5 * N21y * d1nd1n * clough_raw_shape_deriv(9, j, p); 02061 default: 02062 libmesh_error_msg("Invalid shape function index i = " << i); 02063 } 02064 } 02065 default: 02066 libmesh_error_msg("ERROR: Unsupported element type = " << type); 02067 } 02068 } 02069 // 3rd-order Clough-Tocher element 02070 case THIRD: 02071 { 02072 switch (type) 02073 { 02074 // C1 functions on the Clough-Tocher triangle. 02075 case TRI6: 02076 { 02077 libmesh_assert_less (i, 12); 02078 02079 // FIXME: it would be nice to calculate (and cache) 02080 // clough_raw_shape(j,p) only once per triangle, not 1-7 02081 // times 02082 switch (i) 02083 { 02084 // Note: these DoF numbers are "scrambled" because my 02085 // initial numbering conventions didn't match libMesh 02086 case 0: 02087 return clough_raw_shape_deriv(0, j, p) 02088 + d1d2n * clough_raw_shape_deriv(10, j, p) 02089 + d1d3n * clough_raw_shape_deriv(11, j, p); 02090 case 3: 02091 return clough_raw_shape_deriv(1, j, p) 02092 + d2d3n * clough_raw_shape_deriv(11, j, p) 02093 + d2d1n * clough_raw_shape_deriv(9, j, p); 02094 case 6: 02095 return clough_raw_shape_deriv(2, j, p) 02096 + d3d1n * clough_raw_shape_deriv(9, j, p) 02097 + d3d2n * clough_raw_shape_deriv(10, j, p); 02098 case 1: 02099 return d1xd1x * clough_raw_shape_deriv(3, j, p) 02100 + d1xd1y * clough_raw_shape_deriv(4, j, p) 02101 + d1xd2n * clough_raw_shape_deriv(10, j, p) 02102 + d1xd3n * clough_raw_shape_deriv(11, j, p); 02103 case 2: 02104 return d1yd1y * clough_raw_shape_deriv(4, j, p) 02105 + d1yd1x * clough_raw_shape_deriv(3, j, p) 02106 + d1yd2n * clough_raw_shape_deriv(10, j, p) 02107 + d1yd3n * clough_raw_shape_deriv(11, j, p); 02108 case 4: 02109 return d2xd2x * clough_raw_shape_deriv(5, j, p) 02110 + d2xd2y * clough_raw_shape_deriv(6, j, p) 02111 + d2xd3n * clough_raw_shape_deriv(11, j, p) 02112 + d2xd1n * clough_raw_shape_deriv(9, j, p); 02113 case 5: 02114 return d2yd2y * clough_raw_shape_deriv(6, j, p) 02115 + d2yd2x * clough_raw_shape_deriv(5, j, p) 02116 + d2yd3n * clough_raw_shape_deriv(11, j, p) 02117 + d2yd1n * clough_raw_shape_deriv(9, j, p); 02118 case 7: 02119 return d3xd3x * clough_raw_shape_deriv(7, j, p) 02120 + d3xd3y * clough_raw_shape_deriv(8, j, p) 02121 + d3xd1n * clough_raw_shape_deriv(9, j, p) 02122 + d3xd2n * clough_raw_shape_deriv(10, j, p); 02123 case 8: 02124 return d3yd3y * clough_raw_shape_deriv(8, j, p) 02125 + d3yd3x * clough_raw_shape_deriv(7, j, p) 02126 + d3yd1n * clough_raw_shape_deriv(9, j, p) 02127 + d3yd2n * clough_raw_shape_deriv(10, j, p); 02128 case 10: 02129 return d1nd1n * clough_raw_shape_deriv(9, j, p); 02130 case 11: 02131 return d2nd2n * clough_raw_shape_deriv(10, j, p); 02132 case 9: 02133 return d3nd3n * clough_raw_shape_deriv(11, j, p); 02134 02135 default: 02136 libmesh_error_msg("Invalid shape function index i = " << i); 02137 } 02138 } 02139 default: 02140 libmesh_error_msg("ERROR: Unsupported element type = " << type); 02141 } 02142 } 02143 // by default throw an error 02144 default: 02145 libmesh_error_msg("ERROR: Unsupported polynomial order = " << order); 02146 } 02147 02148 libmesh_error_msg("We'll never get here!"); 02149 return 0.; 02150 } 02151 02152 02153 02154 template <> 02155 Real FE<2,CLOUGH>::shape_second_deriv(const Elem* elem, 02156 const Order order, 02157 const unsigned int i, 02158 const unsigned int j, 02159 const Point& p) 02160 { 02161 libmesh_assert(elem); 02162 02163 clough_compute_coefs(elem); 02164 02165 const ElemType type = elem->type(); 02166 02167 const Order totalorder = static_cast<Order>(order + elem->p_level()); 02168 02169 switch (totalorder) 02170 { 02171 // 2nd-order restricted Clough-Tocher element 02172 case SECOND: 02173 { 02174 switch (type) 02175 { 02176 // C1 functions on the Clough-Tocher triangle. 02177 case TRI6: 02178 { 02179 libmesh_assert_less (i, 9); 02180 // FIXME: it would be nice to calculate (and cache) 02181 // clough_raw_shape(j,p) only once per triangle, not 1-7 02182 // times 02183 switch (i) 02184 { 02185 // Note: these DoF numbers are "scrambled" because my 02186 // initial numbering conventions didn't match libMesh 02187 case 0: 02188 return clough_raw_shape_second_deriv(0, j, p) 02189 + d1d2n * clough_raw_shape_second_deriv(10, j, p) 02190 + d1d3n * clough_raw_shape_second_deriv(11, j, p); 02191 case 3: 02192 return clough_raw_shape_second_deriv(1, j, p) 02193 + d2d3n * clough_raw_shape_second_deriv(11, j, p) 02194 + d2d1n * clough_raw_shape_second_deriv(9, j, p); 02195 case 6: 02196 return clough_raw_shape_second_deriv(2, j, p) 02197 + d3d1n * clough_raw_shape_second_deriv(9, j, p) 02198 + d3d2n * clough_raw_shape_second_deriv(10, j, p); 02199 case 1: 02200 return d1xd1x * clough_raw_shape_second_deriv(3, j, p) 02201 + d1xd1y * clough_raw_shape_second_deriv(4, j, p) 02202 + d1xd2n * clough_raw_shape_second_deriv(10, j, p) 02203 + d1xd3n * clough_raw_shape_second_deriv(11, j, p) 02204 + 0.5 * N01x * d3nd3n * clough_raw_shape_second_deriv(11, j, p) 02205 + 0.5 * N02x * d2nd2n * clough_raw_shape_second_deriv(10, j, p); 02206 case 2: 02207 return d1yd1y * clough_raw_shape_second_deriv(4, j, p) 02208 + d1yd1x * clough_raw_shape_second_deriv(3, j, p) 02209 + d1yd2n * clough_raw_shape_second_deriv(10, j, p) 02210 + d1yd3n * clough_raw_shape_second_deriv(11, j, p) 02211 + 0.5 * N01y * d3nd3n * clough_raw_shape_second_deriv(11, j, p) 02212 + 0.5 * N02y * d2nd2n * clough_raw_shape_second_deriv(10, j, p); 02213 case 4: 02214 return d2xd2x * clough_raw_shape_second_deriv(5, j, p) 02215 + d2xd2y * clough_raw_shape_second_deriv(6, j, p) 02216 + d2xd3n * clough_raw_shape_second_deriv(11, j, p) 02217 + d2xd1n * clough_raw_shape_second_deriv(9, j, p) 02218 + 0.5 * N10x * d3nd3n * clough_raw_shape_second_deriv(11, j, p) 02219 + 0.5 * N12x * d1nd1n * clough_raw_shape_second_deriv(9, j, p); 02220 case 5: 02221 return d2yd2y * clough_raw_shape_second_deriv(6, j, p) 02222 + d2yd2x * clough_raw_shape_second_deriv(5, j, p) 02223 + d2yd3n * clough_raw_shape_second_deriv(11, j, p) 02224 + d2yd1n * clough_raw_shape_second_deriv(9, j, p) 02225 + 0.5 * N10y * d3nd3n * clough_raw_shape_second_deriv(11, j, p) 02226 + 0.5 * N12y * d1nd1n * clough_raw_shape_second_deriv(9, j, p); 02227 case 7: 02228 return d3xd3x * clough_raw_shape_second_deriv(7, j, p) 02229 + d3xd3y * clough_raw_shape_second_deriv(8, j, p) 02230 + d3xd1n * clough_raw_shape_second_deriv(9, j, p) 02231 + d3xd2n * clough_raw_shape_second_deriv(10, j, p) 02232 + 0.5 * N20x * d2nd2n * clough_raw_shape_second_deriv(10, j, p) 02233 + 0.5 * N21x * d1nd1n * clough_raw_shape_second_deriv(9, j, p); 02234 case 8: 02235 return d3yd3y * clough_raw_shape_second_deriv(8, j, p) 02236 + d3yd3x * clough_raw_shape_second_deriv(7, j, p) 02237 + d3yd1n * clough_raw_shape_second_deriv(9, j, p) 02238 + d3yd2n * clough_raw_shape_second_deriv(10, j, p) 02239 + 0.5 * N20y * d2nd2n * clough_raw_shape_second_deriv(10, j, p) 02240 + 0.5 * N21y * d1nd1n * clough_raw_shape_second_deriv(9, j, p); 02241 default: 02242 libmesh_error_msg("Invalid shape function index i = " << i); 02243 } 02244 } 02245 default: 02246 libmesh_error_msg("ERROR: Unsupported element type = " << type); 02247 } 02248 } 02249 // 3rd-order Clough-Tocher element 02250 case THIRD: 02251 { 02252 switch (type) 02253 { 02254 // C1 functions on the Clough-Tocher triangle. 02255 case TRI6: 02256 { 02257 libmesh_assert_less (i, 12); 02258 02259 // FIXME: it would be nice to calculate (and cache) 02260 // clough_raw_shape(j,p) only once per triangle, not 1-7 02261 // times 02262 switch (i) 02263 { 02264 // Note: these DoF numbers are "scrambled" because my 02265 // initial numbering conventions didn't match libMesh 02266 case 0: 02267 return clough_raw_shape_second_deriv(0, j, p) 02268 + d1d2n * clough_raw_shape_second_deriv(10, j, p) 02269 + d1d3n * clough_raw_shape_second_deriv(11, j, p); 02270 case 3: 02271 return clough_raw_shape_second_deriv(1, j, p) 02272 + d2d3n * clough_raw_shape_second_deriv(11, j, p) 02273 + d2d1n * clough_raw_shape_second_deriv(9, j, p); 02274 case 6: 02275 return clough_raw_shape_second_deriv(2, j, p) 02276 + d3d1n * clough_raw_shape_second_deriv(9, j, p) 02277 + d3d2n * clough_raw_shape_second_deriv(10, j, p); 02278 case 1: 02279 return d1xd1x * clough_raw_shape_second_deriv(3, j, p) 02280 + d1xd1y * clough_raw_shape_second_deriv(4, j, p) 02281 + d1xd2n * clough_raw_shape_second_deriv(10, j, p) 02282 + d1xd3n * clough_raw_shape_second_deriv(11, j, p); 02283 case 2: 02284 return d1yd1y * clough_raw_shape_second_deriv(4, j, p) 02285 + d1yd1x * clough_raw_shape_second_deriv(3, j, p) 02286 + d1yd2n * clough_raw_shape_second_deriv(10, j, p) 02287 + d1yd3n * clough_raw_shape_second_deriv(11, j, p); 02288 case 4: 02289 return d2xd2x * clough_raw_shape_second_deriv(5, j, p) 02290 + d2xd2y * clough_raw_shape_second_deriv(6, j, p) 02291 + d2xd3n * clough_raw_shape_second_deriv(11, j, p) 02292 + d2xd1n * clough_raw_shape_second_deriv(9, j, p); 02293 case 5: 02294 return d2yd2y * clough_raw_shape_second_deriv(6, j, p) 02295 + d2yd2x * clough_raw_shape_second_deriv(5, j, p) 02296 + d2yd3n * clough_raw_shape_second_deriv(11, j, p) 02297 + d2yd1n * clough_raw_shape_second_deriv(9, j, p); 02298 case 7: 02299 return d3xd3x * clough_raw_shape_second_deriv(7, j, p) 02300 + d3xd3y * clough_raw_shape_second_deriv(8, j, p) 02301 + d3xd1n * clough_raw_shape_second_deriv(9, j, p) 02302 + d3xd2n * clough_raw_shape_second_deriv(10, j, p); 02303 case 8: 02304 return d3yd3y * clough_raw_shape_second_deriv(8, j, p) 02305 + d3yd3x * clough_raw_shape_second_deriv(7, j, p) 02306 + d3yd1n * clough_raw_shape_second_deriv(9, j, p) 02307 + d3yd2n * clough_raw_shape_second_deriv(10, j, p); 02308 case 10: 02309 return d1nd1n * clough_raw_shape_second_deriv(9, j, p); 02310 case 11: 02311 return d2nd2n * clough_raw_shape_second_deriv(10, j, p); 02312 case 9: 02313 return d3nd3n * clough_raw_shape_second_deriv(11, j, p); 02314 02315 default: 02316 libmesh_error_msg("Invalid shape function index i = " << i); 02317 } 02318 } 02319 default: 02320 libmesh_error_msg("ERROR: Unsupported element type = " << type); 02321 } 02322 } 02323 // by default throw an error 02324 default: 02325 libmesh_error_msg("ERROR: Unsupported polynomial order = " << order); 02326 } 02327 02328 libmesh_error_msg("We'll never get here!"); 02329 return 0.; 02330 } 02331 02332 } // namespace libMesh
generated by
