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