MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /* ***************************************************************** 00002 MESQUITE -- The Mesh Quality Improvement Toolkit 00003 00004 Copyright 2006 Lawrence Livermore National Laboratory. Under 00005 the terms of Contract B545069 with the University of Wisconsin -- 00006 Madison, Lawrence Livermore National Laboratory retains certain 00007 rights in this software. 00008 00009 This library is free software; you can redistribute it and/or 00010 modify it under the terms of the GNU Lesser General Public 00011 License as published by the Free Software Foundation; either 00012 version 2.1 of the License, or (at your option) any later version. 00013 00014 This library is distributed in the hope that it will be useful, 00015 but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00017 Lesser General Public License for more details. 00018 00019 You should have received a copy of the GNU Lesser General Public License 00020 (lgpl.txt) along with this library; if not, write to the Free Software 00021 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00022 00023 kraftche@cae.wisc.edu 00024 00025 ***************************************************************** */ 00026 /** \file QuadLagrangeShape.cpp 00027 * \author Jason Kraftcheck 00028 */ 00029 00030 #include "QuadLagrangeShape.hpp" 00031 #include "LinearQuadrilateral.hpp" 00032 #include "MsqError.hpp" 00033 00034 namespace MBMesquite 00035 { 00036 00037 EntityTopology QuadLagrangeShape::element_topology() const 00038 { 00039 return QUADRILATERAL; 00040 } 00041 00042 int QuadLagrangeShape::num_nodes() const 00043 { 00044 return 9; 00045 } 00046 00047 void QuadLagrangeShape::coefficients( Sample loc, NodeSet nodeset, double* coeff_out, size_t* indices_out, 00048 size_t& num_coeff, MsqError& err ) const 00049 { 00050 if( !nodeset.have_any_mid_node() ) 00051 { 00052 LinearQuadrilateral::coefficients( loc, nodeset, coeff_out, indices_out, num_coeff ); 00053 return; 00054 } 00055 00056 switch( loc.dimension ) 00057 { 00058 case 0: 00059 num_coeff = 1; 00060 indices_out[0] = loc.number; 00061 coeff_out[0] = 1.0; 00062 break; 00063 case 1: 00064 coeff_out[0] = coeff_out[1] = coeff_out[2] = coeff_out[3] = coeff_out[4] = coeff_out[5] = coeff_out[6] = 00065 coeff_out[7] = coeff_out[8] = 0.0; 00066 if( nodeset.mid_edge_node( loc.number ) ) 00067 { 00068 // if mid-edge node is present 00069 num_coeff = 1; 00070 indices_out[0] = loc.number + 4; 00071 coeff_out[0] = 1.0; 00072 } 00073 else 00074 { 00075 // If mid-edge node is not present, mapping function value 00076 // for linear edge is even weight of adjacent vertices. 00077 num_coeff = 2; 00078 indices_out[0] = loc.number; 00079 indices_out[1] = ( loc.number + 1 ) % 4; 00080 coeff_out[0] = 0.5; 00081 coeff_out[1] = 0.5; 00082 } 00083 break; 00084 case 2: 00085 if( nodeset.mid_face_node( 0 ) ) 00086 { // if quad center node is present 00087 num_coeff = 1; 00088 indices_out[0] = 8; 00089 coeff_out[0] = 1.0; 00090 } 00091 else 00092 { 00093 // for linear element, (no mid-edge nodes), all corners contribute 1/4. 00094 num_coeff = 4; 00095 indices_out[0] = 0; 00096 indices_out[1] = 1; 00097 indices_out[2] = 2; 00098 indices_out[3] = 3; 00099 coeff_out[0] = 0.25; 00100 coeff_out[1] = 0.25; 00101 coeff_out[2] = 0.25; 00102 coeff_out[3] = 0.25; 00103 // add in contribution for any mid-edge nodes present 00104 for( int i = 0; i < 4; ++i ) 00105 { // for each edge 00106 if( nodeset.mid_edge_node( i ) ) 00107 { 00108 indices_out[num_coeff] = i + 4; 00109 coeff_out[num_coeff] = 0.5; 00110 coeff_out[i] -= 0.25; 00111 coeff_out[( i + 1 ) % 4] -= 0.25; 00112 ++num_coeff; 00113 } 00114 } 00115 } 00116 break; 00117 default: 00118 MSQ_SETERR( err ) 00119 ( MsqError::UNSUPPORTED_ELEMENT, 00120 "Request for dimension %d mapping function value" 00121 "for a quadrilateral element", 00122 loc.dimension ); 00123 } 00124 } 00125 00126 static void derivatives_at_corner( unsigned corner, NodeSet nodeset, size_t* vertices, MsqVector< 2 >* derivs, 00127 size_t& num_vtx ) 00128 { 00129 static const unsigned xi_adj_corners[] = { 1, 0, 3, 2 }; 00130 static const unsigned xi_adj_edges[] = { 0, 0, 2, 2 }; 00131 static const unsigned eta_adj_corners[] = { 3, 2, 1, 0 }; 00132 static const unsigned eta_adj_edges[] = { 3, 1, 1, 3 }; 00133 00134 static const double corner_xi[] = { -1, 1, 1, -1 }; // xi values by corner 00135 static const double corner_eta[] = { -1, -1, 1, 1 }; // eta values by corner 00136 static const double other_xi[] = { 1, -1, -1, 1 }; // xi values for adjacent corner in xi direction 00137 static const double other_eta[] = { 1, 1, -1, -1 }; // eta values for adjcent corner in eta direction 00138 static const double mid_xi[] = { 2, -2, -2, 2 }; // xi values for mid-node in xi direction 00139 static const double mid_eta[] = { 2, 2, -2, -2 }; // eta values for mid-node in eta direction 00140 00141 num_vtx = 3; 00142 vertices[0] = corner; 00143 vertices[1] = xi_adj_corners[corner]; 00144 vertices[2] = eta_adj_corners[corner]; 00145 00146 derivs[0][0] = corner_xi[corner]; 00147 derivs[0][1] = corner_eta[corner]; 00148 derivs[1][0] = other_xi[corner]; 00149 derivs[1][1] = 0.0; 00150 derivs[2][0] = 0.0; 00151 derivs[2][1] = other_eta[corner]; 00152 00153 if( nodeset.mid_edge_node( xi_adj_edges[corner] ) ) 00154 { 00155 vertices[num_vtx] = 4 + xi_adj_edges[corner]; 00156 derivs[num_vtx][0] = 2.0 * mid_xi[corner]; 00157 derivs[num_vtx][1] = 0.0; 00158 derivs[0][0] -= mid_xi[corner]; 00159 derivs[1][0] -= mid_xi[corner]; 00160 ++num_vtx; 00161 } 00162 00163 if( nodeset.mid_edge_node( eta_adj_edges[corner] ) ) 00164 { 00165 vertices[num_vtx] = 4 + eta_adj_edges[corner]; 00166 derivs[num_vtx][0] = 0.0; 00167 derivs[num_vtx][1] = 2.0 * mid_eta[corner]; 00168 derivs[0][1] -= mid_eta[corner]; 00169 derivs[2][1] -= mid_eta[corner]; 00170 ++num_vtx; 00171 } 00172 } 00173 00174 static void derivatives_at_mid_edge( unsigned edge, NodeSet nodeset, size_t* vertices, MsqVector< 2 >* derivs, 00175 size_t& num_vtx ) 00176 { 00177 static const double values[][9] = { { -0.5, -0.5, 0.5, 0.5, -1.0, 2.0, 1.0, 2.0, 4.0 }, 00178 { -0.5, 0.5, 0.5, -0.5, -2.0, 1.0, -2.0, -1.0, -4.0 }, 00179 { -0.5, -0.5, 0.5, 0.5, -1.0, -2.0, 1.0, -2.0, -4.0 }, 00180 { -0.5, 0.5, 0.5, -0.5, 2.0, 1.0, 2.0, -1.0, 4.0 } }; 00181 static const double edge_values[][2] = { { -1, 1 }, { -1, 1 }, { 1, -1 }, { 1, -1 } }; 00182 const unsigned prev_corner = edge; // index of start vertex of edge 00183 const unsigned next_corner = ( edge + 1 ) % 4; // index of end vertex of edge 00184 const unsigned is_eta_edge = edge % 2; // edge is xi = +/- 0 00185 const unsigned is_xi_edge = 1 - is_eta_edge; // edge is eta = +/- 0 00186 // const unsigned mid_edge_node = edge + 4; // mid-edge node index 00187 const unsigned prev_opposite = ( prev_corner + 3 ) % 4; // index of corner adjacent to prev_corner 00188 const unsigned next_opposite = ( next_corner + 1 ) % 4; // index of corner adjacent to next_corner; 00189 00190 // First do derivatives along edge (e.g. wrt xi if edge is eta = +/-1) 00191 num_vtx = 2; 00192 vertices[0] = prev_corner; 00193 vertices[1] = next_corner; 00194 derivs[0][is_eta_edge] = edge_values[edge][0]; 00195 derivs[0][is_xi_edge] = 0.0; 00196 derivs[1][is_eta_edge] = edge_values[edge][1]; 00197 derivs[1][is_xi_edge] = 0.0; 00198 // That's it for the edge-direction derivatives. No other vertices contribute. 00199 00200 // Next handle the linear element case. Handle this as a special case first, 00201 // so the generalized solution doesn't impact performance for linear elements 00202 // too much. 00203 if( !nodeset.have_any_mid_node() ) 00204 { 00205 num_vtx = 4; 00206 vertices[2] = prev_opposite; 00207 vertices[3] = next_opposite; 00208 derivs[0][is_xi_edge] = values[edge][prev_corner]; 00209 derivs[1][is_xi_edge] = values[edge][next_corner]; 00210 derivs[2][is_xi_edge] = values[edge][prev_opposite]; 00211 derivs[2][is_eta_edge] = 0.0; 00212 derivs[3][is_xi_edge] = values[edge][next_opposite]; 00213 derivs[3][is_eta_edge] = 0.0; 00214 return; 00215 } 00216 00217 // Initial (linear) contribution for each corner 00218 double v[4] = { values[edge][0], values[edge][1], values[edge][2], values[edge][3] }; 00219 00220 // If mid-face node is present 00221 double v8 = 0.0; 00222 if( nodeset.mid_face_node( 0 ) ) 00223 { 00224 v8 = values[edge][8]; 00225 vertices[num_vtx] = 8; 00226 derivs[num_vtx][is_eta_edge] = 0.0; 00227 derivs[num_vtx][is_xi_edge] = v8; 00228 v[0] -= 0.25 * v8; 00229 v[1] -= 0.25 * v8; 00230 v[2] -= 0.25 * v8; 00231 v[3] -= 0.25 * v8; 00232 ++num_vtx; 00233 } 00234 00235 // If mid-edge nodes are present 00236 for( unsigned i = 0; i < 4; ++i ) 00237 { 00238 if( nodeset.mid_edge_node( i ) ) 00239 { 00240 const double value = values[edge][i + 4] - 0.5 * v8; 00241 if( fabs( value ) > 0.125 ) 00242 { 00243 v[i] -= 0.5 * value; 00244 v[( i + 1 ) % 4] -= 0.5 * value; 00245 vertices[num_vtx] = i + 4; 00246 derivs[num_vtx][is_eta_edge] = 0.0; 00247 derivs[num_vtx][is_xi_edge] = value; 00248 ++num_vtx; 00249 } 00250 } 00251 } 00252 00253 // update values for adjacent corners 00254 derivs[0][is_xi_edge] = v[prev_corner]; 00255 derivs[1][is_xi_edge] = v[next_corner]; 00256 // do other two corners 00257 if( fabs( v[prev_opposite] ) > 0.125 ) 00258 { 00259 vertices[num_vtx] = prev_opposite; 00260 derivs[num_vtx][is_eta_edge] = 0.0; 00261 derivs[num_vtx][is_xi_edge] = v[prev_opposite]; 00262 ++num_vtx; 00263 } 00264 if( fabs( v[next_opposite] ) > 0.125 ) 00265 { 00266 vertices[num_vtx] = next_opposite; 00267 derivs[num_vtx][is_eta_edge] = 0.0; 00268 derivs[num_vtx][is_xi_edge] = v[next_opposite]; 00269 ++num_vtx; 00270 } 00271 } 00272 00273 static void derivatives_at_mid_elem( NodeSet nodeset, size_t* vertices, MsqVector< 2 >* derivs, size_t& num_vtx ) 00274 { 00275 // fast linear case 00276 // This is provided as an optimization for linear elements. 00277 // If this block of code were removed, the general-case code 00278 // below should produce the same result. 00279 if( !nodeset.have_any_mid_node() ) 00280 { 00281 num_vtx = 4; 00282 vertices[0] = 0; 00283 derivs[0][0] = -0.5; 00284 derivs[0][1] = -0.5; 00285 vertices[1] = 1; 00286 derivs[1][0] = 0.5; 00287 derivs[1][1] = -0.5; 00288 vertices[2] = 2; 00289 derivs[2][0] = 0.5; 00290 derivs[2][1] = 0.5; 00291 vertices[3] = 3; 00292 derivs[3][0] = -0.5; 00293 derivs[3][1] = 0.5; 00294 return; 00295 } 00296 00297 num_vtx = 0; 00298 00299 // N_0 00300 if( !nodeset.both_edge_nodes( 0, 3 ) ) 00301 { // if eiter adjacent mid-edge node is missing 00302 vertices[num_vtx] = 0; 00303 derivs[num_vtx][0] = nodeset.mid_edge_node( 3 ) ? 0.0 : -0.5; 00304 derivs[num_vtx][1] = nodeset.mid_edge_node( 0 ) ? 0.0 : -0.5; 00305 ++num_vtx; 00306 } 00307 00308 // N_1 00309 if( !nodeset.both_edge_nodes( 0, 1 ) ) 00310 { // if eiter adjacent mid-edge node is missing 00311 vertices[num_vtx] = 1; 00312 derivs[num_vtx][0] = nodeset.mid_edge_node( 1 ) ? 0.0 : 0.5; 00313 derivs[num_vtx][1] = nodeset.mid_edge_node( 0 ) ? 0.0 : -0.5; 00314 ++num_vtx; 00315 } 00316 00317 // N_2 00318 if( !nodeset.both_edge_nodes( 1, 2 ) ) 00319 { // if eiter adjacent mid-edge node is missing 00320 vertices[num_vtx] = 2; 00321 derivs[num_vtx][0] = nodeset.mid_edge_node( 1 ) ? 0.0 : 0.5; 00322 derivs[num_vtx][1] = nodeset.mid_edge_node( 2 ) ? 0.0 : 0.5; 00323 ++num_vtx; 00324 } 00325 00326 // N_3 00327 if( !nodeset.both_edge_nodes( 2, 3 ) ) 00328 { // if eiter adjacent mid-edge node is missing 00329 vertices[num_vtx] = 3; 00330 derivs[num_vtx][0] = nodeset.mid_edge_node( 3 ) ? 0.0 : -0.5; 00331 derivs[num_vtx][1] = nodeset.mid_edge_node( 2 ) ? 0.0 : 0.5; 00332 ++num_vtx; 00333 } 00334 00335 // N_4 00336 if( nodeset.mid_edge_node( 0 ) ) 00337 { 00338 vertices[num_vtx] = 4; 00339 derivs[num_vtx][0] = 0.0; 00340 derivs[num_vtx][1] = -1.0; 00341 ++num_vtx; 00342 } 00343 00344 // N_5 00345 if( nodeset.mid_edge_node( 1 ) ) 00346 { 00347 vertices[num_vtx] = 5; 00348 derivs[num_vtx][0] = 1.0; 00349 derivs[num_vtx][1] = 0.0; 00350 ++num_vtx; 00351 } 00352 00353 // N_6 00354 if( nodeset.mid_edge_node( 2 ) ) 00355 { 00356 vertices[num_vtx] = 6; 00357 derivs[num_vtx][0] = 0.0; 00358 derivs[num_vtx][1] = 1.0; 00359 ++num_vtx; 00360 } 00361 00362 // N_7 00363 if( nodeset.mid_edge_node( 3 ) ) 00364 { 00365 vertices[num_vtx] = 7; 00366 derivs[num_vtx][0] = -1.0; 00367 derivs[num_vtx][1] = 0.0; 00368 ++num_vtx; 00369 } 00370 00371 // N_8 (mid-quad node) never contributes to Jacobian at element center!!! 00372 } 00373 00374 void QuadLagrangeShape::derivatives( Sample loc, NodeSet nodeset, size_t* vertex_indices_out, 00375 MsqVector< 2 >* d_coeff_d_xi_out, size_t& num_vtx, MsqError& err ) const 00376 { 00377 if( !nodeset.have_any_mid_node() ) 00378 { 00379 LinearQuadrilateral::derivatives( loc, nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx ); 00380 return; 00381 } 00382 00383 switch( loc.dimension ) 00384 { 00385 case 0: 00386 derivatives_at_corner( loc.number, nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx ); 00387 break; 00388 case 1: 00389 derivatives_at_mid_edge( loc.number, nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx ); 00390 break; 00391 case 2: 00392 derivatives_at_mid_elem( nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx ); 00393 break; 00394 default: 00395 MSQ_SETERR( err )( "Invalid/unsupported logical dimension", MsqError::INVALID_ARG ); 00396 } 00397 } 00398 00399 void QuadLagrangeShape::ideal( Sample, MsqMatrix< 3, 2 >& J, MsqError& ) const 00400 { 00401 J( 0, 0 ) = J( 1, 1 ) = 1.0; 00402 J( 0, 1 ) = J( 1, 0 ) = 0.0; 00403 J( 2, 0 ) = J( 2, 1 ) = 0.0; 00404 } 00405 00406 } // namespace MBMesquite