MOAB: Mesh Oriented datABase  (version 5.4.0)
LinearPrism.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 LinearPrism.cpp
00028  *  \brief mapping function for linear prism
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "MsqError.hpp"
00034 #include "LinearPrism.hpp"
00035 
00036 namespace MBMesquite
00037 {
00038 
00039 static const char* nonlinear_error = "Attempt to use LinearPrism mapping function for a nonlinear element\n";
00040 
00041 EntityTopology LinearPrism::element_topology() const
00042 {
00043     return PRISM;
00044 }
00045 
00046 int LinearPrism::num_nodes() const
00047 {
00048     return 6;
00049 }
00050 
00051 static const int edge_beg[]  = { 0, 1, 2, 0, 1, 2, 3, 4, 5 };
00052 static const int edge_end[]  = { 1, 2, 0, 3, 4, 5, 4, 5, 3 };
00053 static const int faces[5][5] = { { 4, 0, 1, 4, 3 },
00054                                  { 4, 1, 2, 5, 4 },
00055                                  { 4, 2, 0, 3, 5 },
00056                                  { 3, 0, 1, 2, -1 },
00057                                  { 3, 3, 4, 5, -1 } };
00058 
00059 static void coefficients_at_corner( unsigned corner, double* coeff_out, size_t* indices_out, size_t& num_coeff )
00060 {
00061     num_coeff      = 1;
00062     indices_out[0] = corner;
00063     coeff_out[0]   = 1.0;
00064 }
00065 
00066 static void coefficients_at_mid_edge( unsigned edge, double* coeff_out, size_t* indices_out, size_t& num_coeff )
00067 {
00068     num_coeff      = 2;
00069     indices_out[0] = edge_beg[edge];
00070     indices_out[1] = edge_end[edge];
00071     coeff_out[0]   = 0.5;
00072     coeff_out[1]   = 0.5;
00073 }
00074 
00075 static void coefficients_at_mid_face( unsigned face, double* coeff_out, size_t* indices_out, size_t& num_coeff )
00076 {
00077     double f;
00078     if( faces[face][0] == 4 )
00079     {
00080         num_coeff      = 4;
00081         f              = 0.25;
00082         indices_out[3] = faces[face][4];
00083         coeff_out[3]   = f;
00084     }
00085     else
00086     {
00087         num_coeff = 3;
00088         f         = MSQ_ONE_THIRD;
00089     }
00090 
00091     coeff_out[0]   = f;
00092     coeff_out[1]   = f;
00093     coeff_out[2]   = f;
00094     indices_out[0] = faces[face][1];
00095     indices_out[1] = faces[face][2];
00096     indices_out[2] = faces[face][3];
00097 }
00098 
00099 static void coefficients_at_mid_elem( double* coeff_out, size_t* indices_out, size_t& num_coeff )
00100 {
00101     num_coeff          = 6;
00102     const double sixth = 1.0 / 6.0;
00103     coeff_out[0]       = sixth;
00104     coeff_out[1]       = sixth;
00105     coeff_out[2]       = sixth;
00106     coeff_out[3]       = sixth;
00107     coeff_out[4]       = sixth;
00108     coeff_out[5]       = sixth;
00109     indices_out[0]     = 0;
00110     indices_out[1]     = 1;
00111     indices_out[2]     = 2;
00112     indices_out[3]     = 3;
00113     indices_out[4]     = 4;
00114     indices_out[5]     = 5;
00115 }
00116 
00117 void LinearPrism::coefficients( Sample loc,
00118                                 NodeSet nodeset,
00119                                 double* coeff_out,
00120                                 size_t* indices_out,
00121                                 size_t& num_coeff,
00122                                 MsqError& err ) const
00123 {
00124     if( nodeset.have_any_mid_node() )
00125     {
00126         MSQ_SETERR( err )( nonlinear_error, MsqError::UNSUPPORTED_ELEMENT );
00127         return;
00128     }
00129 
00130     switch( loc.dimension )
00131     {
00132         case 0:
00133             coefficients_at_corner( loc.number, coeff_out, indices_out, num_coeff );
00134             break;
00135         case 1:
00136             coefficients_at_mid_edge( loc.number, coeff_out, indices_out, num_coeff );
00137             break;
00138         case 2:
00139             coefficients_at_mid_face( loc.number, coeff_out, indices_out, num_coeff );
00140             break;
00141         case 3:
00142             coefficients_at_mid_elem( coeff_out, indices_out, num_coeff );
00143             break;
00144         default:
00145             MSQ_SETERR( err )( "Invalid/unsupported logical dimension", MsqError::INVALID_ARG );
00146     }
00147 }
00148 
00149 static void derivatives_at_corner( unsigned corner,
00150                                    size_t* vertex_indices_out,
00151                                    MsqVector< 3 >* d_coeff_d_xi_out,
00152                                    size_t& num_vtx )
00153 {
00154     int tri = ( corner / 3 );  // 0 for xi=0, 1 for xi=1
00155     int tv  = corner % 3;      // index of corner with xi=constant triangle
00156 
00157     num_vtx = 4;
00158     // three vertices within the xi=constant triangle
00159     vertex_indices_out[0] = 3 * tri;
00160     vertex_indices_out[1] = 3 * tri + 1;
00161     vertex_indices_out[2] = 3 * tri + 2;
00162     // vertex adjacent to corner in other triangle
00163     vertex_indices_out[3] = 3 - 6 * tri + corner;
00164 
00165     // three vertices within the xi=constant triangle
00166     d_coeff_d_xi_out[0][0] = 0.0;
00167     d_coeff_d_xi_out[0][1] = -1.0;
00168     d_coeff_d_xi_out[0][2] = -1.0;
00169     d_coeff_d_xi_out[1][0] = 0.0;
00170     d_coeff_d_xi_out[1][1] = 1.0;
00171     d_coeff_d_xi_out[1][2] = 0.0;
00172     d_coeff_d_xi_out[2][0] = 0.0;
00173     d_coeff_d_xi_out[2][1] = 0.0;
00174     d_coeff_d_xi_out[2][2] = 1.0;
00175     // fix dxi value for input corner
00176     d_coeff_d_xi_out[tv][0] = 2 * tri - 1;
00177     // vertex adjacent to corner in other triangle
00178     d_coeff_d_xi_out[3][0] = 1 - 2 * tri;
00179     d_coeff_d_xi_out[3][1] = 0.0;
00180     d_coeff_d_xi_out[3][2] = 0.0;
00181 }
00182 
00183 static void derivatives_at_mid_edge( unsigned edge,
00184                                      size_t* vertex_indices_out,
00185                                      MsqVector< 3 >* d_coeff_d_xi_out,
00186                                      size_t& num_vtx )
00187 {
00188     int opp;  // vertex opposite edge in same triagle
00189 
00190     switch( edge / 3 )
00191     {
00192         case 0:  // triangle at xi = 0
00193             opp = ( edge + 2 ) % 3;
00194 
00195             num_vtx = 5;
00196             // vertices in this xi = 0 triagnle
00197             vertex_indices_out[0] = 0;
00198             vertex_indices_out[1] = 1;
00199             vertex_indices_out[2] = 2;
00200             // adjacent vertices in xi = 1 triangle
00201             vertex_indices_out[3] = 3 + edge;
00202             vertex_indices_out[4] = 3 + ( edge + 1 ) % 3;
00203 
00204             // vertices in this xi = 0 triagnle
00205             d_coeff_d_xi_out[0][0] = -0.5;
00206             d_coeff_d_xi_out[0][1] = -1.0;
00207             d_coeff_d_xi_out[0][2] = -1.0;
00208             d_coeff_d_xi_out[1][0] = -0.5;
00209             d_coeff_d_xi_out[1][1] = 1.0;
00210             d_coeff_d_xi_out[1][2] = 0.0;
00211             d_coeff_d_xi_out[2][0] = -0.5;
00212             d_coeff_d_xi_out[2][1] = 0.0;
00213             d_coeff_d_xi_out[2][2] = 1.0;
00214             // clear dxi for vertex opposite edge in xi = 0 triangle
00215             d_coeff_d_xi_out[opp][0] = 0.0;
00216             // adjacent vertices in xi = 1 triangle
00217             d_coeff_d_xi_out[3][0] = 0.5;
00218             d_coeff_d_xi_out[3][1] = 0.0;
00219             d_coeff_d_xi_out[3][2] = 0.0;
00220             d_coeff_d_xi_out[4][0] = 0.5;
00221             d_coeff_d_xi_out[4][1] = 0.0;
00222             d_coeff_d_xi_out[4][2] = 0.0;
00223             break;
00224 
00225         case 1:  // lateral edges (not in either triangle)
00226             num_vtx               = 6;
00227             vertex_indices_out[0] = 0;
00228             vertex_indices_out[1] = 1;
00229             vertex_indices_out[2] = 2;
00230             vertex_indices_out[3] = 3;
00231             vertex_indices_out[4] = 4;
00232             vertex_indices_out[5] = 5;
00233 
00234             // set all deta & dzeta values, zero all dxi values
00235             d_coeff_d_xi_out[0][0] = 0.0;
00236             d_coeff_d_xi_out[0][1] = -0.5;
00237             d_coeff_d_xi_out[0][2] = -0.5;
00238             d_coeff_d_xi_out[1][0] = 0.0;
00239             d_coeff_d_xi_out[1][1] = 0.5;
00240             d_coeff_d_xi_out[1][2] = 0.0;
00241             d_coeff_d_xi_out[2][0] = 0.0;
00242             d_coeff_d_xi_out[2][1] = 0.0;
00243             d_coeff_d_xi_out[2][2] = 0.5;
00244             d_coeff_d_xi_out[3][0] = 0.0;
00245             d_coeff_d_xi_out[3][1] = -0.5;
00246             d_coeff_d_xi_out[3][2] = -0.5;
00247             d_coeff_d_xi_out[4][0] = 0.0;
00248             d_coeff_d_xi_out[4][1] = 0.5;
00249             d_coeff_d_xi_out[4][2] = 0.0;
00250             d_coeff_d_xi_out[5][0] = 0.0;
00251             d_coeff_d_xi_out[5][1] = 0.0;
00252             d_coeff_d_xi_out[5][2] = 0.5;
00253 
00254             // set dxi values for end points of edge
00255             d_coeff_d_xi_out[( edge - 3 )][0] = -1;
00256             d_coeff_d_xi_out[edge][0]         = 1;
00257             break;
00258 
00259         case 2:  // triangle at xi = 1
00260             opp = ( edge + 2 ) % 3;
00261 
00262             num_vtx = 5;
00263             // vertices in this xi = 1 triagnle
00264             vertex_indices_out[0] = 3;
00265             vertex_indices_out[1] = 4;
00266             vertex_indices_out[2] = 5;
00267             // adjacent vertices in xi = 1 triangle
00268             vertex_indices_out[3] = edge - 6;
00269             vertex_indices_out[4] = ( edge - 5 ) % 3;
00270 
00271             // vertices in this xi = 1 triagnle
00272             d_coeff_d_xi_out[0][0] = 0.5;
00273             d_coeff_d_xi_out[0][1] = -1.0;
00274             d_coeff_d_xi_out[0][2] = -1.0;
00275             d_coeff_d_xi_out[1][0] = 0.5;
00276             d_coeff_d_xi_out[1][1] = 1.0;
00277             d_coeff_d_xi_out[1][2] = 0.0;
00278             d_coeff_d_xi_out[2][0] = 0.5;
00279             d_coeff_d_xi_out[2][1] = 0.0;
00280             d_coeff_d_xi_out[2][2] = 1.0;
00281             // clear dxi for vertex opposite edge in xi = 1 triangle
00282             d_coeff_d_xi_out[opp][0] = 0.0;
00283             // adjacent vertices in xi = 0 triangle
00284             d_coeff_d_xi_out[3][0] = -0.5;
00285             d_coeff_d_xi_out[3][1] = 0.0;
00286             d_coeff_d_xi_out[3][2] = 0.0;
00287             d_coeff_d_xi_out[4][0] = -0.5;
00288             d_coeff_d_xi_out[4][1] = 0.0;
00289             d_coeff_d_xi_out[4][2] = 0.0;
00290             break;
00291     }
00292 }
00293 static void derivatives_at_mid_face( unsigned face,
00294                                      size_t* vertex_indices_out,
00295                                      MsqVector< 3 >* d_coeff_d_xi_out,
00296                                      size_t& num_vtx )
00297 {
00298     num_vtx               = 6;
00299     vertex_indices_out[0] = 0;
00300     vertex_indices_out[1] = 1;
00301     vertex_indices_out[2] = 2;
00302     vertex_indices_out[3] = 3;
00303     vertex_indices_out[4] = 4;
00304     vertex_indices_out[5] = 5;
00305 
00306     int opp;         // start vtx of edge opposite from quad face
00307     int tri_offset;  // offset in d_coeff_d_xi_out for triangle containing edge
00308 
00309     if( face < 3 )
00310     {  // quad face
00311         // set all values
00312         d_coeff_d_xi_out[0][0] = -0.5;
00313         d_coeff_d_xi_out[0][1] = -0.5;
00314         d_coeff_d_xi_out[0][2] = -0.5;
00315         d_coeff_d_xi_out[1][0] = -0.5;
00316         d_coeff_d_xi_out[1][1] = 0.5;
00317         d_coeff_d_xi_out[1][2] = 0.0;
00318         d_coeff_d_xi_out[2][0] = -0.5;
00319         d_coeff_d_xi_out[2][1] = 0.0;
00320         d_coeff_d_xi_out[2][2] = 0.5;
00321         d_coeff_d_xi_out[3][0] = 0.5;
00322         d_coeff_d_xi_out[3][1] = -0.5;
00323         d_coeff_d_xi_out[3][2] = -0.5;
00324         d_coeff_d_xi_out[4][0] = 0.5;
00325         d_coeff_d_xi_out[4][1] = 0.5;
00326         d_coeff_d_xi_out[4][2] = 0.0;
00327         d_coeff_d_xi_out[5][0] = 0.5;
00328         d_coeff_d_xi_out[5][1] = 0.0;
00329         d_coeff_d_xi_out[5][2] = 0.5;
00330         // clear dxi for ends of edge opposite from face
00331         opp                              = ( face + 2 ) % 3;
00332         d_coeff_d_xi_out[opp][0]         = 0.0;
00333         d_coeff_d_xi_out[( opp + 3 )][0] = 0.0;
00334     }
00335     else
00336     {  // triangular faces
00337         // set all xi values, zero all other values
00338         const double third     = 1. / 3;
00339         d_coeff_d_xi_out[0][0] = -third;
00340         d_coeff_d_xi_out[0][1] = 0;
00341         d_coeff_d_xi_out[0][2] = 0;
00342         d_coeff_d_xi_out[1][0] = -third;
00343         d_coeff_d_xi_out[1][1] = 0;
00344         d_coeff_d_xi_out[1][2] = 0;
00345         d_coeff_d_xi_out[2][0] = -third;
00346         d_coeff_d_xi_out[2][1] = 0;
00347         d_coeff_d_xi_out[2][2] = 0;
00348         d_coeff_d_xi_out[3][0] = third;
00349         d_coeff_d_xi_out[3][1] = 0;
00350         d_coeff_d_xi_out[3][2] = 0;
00351         d_coeff_d_xi_out[4][0] = third;
00352         d_coeff_d_xi_out[4][1] = 0;
00353         d_coeff_d_xi_out[4][2] = 0;
00354         d_coeff_d_xi_out[5][0] = third;
00355         d_coeff_d_xi_out[5][1] = 0;
00356         d_coeff_d_xi_out[5][2] = 0;
00357         // set deta and dzeta values for vertices in same triangle as edge
00358         tri_offset                          = 3 * ( face - 3 );  // either 0 or 3
00359         d_coeff_d_xi_out[tri_offset][1]     = -1.0;
00360         d_coeff_d_xi_out[tri_offset][2]     = -1.0;
00361         d_coeff_d_xi_out[tri_offset + 1][1] = 1.0;
00362         d_coeff_d_xi_out[tri_offset + 2][2] = 1.0;
00363     }
00364 }
00365 static void derivatives_at_mid_elem( size_t* vertex_indices_out, MsqVector< 3 >* d_coeff_d_xi_out, size_t& num_vtx )
00366 {
00367     const double third = 1. / 3;
00368 
00369     num_vtx = 6;
00370     ;
00371     vertex_indices_out[0] = 0;
00372     vertex_indices_out[1] = 1;
00373     vertex_indices_out[2] = 2;
00374     vertex_indices_out[3] = 3;
00375     vertex_indices_out[4] = 4;
00376     vertex_indices_out[5] = 5;
00377 
00378     d_coeff_d_xi_out[0][0] = -third;
00379     d_coeff_d_xi_out[0][1] = -0.5;
00380     d_coeff_d_xi_out[0][2] = -0.5;
00381     d_coeff_d_xi_out[1][0] = -third;
00382     d_coeff_d_xi_out[1][1] = 0.5;
00383     d_coeff_d_xi_out[1][2] = 0.0;
00384     d_coeff_d_xi_out[2][0] = -third;
00385     d_coeff_d_xi_out[2][1] = 0.0;
00386     d_coeff_d_xi_out[2][2] = 0.5;
00387     d_coeff_d_xi_out[3][0] = third;
00388     d_coeff_d_xi_out[3][1] = -0.5;
00389     d_coeff_d_xi_out[3][2] = -0.5;
00390     d_coeff_d_xi_out[4][0] = third;
00391     d_coeff_d_xi_out[4][1] = 0.5;
00392     d_coeff_d_xi_out[4][2] = 0.0;
00393     d_coeff_d_xi_out[5][0] = third;
00394     d_coeff_d_xi_out[5][1] = 0.0;
00395     d_coeff_d_xi_out[5][2] = 0.5;
00396 }
00397 
00398 void LinearPrism::derivatives( Sample loc,
00399                                NodeSet nodeset,
00400                                size_t* vertex_indices_out,
00401                                MsqVector< 3 >* d_coeff_d_xi_out,
00402                                size_t& num_vtx,
00403                                MsqError& err ) const
00404 {
00405     if( nodeset.have_any_mid_node() )
00406     {
00407         MSQ_SETERR( err )( nonlinear_error, MsqError::UNSUPPORTED_ELEMENT );
00408         return;
00409     }
00410 
00411     switch( loc.dimension )
00412     {
00413         case 0:
00414             derivatives_at_corner( loc.number, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00415             break;
00416         case 1:
00417             derivatives_at_mid_edge( loc.number, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00418             break;
00419         case 2:
00420             derivatives_at_mid_face( loc.number, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00421             break;
00422         case 3:
00423             derivatives_at_mid_elem( vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00424             break;
00425         default:
00426             MSQ_SETERR( err )( "Invalid/unsupported logical dimension", MsqError::INVALID_ARG );
00427     }
00428 }
00429 
00430 void LinearPrism::ideal( Sample, MsqMatrix< 3, 3 >& J, MsqError& ) const
00431 {
00432     const double a = 0.52455753171082409;  // 2^(-2/3) * 3^(-1/6)
00433     const double b = 0.90856029641606983;  // a * sqrt(3) = 1/2 cbrt(6)
00434 
00435     J( 0, 0 ) = 2 * a;
00436     J( 0, 1 ) = 0.0;
00437     J( 0, 2 ) = 0.0;
00438     J( 1, 0 ) = 0.0;
00439     J( 1, 1 ) = 2 * a;
00440     J( 1, 2 ) = a;
00441     J( 2, 0 ) = 0.0;
00442     J( 2, 1 ) = 0.0;
00443     J( 2, 2 ) = b;
00444 }
00445 
00446 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines