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