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 (2006) kraftche@cae.wisc.edu 00024 00025 ***************************************************************** */ 00026 00027 /** \file LinearHexahedron.cpp 00028 * \brief Implement shape function for linear hexahedron 00029 * \author Jason Kraftcheck 00030 */ 00031 00032 #include "Mesquite.hpp" 00033 #include "MsqError.hpp" 00034 #include "LinearHexahedron.hpp" 00035 00036 namespace MBMesquite 00037 { 00038 00039 static const char* nonlinear_error = "Attempt to use LinearHexahedron mapping function for a nonlinear element\n"; 00040 00041 EntityTopology LinearHexahedron::element_topology() const 00042 { 00043 return HEXAHEDRON; 00044 } 00045 00046 int LinearHexahedron::num_nodes() const 00047 { 00048 return 8; 00049 } 00050 00051 static inline int coeff_xi_sign( unsigned coeff ) 00052 { 00053 return 2 * ( ( ( coeff + 1 ) / 2 ) % 2 ) - 1; 00054 } 00055 static inline int coeff_eta_sign( unsigned coeff ) 00056 { 00057 return 2 * ( ( coeff / 2 ) % 2 ) - 1; 00058 } 00059 static inline int coeff_zeta_sign( unsigned coeff ) 00060 { 00061 return 2 * ( coeff / 4 ) - 1; 00062 } 00063 00064 static void coefficients_at_corner( unsigned corner, double* coeff_out, size_t* indices_out, size_t& num_coeff ) 00065 { 00066 num_coeff = 1; 00067 coeff_out[0] = 1.0; 00068 indices_out[0] = corner; 00069 } 00070 00071 const unsigned xi = 0, eta = 1, zeta = 2; 00072 const int edge_dir[] = { xi, eta, xi, eta, zeta, zeta, zeta, zeta, xi, eta, xi, eta }; 00073 const int edge_beg[] = { 0, 1, 2, 3, 0, 1, 2, 3, 4, 5, 6, 7 }; // start vertex by edge number 00074 const int edge_end[] = { 1, 2, 3, 0, 4, 5, 6, 7, 5, 6, 7, 4 }; // end vetex by edge number 00075 const int edge_opposite[] = { 10, 11, 8, 9, 6, 7, 4, 5, 2, 3, 0, 1 }; // opposite edge 00076 const int edge_beg_orth1[] = { 00077 3, 5, 1, 7, 1, 0, 3, 2, 7, 1, 5, 3 00078 }; // vtx adjacent to edge start in edge_dir[e]+1 direction 00079 const int edge_beg_orth2[] = { 00080 4, 0, 6, 2, 3, 2, 1, 0, 0, 4, 2, 6 00081 }; // vtx adjacent to edge start in edge_dir[e]+2 direction 00082 const int edge_end_orth1[] = { 00083 2, 6, 0, 4, 5, 4, 7, 6, 6, 2, 4, 0 00084 }; // vtx adjacent to edge end in edge_dir[e]+1 direction 00085 const int edge_end_orth2[] = { 00086 5, 3, 7, 1, 7, 6, 5, 4, 1, 7, 3, 5 00087 }; // vtx adjacent to edge end in edge_dir[e]+2 direction 00088 00089 static void coefficients_at_mid_edge( unsigned edge, double* coeff_out, size_t* indices_out, size_t& num_coeff ) 00090 { 00091 num_coeff = 2; 00092 coeff_out[0] = 0.5; 00093 coeff_out[1] = 0.5; 00094 indices_out[0] = edge_beg[edge]; 00095 indices_out[1] = edge_end[edge]; 00096 } 00097 00098 const int face_vtx[6][4] = { { 0, 1, 4, 5 }, // face 0 vertices 00099 { 1, 2, 5, 6 }, // face 1 vertices 00100 { 2, 3, 6, 7 }, // face 2 00101 { 0, 3, 4, 7 }, // face 3 00102 { 0, 1, 2, 3 }, // face 4 00103 { 4, 5, 6, 7 } }; // face 5 00104 const int face_opp[6] = { 2, 3, 0, 1, 5, 4 }; // opposite faces on hex 00105 const int face_dir[6] = { eta, xi, eta, xi, zeta, zeta }; // normal direction 00106 00107 static void coefficients_at_mid_face( unsigned face, double* coeff_out, size_t* indices_out, size_t& num_coeff ) 00108 { 00109 num_coeff = 4; 00110 coeff_out[0] = 0.25; 00111 coeff_out[1] = 0.25; 00112 coeff_out[2] = 0.25; 00113 coeff_out[3] = 0.25; 00114 indices_out[0] = face_vtx[face][0]; 00115 indices_out[1] = face_vtx[face][1]; 00116 indices_out[2] = face_vtx[face][2]; 00117 indices_out[3] = face_vtx[face][3]; 00118 } 00119 00120 static void coefficients_at_mid_elem( double* coeff_out, size_t* indices_out, size_t& num_coeff ) 00121 { 00122 num_coeff = 8; 00123 coeff_out[0] = 0.125; 00124 coeff_out[1] = 0.125; 00125 coeff_out[2] = 0.125; 00126 coeff_out[3] = 0.125; 00127 coeff_out[4] = 0.125; 00128 coeff_out[5] = 0.125; 00129 coeff_out[6] = 0.125; 00130 coeff_out[7] = 0.125; 00131 indices_out[0] = 0; 00132 indices_out[1] = 1; 00133 indices_out[2] = 2; 00134 indices_out[3] = 3; 00135 indices_out[4] = 4; 00136 indices_out[5] = 5; 00137 indices_out[6] = 6; 00138 indices_out[7] = 7; 00139 } 00140 00141 void LinearHexahedron::coefficients( Sample loc, NodeSet nodeset, double* coeff_out, size_t* indices_out, 00142 size_t& num_coeff, MsqError& err ) const 00143 { 00144 if( nodeset.have_any_mid_node() ) 00145 { 00146 MSQ_SETERR( err )( nonlinear_error, MsqError::UNSUPPORTED_ELEMENT ); 00147 return; 00148 } 00149 00150 switch( loc.dimension ) 00151 { 00152 case 0: 00153 coefficients_at_corner( loc.number, coeff_out, indices_out, num_coeff ); 00154 break; 00155 case 1: 00156 coefficients_at_mid_edge( loc.number, coeff_out, indices_out, num_coeff ); 00157 break; 00158 case 2: 00159 coefficients_at_mid_face( loc.number, coeff_out, indices_out, num_coeff ); 00160 break; 00161 case 3: 00162 coefficients_at_mid_elem( coeff_out, indices_out, num_coeff ); 00163 break; 00164 default: 00165 MSQ_SETERR( err )( "Invalid/unsupported logical dimension", MsqError::INVALID_ARG ); 00166 } 00167 } 00168 00169 static void derivatives_at_corner( unsigned corner, size_t* vertex_indices_out, MsqVector< 3 >* d_coeff_d_xi_out, 00170 size_t& num_vtx ) 00171 { 00172 const int xi_sign = coeff_xi_sign( corner ); 00173 const int eta_sign = coeff_eta_sign( corner ); 00174 const int zeta_sign = coeff_zeta_sign( corner ); 00175 const int offset = 4 * ( corner / 4 ); 00176 // const int adj_in_xi = offset + (9 - corner)%4; 00177 const int adj_in_xi = corner + 1 - 2 * ( corner % 2 ); 00178 const int adj_in_eta = 3 + 2 * offset - (int)corner; 00179 const int adj_in_zeta = corner - zeta_sign * 4; 00180 00181 num_vtx = 4; 00182 vertex_indices_out[0] = corner; 00183 vertex_indices_out[1] = adj_in_xi; 00184 vertex_indices_out[2] = adj_in_eta; 00185 vertex_indices_out[3] = adj_in_zeta; 00186 00187 d_coeff_d_xi_out[0][0] = xi_sign; 00188 d_coeff_d_xi_out[0][1] = eta_sign; 00189 d_coeff_d_xi_out[0][2] = zeta_sign; 00190 00191 d_coeff_d_xi_out[1][0] = -xi_sign; 00192 d_coeff_d_xi_out[1][1] = 0.0; 00193 d_coeff_d_xi_out[1][2] = 0.0; 00194 00195 d_coeff_d_xi_out[2][0] = 0.0; 00196 d_coeff_d_xi_out[2][1] = -eta_sign; 00197 d_coeff_d_xi_out[2][2] = 0.0; 00198 00199 d_coeff_d_xi_out[3][0] = 0.0; 00200 d_coeff_d_xi_out[3][1] = 0.0; 00201 d_coeff_d_xi_out[3][2] = -zeta_sign; 00202 } 00203 00204 static void derivatives_at_mid_edge( unsigned edge, size_t* vertex_indices_out, MsqVector< 3 >* d_coeff_d_xi_out, 00205 size_t& num_vtx ) 00206 { 00207 const int direction = edge_dir[edge]; 00208 const int ortho1 = ( direction + 1 ) % 3; 00209 const int ortho2 = ( direction + 2 ) % 3; 00210 const int vtx = edge_beg[edge]; 00211 const int signs[] = { coeff_xi_sign( vtx ), coeff_eta_sign( vtx ), coeff_zeta_sign( vtx ) }; 00212 const int sign_dir = signs[direction]; 00213 const int sign_or1 = signs[ortho1]; 00214 const int sign_or2 = signs[ortho2]; 00215 00216 num_vtx = 6; 00217 vertex_indices_out[0] = edge_beg[edge]; 00218 vertex_indices_out[1] = edge_end[edge]; 00219 vertex_indices_out[2] = edge_beg_orth1[edge]; 00220 vertex_indices_out[3] = edge_end_orth1[edge]; 00221 vertex_indices_out[4] = edge_beg_orth2[edge]; 00222 vertex_indices_out[5] = edge_end_orth2[edge]; 00223 00224 d_coeff_d_xi_out[0][direction] = sign_dir; 00225 d_coeff_d_xi_out[0][ortho1] = sign_or1 * 0.5; 00226 d_coeff_d_xi_out[0][ortho2] = sign_or2 * 0.5; 00227 00228 d_coeff_d_xi_out[1][direction] = -sign_dir; 00229 d_coeff_d_xi_out[1][ortho1] = sign_or1 * 0.5; 00230 d_coeff_d_xi_out[1][ortho2] = sign_or2 * 0.5; 00231 00232 d_coeff_d_xi_out[2][direction] = 0.0; 00233 d_coeff_d_xi_out[2][ortho1] = -sign_or1 * 0.5; 00234 d_coeff_d_xi_out[2][ortho2] = 0.0; 00235 00236 d_coeff_d_xi_out[3][direction] = 0.0; 00237 d_coeff_d_xi_out[3][ortho1] = -sign_or1 * 0.5; 00238 d_coeff_d_xi_out[3][ortho2] = 0.0; 00239 00240 d_coeff_d_xi_out[4][direction] = 0.0; 00241 d_coeff_d_xi_out[4][ortho1] = 0.0; 00242 d_coeff_d_xi_out[4][ortho2] = -sign_or2 * 0.5; 00243 00244 d_coeff_d_xi_out[5][direction] = 0.0; 00245 d_coeff_d_xi_out[5][ortho1] = 0.0; 00246 d_coeff_d_xi_out[5][ortho2] = -sign_or2 * 0.5; 00247 } 00248 00249 static void derivatives_at_mid_face( unsigned face, size_t* vertex_indices_out, MsqVector< 3 >* d_coeff_d_xi_out, 00250 size_t& num_vtx ) 00251 { 00252 const int vtx_signs[4][3] = { { coeff_xi_sign( face_vtx[face][0] ), coeff_eta_sign( face_vtx[face][0] ), 00253 coeff_zeta_sign( face_vtx[face][0] ) }, 00254 { coeff_xi_sign( face_vtx[face][1] ), coeff_eta_sign( face_vtx[face][1] ), 00255 coeff_zeta_sign( face_vtx[face][1] ) }, 00256 { coeff_xi_sign( face_vtx[face][2] ), coeff_eta_sign( face_vtx[face][2] ), 00257 coeff_zeta_sign( face_vtx[face][2] ) }, 00258 { coeff_xi_sign( face_vtx[face][3] ), coeff_eta_sign( face_vtx[face][3] ), 00259 coeff_zeta_sign( face_vtx[face][3] ) } }; 00260 const int ortho_dir = face_dir[face]; 00261 const int face_dir1 = ( ortho_dir + 1 ) % 3; 00262 const int face_dir2 = ( ortho_dir + 2 ) % 3; 00263 const int ortho_sign = vtx_signs[0][ortho_dir]; // same for all 4 verts 00264 00265 num_vtx = 8; 00266 vertex_indices_out[0] = face_vtx[face][0]; 00267 vertex_indices_out[1] = face_vtx[face][1]; 00268 vertex_indices_out[2] = face_vtx[face][2]; 00269 vertex_indices_out[3] = face_vtx[face][3]; 00270 vertex_indices_out[4] = face_vtx[face_opp[face]][0]; 00271 vertex_indices_out[5] = face_vtx[face_opp[face]][1]; 00272 vertex_indices_out[6] = face_vtx[face_opp[face]][2]; 00273 vertex_indices_out[7] = face_vtx[face_opp[face]][3]; 00274 00275 d_coeff_d_xi_out[0][ortho_dir] = ortho_sign * 0.25; 00276 d_coeff_d_xi_out[0][face_dir1] = vtx_signs[0][face_dir1] * 0.5; 00277 d_coeff_d_xi_out[0][face_dir2] = vtx_signs[0][face_dir2] * 0.5; 00278 00279 d_coeff_d_xi_out[1][ortho_dir] = ortho_sign * 0.25; 00280 d_coeff_d_xi_out[1][face_dir1] = vtx_signs[1][face_dir1] * 0.5; 00281 d_coeff_d_xi_out[1][face_dir2] = vtx_signs[1][face_dir2] * 0.5; 00282 00283 d_coeff_d_xi_out[2][ortho_dir] = ortho_sign * 0.25; 00284 d_coeff_d_xi_out[2][face_dir1] = vtx_signs[2][face_dir1] * 0.5; 00285 d_coeff_d_xi_out[2][face_dir2] = vtx_signs[2][face_dir2] * 0.5; 00286 00287 d_coeff_d_xi_out[3][ortho_dir] = ortho_sign * 0.25; 00288 d_coeff_d_xi_out[3][face_dir1] = vtx_signs[3][face_dir1] * 0.5; 00289 d_coeff_d_xi_out[3][face_dir2] = vtx_signs[3][face_dir2] * 0.5; 00290 00291 d_coeff_d_xi_out[4][ortho_dir] = -ortho_sign * 0.25; 00292 d_coeff_d_xi_out[4][face_dir1] = 0.0; 00293 d_coeff_d_xi_out[4][face_dir2] = 0.0; 00294 00295 d_coeff_d_xi_out[5][ortho_dir] = -ortho_sign * 0.25; 00296 d_coeff_d_xi_out[5][face_dir1] = 0.0; 00297 d_coeff_d_xi_out[5][face_dir2] = 0.0; 00298 00299 d_coeff_d_xi_out[6][ortho_dir] = -ortho_sign * 0.25; 00300 d_coeff_d_xi_out[6][face_dir1] = 0.0; 00301 d_coeff_d_xi_out[6][face_dir2] = 0.0; 00302 00303 d_coeff_d_xi_out[7][ortho_dir] = -ortho_sign * 0.25; 00304 d_coeff_d_xi_out[7][face_dir1] = 0.0; 00305 d_coeff_d_xi_out[7][face_dir2] = 0.0; 00306 } 00307 00308 static void derivatives_at_mid_elem( size_t* vertex_indices_out, MsqVector< 3 >* d_coeff_d_xi_out, size_t& num_vtx ) 00309 00310 { 00311 num_vtx = 8; 00312 vertex_indices_out[0] = 0; 00313 vertex_indices_out[1] = 1; 00314 vertex_indices_out[2] = 2; 00315 vertex_indices_out[3] = 3; 00316 vertex_indices_out[4] = 4; 00317 vertex_indices_out[5] = 5; 00318 vertex_indices_out[6] = 6; 00319 vertex_indices_out[7] = 7; 00320 00321 d_coeff_d_xi_out[0][0] = -0.25; 00322 d_coeff_d_xi_out[0][1] = -0.25; 00323 d_coeff_d_xi_out[0][2] = -0.25; 00324 00325 d_coeff_d_xi_out[1][0] = 0.25; 00326 d_coeff_d_xi_out[1][1] = -0.25; 00327 d_coeff_d_xi_out[1][2] = -0.25; 00328 00329 d_coeff_d_xi_out[2][0] = 0.25; 00330 d_coeff_d_xi_out[2][1] = 0.25; 00331 d_coeff_d_xi_out[2][2] = -0.25; 00332 00333 d_coeff_d_xi_out[3][0] = -0.25; 00334 d_coeff_d_xi_out[3][1] = 0.25; 00335 d_coeff_d_xi_out[3][2] = -0.25; 00336 00337 d_coeff_d_xi_out[4][0] = -0.25; 00338 d_coeff_d_xi_out[4][1] = -0.25; 00339 d_coeff_d_xi_out[4][2] = 0.25; 00340 00341 d_coeff_d_xi_out[5][0] = 0.25; 00342 d_coeff_d_xi_out[5][1] = -0.25; 00343 d_coeff_d_xi_out[5][2] = 0.25; 00344 00345 d_coeff_d_xi_out[6][0] = 0.25; 00346 d_coeff_d_xi_out[6][1] = 0.25; 00347 d_coeff_d_xi_out[6][2] = 0.25; 00348 00349 d_coeff_d_xi_out[7][0] = -0.25; 00350 d_coeff_d_xi_out[7][1] = 0.25; 00351 d_coeff_d_xi_out[7][2] = 0.25; 00352 } 00353 00354 void LinearHexahedron::derivatives( Sample loc, NodeSet nodeset, size_t* vertex_indices_out, 00355 MsqVector< 3 >* d_coeff_d_xi_out, size_t& num_vtx, MsqError& err ) const 00356 { 00357 if( nodeset.have_any_mid_node() ) 00358 { 00359 MSQ_SETERR( err )( nonlinear_error, MsqError::UNSUPPORTED_ELEMENT ); 00360 return; 00361 } 00362 00363 switch( loc.dimension ) 00364 { 00365 case 0: 00366 derivatives_at_corner( loc.number, vertex_indices_out, d_coeff_d_xi_out, num_vtx ); 00367 break; 00368 case 1: 00369 derivatives_at_mid_edge( loc.number, vertex_indices_out, d_coeff_d_xi_out, num_vtx ); 00370 break; 00371 case 2: 00372 derivatives_at_mid_face( loc.number, vertex_indices_out, d_coeff_d_xi_out, num_vtx ); 00373 break; 00374 case 3: 00375 derivatives_at_mid_elem( vertex_indices_out, d_coeff_d_xi_out, num_vtx ); 00376 break; 00377 default: 00378 MSQ_SETERR( err )( "Invalid/unsupported logical dimension", MsqError::INVALID_ARG ); 00379 } 00380 } 00381 00382 void LinearHexahedron::ideal( Sample, MsqMatrix< 3, 3 >& J, MsqError& ) const 00383 { 00384 J( 0, 0 ) = J( 1, 1 ) = J( 2, 2 ) = 1.0; 00385 J( 1, 0 ) = J( 0, 1 ) = J( 0, 2 ) = 0.0; 00386 J( 2, 0 ) = J( 2, 1 ) = J( 1, 2 ) = 0.0; 00387 } 00388 00389 } // namespace MBMesquite