MOAB: Mesh Oriented datABase
(version 5.3.1)
|
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 <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() ) { 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