MOAB: Mesh Oriented datABase  (version 5.4.1)
TriLagrangeShape.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2006 Sandia National Laboratories.  Developed at the
00005     University of Wisconsin--Madison under SNL contract number
00006     624796.  The U.S. Government and the University of Wisconsin
00007     retain certain rights to 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) [email protected]
00024 
00025   ***************************************************************** */
00026 
00027 /** \file TriLagrangeShape.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "TriLagrangeShape.hpp"
00034 #include "MsqError.hpp"
00035 #include <cassert>
00036 
00037 namespace MBMesquite
00038 {
00039 
00040 EntityTopology TriLagrangeShape::element_topology() const
00041 {
00042     return TRIANGLE;
00043 }
00044 
00045 int TriLagrangeShape::num_nodes() const
00046 {
00047     return 6;
00048 }
00049 
00050 NodeSet TriLagrangeShape::sample_points( NodeSet ns ) const
00051 {
00052     if( ns.have_any_mid_node() )
00053     {
00054         ns.set_all_corner_nodes( TRIANGLE );
00055     }
00056     else
00057     {
00058         ns.clear();
00059         ns.set_mid_face_node( 0 );
00060     }
00061     return ns;
00062 }
00063 
00064 void TriLagrangeShape::coefficients( Sample loc,
00065                                      NodeSet nodeset,
00066                                      double* coeff_out,
00067                                      size_t* indices_out,
00068                                      size_t& num_coeff,
00069                                      MsqError& err ) const
00070 {
00071     if( nodeset.have_any_mid_face_node() )
00072     {
00073         MSQ_SETERR( err )
00074         ( "TriLagrangeShape does not support mid-element nodes", MsqError::UNSUPPORTED_ELEMENT );
00075         return;
00076     }
00077 
00078     switch( loc.dimension )
00079     {
00080         case 0:
00081             num_coeff      = 1;
00082             indices_out[0] = loc.number;
00083             coeff_out[0]   = 1.0;
00084             break;
00085         case 1:
00086             if( nodeset.mid_edge_node( loc.number ) )
00087             {  // if mid-edge node is present
00088                 num_coeff      = 1;
00089                 indices_out[0] = 3 + loc.number;
00090                 coeff_out[0]   = 1.0;
00091             }
00092             else
00093             {  // no mid node on edge
00094                 num_coeff      = 2;
00095                 indices_out[0] = loc.number;
00096                 indices_out[1] = ( loc.number + 1 ) % 3;
00097                 coeff_out[0]   = 0.5;
00098                 coeff_out[1]   = 0.5;
00099             }
00100             break;
00101         case 2:
00102             num_coeff      = 3;
00103             indices_out[0] = 0;
00104             indices_out[1] = 1;
00105             indices_out[2] = 2;
00106             coeff_out[0]   = 1.0 / 3.0;
00107             coeff_out[1]   = 1.0 / 3.0;
00108             coeff_out[2]   = 1.0 / 3.0;
00109             for( int i = 0; i < 3; ++i )
00110             {  // for each mid-edge node
00111                 if( nodeset.mid_edge_node( i ) )
00112                 {  // if node is present
00113                     indices_out[num_coeff] = i + 3;
00114                     // coeff for mid-edge node
00115                     coeff_out[num_coeff] = 4.0 / 9.0;
00116                     // adjust coeff for adj corner nodes
00117                     coeff_out[i] -= 2.0 / 9.0;
00118                     coeff_out[( i + 1 ) % 3] -= 2.0 / 9.0;
00119                     // update count
00120                     ++num_coeff;
00121                 }
00122             }
00123             break;
00124         default:
00125             MSQ_SETERR( err )
00126             ( MsqError::UNSUPPORTED_ELEMENT,
00127               "Request for dimension %d mapping function value"
00128               "for a triangular element",
00129               loc.dimension );
00130     }
00131 }
00132 
00133 static inline void get_linear_derivatives( size_t* vertices, MsqVector< 2 >* derivs )
00134 {
00135     vertices[0]  = 0;
00136     vertices[1]  = 1;
00137     vertices[2]  = 2;
00138     derivs[0][0] = -1.0;
00139     derivs[0][1] = -1.0;
00140     derivs[1][0] = 1.0;
00141     derivs[1][1] = 0.0;
00142     derivs[2][0] = 0.0;
00143     derivs[2][1] = 1.0;
00144 }
00145 
00146 static void derivatives_at_corner( unsigned corner,
00147                                    NodeSet nodeset,
00148                                    size_t* vertices,
00149                                    MsqVector< 2 >* derivs,
00150                                    size_t& num_vtx )
00151 {
00152     num_vtx = 3;
00153     get_linear_derivatives( vertices, derivs );
00154     switch( corner )
00155     {
00156         case 0:
00157             if( nodeset.mid_edge_node( 0 ) )
00158             {
00159                 vertices[num_vtx]  = 3;
00160                 derivs[num_vtx][0] = 4.0;
00161                 derivs[num_vtx][1] = 0.0;
00162                 ++num_vtx;
00163                 derivs[0][0] -= 2.0;
00164                 derivs[1][0] -= 2.0;
00165             }
00166             if( nodeset.mid_edge_node( 2 ) )
00167             {
00168                 vertices[num_vtx]  = 5;
00169                 derivs[num_vtx][0] = 0.0;
00170                 derivs[num_vtx][1] = 4.0;
00171                 ++num_vtx;
00172                 derivs[0][1] -= 2.0;
00173                 derivs[2][1] -= 2.0;
00174             }
00175             break;
00176 
00177         case 1:
00178             if( nodeset.mid_edge_node( 0 ) )
00179             {
00180                 vertices[num_vtx]  = 3;
00181                 derivs[num_vtx][0] = -4.0;
00182                 derivs[num_vtx][1] = -4.0;
00183                 ++num_vtx;
00184                 derivs[0][0] += 2.0;
00185                 derivs[0][1] += 2.0;
00186                 derivs[1][0] += 2.0;
00187                 derivs[1][1] += 2.0;
00188             }
00189             if( nodeset.mid_edge_node( 1 ) )
00190             {
00191                 vertices[num_vtx]  = 4;
00192                 derivs[num_vtx][0] = 0.0;
00193                 derivs[num_vtx][1] = 4.0;
00194                 ++num_vtx;
00195                 derivs[1][1] -= 2.0;
00196                 derivs[2][1] -= 2.0;
00197             }
00198             break;
00199 
00200         case 2:
00201             if( nodeset.mid_edge_node( 1 ) )
00202             {
00203                 vertices[num_vtx]  = 4;
00204                 derivs[num_vtx][0] = 4.0;
00205                 derivs[num_vtx][1] = 0.0;
00206                 ++num_vtx;
00207                 derivs[1][0] -= 2.0;
00208                 derivs[2][0] -= 2.0;
00209             }
00210             if( nodeset.mid_edge_node( 2 ) )
00211             {
00212                 vertices[num_vtx]  = 5;
00213                 derivs[num_vtx][0] = -4.0;
00214                 derivs[num_vtx][1] = -4.0;
00215                 ++num_vtx;
00216                 derivs[0][0] += 2.0;
00217                 derivs[0][1] += 2.0;
00218                 derivs[2][0] += 2.0;
00219                 derivs[2][1] += 2.0;
00220             }
00221             break;
00222     }
00223 }
00224 
00225 static const double edr[3][3] = { { 0.0, -2.0, 2.0 }, { 0.0, 2.0, 2.0 }, { 0.0, -2.0, -2.0 } };
00226 static const double eds[3][3] = { { -2.0, -2.0, 0.0 }, { 2.0, 2.0, 0.0 }, { 2.0, -2.0, 0.0 } };
00227 
00228 static void derivatives_at_mid_edge( unsigned edge,
00229                                      NodeSet nodeset,
00230                                      size_t* vertices,
00231                                      MsqVector< 2 >* derivs,
00232                                      size_t& num_vtx )
00233 {
00234     // The mid-edge behavior is rather strange.
00235     // A corner vertex contributes to the jacobian
00236     // at the mid-point of the opposite edge unless
00237     // one, but *not* both of the adjacent mid-edge
00238     // nodes is present.
00239 
00240     // The value for each corner is incremented when
00241     // a mid-side node is present.  If the value is
00242     // 0 when finished, the corner doesn't contribute.
00243     // Initialize values to 0 for corners adjacent to
00244     // edge so they are never zero.
00245     int corner_count[3]            = { 1, 1, 1 };
00246     corner_count[( edge + 2 ) % 3] = -1;
00247 
00248     // begin with derivatives for linear tri
00249     double corner_derivs[3][2] = { { -1.0, -1.0 }, { 1.0, 0.0 }, { 0.0, 1.0 } };
00250 
00251     // do mid-side nodes first
00252     num_vtx = 0;
00253     for( unsigned i = 0; i < 3; ++i )
00254     {
00255         if( nodeset.mid_edge_node( i ) )
00256         {
00257             vertices[num_vtx]  = i + 3;
00258             derivs[num_vtx][0] = edr[i][edge];
00259             derivs[num_vtx][1] = eds[i][edge];
00260             ++num_vtx;
00261 
00262             int a = ( i + 1 ) % 3;
00263             corner_derivs[i][0] -= 0.5 * edr[i][edge];
00264             corner_derivs[i][1] -= 0.5 * eds[i][edge];
00265             corner_derivs[a][0] -= 0.5 * edr[i][edge];
00266             corner_derivs[a][1] -= 0.5 * eds[i][edge];
00267             ++corner_count[i];
00268             ++corner_count[a];
00269         }
00270     }
00271 
00272     // now add corner nodes to list
00273     for( unsigned i = 0; i < 3; ++i )
00274     {
00275         if( corner_count[i] )
00276         {
00277             vertices[num_vtx]  = i;
00278             derivs[num_vtx][0] = corner_derivs[i][0];
00279             derivs[num_vtx][1] = corner_derivs[i][1];
00280             ++num_vtx;
00281         }
00282     }
00283 }
00284 
00285 static const double fdr[] = { 0.0, 4.0 / 3.0, -4.0 / 3.0 };
00286 static const double fds[] = { -4.0 / 3.0, 4.0 / 3.0, 0.0 };
00287 
00288 static void derivatives_at_mid_elem( NodeSet nodeset, size_t* vertices, MsqVector< 2 >* derivs, size_t& num_vtx )
00289 {
00290     get_linear_derivatives( vertices, derivs );
00291     num_vtx = 3;
00292     for( unsigned i = 0; i < 3; ++i )
00293     {
00294         if( nodeset.mid_edge_node( i ) )
00295         {
00296             vertices[num_vtx]  = i + 3;
00297             derivs[num_vtx][0] = fdr[i];
00298             derivs[num_vtx][1] = fds[i];
00299             ++num_vtx;
00300 
00301             int a = ( i + 1 ) % 3;
00302             derivs[i][0] -= 0.5 * fdr[i];
00303             derivs[i][1] -= 0.5 * fds[i];
00304             derivs[a][0] -= 0.5 * fdr[i];
00305             derivs[a][1] -= 0.5 * fds[i];
00306         }
00307     }
00308 }
00309 
00310 void TriLagrangeShape::derivatives( Sample loc,
00311                                     NodeSet nodeset,
00312                                     size_t* vertex_indices_out,
00313                                     MsqVector< 2 >* d_coeff_d_xi_out,
00314                                     size_t& num_vtx,
00315                                     MsqError& err ) const
00316 {
00317     if( !nodeset.have_any_mid_node() )
00318     {
00319         num_vtx = 3;
00320         get_linear_derivatives( vertex_indices_out, d_coeff_d_xi_out );
00321         return;
00322     }
00323 
00324     if( nodeset.have_any_mid_face_node() )
00325     {
00326         MSQ_SETERR( err )
00327         ( "TriLagrangeShape does not support mid-element nodes", MsqError::UNSUPPORTED_ELEMENT );
00328         return;
00329     }
00330 
00331     switch( loc.dimension )
00332     {
00333         case 0:
00334             derivatives_at_corner( loc.number, nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00335             break;
00336         case 1:
00337             derivatives_at_mid_edge( loc.number, nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00338             break;
00339         case 2:
00340             derivatives_at_mid_elem( nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00341             break;
00342         default:
00343             MSQ_SETERR( err )( "Invalid/unsupported logical dimension", MsqError::INVALID_ARG );
00344     }
00345 }
00346 
00347 void TriLagrangeShape::ideal( Sample, MsqMatrix< 3, 2 >& J, MsqError& ) const
00348 {
00349     //  const double g = 1.074569931823542; // sqrt(2/sqrt(3));
00350     //  J(0,0) = g;    J(0,1) = 0.5*g;
00351     //  J(1,0) = 0.0;  J(1,1) = 1.0/g;
00352     //  J(2,0) = 0.0;  J(2,1) = 0.0;
00353     const double a = 0.70710678118654746;  // 1/sqrt(2)
00354     const double b = 1.3160740129524924;   // 4th root of 3
00355     J( 0, 0 )      = -a / b;
00356     J( 0, 1 )      = a / b;
00357     J( 1, 0 )      = -a * b;
00358     J( 1, 1 )      = -a * b;
00359     J( 2, 0 )      = 0.0;
00360     J( 2, 1 )      = 0.0;
00361 }
00362 
00363 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines