MOAB: Mesh Oriented datABase  (version 5.4.0)
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,
00048                                       NodeSet nodeset,
00049                                       double* coeff_out,
00050                                       size_t* indices_out,
00051                                       size_t& num_coeff,
00052                                       MsqError& err ) const
00053 {
00054     if( !nodeset.have_any_mid_node() )
00055     {
00056         LinearQuadrilateral::coefficients( loc, nodeset, coeff_out, indices_out, num_coeff );
00057         return;
00058     }
00059 
00060     switch( loc.dimension )
00061     {
00062         case 0:
00063             num_coeff      = 1;
00064             indices_out[0] = loc.number;
00065             coeff_out[0]   = 1.0;
00066             break;
00067         case 1:
00068             coeff_out[0] = coeff_out[1] = coeff_out[2] = coeff_out[3] = coeff_out[4] = coeff_out[5] = coeff_out[6] =
00069                 coeff_out[7] = coeff_out[8] = 0.0;
00070             if( nodeset.mid_edge_node( loc.number ) )
00071             {
00072                 // if mid-edge node is present
00073                 num_coeff      = 1;
00074                 indices_out[0] = loc.number + 4;
00075                 coeff_out[0]   = 1.0;
00076             }
00077             else
00078             {
00079                 // If mid-edge node is not present, mapping function value
00080                 // for linear edge is even weight of adjacent vertices.
00081                 num_coeff      = 2;
00082                 indices_out[0] = loc.number;
00083                 indices_out[1] = ( loc.number + 1 ) % 4;
00084                 coeff_out[0]   = 0.5;
00085                 coeff_out[1]   = 0.5;
00086             }
00087             break;
00088         case 2:
00089             if( nodeset.mid_face_node( 0 ) )
00090             {  // if quad center node is present
00091                 num_coeff      = 1;
00092                 indices_out[0] = 8;
00093                 coeff_out[0]   = 1.0;
00094             }
00095             else
00096             {
00097                 // for linear element, (no mid-edge nodes), all corners contribute 1/4.
00098                 num_coeff      = 4;
00099                 indices_out[0] = 0;
00100                 indices_out[1] = 1;
00101                 indices_out[2] = 2;
00102                 indices_out[3] = 3;
00103                 coeff_out[0]   = 0.25;
00104                 coeff_out[1]   = 0.25;
00105                 coeff_out[2]   = 0.25;
00106                 coeff_out[3]   = 0.25;
00107                 // add in contribution for any mid-edge nodes present
00108                 for( int i = 0; i < 4; ++i )
00109                 {  // for each edge
00110                     if( nodeset.mid_edge_node( i ) )
00111                     {
00112                         indices_out[num_coeff] = i + 4;
00113                         coeff_out[num_coeff]   = 0.5;
00114                         coeff_out[i] -= 0.25;
00115                         coeff_out[( i + 1 ) % 4] -= 0.25;
00116                         ++num_coeff;
00117                     }
00118                 }
00119             }
00120             break;
00121         default:
00122             MSQ_SETERR( err )
00123             ( MsqError::UNSUPPORTED_ELEMENT,
00124               "Request for dimension %d mapping function value"
00125               "for a quadrilateral element",
00126               loc.dimension );
00127     }
00128 }
00129 
00130 static void derivatives_at_corner( unsigned corner,
00131                                    NodeSet nodeset,
00132                                    size_t* vertices,
00133                                    MsqVector< 2 >* derivs,
00134                                    size_t& num_vtx )
00135 {
00136     static const unsigned xi_adj_corners[]  = { 1, 0, 3, 2 };
00137     static const unsigned xi_adj_edges[]    = { 0, 0, 2, 2 };
00138     static const unsigned eta_adj_corners[] = { 3, 2, 1, 0 };
00139     static const unsigned eta_adj_edges[]   = { 3, 1, 1, 3 };
00140 
00141     static const double corner_xi[]  = { -1, 1, 1, -1 };  // xi values by corner
00142     static const double corner_eta[] = { -1, -1, 1, 1 };  // eta values by corner
00143     static const double other_xi[]   = { 1, -1, -1, 1 };  // xi values for adjacent corner in xi direction
00144     static const double other_eta[]  = { 1, 1, -1, -1 };  // eta values for adjcent corner in eta direction
00145     static const double mid_xi[]     = { 2, -2, -2, 2 };  // xi values for mid-node in xi direction
00146     static const double mid_eta[]    = { 2, 2, -2, -2 };  // eta values for mid-node in eta direction
00147 
00148     num_vtx     = 3;
00149     vertices[0] = corner;
00150     vertices[1] = xi_adj_corners[corner];
00151     vertices[2] = eta_adj_corners[corner];
00152 
00153     derivs[0][0] = corner_xi[corner];
00154     derivs[0][1] = corner_eta[corner];
00155     derivs[1][0] = other_xi[corner];
00156     derivs[1][1] = 0.0;
00157     derivs[2][0] = 0.0;
00158     derivs[2][1] = other_eta[corner];
00159 
00160     if( nodeset.mid_edge_node( xi_adj_edges[corner] ) )
00161     {
00162         vertices[num_vtx]  = 4 + xi_adj_edges[corner];
00163         derivs[num_vtx][0] = 2.0 * mid_xi[corner];
00164         derivs[num_vtx][1] = 0.0;
00165         derivs[0][0] -= mid_xi[corner];
00166         derivs[1][0] -= mid_xi[corner];
00167         ++num_vtx;
00168     }
00169 
00170     if( nodeset.mid_edge_node( eta_adj_edges[corner] ) )
00171     {
00172         vertices[num_vtx]  = 4 + eta_adj_edges[corner];
00173         derivs[num_vtx][0] = 0.0;
00174         derivs[num_vtx][1] = 2.0 * mid_eta[corner];
00175         derivs[0][1] -= mid_eta[corner];
00176         derivs[2][1] -= mid_eta[corner];
00177         ++num_vtx;
00178     }
00179 }
00180 
00181 static void derivatives_at_mid_edge( unsigned edge,
00182                                      NodeSet nodeset,
00183                                      size_t* vertices,
00184                                      MsqVector< 2 >* derivs,
00185                                      size_t& num_vtx )
00186 {
00187     static const double values[][9]      = { { -0.5, -0.5, 0.5, 0.5, -1.0, 2.0, 1.0, 2.0, 4.0 },
00188                                         { -0.5, 0.5, 0.5, -0.5, -2.0, 1.0, -2.0, -1.0, -4.0 },
00189                                         { -0.5, -0.5, 0.5, 0.5, -1.0, -2.0, 1.0, -2.0, -4.0 },
00190                                         { -0.5, 0.5, 0.5, -0.5, 2.0, 1.0, 2.0, -1.0, 4.0 } };
00191     static const double edge_values[][2] = { { -1, 1 }, { -1, 1 }, { 1, -1 }, { 1, -1 } };
00192     const unsigned prev_corner           = edge;              // index of start vertex of edge
00193     const unsigned next_corner           = ( edge + 1 ) % 4;  // index of end vertex of edge
00194     const unsigned is_eta_edge           = edge % 2;          // edge is xi = +/- 0
00195     const unsigned is_xi_edge            = 1 - is_eta_edge;   // edge is eta = +/- 0
00196     // const unsigned mid_edge_node = edge + 4;     // mid-edge node index
00197     const unsigned prev_opposite = ( prev_corner + 3 ) % 4;  // index of corner adjacent to prev_corner
00198     const unsigned next_opposite = ( next_corner + 1 ) % 4;  // index of corner adjacent to next_corner;
00199 
00200     // First do derivatives along edge (e.g. wrt xi if edge is eta = +/-1)
00201     num_vtx                = 2;
00202     vertices[0]            = prev_corner;
00203     vertices[1]            = next_corner;
00204     derivs[0][is_eta_edge] = edge_values[edge][0];
00205     derivs[0][is_xi_edge]  = 0.0;
00206     derivs[1][is_eta_edge] = edge_values[edge][1];
00207     derivs[1][is_xi_edge]  = 0.0;
00208     // That's it for the edge-direction derivatives.  No other vertices contribute.
00209 
00210     // Next handle the linear element case.  Handle this as a special case first,
00211     // so the generalized solution doesn't impact performance for linear elements
00212     // too much.
00213     if( !nodeset.have_any_mid_node() )
00214     {
00215         num_vtx                = 4;
00216         vertices[2]            = prev_opposite;
00217         vertices[3]            = next_opposite;
00218         derivs[0][is_xi_edge]  = values[edge][prev_corner];
00219         derivs[1][is_xi_edge]  = values[edge][next_corner];
00220         derivs[2][is_xi_edge]  = values[edge][prev_opposite];
00221         derivs[2][is_eta_edge] = 0.0;
00222         derivs[3][is_xi_edge]  = values[edge][next_opposite];
00223         derivs[3][is_eta_edge] = 0.0;
00224         return;
00225     }
00226 
00227     // Initial (linear) contribution for each corner
00228     double v[4] = { values[edge][0], values[edge][1], values[edge][2], values[edge][3] };
00229 
00230     // If mid-face node is present
00231     double v8 = 0.0;
00232     if( nodeset.mid_face_node( 0 ) )
00233     {
00234         v8                           = values[edge][8];
00235         vertices[num_vtx]            = 8;
00236         derivs[num_vtx][is_eta_edge] = 0.0;
00237         derivs[num_vtx][is_xi_edge]  = v8;
00238         v[0] -= 0.25 * v8;
00239         v[1] -= 0.25 * v8;
00240         v[2] -= 0.25 * v8;
00241         v[3] -= 0.25 * v8;
00242         ++num_vtx;
00243     }
00244 
00245     // If mid-edge nodes are present
00246     for( unsigned i = 0; i < 4; ++i )
00247     {
00248         if( nodeset.mid_edge_node( i ) )
00249         {
00250             const double value = values[edge][i + 4] - 0.5 * v8;
00251             if( fabs( value ) > 0.125 )
00252             {
00253                 v[i] -= 0.5 * value;
00254                 v[( i + 1 ) % 4] -= 0.5 * value;
00255                 vertices[num_vtx]            = i + 4;
00256                 derivs[num_vtx][is_eta_edge] = 0.0;
00257                 derivs[num_vtx][is_xi_edge]  = value;
00258                 ++num_vtx;
00259             }
00260         }
00261     }
00262 
00263     // update values for adjacent corners
00264     derivs[0][is_xi_edge] = v[prev_corner];
00265     derivs[1][is_xi_edge] = v[next_corner];
00266     // do other two corners
00267     if( fabs( v[prev_opposite] ) > 0.125 )
00268     {
00269         vertices[num_vtx]            = prev_opposite;
00270         derivs[num_vtx][is_eta_edge] = 0.0;
00271         derivs[num_vtx][is_xi_edge]  = v[prev_opposite];
00272         ++num_vtx;
00273     }
00274     if( fabs( v[next_opposite] ) > 0.125 )
00275     {
00276         vertices[num_vtx]            = next_opposite;
00277         derivs[num_vtx][is_eta_edge] = 0.0;
00278         derivs[num_vtx][is_xi_edge]  = v[next_opposite];
00279         ++num_vtx;
00280     }
00281 }
00282 
00283 static void derivatives_at_mid_elem( NodeSet nodeset, size_t* vertices, MsqVector< 2 >* derivs, size_t& num_vtx )
00284 {
00285     // fast linear case
00286     // This is provided as an optimization for linear elements.
00287     // If this block of code were removed, the general-case code
00288     // below should produce the same result.
00289     if( !nodeset.have_any_mid_node() )
00290     {
00291         num_vtx      = 4;
00292         vertices[0]  = 0;
00293         derivs[0][0] = -0.5;
00294         derivs[0][1] = -0.5;
00295         vertices[1]  = 1;
00296         derivs[1][0] = 0.5;
00297         derivs[1][1] = -0.5;
00298         vertices[2]  = 2;
00299         derivs[2][0] = 0.5;
00300         derivs[2][1] = 0.5;
00301         vertices[3]  = 3;
00302         derivs[3][0] = -0.5;
00303         derivs[3][1] = 0.5;
00304         return;
00305     }
00306 
00307     num_vtx = 0;
00308 
00309     // N_0
00310     if( !nodeset.both_edge_nodes( 0, 3 ) )
00311     {  // if eiter adjacent mid-edge node is missing
00312         vertices[num_vtx]  = 0;
00313         derivs[num_vtx][0] = nodeset.mid_edge_node( 3 ) ? 0.0 : -0.5;
00314         derivs[num_vtx][1] = nodeset.mid_edge_node( 0 ) ? 0.0 : -0.5;
00315         ++num_vtx;
00316     }
00317 
00318     // N_1
00319     if( !nodeset.both_edge_nodes( 0, 1 ) )
00320     {  // if eiter adjacent mid-edge node is missing
00321         vertices[num_vtx]  = 1;
00322         derivs[num_vtx][0] = nodeset.mid_edge_node( 1 ) ? 0.0 : 0.5;
00323         derivs[num_vtx][1] = nodeset.mid_edge_node( 0 ) ? 0.0 : -0.5;
00324         ++num_vtx;
00325     }
00326 
00327     // N_2
00328     if( !nodeset.both_edge_nodes( 1, 2 ) )
00329     {  // if eiter adjacent mid-edge node is missing
00330         vertices[num_vtx]  = 2;
00331         derivs[num_vtx][0] = nodeset.mid_edge_node( 1 ) ? 0.0 : 0.5;
00332         derivs[num_vtx][1] = nodeset.mid_edge_node( 2 ) ? 0.0 : 0.5;
00333         ++num_vtx;
00334     }
00335 
00336     // N_3
00337     if( !nodeset.both_edge_nodes( 2, 3 ) )
00338     {  // if eiter adjacent mid-edge node is missing
00339         vertices[num_vtx]  = 3;
00340         derivs[num_vtx][0] = nodeset.mid_edge_node( 3 ) ? 0.0 : -0.5;
00341         derivs[num_vtx][1] = nodeset.mid_edge_node( 2 ) ? 0.0 : 0.5;
00342         ++num_vtx;
00343     }
00344 
00345     // N_4
00346     if( nodeset.mid_edge_node( 0 ) )
00347     {
00348         vertices[num_vtx]  = 4;
00349         derivs[num_vtx][0] = 0.0;
00350         derivs[num_vtx][1] = -1.0;
00351         ++num_vtx;
00352     }
00353 
00354     // N_5
00355     if( nodeset.mid_edge_node( 1 ) )
00356     {
00357         vertices[num_vtx]  = 5;
00358         derivs[num_vtx][0] = 1.0;
00359         derivs[num_vtx][1] = 0.0;
00360         ++num_vtx;
00361     }
00362 
00363     // N_6
00364     if( nodeset.mid_edge_node( 2 ) )
00365     {
00366         vertices[num_vtx]  = 6;
00367         derivs[num_vtx][0] = 0.0;
00368         derivs[num_vtx][1] = 1.0;
00369         ++num_vtx;
00370     }
00371 
00372     // N_7
00373     if( nodeset.mid_edge_node( 3 ) )
00374     {
00375         vertices[num_vtx]  = 7;
00376         derivs[num_vtx][0] = -1.0;
00377         derivs[num_vtx][1] = 0.0;
00378         ++num_vtx;
00379     }
00380 
00381     // N_8 (mid-quad node) never contributes to Jacobian at element center!!!
00382 }
00383 
00384 void QuadLagrangeShape::derivatives( Sample loc,
00385                                      NodeSet nodeset,
00386                                      size_t* vertex_indices_out,
00387                                      MsqVector< 2 >* d_coeff_d_xi_out,
00388                                      size_t& num_vtx,
00389                                      MsqError& err ) const
00390 {
00391     if( !nodeset.have_any_mid_node() )
00392     {
00393         LinearQuadrilateral::derivatives( loc, nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00394         return;
00395     }
00396 
00397     switch( loc.dimension )
00398     {
00399         case 0:
00400             derivatives_at_corner( loc.number, nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00401             break;
00402         case 1:
00403             derivatives_at_mid_edge( loc.number, nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00404             break;
00405         case 2:
00406             derivatives_at_mid_elem( nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00407             break;
00408         default:
00409             MSQ_SETERR( err )( "Invalid/unsupported logical dimension", MsqError::INVALID_ARG );
00410     }
00411 }
00412 
00413 void QuadLagrangeShape::ideal( Sample, MsqMatrix< 3, 2 >& J, MsqError& ) const
00414 {
00415     J( 0, 0 ) = J( 1, 1 ) = 1.0;
00416     J( 0, 1 ) = J( 1, 0 ) = 0.0;
00417     J( 2, 0 ) = J( 2, 1 ) = 0.0;
00418 }
00419 
00420 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines