MOAB: Mesh Oriented datABase  (version 5.2.1)
LinearPyramid.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 LinearPyramid.cpp
00028  *  \brief LinearPyramid implementation
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "LinearPyramid.hpp"
00034 #include "MsqError.hpp"
00035 
00036 namespace MBMesquite
00037 {
00038 
00039 static const char* nonlinear_error = "Attempt to use LinearTriangle mapping function for a nonlinear element\n";
00040 
00041 static inline void set_equal_derivatives( double value, size_t* indices, MsqVector< 3 >* derivs, size_t& num_vtx )
00042 {
00043     num_vtx    = 5;
00044     indices[0] = 0;
00045     indices[1] = 1;
00046     indices[2] = 2;
00047     indices[3] = 3;
00048     indices[4] = 4;
00049 
00050     derivs[0][0] = -value;
00051     derivs[0][1] = -value;
00052     derivs[0][2] = -0.25;
00053 
00054     derivs[1][0] = value;
00055     derivs[1][1] = -value;
00056     derivs[1][2] = -0.25;
00057 
00058     derivs[2][0] = value;
00059     derivs[2][1] = value;
00060     derivs[2][2] = -0.25;
00061 
00062     derivs[3][0] = -value;
00063     derivs[3][1] = value;
00064     derivs[3][2] = -0.25;
00065 
00066     derivs[4][0] = 0.0;
00067     derivs[4][1] = 0.0;
00068     derivs[4][2] = 1.0;
00069 }
00070 
00071 static inline void set_edge_derivatives( unsigned base_corner, double value, size_t* indices, MsqVector< 3 >* derivs,
00072                                          size_t& num_vtx )
00073 {
00074     const int direction = base_corner % 2;
00075     const int edge_beg  = base_corner;
00076     const int edge_end  = ( base_corner + 1 ) % 4;
00077     const int adj_end   = ( base_corner + 2 ) % 4;
00078     const int adj_beg   = ( base_corner + 3 ) % 4;
00079     const int dir_sign  = 2 * ( edge_beg / 2 ) - 1;
00080     const int oth_sign  = 2 * ( ( edge_beg + 1 ) / 2 % 2 ) - 1;
00081 
00082     num_vtx    = 5;
00083     indices[0] = edge_beg;
00084     indices[1] = edge_end;
00085     indices[2] = adj_end;
00086     indices[3] = adj_beg;
00087     indices[4] = 4;
00088 
00089     derivs[0][direction]     = 2 * dir_sign * value;
00090     derivs[0][1 - direction] = oth_sign * value;
00091     derivs[0][2]             = -0.5;
00092 
00093     derivs[1][direction]     = -2 * dir_sign * value;
00094     derivs[1][1 - direction] = oth_sign * value;
00095     derivs[1][2]             = -0.5;
00096 
00097     derivs[2][direction]     = 0.0;
00098     derivs[2][1 - direction] = -oth_sign * value;
00099     derivs[2][2]             = 0.0;
00100 
00101     derivs[3][direction]     = 0.0;
00102     derivs[3][1 - direction] = -oth_sign * value;
00103     derivs[3][2]             = 0.0;
00104 
00105     derivs[4][0] = 0.0;
00106     derivs[4][1] = 0.0;
00107     derivs[4][2] = 1.0;
00108 }
00109 
00110 static inline void set_corner_derivatives( unsigned corner, double value, size_t* indices, MsqVector< 3 >* derivs,
00111                                            size_t& num_vtx )
00112 {
00113     const unsigned adj_in_xi  = ( 5 - corner ) % 4;
00114     const unsigned adj_in_eta = 3 - corner;
00115 
00116     const int dxi_sign      = 2 * ( ( corner + 1 ) / 2 % 2 ) - 1;
00117     const int deta_sign     = 2 * ( corner / 2 ) - 1;
00118     const double dxi_value  = dxi_sign * value;
00119     const double deta_value = deta_sign * value;
00120 
00121     num_vtx    = 4;
00122     indices[0] = corner;
00123     indices[1] = adj_in_xi;
00124     indices[2] = adj_in_eta;
00125     indices[3] = 4;
00126 
00127     derivs[0][0] = dxi_value;
00128     derivs[0][1] = deta_value;
00129     derivs[0][2] = -1.0;
00130 
00131     derivs[1][0] = -dxi_value;
00132     derivs[1][1] = 0.0;
00133     derivs[1][2] = 0.0;
00134 
00135     derivs[2][0] = 0.0;
00136     derivs[2][1] = -deta_value;
00137     derivs[2][2] = 0.0;
00138 
00139     derivs[3][0] = 0.0;
00140     derivs[3][1] = 0.0;
00141     derivs[3][2] = 1.0;
00142 }
00143 
00144 EntityTopology LinearPyramid::element_topology() const
00145 {
00146     return PYRAMID;
00147 }
00148 
00149 int LinearPyramid::num_nodes() const
00150 {
00151     return 5;
00152 }
00153 
00154 NodeSet LinearPyramid::sample_points( NodeSet ) const
00155 {
00156     NodeSet result;
00157     result.set_all_corner_nodes( PYRAMID );
00158     result.clear_corner_node( 4 );
00159     return result;
00160 }
00161 
00162 static void coefficients_at_corner( unsigned corner, double* coeff_out, size_t* indices_out, size_t& num_coeff )
00163 {
00164     num_coeff      = 1;
00165     indices_out[0] = corner;
00166     coeff_out[0]   = 1.0;
00167 }
00168 
00169 static void coefficients_at_mid_edge( unsigned edge, double* coeff_out, size_t* indices_out, size_t& num_coeff )
00170 {
00171     num_coeff    = 2;
00172     coeff_out[0] = 0.5;
00173     coeff_out[1] = 0.5;
00174 
00175     if( edge < 4 )
00176     {
00177         indices_out[0] = edge;
00178         indices_out[1] = ( edge + 1 ) % 4;
00179     }
00180     else
00181     {
00182         indices_out[0] = edge - 4;
00183         indices_out[1] = 4;
00184     }
00185 }
00186 
00187 static void coefficients_at_mid_face( unsigned face, double* coeff_out, size_t* indices_out, size_t& num_coeff )
00188 {
00189     if( face == 4 )
00190     {
00191         num_coeff      = 4;
00192         coeff_out[0]   = 0.25;
00193         coeff_out[1]   = 0.25;
00194         coeff_out[2]   = 0.25;
00195         coeff_out[3]   = 0.25;
00196         indices_out[0] = 0;
00197         indices_out[1] = 1;
00198         indices_out[2] = 2;
00199         indices_out[3] = 3;
00200     }
00201     else
00202     {
00203         num_coeff      = 3;
00204         indices_out[0] = face;
00205         indices_out[1] = ( face + 1 ) % 4;
00206         indices_out[2] = 4;
00207         coeff_out[0]   = 0.25;
00208         coeff_out[1]   = 0.25;
00209         coeff_out[2]   = 0.50;
00210     }
00211 }
00212 
00213 static void coefficients_at_mid_elem( double* coeff_out, size_t* indices_out, size_t& num_coeff )
00214 {
00215     num_coeff      = 5;
00216     coeff_out[0]   = 0.125;
00217     coeff_out[1]   = 0.125;
00218     coeff_out[2]   = 0.125;
00219     coeff_out[3]   = 0.125;
00220     coeff_out[4]   = 0.500;
00221     indices_out[0] = 0;
00222     indices_out[1] = 1;
00223     indices_out[2] = 2;
00224     indices_out[3] = 3;
00225     indices_out[4] = 4;
00226 }
00227 
00228 void LinearPyramid::coefficients( Sample loc, NodeSet nodeset, double* coeff_out, size_t* indices_out,
00229                                   size_t& num_coeff, MsqError& err ) const
00230 {
00231     if( nodeset.have_any_mid_node() )
00232     {
00233         MSQ_SETERR( err )( nonlinear_error, MsqError::UNSUPPORTED_ELEMENT );
00234         return;
00235     }
00236 
00237     switch( loc.dimension )
00238     {
00239         case 0:
00240             coefficients_at_corner( loc.number, coeff_out, indices_out, num_coeff );
00241             break;
00242         case 1:
00243             coefficients_at_mid_edge( loc.number, coeff_out, indices_out, num_coeff );
00244             break;
00245         case 2:
00246             coefficients_at_mid_face( loc.number, coeff_out, indices_out, num_coeff );
00247             break;
00248         case 3:
00249             coefficients_at_mid_elem( coeff_out, indices_out, num_coeff );
00250             break;
00251         default:
00252             MSQ_SETERR( err )( "Invalid/unsupported logical dimension", MsqError::INVALID_ARG );
00253     }
00254 }
00255 
00256 void LinearPyramid::derivatives( Sample loc, NodeSet nodeset, size_t* vertex_indices_out,
00257                                  MsqVector< 3 >* d_coeff_d_xi_out, size_t& num_vtx, MsqError& err ) const
00258 {
00259     if( nodeset.have_any_mid_node() )
00260     {
00261         MSQ_SETERR( err )( nonlinear_error, MsqError::UNSUPPORTED_ELEMENT );
00262         return;
00263     }
00264 
00265     switch( loc.dimension )
00266     {
00267         case 0:
00268             if( loc.number == 4 ) { set_equal_derivatives( 0.0, vertex_indices_out, d_coeff_d_xi_out, num_vtx ); }
00269             else
00270             {
00271                 set_corner_derivatives( loc.number, 1.0, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00272             }
00273             break;
00274         case 1:
00275             if( loc.number < 4 )
00276             { set_edge_derivatives( loc.number, 0.50, vertex_indices_out, d_coeff_d_xi_out, num_vtx ); }
00277             else
00278             {
00279                 set_corner_derivatives( loc.number - 4, 0.50, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00280             }
00281             break;
00282         case 2:
00283             if( loc.number == 4 ) { set_equal_derivatives( 0.5, vertex_indices_out, d_coeff_d_xi_out, num_vtx ); }
00284             else
00285             {
00286                 set_edge_derivatives( loc.number, 0.25, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00287             }
00288             break;
00289         case 3:
00290             set_equal_derivatives( 0.25, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00291             break;
00292         default:
00293             MSQ_SETERR( err )( "Invalid/unsupported logical dimension", MsqError::INVALID_ARG );
00294     }
00295 }
00296 
00297 void LinearPyramid::ideal( Sample location, MsqMatrix< 3, 3 >& J, MsqError& ) const
00298 {
00299     // For an ideal element with unit edge length at the base and unit
00300     // height, the Jacobian matrix is:
00301     // | 1-zeta    0      1/2 - xi  |
00302     // |  0       1-zeta  1/2 - eta |
00303     // |  0        0       1        |
00304     //
00305     // The coefficient to produce a unit determinant
00306     // is therefore (1-zeta)^(-2/3).
00307     //
00308     // Thus the unit-determinate ideal element Jacobian
00309     // is, given alpha = (1-zeta)^(-1/3):
00310     //
00311     // | 1/alpha  0      alpha^2 (1/2 - xi) |
00312     // | 0       1/alpha alpha^2 (1/2 - eta)|
00313     // | 0        0      alpha^2            |
00314     //
00315     // There are only three zeta values of interest:
00316     //  zeta = 1 : the degenerate case
00317     //  zeta = 0 : both 1/alpha and alpha^2 are 1.0
00318     //  zeta = 1/2 : 1/alpha = 1/cbrt(2.0) and alpha^2 = 2*(1/alpha)
00319 
00320     // special case for apex
00321     if( location.dimension == 0 && location.number == 4 )
00322     {
00323         J = MsqMatrix< 3, 3 >( 0.0 );
00324         return;
00325     }
00326 
00327     // These are always zero
00328     J( 0, 1 ) = J( 1, 0 ) = J( 2, 0 ) = J( 2, 1 ) = 0.0;
00329 
00330     // Set diagonal terms and magnitude of terms in 3rd column based on zeta
00331 
00332     // All of the zeta=0 locations
00333     double f;
00334     if( location.dimension == 0 || ( location.dimension == 1 && location.number < 4 ) ||
00335         ( location.dimension == 2 && location.number == 4 ) )
00336     {
00337         J( 0, 0 ) = J( 1, 1 ) = J( 2, 2 ) = 1.0;
00338         f                                 = 0.5;
00339     }
00340     // all of the zeta=1/2 locations
00341     else
00342     {
00343         f = J( 0, 0 ) = J( 1, 1 ) = 0.79370052598409979;
00344         J( 2, 2 )                 = 2.0 * f;
00345     }
00346 
00347     // Set terms in 3rd column based on xi,eta
00348 
00349     // The xi = eta = 0.5 locations (mid-element in xi and eta)
00350     if( location.dimension == 3 || ( location.dimension == 2 && location.number == 4 ) )
00351     { J( 0, 2 ) = J( 1, 2 ) = 0.0; }
00352     // The corner locations
00353     else if( location.dimension == 0 || ( location.dimension == 1 && location.number >= 4 ) )
00354     {
00355         switch( location.number % 4 )
00356         {
00357             case 0:
00358                 J( 0, 2 ) = f;
00359                 J( 1, 2 ) = f;
00360                 break;
00361             case 1:
00362                 J( 0, 2 ) = -f;
00363                 J( 1, 2 ) = f;
00364                 break;
00365             case 2:
00366                 J( 0, 2 ) = -f;
00367                 J( 1, 2 ) = -f;
00368                 break;
00369             case 3:
00370                 J( 0, 2 ) = f;
00371                 J( 1, 2 ) = -f;
00372                 break;
00373         }
00374     }
00375     // The mid-edge locations
00376     else
00377     {
00378         switch( location.number )
00379         {
00380             case 0:
00381                 J( 0, 2 ) = 0;
00382                 J( 1, 2 ) = f;
00383                 break;
00384             case 1:
00385                 J( 0, 2 ) = -f;
00386                 J( 1, 2 ) = 0;
00387                 break;
00388             case 2:
00389                 J( 0, 2 ) = 0;
00390                 J( 1, 2 ) = -f;
00391                 break;
00392             case 3:
00393                 J( 0, 2 ) = f;
00394                 J( 1, 2 ) = 0;
00395                 break;
00396         }
00397     }
00398 }
00399 
00400 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines