MOAB: Mesh Oriented datABase  (version 5.4.0)
LinearHexahedron.cpp
Go to the documentation of this file.
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[] = { 3, 5, 1, 7, 1, 0,
00077                                3, 2, 7, 1, 5, 3 };  // vtx adjacent to edge start in edge_dir[e]+1 direction
00078 const int edge_beg_orth2[] = { 4, 0, 6, 2, 3, 2,
00079                                1, 0, 0, 4, 2, 6 };  // vtx adjacent to edge start in edge_dir[e]+2 direction
00080 const int edge_end_orth1[] = { 2, 6, 0, 4, 5, 4,
00081                                7, 6, 6, 2, 4, 0 };  // vtx adjacent to edge end in edge_dir[e]+1 direction
00082 const int edge_end_orth2[] = { 5, 3, 7, 1, 7, 6,
00083                                5, 4, 1, 7, 3, 5 };  // vtx adjacent to edge end in edge_dir[e]+2 direction
00084 
00085 static void coefficients_at_mid_edge( unsigned edge, double* coeff_out, size_t* indices_out, size_t& num_coeff )
00086 {
00087     num_coeff      = 2;
00088     coeff_out[0]   = 0.5;
00089     coeff_out[1]   = 0.5;
00090     indices_out[0] = edge_beg[edge];
00091     indices_out[1] = edge_end[edge];
00092 }
00093 
00094 const int face_vtx[6][4] = { { 0, 1, 4, 5 },                  // face 0 vertices
00095                              { 1, 2, 5, 6 },                  // face 1 vertices
00096                              { 2, 3, 6, 7 },                  // face 2
00097                              { 0, 3, 4, 7 },                  // face 3
00098                              { 0, 1, 2, 3 },                  // face 4
00099                              { 4, 5, 6, 7 } };                // face 5
00100 const int face_opp[6]    = { 2, 3, 0, 1, 5, 4 };              // opposite faces on hex
00101 const int face_dir[6]    = { eta, xi, eta, xi, zeta, zeta };  // normal direction
00102 
00103 static void coefficients_at_mid_face( unsigned face, double* coeff_out, size_t* indices_out, size_t& num_coeff )
00104 {
00105     num_coeff      = 4;
00106     coeff_out[0]   = 0.25;
00107     coeff_out[1]   = 0.25;
00108     coeff_out[2]   = 0.25;
00109     coeff_out[3]   = 0.25;
00110     indices_out[0] = face_vtx[face][0];
00111     indices_out[1] = face_vtx[face][1];
00112     indices_out[2] = face_vtx[face][2];
00113     indices_out[3] = face_vtx[face][3];
00114 }
00115 
00116 static void coefficients_at_mid_elem( double* coeff_out, size_t* indices_out, size_t& num_coeff )
00117 {
00118     num_coeff      = 8;
00119     coeff_out[0]   = 0.125;
00120     coeff_out[1]   = 0.125;
00121     coeff_out[2]   = 0.125;
00122     coeff_out[3]   = 0.125;
00123     coeff_out[4]   = 0.125;
00124     coeff_out[5]   = 0.125;
00125     coeff_out[6]   = 0.125;
00126     coeff_out[7]   = 0.125;
00127     indices_out[0] = 0;
00128     indices_out[1] = 1;
00129     indices_out[2] = 2;
00130     indices_out[3] = 3;
00131     indices_out[4] = 4;
00132     indices_out[5] = 5;
00133     indices_out[6] = 6;
00134     indices_out[7] = 7;
00135 }
00136 
00137 void LinearHexahedron::coefficients( Sample loc,
00138                                      NodeSet nodeset,
00139                                      double* coeff_out,
00140                                      size_t* indices_out,
00141                                      size_t& num_coeff,
00142                                      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,
00170                                    size_t* vertex_indices_out,
00171                                    MsqVector< 3 >* d_coeff_d_xi_out,
00172                                    size_t& num_vtx )
00173 {
00174     const int xi_sign   = coeff_xi_sign( corner );
00175     const int eta_sign  = coeff_eta_sign( corner );
00176     const int zeta_sign = coeff_zeta_sign( corner );
00177     const int offset    = 4 * ( corner / 4 );
00178     // const int adj_in_xi   = offset + (9 - corner)%4;
00179     const int adj_in_xi   = corner + 1 - 2 * ( corner % 2 );
00180     const int adj_in_eta  = 3 + 2 * offset - (int)corner;
00181     const int adj_in_zeta = corner - zeta_sign * 4;
00182 
00183     num_vtx               = 4;
00184     vertex_indices_out[0] = corner;
00185     vertex_indices_out[1] = adj_in_xi;
00186     vertex_indices_out[2] = adj_in_eta;
00187     vertex_indices_out[3] = adj_in_zeta;
00188 
00189     d_coeff_d_xi_out[0][0] = xi_sign;
00190     d_coeff_d_xi_out[0][1] = eta_sign;
00191     d_coeff_d_xi_out[0][2] = zeta_sign;
00192 
00193     d_coeff_d_xi_out[1][0] = -xi_sign;
00194     d_coeff_d_xi_out[1][1] = 0.0;
00195     d_coeff_d_xi_out[1][2] = 0.0;
00196 
00197     d_coeff_d_xi_out[2][0] = 0.0;
00198     d_coeff_d_xi_out[2][1] = -eta_sign;
00199     d_coeff_d_xi_out[2][2] = 0.0;
00200 
00201     d_coeff_d_xi_out[3][0] = 0.0;
00202     d_coeff_d_xi_out[3][1] = 0.0;
00203     d_coeff_d_xi_out[3][2] = -zeta_sign;
00204 }
00205 
00206 static void derivatives_at_mid_edge( unsigned edge,
00207                                      size_t* vertex_indices_out,
00208                                      MsqVector< 3 >* d_coeff_d_xi_out,
00209                                      size_t& num_vtx )
00210 {
00211     const int direction = edge_dir[edge];
00212     const int ortho1    = ( direction + 1 ) % 3;
00213     const int ortho2    = ( direction + 2 ) % 3;
00214     const int vtx       = edge_beg[edge];
00215     const int signs[]   = { coeff_xi_sign( vtx ), coeff_eta_sign( vtx ), coeff_zeta_sign( vtx ) };
00216     const int sign_dir  = signs[direction];
00217     const int sign_or1  = signs[ortho1];
00218     const int sign_or2  = signs[ortho2];
00219 
00220     num_vtx               = 6;
00221     vertex_indices_out[0] = edge_beg[edge];
00222     vertex_indices_out[1] = edge_end[edge];
00223     vertex_indices_out[2] = edge_beg_orth1[edge];
00224     vertex_indices_out[3] = edge_end_orth1[edge];
00225     vertex_indices_out[4] = edge_beg_orth2[edge];
00226     vertex_indices_out[5] = edge_end_orth2[edge];
00227 
00228     d_coeff_d_xi_out[0][direction] = sign_dir;
00229     d_coeff_d_xi_out[0][ortho1]    = sign_or1 * 0.5;
00230     d_coeff_d_xi_out[0][ortho2]    = sign_or2 * 0.5;
00231 
00232     d_coeff_d_xi_out[1][direction] = -sign_dir;
00233     d_coeff_d_xi_out[1][ortho1]    = sign_or1 * 0.5;
00234     d_coeff_d_xi_out[1][ortho2]    = sign_or2 * 0.5;
00235 
00236     d_coeff_d_xi_out[2][direction] = 0.0;
00237     d_coeff_d_xi_out[2][ortho1]    = -sign_or1 * 0.5;
00238     d_coeff_d_xi_out[2][ortho2]    = 0.0;
00239 
00240     d_coeff_d_xi_out[3][direction] = 0.0;
00241     d_coeff_d_xi_out[3][ortho1]    = -sign_or1 * 0.5;
00242     d_coeff_d_xi_out[3][ortho2]    = 0.0;
00243 
00244     d_coeff_d_xi_out[4][direction] = 0.0;
00245     d_coeff_d_xi_out[4][ortho1]    = 0.0;
00246     d_coeff_d_xi_out[4][ortho2]    = -sign_or2 * 0.5;
00247 
00248     d_coeff_d_xi_out[5][direction] = 0.0;
00249     d_coeff_d_xi_out[5][ortho1]    = 0.0;
00250     d_coeff_d_xi_out[5][ortho2]    = -sign_or2 * 0.5;
00251 }
00252 
00253 static void derivatives_at_mid_face( unsigned face,
00254                                      size_t* vertex_indices_out,
00255                                      MsqVector< 3 >* d_coeff_d_xi_out,
00256                                      size_t& num_vtx )
00257 {
00258     const int vtx_signs[4][3] = { { coeff_xi_sign( face_vtx[face][0] ), coeff_eta_sign( face_vtx[face][0] ),
00259                                     coeff_zeta_sign( face_vtx[face][0] ) },
00260                                   { coeff_xi_sign( face_vtx[face][1] ), coeff_eta_sign( face_vtx[face][1] ),
00261                                     coeff_zeta_sign( face_vtx[face][1] ) },
00262                                   { coeff_xi_sign( face_vtx[face][2] ), coeff_eta_sign( face_vtx[face][2] ),
00263                                     coeff_zeta_sign( face_vtx[face][2] ) },
00264                                   { coeff_xi_sign( face_vtx[face][3] ), coeff_eta_sign( face_vtx[face][3] ),
00265                                     coeff_zeta_sign( face_vtx[face][3] ) } };
00266     const int ortho_dir       = face_dir[face];
00267     const int face_dir1       = ( ortho_dir + 1 ) % 3;
00268     const int face_dir2       = ( ortho_dir + 2 ) % 3;
00269     const int ortho_sign      = vtx_signs[0][ortho_dir];  // same for all 4 verts
00270 
00271     num_vtx               = 8;
00272     vertex_indices_out[0] = face_vtx[face][0];
00273     vertex_indices_out[1] = face_vtx[face][1];
00274     vertex_indices_out[2] = face_vtx[face][2];
00275     vertex_indices_out[3] = face_vtx[face][3];
00276     vertex_indices_out[4] = face_vtx[face_opp[face]][0];
00277     vertex_indices_out[5] = face_vtx[face_opp[face]][1];
00278     vertex_indices_out[6] = face_vtx[face_opp[face]][2];
00279     vertex_indices_out[7] = face_vtx[face_opp[face]][3];
00280 
00281     d_coeff_d_xi_out[0][ortho_dir] = ortho_sign * 0.25;
00282     d_coeff_d_xi_out[0][face_dir1] = vtx_signs[0][face_dir1] * 0.5;
00283     d_coeff_d_xi_out[0][face_dir2] = vtx_signs[0][face_dir2] * 0.5;
00284 
00285     d_coeff_d_xi_out[1][ortho_dir] = ortho_sign * 0.25;
00286     d_coeff_d_xi_out[1][face_dir1] = vtx_signs[1][face_dir1] * 0.5;
00287     d_coeff_d_xi_out[1][face_dir2] = vtx_signs[1][face_dir2] * 0.5;
00288 
00289     d_coeff_d_xi_out[2][ortho_dir] = ortho_sign * 0.25;
00290     d_coeff_d_xi_out[2][face_dir1] = vtx_signs[2][face_dir1] * 0.5;
00291     d_coeff_d_xi_out[2][face_dir2] = vtx_signs[2][face_dir2] * 0.5;
00292 
00293     d_coeff_d_xi_out[3][ortho_dir] = ortho_sign * 0.25;
00294     d_coeff_d_xi_out[3][face_dir1] = vtx_signs[3][face_dir1] * 0.5;
00295     d_coeff_d_xi_out[3][face_dir2] = vtx_signs[3][face_dir2] * 0.5;
00296 
00297     d_coeff_d_xi_out[4][ortho_dir] = -ortho_sign * 0.25;
00298     d_coeff_d_xi_out[4][face_dir1] = 0.0;
00299     d_coeff_d_xi_out[4][face_dir2] = 0.0;
00300 
00301     d_coeff_d_xi_out[5][ortho_dir] = -ortho_sign * 0.25;
00302     d_coeff_d_xi_out[5][face_dir1] = 0.0;
00303     d_coeff_d_xi_out[5][face_dir2] = 0.0;
00304 
00305     d_coeff_d_xi_out[6][ortho_dir] = -ortho_sign * 0.25;
00306     d_coeff_d_xi_out[6][face_dir1] = 0.0;
00307     d_coeff_d_xi_out[6][face_dir2] = 0.0;
00308 
00309     d_coeff_d_xi_out[7][ortho_dir] = -ortho_sign * 0.25;
00310     d_coeff_d_xi_out[7][face_dir1] = 0.0;
00311     d_coeff_d_xi_out[7][face_dir2] = 0.0;
00312 }
00313 
00314 static void derivatives_at_mid_elem( size_t* vertex_indices_out, MsqVector< 3 >* d_coeff_d_xi_out, size_t& num_vtx )
00315 
00316 {
00317     num_vtx               = 8;
00318     vertex_indices_out[0] = 0;
00319     vertex_indices_out[1] = 1;
00320     vertex_indices_out[2] = 2;
00321     vertex_indices_out[3] = 3;
00322     vertex_indices_out[4] = 4;
00323     vertex_indices_out[5] = 5;
00324     vertex_indices_out[6] = 6;
00325     vertex_indices_out[7] = 7;
00326 
00327     d_coeff_d_xi_out[0][0] = -0.25;
00328     d_coeff_d_xi_out[0][1] = -0.25;
00329     d_coeff_d_xi_out[0][2] = -0.25;
00330 
00331     d_coeff_d_xi_out[1][0] = 0.25;
00332     d_coeff_d_xi_out[1][1] = -0.25;
00333     d_coeff_d_xi_out[1][2] = -0.25;
00334 
00335     d_coeff_d_xi_out[2][0] = 0.25;
00336     d_coeff_d_xi_out[2][1] = 0.25;
00337     d_coeff_d_xi_out[2][2] = -0.25;
00338 
00339     d_coeff_d_xi_out[3][0] = -0.25;
00340     d_coeff_d_xi_out[3][1] = 0.25;
00341     d_coeff_d_xi_out[3][2] = -0.25;
00342 
00343     d_coeff_d_xi_out[4][0] = -0.25;
00344     d_coeff_d_xi_out[4][1] = -0.25;
00345     d_coeff_d_xi_out[4][2] = 0.25;
00346 
00347     d_coeff_d_xi_out[5][0] = 0.25;
00348     d_coeff_d_xi_out[5][1] = -0.25;
00349     d_coeff_d_xi_out[5][2] = 0.25;
00350 
00351     d_coeff_d_xi_out[6][0] = 0.25;
00352     d_coeff_d_xi_out[6][1] = 0.25;
00353     d_coeff_d_xi_out[6][2] = 0.25;
00354 
00355     d_coeff_d_xi_out[7][0] = -0.25;
00356     d_coeff_d_xi_out[7][1] = 0.25;
00357     d_coeff_d_xi_out[7][2] = 0.25;
00358 }
00359 
00360 void LinearHexahedron::derivatives( Sample loc,
00361                                     NodeSet nodeset,
00362                                     size_t* vertex_indices_out,
00363                                     MsqVector< 3 >* d_coeff_d_xi_out,
00364                                     size_t& num_vtx,
00365                                     MsqError& err ) const
00366 {
00367     if( nodeset.have_any_mid_node() )
00368     {
00369         MSQ_SETERR( err )( nonlinear_error, MsqError::UNSUPPORTED_ELEMENT );
00370         return;
00371     }
00372 
00373     switch( loc.dimension )
00374     {
00375         case 0:
00376             derivatives_at_corner( loc.number, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00377             break;
00378         case 1:
00379             derivatives_at_mid_edge( loc.number, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00380             break;
00381         case 2:
00382             derivatives_at_mid_face( loc.number, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00383             break;
00384         case 3:
00385             derivatives_at_mid_elem( vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00386             break;
00387         default:
00388             MSQ_SETERR( err )( "Invalid/unsupported logical dimension", MsqError::INVALID_ARG );
00389     }
00390 }
00391 
00392 void LinearHexahedron::ideal( Sample, MsqMatrix< 3, 3 >& J, MsqError& ) const
00393 {
00394     J( 0, 0 ) = J( 1, 1 ) = J( 2, 2 ) = 1.0;
00395     J( 1, 0 ) = J( 0, 1 ) = J( 0, 2 ) = 0.0;
00396     J( 2, 0 ) = J( 2, 1 ) = J( 1, 2 ) = 0.0;
00397 }
00398 
00399 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines