MOAB: Mesh Oriented datABase  (version 5.2.1)
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[] = {
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines