fe_clough_shape_2D.C
Go to the documentation of this file.00001 // The libMesh Finite Element Library. 00002 // Copyright (C) 2002-2012 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 = libMesh::invalid_uint; 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 case 1: 00552 switch (subtriangle_lookup(p)) 00553 { 00554 case 0: 00555 return 6 - 12*xi; 00556 case 1: 00557 return 18 - 27*xi - 21*eta; 00558 case 2: 00559 return 6 - 15*xi + 3*eta; 00560 } 00561 case 2: 00562 switch (subtriangle_lookup(p)) 00563 { 00564 case 0: 00565 return 0; 00566 case 1: 00567 return 12 - 15*xi - 21*eta; 00568 case 2: 00569 return -3*xi + 3*eta; 00570 } 00571 case 3: 00572 switch (subtriangle_lookup(p)) 00573 { 00574 case 0: 00575 return -4 + 6*xi; 00576 case 1: 00577 return -9 + 13*xi + 8*eta; 00578 case 2: 00579 return -1 - 7*xi + 4*eta; 00580 } 00581 case 4: 00582 switch (subtriangle_lookup(p)) 00583 { 00584 case 0: 00585 return 4*eta; 00586 case 1: 00587 return 1 - 2*xi + 3*eta; 00588 case 2: 00589 return -3 + 14*xi - eta; 00590 } 00591 case 5: 00592 switch (subtriangle_lookup(p)) 00593 { 00594 case 0: 00595 return -2 + 6*xi; 00596 case 1: 00597 return -4 + 17./2.*xi + 7./2.*eta; 00598 case 2: 00599 return -2 + 13./2.*xi - 1./2.*eta; 00600 } 00601 case 6: 00602 switch (subtriangle_lookup(p)) 00603 { 00604 case 0: 00605 return 4*eta; 00606 case 1: 00607 return 9 - 23./2.*xi - 23./2.*eta; 00608 case 2: 00609 return -1 + 5./2.*xi + 9./2.*eta; 00610 } 00611 case 7: 00612 switch (subtriangle_lookup(p)) 00613 { 00614 case 0: 00615 return 0; 00616 case 1: 00617 return 7 - 17./2.*xi - 25./2.*eta; 00618 case 2: 00619 return 1 - 13./2.*xi + 7./2.*eta; 00620 } 00621 case 8: 00622 switch (subtriangle_lookup(p)) 00623 { 00624 case 0: 00625 return 0; 00626 case 1: 00627 return -2 + 5./2.*xi + 7./2.*eta; 00628 case 2: 00629 return 1./2.*xi - 1./2.*eta; 00630 } 00631 case 9: 00632 switch (subtriangle_lookup(p)) 00633 { 00634 case 0: 00635 return 0; 00636 case 1: 00637 return std::sqrt(2.) * (8 - 10*xi - 14*eta); 00638 case 2: 00639 return std::sqrt(2.) * (-2*xi + 2*eta); 00640 } 00641 case 10: 00642 switch (subtriangle_lookup(p)) 00643 { 00644 case 0: 00645 return 0; 00646 case 1: 00647 return -4 + 4*xi + 8*eta; 00648 case 2: 00649 return -4 + 20*xi - 8*eta; 00650 } 00651 case 11: 00652 switch (subtriangle_lookup(p)) 00653 { 00654 case 0: 00655 return -8*eta; 00656 case 1: 00657 return -12 + 16*xi + 12*eta; 00658 case 2: 00659 return 4 - 16*xi - 4*eta; 00660 } 00661 } 00662 00663 // second derivative in xi-eta direction 00664 case 1: 00665 switch (basis_num) 00666 { 00667 case 0: 00668 switch (subtriangle_lookup(p)) 00669 { 00670 case 0: 00671 return -6*eta; 00672 case 1: 00673 return -30 +42*xi 00674 + 42*eta; 00675 case 2: 00676 return -6*xi; 00677 } 00678 case 1: 00679 switch (subtriangle_lookup(p)) 00680 { 00681 case 0: 00682 return + 3*eta; 00683 case 1: 00684 return 15 - 21*xi - 21*eta; 00685 case 2: 00686 return 3*xi; 00687 } 00688 case 2: 00689 switch (subtriangle_lookup(p)) 00690 { 00691 case 0: 00692 return 3*eta; 00693 case 1: 00694 return 15 - 21*xi - 21*eta; 00695 case 2: 00696 return 3*xi; 00697 } 00698 case 3: 00699 switch (subtriangle_lookup(p)) 00700 { 00701 case 0: 00702 return -eta; 00703 case 1: 00704 return -4 + 8*xi + 3*eta; 00705 case 2: 00706 return -3 + 4*xi + 4*eta; 00707 } 00708 case 4: 00709 switch (subtriangle_lookup(p)) 00710 { 00711 case 0: 00712 return -3 + 4*xi + 4*eta; 00713 case 1: 00714 return - 4 + 3*xi + 8*eta; 00715 case 2: 00716 return -xi; 00717 } 00718 case 5: 00719 switch (subtriangle_lookup(p)) 00720 { 00721 case 0: 00722 return - 1./2.*eta; 00723 case 1: 00724 return -5./2. + 7./2.*xi + 7./2.*eta; 00725 case 2: 00726 return - 1./2.*xi; 00727 } 00728 case 6: 00729 switch (subtriangle_lookup(p)) 00730 { 00731 case 0: 00732 return -1 + 4*xi + 7./2.*eta; 00733 case 1: 00734 return 19./2. - 23./2.*xi - 25./2.*eta; 00735 case 2: 00736 return 9./2.*xi; 00737 } 00738 case 7: 00739 switch (subtriangle_lookup(p)) 00740 { 00741 case 0: 00742 return 9./2.*eta; 00743 case 1: 00744 return 19./2. - 25./2.*xi - 23./2.*eta; 00745 case 2: 00746 return -1 + 7./2.*xi + 4*eta; 00747 } 00748 case 8: 00749 switch (subtriangle_lookup(p)) 00750 { 00751 case 0: 00752 return -1./2.*eta; 00753 case 1: 00754 return -5./2. + 7./2.*xi + 7./2.*eta; 00755 case 2: 00756 return -1./2.*xi; 00757 } 00758 case 9: 00759 switch (subtriangle_lookup(p)) 00760 { 00761 case 0: 00762 return std::sqrt(2.) * (2*eta); 00763 case 1: 00764 return std::sqrt(2.) * (10 - 14*xi - 14*eta); 00765 case 2: 00766 return std::sqrt(2.) * (2*xi); 00767 } 00768 case 10: 00769 switch (subtriangle_lookup(p)) 00770 { 00771 case 0: 00772 return -4*eta; 00773 case 1: 00774 return - 8 + 8*xi + 12*eta; 00775 case 2: 00776 return 4 - 8*xi - 8*eta; 00777 } 00778 case 11: 00779 switch (subtriangle_lookup(p)) 00780 { 00781 case 0: 00782 return 4 - 8*xi - 8*eta; 00783 case 1: 00784 return -8 + 12*xi + 8*eta; 00785 case 2: 00786 return -4*xi; 00787 } 00788 } 00789 00790 // second derivative in eta-eta direction 00791 case 2: 00792 switch (basis_num) 00793 { 00794 case 0: 00795 switch (subtriangle_lookup(p)) 00796 { 00797 case 0: 00798 return -6 - 6*xi + 18*eta; 00799 case 1: 00800 return -30 + 42*xi + 42*eta; 00801 case 2: 00802 return -6 + 12*eta; 00803 } 00804 case 1: 00805 switch (subtriangle_lookup(p)) 00806 { 00807 case 0: 00808 return 3*xi - 3*eta; 00809 case 1: 00810 return 12 - 21*xi - 15*eta; 00811 case 2: 00812 return 0; 00813 } 00814 case 2: 00815 switch (subtriangle_lookup(p)) 00816 { 00817 case 0: 00818 return 6 + 3*xi - 15*eta; 00819 case 1: 00820 return 18 - 21.*xi - 27*eta; 00821 case 2: 00822 return 6 - 12*eta; 00823 } 00824 case 3: 00825 switch (subtriangle_lookup(p)) 00826 { 00827 case 0: 00828 return -3 - xi + 14*eta; 00829 case 1: 00830 return 1 + 3*xi - 2*eta; 00831 case 2: 00832 return 4*xi; 00833 } 00834 case 4: 00835 switch (subtriangle_lookup(p)) 00836 { 00837 case 0: 00838 return -1 + 4*xi - 7*eta; 00839 case 1: 00840 return -9 + 8*xi + 13*eta; 00841 case 2: 00842 return -4 + 6*eta; 00843 } 00844 case 5: 00845 switch (subtriangle_lookup(p)) 00846 { 00847 case 0: 00848 return - 1./2.*xi + 1./2.*eta; 00849 case 1: 00850 return -2 + 7./2.*xi + 5./2.*eta; 00851 case 2: 00852 return 0; 00853 } 00854 case 6: 00855 switch (subtriangle_lookup(p)) 00856 { 00857 case 0: 00858 return 1 + 7./2.*xi - 13./2.*eta; 00859 case 1: 00860 return 7 - 25./2.*xi - 17./2.*eta; 00861 case 2: 00862 return 0; 00863 } 00864 case 7: 00865 switch (subtriangle_lookup(p)) 00866 { 00867 case 0: 00868 return -1 + 9./2.*xi + 5./2.*eta; 00869 case 1: 00870 return 9 - 23./2.*xi - 23./2.*eta; 00871 case 2: 00872 return 4*xi; 00873 } 00874 case 8: 00875 switch (subtriangle_lookup(p)) 00876 { 00877 case 0: 00878 return -2 - 1./2.*xi + 13./2.*eta; 00879 case 1: 00880 return -4 + 7./2.*xi + 17./2.*eta; 00881 case 2: 00882 return -2 + 6*eta; 00883 } 00884 case 9: 00885 switch (subtriangle_lookup(p)) 00886 { 00887 case 0: 00888 return std::sqrt(2.) * (2*xi - 2*eta); 00889 case 1: 00890 return std::sqrt(2.) * (8 - 14*xi - 10*eta); 00891 case 2: 00892 return 0; 00893 } 00894 case 10: 00895 switch (subtriangle_lookup(p)) 00896 { 00897 case 0: 00898 return 4 - 4*xi - 16*eta; 00899 case 1: 00900 return -12 + 12*xi + 16*eta; 00901 case 2: 00902 return -8*xi; 00903 } 00904 case 11: 00905 switch (subtriangle_lookup(p)) 00906 { 00907 case 0: 00908 return -4 - 8*xi + 20*eta; 00909 case 1: 00910 return -4 + 8*xi + 4*eta; 00911 case 2: 00912 return 0; 00913 } 00914 } 00915 } 00916 00917 libmesh_error(); 00918 return 0.; 00919 } 00920 00921 00922 00923 Real clough_raw_shape_deriv(const unsigned int basis_num, 00924 const unsigned int deriv_type, 00925 const Point& p) 00926 { 00927 Real xi = p(0), eta = p(1); 00928 00929 switch (deriv_type) 00930 { 00931 // derivative in xi direction 00932 case 0: 00933 switch (basis_num) 00934 { 00935 case 0: 00936 switch (subtriangle_lookup(p)) 00937 { 00938 case 0: 00939 return -6*xi + 6*xi*xi 00940 - 3*eta*eta; 00941 case 1: 00942 return 9 - 30*xi + 21*xi*xi 00943 - 30*eta + 42*xi*eta 00944 + 21*eta*eta; 00945 case 2: 00946 return -6*xi + 9*xi*xi 00947 - 6*xi*eta; 00948 } 00949 case 1: 00950 switch (subtriangle_lookup(p)) 00951 { 00952 case 0: 00953 return 6*xi - 6*xi*xi 00954 + 3./2.*eta*eta; 00955 case 1: 00956 return - 9./2. + 18*xi - 27./2.*xi*xi 00957 + 15*eta - 21*xi*eta 00958 - 21./2.*eta*eta; 00959 case 2: 00960 return 6*xi - 15./2.*xi*xi 00961 + 3*xi*eta; 00962 } 00963 case 2: 00964 switch (subtriangle_lookup(p)) 00965 { 00966 case 0: 00967 return 3./2.*eta*eta; 00968 case 1: 00969 return - 9./2. + 12*xi - 15./2.*xi*xi 00970 + 15*eta - 21*xi*eta 00971 - 21./2.*eta*eta; 00972 case 2: 00973 return -3./2.*xi*xi 00974 + 3*xi*eta; 00975 } 00976 case 3: 00977 switch (subtriangle_lookup(p)) 00978 { 00979 case 0: 00980 return 1 - 4*xi + 3*xi*xi 00981 - 1./2.*eta*eta; 00982 case 1: 00983 return 5./2. - 9*xi + 13./2.*xi*xi 00984 - 4*eta + 8*xi*eta 00985 + 3./2.*eta*eta; 00986 case 2: 00987 return 1 - xi - 7./2.*xi*xi 00988 - 3*eta + 4*xi*eta 00989 + 2*eta*eta; 00990 } 00991 case 4: 00992 switch (subtriangle_lookup(p)) 00993 { 00994 case 0: 00995 return - 3*eta + 4*xi*eta 00996 + 2*eta*eta; 00997 case 1: 00998 return xi - xi*xi 00999 - 4*eta + 3*xi*eta 01000 + 4*eta*eta; 01001 case 2: 01002 return -3*xi + 7*xi*xi 01003 - xi*eta; 01004 } 01005 case 5: 01006 switch (subtriangle_lookup(p)) 01007 { 01008 case 0: 01009 return -2*xi + 3*xi*xi 01010 - 1./4.*eta*eta; 01011 case 1: 01012 return 3./4. - 4*xi + 17./4.*xi*xi 01013 - 5./2.*eta + 7./2.*xi*eta 01014 + 7./4.*eta*eta; 01015 case 2: 01016 return -2*xi + 13./4.*xi*xi 01017 - 1./2.*xi*eta; 01018 } 01019 case 6: 01020 switch (subtriangle_lookup(p)) 01021 { 01022 case 0: 01023 return -eta + 4*xi*eta 01024 + 7./4.*eta*eta; 01025 case 1: 01026 return -13./4. + 9*xi - 23./4.*xi*xi 01027 + 19./2.*eta - 23./2.*xi*eta 01028 - 25./4.*eta*eta; 01029 case 2: 01030 return -xi + 5./4.*xi*xi 01031 + 9./2.*xi*eta; 01032 } 01033 case 7: 01034 switch (subtriangle_lookup(p)) 01035 { 01036 case 0: 01037 return 9./4.*eta*eta; 01038 case 1: 01039 return - 11./4. + 7*xi - 17./4.*xi*xi 01040 + 19./2.*eta - 25./2.*xi*eta 01041 - 23./4.*eta*eta; 01042 case 2: 01043 return xi - 13./4.*xi*xi 01044 - eta + 7./2.*xi*eta + 2*eta*eta; 01045 } 01046 case 8: 01047 switch (subtriangle_lookup(p)) 01048 { 01049 case 0: 01050 return - 1./4.*eta*eta; 01051 case 1: 01052 return 3./4. - 2*xi + 5./4.*xi*xi 01053 - 5./2.*eta + 7./2.*xi*eta 01054 + 7./4.*eta*eta; 01055 case 2: 01056 return 1./4.*xi*xi 01057 - 1./2.*xi*eta; 01058 } 01059 case 9: 01060 switch (subtriangle_lookup(p)) 01061 { 01062 case 0: 01063 return std::sqrt(2.) * eta*eta; 01064 case 1: 01065 return std::sqrt(2.) * (-3 + 8*xi - 5*xi*xi 01066 + 10*eta - 14*xi*eta 01067 - 7*eta*eta); 01068 case 2: 01069 return std::sqrt(2.) * (-xi*xi + 2*xi*eta); 01070 } 01071 case 10: 01072 switch (subtriangle_lookup(p)) 01073 { 01074 case 0: 01075 return -2*eta*eta; 01076 case 1: 01077 return 2 - 4*xi + 2*xi*xi 01078 - 8*eta + 8*xi*eta 01079 + 6*eta*eta; 01080 case 2: 01081 return -4*xi + 10*xi*xi 01082 + 4*eta - 8*xi*eta 01083 - 4*eta*eta; 01084 } 01085 case 11: 01086 switch (subtriangle_lookup(p)) 01087 { 01088 case 0: 01089 return 4*eta - 8*xi*eta 01090 - 4*eta*eta; 01091 case 1: 01092 return 4 - 12*xi + 8*xi*xi 01093 - 8*eta + 12*xi*eta 01094 + 4*eta*eta; 01095 case 2: 01096 return 4*xi - 8*xi*xi 01097 - 4*xi*eta; 01098 } 01099 } 01100 01101 // derivative in eta direction 01102 case 1: 01103 switch (basis_num) 01104 { 01105 case 0: 01106 switch (subtriangle_lookup(p)) 01107 { 01108 case 0: 01109 return - 6*eta - 6*xi*eta + 9*eta*eta; 01110 case 1: 01111 return 9 - 30*xi + 21*xi*xi 01112 - 30*eta + 42*xi*eta + 21*eta*eta; 01113 case 2: 01114 return - 3*xi*xi 01115 - 6*eta + 6*eta*eta; 01116 } 01117 case 1: 01118 switch (subtriangle_lookup(p)) 01119 { 01120 case 0: 01121 return + 3*xi*eta 01122 - 3./2.*eta*eta; 01123 case 1: 01124 return - 9./2. + 15*xi - 21./2.*xi*xi 01125 + 12*eta - 21*xi*eta - 15./2.*eta*eta; 01126 case 2: 01127 return + 3./2.*xi*xi; 01128 } 01129 case 2: 01130 switch (subtriangle_lookup(p)) 01131 { 01132 case 0: 01133 return 6*eta + 3*xi*eta - 15./2.*eta*eta; 01134 case 1: 01135 return - 9./2. + 15*xi - 21./2.*xi*xi 01136 + 18*eta - 21.*xi*eta - 27./2.*eta*eta; 01137 case 2: 01138 return 3./2.*xi*xi 01139 + 6*eta - 6*eta*eta; 01140 } 01141 case 3: 01142 switch (subtriangle_lookup(p)) 01143 { 01144 case 0: 01145 return - 3*eta - xi*eta + 7*eta*eta; 01146 case 1: 01147 return - 4*xi + 4*xi*xi 01148 + eta + 3*xi*eta - eta*eta; 01149 case 2: 01150 return - 3*xi + 2*xi*xi 01151 + 4*xi*eta; 01152 } 01153 case 4: 01154 switch (subtriangle_lookup(p)) 01155 { 01156 case 0: 01157 return 1 - 3*xi + 2*xi*xi 01158 - eta + 4*xi*eta - 7./2.*eta*eta; 01159 case 1: 01160 return 5./2. - 4*xi + 3./2.*xi*xi 01161 - 9.*eta + 8*xi*eta + 13./2.*eta*eta; 01162 case 2: 01163 return 1 - 1./2.*xi*xi - 4*eta + 3*eta*eta; 01164 } 01165 case 5: 01166 switch (subtriangle_lookup(p)) 01167 { 01168 case 0: 01169 return - 1./2.*xi*eta + 1./4.*eta*eta; 01170 case 1: 01171 return 3./4. - 5./2.*xi + 7./4.*xi*xi 01172 - 2*eta + 7./2.*xi*eta + 5./4.*eta*eta; 01173 case 2: 01174 return - 1./4.*xi*xi; 01175 } 01176 case 6: 01177 switch (subtriangle_lookup(p)) 01178 { 01179 case 0: 01180 return -xi + 2*xi*xi 01181 + eta + 7./2.*xi*eta - 13./4.*eta*eta; 01182 case 1: 01183 return - 11./4. + 19./2.*xi - 23./4.*xi*xi 01184 + 7*eta - 25./2.*xi*eta - 17./4.*eta*eta; 01185 case 2: 01186 return 9./4.*xi*xi; 01187 } 01188 case 7: 01189 switch (subtriangle_lookup(p)) 01190 { 01191 case 0: 01192 return -eta + 9./2.*xi*eta + 5./4.*eta*eta; 01193 case 1: 01194 return - 13./4. + 19./2.*xi - 25./4.*xi*xi 01195 + 9*eta - 23./2.*xi*eta - 23./4.*eta*eta; 01196 case 2: 01197 return - xi + 7./4.*xi*xi + 4*xi*eta; 01198 } 01199 case 8: 01200 switch (subtriangle_lookup(p)) 01201 { 01202 case 0: 01203 return -2*eta - 1./2.*xi*eta + 13./4.*eta*eta; 01204 case 1: 01205 return 3./4. - 5./2.*xi + 7./4.*xi*xi 01206 - 4*eta + 7./2.*xi*eta + 17./4.*eta*eta; 01207 case 2: 01208 return - 1./4.*xi*xi 01209 - 2*eta + 3*eta*eta; 01210 } 01211 case 9: 01212 switch (subtriangle_lookup(p)) 01213 { 01214 case 0: 01215 return std::sqrt(2.) * (2*xi*eta - eta*eta); 01216 case 1: 01217 return std::sqrt(2.) * (- 3 + 10*xi - 7*xi*xi 01218 + 8*eta - 14*xi*eta - 5*eta*eta); 01219 case 2: 01220 return std::sqrt(2.) * (xi*xi); 01221 } 01222 case 10: 01223 switch (subtriangle_lookup(p)) 01224 { 01225 case 0: 01226 return 4*eta - 4*xi*eta - 8*eta*eta; 01227 case 1: 01228 return 4 - 8*xi + 4*xi*xi 01229 - 12*eta + 12*xi*eta + 8*eta*eta; 01230 case 2: 01231 return 4*xi - 4*xi*xi 01232 - 8*xi*eta; 01233 } 01234 case 11: 01235 switch (subtriangle_lookup(p)) 01236 { 01237 case 0: 01238 return 4*xi - 4*xi*xi 01239 - 4*eta - 8*xi*eta + 10.*eta*eta; 01240 case 1: 01241 return 2 - 8*xi + 6*xi*xi 01242 - 4*eta + 8*xi*eta + 2*eta*eta; 01243 case 2: 01244 return - 2*xi*xi; 01245 } 01246 } 01247 } 01248 01249 libmesh_error(); 01250 return 0.; 01251 } 01252 01253 Real clough_raw_shape(const unsigned int basis_num, 01254 const Point& p) 01255 { 01256 Real xi = p(0), eta = p(1); 01257 01258 switch (basis_num) 01259 { 01260 case 0: 01261 switch (subtriangle_lookup(p)) 01262 { 01263 case 0: 01264 return 1 - 3*xi*xi + 2*xi*xi*xi 01265 - 3*eta*eta - 3*xi*eta*eta + 3*eta*eta*eta; 01266 case 1: 01267 return -1 + 9*xi - 15*xi*xi + 7*xi*xi*xi 01268 + 9*eta - 30*xi*eta +21*xi*xi*eta 01269 - 15*eta*eta + 21*xi*eta*eta + 7*eta*eta*eta; 01270 case 2: 01271 return 1 - 3*xi*xi + 3*xi*xi*xi 01272 - 3*xi*xi*eta 01273 - 3*eta*eta + 2*eta*eta*eta; 01274 } 01275 case 1: 01276 switch (subtriangle_lookup(p)) 01277 { 01278 case 0: 01279 return 3*xi*xi - 2*xi*xi*xi 01280 + 3./2.*xi*eta*eta 01281 - 1./2.*eta*eta*eta; 01282 case 1: 01283 return 1 - 9./2.*xi + 9*xi*xi - 9./2.*xi*xi*xi 01284 - 9./2.*eta + 15*xi*eta - 21./2.*xi*xi*eta 01285 + 6*eta*eta - 21./2.*xi*eta*eta - 5./2.*eta*eta*eta; 01286 case 2: 01287 return 3*xi*xi - 5./2.*xi*xi*xi 01288 + 3./2.*xi*xi*eta; 01289 } 01290 case 2: 01291 switch (subtriangle_lookup(p)) 01292 { 01293 case 0: 01294 return 3*eta*eta + 3./2.*xi*eta*eta - 5./2.*eta*eta*eta; 01295 case 1: 01296 return 1 - 9./2.*xi + 6*xi*xi - 5./2.*xi*xi*xi 01297 - 9./2.*eta + 15*xi*eta - 21./2.*xi*xi*eta 01298 + 9*eta*eta - 21./2.*xi*eta*eta - 9./2.*eta*eta*eta; 01299 case 2: 01300 return -1./2.*xi*xi*xi 01301 + 3./2.*xi*xi*eta 01302 + 3*eta*eta - 2*eta*eta*eta; 01303 } 01304 case 3: 01305 switch (subtriangle_lookup(p)) 01306 { 01307 case 0: 01308 return xi - 2*xi*xi + xi*xi*xi 01309 - 3./2.*eta*eta - 1./2.*xi*eta*eta + 7./3.*eta*eta*eta; 01310 case 1: 01311 return -1./6. + 5./2.*xi - 9./2.*xi*xi + 13./6.*xi*xi*xi 01312 - 4*xi*eta + 4*xi*xi*eta 01313 + 1./2.*eta*eta + 3./2.*xi*eta*eta - 1./3.*eta*eta*eta; 01314 case 2: 01315 return xi - 1./2.*xi*xi - 7./6.*xi*xi*xi 01316 - 3*xi*eta + 2*xi*xi*eta 01317 + 2*xi*eta*eta; 01318 } 01319 case 4: 01320 switch (subtriangle_lookup(p)) 01321 { 01322 case 0: 01323 return eta - 3*xi*eta + 2*xi*xi*eta 01324 - 1./2.*eta*eta + 2*xi*eta*eta - 7./6.*eta*eta*eta; 01325 case 1: 01326 return -1./6. + 1./2.*xi*xi - 1./3.*xi*xi*xi 01327 + 5./2.*eta - 4*xi*eta + 3./2.*xi*xi*eta 01328 - 9./2.*eta*eta + 4*xi*eta*eta + 13./6.*eta*eta*eta; 01329 case 2: 01330 return -3./2.*xi*xi + 7./3.*xi*xi*xi 01331 + eta - 1./2.*xi*xi*eta - 2*eta*eta + eta*eta*eta; 01332 } 01333 case 5: 01334 switch (subtriangle_lookup(p)) 01335 { 01336 case 0: 01337 return -xi*xi + xi*xi*xi 01338 - 1./4.*xi*eta*eta + 1./12.*eta*eta*eta; 01339 case 1: 01340 return -1./6. + 3./4.*xi - 2*xi*xi + 17./12.*xi*xi*xi 01341 + 3./4.*eta - 5./2.*xi*eta + 7./4.*xi*xi*eta 01342 - eta*eta + 7./4.*xi*eta*eta + 5./12.*eta*eta*eta; 01343 case 2: 01344 return -xi*xi + 13./12.*xi*xi*xi 01345 - 1./4.*xi*xi*eta; 01346 } 01347 case 6: 01348 switch (subtriangle_lookup(p)) 01349 { 01350 case 0: 01351 return -xi*eta + 2*xi*xi*eta 01352 + 1./2.*eta*eta + 7./4.*xi*eta*eta - 13./12.*eta*eta*eta; 01353 case 1: 01354 return 2./3. - 13./4.*xi + 9./2.*xi*xi - 23./12.*xi*xi*xi 01355 - 11./4.*eta + 19./2.*xi*eta - 23./4.*xi*xi*eta 01356 + 7./2.*eta*eta - 25./4.*xi*eta*eta - 17./12.*eta*eta*eta; 01357 case 2: 01358 return -1./2.*xi*xi + 5./12.*xi*xi*xi 01359 + 9./4.*xi*xi*eta; 01360 } 01361 case 7: 01362 switch (subtriangle_lookup(p)) 01363 { 01364 case 0: 01365 return -1./2.*eta*eta + 9./4.*xi*eta*eta + 5./12.*eta*eta*eta; 01366 case 1: 01367 return 2./3. - 11./4.*xi + 7./2.*xi*xi - 17./12.*xi*xi*xi 01368 - 13./4.*eta + 19./2.*xi*eta - 25./4.*xi*xi*eta 01369 + 9./2.*eta*eta - 23./4.*xi*eta*eta - 23./12.*eta*eta*eta; 01370 case 2: 01371 return 1./2.*xi*xi - 13./12.*xi*xi*xi 01372 - xi*eta + 7./4.*xi*xi*eta + 2*xi*eta*eta; 01373 } 01374 case 8: 01375 switch (subtriangle_lookup(p)) 01376 { 01377 case 0: 01378 return -eta*eta - 1./4.*xi*eta*eta + 13./12.*eta*eta*eta; 01379 case 1: 01380 return -1./6. + 3./4.*xi - xi*xi + 5./12.*xi*xi*xi 01381 + 3./4.*eta - 5./2.*xi*eta + 7./4.*xi*xi*eta 01382 - 2*eta*eta + 7./4.*xi*eta*eta + 17./12.*eta*eta*eta; 01383 case 2: 01384 return 1./12.*xi*xi*xi 01385 - 1./4.*xi*xi*eta 01386 - eta*eta + eta*eta*eta; 01387 } 01388 case 9: 01389 switch (subtriangle_lookup(p)) 01390 { 01391 case 0: 01392 return std::sqrt(2.) * (xi*eta*eta - 1./3.*eta*eta*eta); 01393 case 1: 01394 return std::sqrt(2.) * (2./3. - 3*xi + 4*xi*xi - 5./3.*xi*xi*xi 01395 - 3*eta + 10*xi*eta - 7*xi*xi*eta 01396 + 4*eta*eta - 7*xi*eta*eta - 5./3.*eta*eta*eta); 01397 case 2: 01398 return std::sqrt(2.) * (-1./3.*xi*xi*xi + xi*xi*eta); 01399 } 01400 case 10: 01401 switch (subtriangle_lookup(p)) 01402 { 01403 case 0: 01404 return 2*eta*eta - 2*xi*eta*eta - 8./3.*eta*eta*eta; 01405 case 1: 01406 return -2./3. + 2*xi - 2*xi*xi + 2./3.*xi*xi*xi 01407 + 4*eta - 8*xi*eta + 4*xi*xi*eta 01408 - 6*eta*eta + 6*xi*eta*eta + 8./3.*eta*eta*eta; 01409 case 2: 01410 return -2*xi*xi + 10./3.*xi*xi*xi 01411 + 4*xi*eta - 4*xi*xi*eta 01412 - 4*xi*eta*eta; 01413 } 01414 case 11: 01415 switch (subtriangle_lookup(p)) 01416 { 01417 case 0: 01418 return 4*xi*eta - 4*xi*xi*eta 01419 - 2*eta*eta - 4*xi*eta*eta + 10./3.*eta*eta*eta; 01420 case 1: 01421 return -2./3. + 4*xi - 6*xi*xi + 8./3.*xi*xi*xi 01422 + 2*eta - 8*xi*eta + 6*xi*xi*eta 01423 - 2*eta*eta + 4*xi*eta*eta + 2./3.*eta*eta*eta; 01424 case 2: 01425 return 2*xi*xi - 8./3.*xi*xi*xi 01426 - 2*xi*xi*eta; 01427 } 01428 } 01429 01430 libmesh_error(); 01431 return 0.; 01432 } 01433 01434 01435 } // end anonymous namespace 01436 01437 01438 namespace libMesh 01439 { 01440 01441 01442 template <> 01443 Real FE<2,CLOUGH>::shape(const ElemType, 01444 const Order, 01445 const unsigned int, 01446 const Point&) 01447 { 01448 libMesh::err << "Clough-Tocher elements require the real element\n" 01449 << "to construct gradient-based degrees of freedom." 01450 << std::endl; 01451 01452 libmesh_error(); 01453 return 0.; 01454 } 01455 01456 01457 01458 template <> 01459 Real FE<2,CLOUGH>::shape(const Elem* elem, 01460 const Order order, 01461 const unsigned int i, 01462 const Point& p) 01463 { 01464 libmesh_assert(elem); 01465 01466 clough_compute_coefs(elem); 01467 01468 const ElemType type = elem->type(); 01469 01470 const Order totalorder = static_cast<Order>(order + elem->p_level()); 01471 01472 switch (totalorder) 01473 { 01474 // 2nd-order restricted Clough-Tocher element 01475 case SECOND: 01476 { 01477 // There may be a bug in the 2nd order case; the 3rd order 01478 // Clough-Tocher elements are pretty uniformly better anyways 01479 // so use those instead. 01480 libmesh_experimental(); 01481 switch (type) 01482 { 01483 // C1 functions on the Clough-Tocher triangle. 01484 case TRI6: 01485 { 01486 libmesh_assert_less (i, 9); 01487 // FIXME: it would be nice to calculate (and cache) 01488 // clough_raw_shape(j,p) only once per triangle, not 1-7 01489 // times 01490 switch (i) 01491 { 01492 // Note: these DoF numbers are "scrambled" because my 01493 // initial numbering conventions didn't match libMesh 01494 case 0: 01495 return clough_raw_shape(0, p) 01496 + d1d2n * clough_raw_shape(10, p) 01497 + d1d3n * clough_raw_shape(11, p); 01498 case 3: 01499 return clough_raw_shape(1, p) 01500 + d2d3n * clough_raw_shape(11, p) 01501 + d2d1n * clough_raw_shape(9, p); 01502 case 6: 01503 return clough_raw_shape(2, p) 01504 + d3d1n * clough_raw_shape(9, p) 01505 + d3d2n * clough_raw_shape(10, p); 01506 case 1: 01507 return d1xd1x * clough_raw_shape(3, p) 01508 + d1xd1y * clough_raw_shape(4, p) 01509 + d1xd2n * clough_raw_shape(10, p) 01510 + d1xd3n * clough_raw_shape(11, p) 01511 + 0.5 * N01x * d3nd3n * clough_raw_shape(11, p) 01512 + 0.5 * N02x * d2nd2n * clough_raw_shape(10, p); 01513 case 2: 01514 return d1yd1y * clough_raw_shape(4, p) 01515 + d1yd1x * clough_raw_shape(3, p) 01516 + d1yd2n * clough_raw_shape(10, p) 01517 + d1yd3n * clough_raw_shape(11, p) 01518 + 0.5 * N01y * d3nd3n * clough_raw_shape(11, p) 01519 + 0.5 * N02y * d2nd2n * clough_raw_shape(10, p); 01520 case 4: 01521 return d2xd2x * clough_raw_shape(5, p) 01522 + d2xd2y * clough_raw_shape(6, p) 01523 + d2xd3n * clough_raw_shape(11, p) 01524 + d2xd1n * clough_raw_shape(9, p) 01525 + 0.5 * N10x * d3nd3n * clough_raw_shape(11, p) 01526 + 0.5 * N12x * d1nd1n * clough_raw_shape(9, p); 01527 case 5: 01528 return d2yd2y * clough_raw_shape(6, p) 01529 + d2yd2x * clough_raw_shape(5, p) 01530 + d2yd3n * clough_raw_shape(11, p) 01531 + d2yd1n * clough_raw_shape(9, p) 01532 + 0.5 * N10y * d3nd3n * clough_raw_shape(11, p) 01533 + 0.5 * N12y * d1nd1n * clough_raw_shape(9, p); 01534 case 7: 01535 return d3xd3x * clough_raw_shape(7, p) 01536 + d3xd3y * clough_raw_shape(8, p) 01537 + d3xd1n * clough_raw_shape(9, p) 01538 + d3xd2n * clough_raw_shape(10, p) 01539 + 0.5 * N20x * d2nd2n * clough_raw_shape(10, p) 01540 + 0.5 * N21x * d1nd1n * clough_raw_shape(9, p); 01541 case 8: 01542 return d3yd3y * clough_raw_shape(8, p) 01543 + d3yd3x * clough_raw_shape(7, p) 01544 + d3yd1n * clough_raw_shape(9, p) 01545 + d3yd2n * clough_raw_shape(10, p) 01546 + 0.5 * N20y * d2nd2n * clough_raw_shape(10, p) 01547 + 0.5 * N21y * d1nd1n * clough_raw_shape(9, p); 01548 default: 01549 libmesh_error(); 01550 } 01551 } 01552 default: 01553 libMesh::err << "ERROR: Unsupported element type!" << std::endl; 01554 libmesh_error(); 01555 } 01556 } 01557 // 3rd-order Clough-Tocher element 01558 case THIRD: 01559 { 01560 switch (type) 01561 { 01562 // C1 functions on the Clough-Tocher triangle. 01563 case TRI6: 01564 { 01565 libmesh_assert_less (i, 12); 01566 01567 // FIXME: it would be nice to calculate (and cache) 01568 // clough_raw_shape(j,p) only once per triangle, not 1-7 01569 // times 01570 switch (i) 01571 { 01572 // Note: these DoF numbers are "scrambled" because my 01573 // initial numbering conventions didn't match libMesh 01574 case 0: 01575 return clough_raw_shape(0, p) 01576 + d1d2n * clough_raw_shape(10, p) 01577 + d1d3n * clough_raw_shape(11, p); 01578 case 3: 01579 return clough_raw_shape(1, p) 01580 + d2d3n * clough_raw_shape(11, p) 01581 + d2d1n * clough_raw_shape(9, p); 01582 case 6: 01583 return clough_raw_shape(2, p) 01584 + d3d1n * clough_raw_shape(9, p) 01585 + d3d2n * clough_raw_shape(10, p); 01586 case 1: 01587 return d1xd1x * clough_raw_shape(3, p) 01588 + d1xd1y * clough_raw_shape(4, p) 01589 + d1xd2n * clough_raw_shape(10, p) 01590 + d1xd3n * clough_raw_shape(11, p); 01591 case 2: 01592 return d1yd1y * clough_raw_shape(4, p) 01593 + d1yd1x * clough_raw_shape(3, p) 01594 + d1yd2n * clough_raw_shape(10, p) 01595 + d1yd3n * clough_raw_shape(11, p); 01596 case 4: 01597 return d2xd2x * clough_raw_shape(5, p) 01598 + d2xd2y * clough_raw_shape(6, p) 01599 + d2xd3n * clough_raw_shape(11, p) 01600 + d2xd1n * clough_raw_shape(9, p); 01601 case 5: 01602 return d2yd2y * clough_raw_shape(6, p) 01603 + d2yd2x * clough_raw_shape(5, p) 01604 + d2yd3n * clough_raw_shape(11, p) 01605 + d2yd1n * clough_raw_shape(9, p); 01606 case 7: 01607 return d3xd3x * clough_raw_shape(7, p) 01608 + d3xd3y * clough_raw_shape(8, p) 01609 + d3xd1n * clough_raw_shape(9, p) 01610 + d3xd2n * clough_raw_shape(10, p); 01611 case 8: 01612 return d3yd3y * clough_raw_shape(8, p) 01613 + d3yd3x * clough_raw_shape(7, p) 01614 + d3yd1n * clough_raw_shape(9, p) 01615 + d3yd2n * clough_raw_shape(10, p); 01616 case 10: 01617 return d1nd1n * clough_raw_shape(9, p); 01618 case 11: 01619 return d2nd2n * clough_raw_shape(10, p); 01620 case 9: 01621 return d3nd3n * clough_raw_shape(11, p); 01622 01623 default: 01624 libmesh_error(); 01625 } 01626 } 01627 default: 01628 libMesh::err << "ERROR: Unsupported element type!" << std::endl; 01629 libmesh_error(); 01630 } 01631 } 01632 // by default throw an error 01633 default: 01634 libMesh::err << "ERROR: Unsupported polynomial order!" << std::endl; 01635 libmesh_error(); 01636 } 01637 01638 libmesh_error(); 01639 return 0.; 01640 } 01641 01642 01643 01644 template <> 01645 Real FE<2,CLOUGH>::shape_deriv(const ElemType, 01646 const Order, 01647 const unsigned int, 01648 const unsigned int, 01649 const Point&) 01650 { 01651 libMesh::err << "Clough-Tocher elements require the real element\n" 01652 << "to construct gradient-based degrees of freedom." 01653 << std::endl; 01654 01655 libmesh_error(); 01656 return 0.; 01657 } 01658 01659 01660 01661 template <> 01662 Real FE<2,CLOUGH>::shape_deriv(const Elem* elem, 01663 const Order order, 01664 const unsigned int i, 01665 const unsigned int j, 01666 const Point& p) 01667 { 01668 libmesh_assert(elem); 01669 01670 clough_compute_coefs(elem); 01671 01672 const ElemType type = elem->type(); 01673 01674 const Order totalorder = static_cast<Order>(order + elem->p_level()); 01675 01676 switch (totalorder) 01677 { 01678 // 2nd-order restricted Clough-Tocher element 01679 case SECOND: 01680 { 01681 // There may be a bug in the 2nd order case; the 3rd order 01682 // Clough-Tocher elements are pretty uniformly better anyways 01683 // so use those instead. 01684 libmesh_experimental(); 01685 switch (type) 01686 { 01687 // C1 functions on the Clough-Tocher triangle. 01688 case TRI6: 01689 { 01690 libmesh_assert_less (i, 9); 01691 // FIXME: it would be nice to calculate (and cache) 01692 // clough_raw_shape(j,p) only once per triangle, not 1-7 01693 // times 01694 switch (i) 01695 { 01696 // Note: these DoF numbers are "scrambled" because my 01697 // initial numbering conventions didn't match libMesh 01698 case 0: 01699 return clough_raw_shape_deriv(0, j, p) 01700 + d1d2n * clough_raw_shape_deriv(10, j, p) 01701 + d1d3n * clough_raw_shape_deriv(11, j, p); 01702 case 3: 01703 return clough_raw_shape_deriv(1, j, p) 01704 + d2d3n * clough_raw_shape_deriv(11, j, p) 01705 + d2d1n * clough_raw_shape_deriv(9, j, p); 01706 case 6: 01707 return clough_raw_shape_deriv(2, j, p) 01708 + d3d1n * clough_raw_shape_deriv(9, j, p) 01709 + d3d2n * clough_raw_shape_deriv(10, j, p); 01710 case 1: 01711 return d1xd1x * clough_raw_shape_deriv(3, j, p) 01712 + d1xd1y * clough_raw_shape_deriv(4, j, p) 01713 + d1xd2n * clough_raw_shape_deriv(10, j, p) 01714 + d1xd3n * clough_raw_shape_deriv(11, j, p) 01715 + 0.5 * N01x * d3nd3n * clough_raw_shape_deriv(11, j, p) 01716 + 0.5 * N02x * d2nd2n * clough_raw_shape_deriv(10, j, p); 01717 case 2: 01718 return d1yd1y * clough_raw_shape_deriv(4, j, p) 01719 + d1yd1x * clough_raw_shape_deriv(3, j, p) 01720 + d1yd2n * clough_raw_shape_deriv(10, j, p) 01721 + d1yd3n * clough_raw_shape_deriv(11, j, p) 01722 + 0.5 * N01y * d3nd3n * clough_raw_shape_deriv(11, j, p) 01723 + 0.5 * N02y * d2nd2n * clough_raw_shape_deriv(10, j, p); 01724 case 4: 01725 return d2xd2x * clough_raw_shape_deriv(5, j, p) 01726 + d2xd2y * clough_raw_shape_deriv(6, j, p) 01727 + d2xd3n * clough_raw_shape_deriv(11, j, p) 01728 + d2xd1n * clough_raw_shape_deriv(9, j, p) 01729 + 0.5 * N10x * d3nd3n * clough_raw_shape_deriv(11, j, p) 01730 + 0.5 * N12x * d1nd1n * clough_raw_shape_deriv(9, j, p); 01731 case 5: 01732 return d2yd2y * clough_raw_shape_deriv(6, j, p) 01733 + d2yd2x * clough_raw_shape_deriv(5, j, p) 01734 + d2yd3n * clough_raw_shape_deriv(11, j, p) 01735 + d2yd1n * clough_raw_shape_deriv(9, j, p) 01736 + 0.5 * N10y * d3nd3n * clough_raw_shape_deriv(11, j, p) 01737 + 0.5 * N12y * d1nd1n * clough_raw_shape_deriv(9, j, p); 01738 case 7: 01739 return d3xd3x * clough_raw_shape_deriv(7, j, p) 01740 + d3xd3y * clough_raw_shape_deriv(8, j, p) 01741 + d3xd1n * clough_raw_shape_deriv(9, j, p) 01742 + d3xd2n * clough_raw_shape_deriv(10, j, p) 01743 + 0.5 * N20x * d2nd2n * clough_raw_shape_deriv(10, j, p) 01744 + 0.5 * N21x * d1nd1n * clough_raw_shape_deriv(9, j, p); 01745 case 8: 01746 return d3yd3y * clough_raw_shape_deriv(8, j, p) 01747 + d3yd3x * clough_raw_shape_deriv(7, j, p) 01748 + d3yd1n * clough_raw_shape_deriv(9, j, p) 01749 + d3yd2n * clough_raw_shape_deriv(10, j, p) 01750 + 0.5 * N20y * d2nd2n * clough_raw_shape_deriv(10, j, p) 01751 + 0.5 * N21y * d1nd1n * clough_raw_shape_deriv(9, j, p); 01752 default: 01753 libmesh_error(); 01754 } 01755 } 01756 default: 01757 libMesh::err << "ERROR: Unsupported element type!" << std::endl; 01758 libmesh_error(); 01759 } 01760 } 01761 // 3rd-order Clough-Tocher element 01762 case THIRD: 01763 { 01764 switch (type) 01765 { 01766 // C1 functions on the Clough-Tocher triangle. 01767 case TRI6: 01768 { 01769 libmesh_assert_less (i, 12); 01770 01771 // FIXME: it would be nice to calculate (and cache) 01772 // clough_raw_shape(j,p) only once per triangle, not 1-7 01773 // times 01774 switch (i) 01775 { 01776 // Note: these DoF numbers are "scrambled" because my 01777 // initial numbering conventions didn't match libMesh 01778 case 0: 01779 return clough_raw_shape_deriv(0, j, p) 01780 + d1d2n * clough_raw_shape_deriv(10, j, p) 01781 + d1d3n * clough_raw_shape_deriv(11, j, p); 01782 case 3: 01783 return clough_raw_shape_deriv(1, j, p) 01784 + d2d3n * clough_raw_shape_deriv(11, j, p) 01785 + d2d1n * clough_raw_shape_deriv(9, j, p); 01786 case 6: 01787 return clough_raw_shape_deriv(2, j, p) 01788 + d3d1n * clough_raw_shape_deriv(9, j, p) 01789 + d3d2n * clough_raw_shape_deriv(10, j, p); 01790 case 1: 01791 return d1xd1x * clough_raw_shape_deriv(3, j, p) 01792 + d1xd1y * clough_raw_shape_deriv(4, j, p) 01793 + d1xd2n * clough_raw_shape_deriv(10, j, p) 01794 + d1xd3n * clough_raw_shape_deriv(11, j, p); 01795 case 2: 01796 return d1yd1y * clough_raw_shape_deriv(4, j, p) 01797 + d1yd1x * clough_raw_shape_deriv(3, j, p) 01798 + d1yd2n * clough_raw_shape_deriv(10, j, p) 01799 + d1yd3n * clough_raw_shape_deriv(11, j, p); 01800 case 4: 01801 return d2xd2x * clough_raw_shape_deriv(5, j, p) 01802 + d2xd2y * clough_raw_shape_deriv(6, j, p) 01803 + d2xd3n * clough_raw_shape_deriv(11, j, p) 01804 + d2xd1n * clough_raw_shape_deriv(9, j, p); 01805 case 5: 01806 return d2yd2y * clough_raw_shape_deriv(6, j, p) 01807 + d2yd2x * clough_raw_shape_deriv(5, j, p) 01808 + d2yd3n * clough_raw_shape_deriv(11, j, p) 01809 + d2yd1n * clough_raw_shape_deriv(9, j, p); 01810 case 7: 01811 return d3xd3x * clough_raw_shape_deriv(7, j, p) 01812 + d3xd3y * clough_raw_shape_deriv(8, j, p) 01813 + d3xd1n * clough_raw_shape_deriv(9, j, p) 01814 + d3xd2n * clough_raw_shape_deriv(10, j, p); 01815 case 8: 01816 return d3yd3y * clough_raw_shape_deriv(8, j, p) 01817 + d3yd3x * clough_raw_shape_deriv(7, j, p) 01818 + d3yd1n * clough_raw_shape_deriv(9, j, p) 01819 + d3yd2n * clough_raw_shape_deriv(10, j, p); 01820 case 10: 01821 return d1nd1n * clough_raw_shape_deriv(9, j, p); 01822 case 11: 01823 return d2nd2n * clough_raw_shape_deriv(10, j, p); 01824 case 9: 01825 return d3nd3n * clough_raw_shape_deriv(11, j, p); 01826 01827 default: 01828 libmesh_error(); 01829 } 01830 } 01831 default: 01832 libMesh::err << "ERROR: Unsupported element type!" << std::endl; 01833 libmesh_error(); 01834 } 01835 } 01836 // by default throw an error 01837 default: 01838 libMesh::err << "ERROR: Unsupported polynomial order!" << std::endl; 01839 libmesh_error(); 01840 } 01841 01842 libmesh_error(); 01843 return 0.; 01844 } 01845 01846 01847 01848 template <> 01849 Real FE<2,CLOUGH>::shape_second_deriv(const Elem* elem, 01850 const Order order, 01851 const unsigned int i, 01852 const unsigned int j, 01853 const Point& p) 01854 { 01855 libmesh_assert(elem); 01856 01857 clough_compute_coefs(elem); 01858 01859 const ElemType type = elem->type(); 01860 01861 const Order totalorder = static_cast<Order>(order + elem->p_level()); 01862 01863 switch (totalorder) 01864 { 01865 // 2nd-order restricted Clough-Tocher element 01866 case SECOND: 01867 { 01868 switch (type) 01869 { 01870 // C1 functions on the Clough-Tocher triangle. 01871 case TRI6: 01872 { 01873 libmesh_assert_less (i, 9); 01874 // FIXME: it would be nice to calculate (and cache) 01875 // clough_raw_shape(j,p) only once per triangle, not 1-7 01876 // times 01877 switch (i) 01878 { 01879 // Note: these DoF numbers are "scrambled" because my 01880 // initial numbering conventions didn't match libMesh 01881 case 0: 01882 return clough_raw_shape_second_deriv(0, j, p) 01883 + d1d2n * clough_raw_shape_second_deriv(10, j, p) 01884 + d1d3n * clough_raw_shape_second_deriv(11, j, p); 01885 case 3: 01886 return clough_raw_shape_second_deriv(1, j, p) 01887 + d2d3n * clough_raw_shape_second_deriv(11, j, p) 01888 + d2d1n * clough_raw_shape_second_deriv(9, j, p); 01889 case 6: 01890 return clough_raw_shape_second_deriv(2, j, p) 01891 + d3d1n * clough_raw_shape_second_deriv(9, j, p) 01892 + d3d2n * clough_raw_shape_second_deriv(10, j, p); 01893 case 1: 01894 return d1xd1x * clough_raw_shape_second_deriv(3, j, p) 01895 + d1xd1y * clough_raw_shape_second_deriv(4, j, p) 01896 + d1xd2n * clough_raw_shape_second_deriv(10, j, p) 01897 + d1xd3n * clough_raw_shape_second_deriv(11, j, p) 01898 + 0.5 * N01x * d3nd3n * clough_raw_shape_second_deriv(11, j, p) 01899 + 0.5 * N02x * d2nd2n * clough_raw_shape_second_deriv(10, j, p); 01900 case 2: 01901 return d1yd1y * clough_raw_shape_second_deriv(4, j, p) 01902 + d1yd1x * clough_raw_shape_second_deriv(3, j, p) 01903 + d1yd2n * clough_raw_shape_second_deriv(10, j, p) 01904 + d1yd3n * clough_raw_shape_second_deriv(11, j, p) 01905 + 0.5 * N01y * d3nd3n * clough_raw_shape_second_deriv(11, j, p) 01906 + 0.5 * N02y * d2nd2n * clough_raw_shape_second_deriv(10, j, p); 01907 case 4: 01908 return d2xd2x * clough_raw_shape_second_deriv(5, j, p) 01909 + d2xd2y * clough_raw_shape_second_deriv(6, j, p) 01910 + d2xd3n * clough_raw_shape_second_deriv(11, j, p) 01911 + d2xd1n * clough_raw_shape_second_deriv(9, j, p) 01912 + 0.5 * N10x * d3nd3n * clough_raw_shape_second_deriv(11, j, p) 01913 + 0.5 * N12x * d1nd1n * clough_raw_shape_second_deriv(9, j, p); 01914 case 5: 01915 return d2yd2y * clough_raw_shape_second_deriv(6, j, p) 01916 + d2yd2x * clough_raw_shape_second_deriv(5, j, p) 01917 + d2yd3n * clough_raw_shape_second_deriv(11, j, p) 01918 + d2yd1n * clough_raw_shape_second_deriv(9, j, p) 01919 + 0.5 * N10y * d3nd3n * clough_raw_shape_second_deriv(11, j, p) 01920 + 0.5 * N12y * d1nd1n * clough_raw_shape_second_deriv(9, j, p); 01921 case 7: 01922 return d3xd3x * clough_raw_shape_second_deriv(7, j, p) 01923 + d3xd3y * clough_raw_shape_second_deriv(8, j, p) 01924 + d3xd1n * clough_raw_shape_second_deriv(9, j, p) 01925 + d3xd2n * clough_raw_shape_second_deriv(10, j, p) 01926 + 0.5 * N20x * d2nd2n * clough_raw_shape_second_deriv(10, j, p) 01927 + 0.5 * N21x * d1nd1n * clough_raw_shape_second_deriv(9, j, p); 01928 case 8: 01929 return d3yd3y * clough_raw_shape_second_deriv(8, j, p) 01930 + d3yd3x * clough_raw_shape_second_deriv(7, j, p) 01931 + d3yd1n * clough_raw_shape_second_deriv(9, j, p) 01932 + d3yd2n * clough_raw_shape_second_deriv(10, j, p) 01933 + 0.5 * N20y * d2nd2n * clough_raw_shape_second_deriv(10, j, p) 01934 + 0.5 * N21y * d1nd1n * clough_raw_shape_second_deriv(9, j, p); 01935 default: 01936 libmesh_error(); 01937 } 01938 } 01939 default: 01940 libMesh::err << "ERROR: Unsupported element type!" << std::endl; 01941 libmesh_error(); 01942 } 01943 } 01944 // 3rd-order Clough-Tocher element 01945 case THIRD: 01946 { 01947 switch (type) 01948 { 01949 // C1 functions on the Clough-Tocher triangle. 01950 case TRI6: 01951 { 01952 libmesh_assert_less (i, 12); 01953 01954 // FIXME: it would be nice to calculate (and cache) 01955 // clough_raw_shape(j,p) only once per triangle, not 1-7 01956 // times 01957 switch (i) 01958 { 01959 // Note: these DoF numbers are "scrambled" because my 01960 // initial numbering conventions didn't match libMesh 01961 case 0: 01962 return clough_raw_shape_second_deriv(0, j, p) 01963 + d1d2n * clough_raw_shape_second_deriv(10, j, p) 01964 + d1d3n * clough_raw_shape_second_deriv(11, j, p); 01965 case 3: 01966 return clough_raw_shape_second_deriv(1, j, p) 01967 + d2d3n * clough_raw_shape_second_deriv(11, j, p) 01968 + d2d1n * clough_raw_shape_second_deriv(9, j, p); 01969 case 6: 01970 return clough_raw_shape_second_deriv(2, j, p) 01971 + d3d1n * clough_raw_shape_second_deriv(9, j, p) 01972 + d3d2n * clough_raw_shape_second_deriv(10, j, p); 01973 case 1: 01974 return d1xd1x * clough_raw_shape_second_deriv(3, j, p) 01975 + d1xd1y * clough_raw_shape_second_deriv(4, j, p) 01976 + d1xd2n * clough_raw_shape_second_deriv(10, j, p) 01977 + d1xd3n * clough_raw_shape_second_deriv(11, j, p); 01978 case 2: 01979 return d1yd1y * clough_raw_shape_second_deriv(4, j, p) 01980 + d1yd1x * clough_raw_shape_second_deriv(3, j, p) 01981 + d1yd2n * clough_raw_shape_second_deriv(10, j, p) 01982 + d1yd3n * clough_raw_shape_second_deriv(11, j, p); 01983 case 4: 01984 return d2xd2x * clough_raw_shape_second_deriv(5, j, p) 01985 + d2xd2y * clough_raw_shape_second_deriv(6, j, p) 01986 + d2xd3n * clough_raw_shape_second_deriv(11, j, p) 01987 + d2xd1n * clough_raw_shape_second_deriv(9, j, p); 01988 case 5: 01989 return d2yd2y * clough_raw_shape_second_deriv(6, j, p) 01990 + d2yd2x * clough_raw_shape_second_deriv(5, j, p) 01991 + d2yd3n * clough_raw_shape_second_deriv(11, j, p) 01992 + d2yd1n * clough_raw_shape_second_deriv(9, j, p); 01993 case 7: 01994 return d3xd3x * clough_raw_shape_second_deriv(7, j, p) 01995 + d3xd3y * clough_raw_shape_second_deriv(8, j, p) 01996 + d3xd1n * clough_raw_shape_second_deriv(9, j, p) 01997 + d3xd2n * clough_raw_shape_second_deriv(10, j, p); 01998 case 8: 01999 return d3yd3y * clough_raw_shape_second_deriv(8, j, p) 02000 + d3yd3x * clough_raw_shape_second_deriv(7, j, p) 02001 + d3yd1n * clough_raw_shape_second_deriv(9, j, p) 02002 + d3yd2n * clough_raw_shape_second_deriv(10, j, p); 02003 case 10: 02004 return d1nd1n * clough_raw_shape_second_deriv(9, j, p); 02005 case 11: 02006 return d2nd2n * clough_raw_shape_second_deriv(10, j, p); 02007 case 9: 02008 return d3nd3n * clough_raw_shape_second_deriv(11, j, p); 02009 02010 default: 02011 libmesh_error(); 02012 } 02013 } 02014 default: 02015 libMesh::err << "ERROR: Unsupported element type!" << std::endl; 02016 libmesh_error(); 02017 } 02018 } 02019 // by default throw an error 02020 default: 02021 libMesh::err << "ERROR: Unsupported polynomial order!" << std::endl; 02022 libmesh_error(); 02023 } 02024 02025 libmesh_error(); 02026 return 0.; 02027 } 02028 02029 } // namespace libMesh
Site Created By: libMesh Developers
Last modified: February 05 2013 19:54:46 UTC
Hosted By: