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