MOAB: Mesh Oriented datABase  (version 5.3.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, NodeSet nodeset, double* coeff_out, size_t* indices_out, size_t& num_coeff,
00118                                 MsqError& err ) const
00119 {
00120     if( nodeset.have_any_mid_node() )
00121     {
00122         MSQ_SETERR( err )( nonlinear_error, MsqError::UNSUPPORTED_ELEMENT );
00123         return;
00124     }
00125 
00126     switch( loc.dimension )
00127     {
00128         case 0:
00129             coefficients_at_corner( loc.number, coeff_out, indices_out, num_coeff );
00130             break;
00131         case 1:
00132             coefficients_at_mid_edge( loc.number, coeff_out, indices_out, num_coeff );
00133             break;
00134         case 2:
00135             coefficients_at_mid_face( loc.number, coeff_out, indices_out, num_coeff );
00136             break;
00137         case 3:
00138             coefficients_at_mid_elem( coeff_out, indices_out, num_coeff );
00139             break;
00140         default:
00141             MSQ_SETERR( err )( "Invalid/unsupported logical dimension", MsqError::INVALID_ARG );
00142     }
00143 }
00144 
00145 static void derivatives_at_corner( unsigned corner, size_t* vertex_indices_out, MsqVector< 3 >* d_coeff_d_xi_out,
00146                                    size_t& num_vtx )
00147 {
00148     int tri = ( corner / 3 );  // 0 for xi=0, 1 for xi=1
00149     int tv  = corner % 3;      // index of corner with xi=constant triangle
00150 
00151     num_vtx = 4;
00152     // three vertices within the xi=constant triangle
00153     vertex_indices_out[0] = 3 * tri;
00154     vertex_indices_out[1] = 3 * tri + 1;
00155     vertex_indices_out[2] = 3 * tri + 2;
00156     // vertex adjacent to corner in other triangle
00157     vertex_indices_out[3] = 3 - 6 * tri + corner;
00158 
00159     // three vertices within the xi=constant triangle
00160     d_coeff_d_xi_out[0][0] = 0.0;
00161     d_coeff_d_xi_out[0][1] = -1.0;
00162     d_coeff_d_xi_out[0][2] = -1.0;
00163     d_coeff_d_xi_out[1][0] = 0.0;
00164     d_coeff_d_xi_out[1][1] = 1.0;
00165     d_coeff_d_xi_out[1][2] = 0.0;
00166     d_coeff_d_xi_out[2][0] = 0.0;
00167     d_coeff_d_xi_out[2][1] = 0.0;
00168     d_coeff_d_xi_out[2][2] = 1.0;
00169     // fix dxi value for input corner
00170     d_coeff_d_xi_out[tv][0] = 2 * tri - 1;
00171     // vertex adjacent to corner in other triangle
00172     d_coeff_d_xi_out[3][0] = 1 - 2 * tri;
00173     d_coeff_d_xi_out[3][1] = 0.0;
00174     d_coeff_d_xi_out[3][2] = 0.0;
00175 }
00176 
00177 static void derivatives_at_mid_edge( unsigned edge, size_t* vertex_indices_out, MsqVector< 3 >* d_coeff_d_xi_out,
00178                                      size_t& num_vtx )
00179 {
00180     int opp;  // vertex opposite edge in same triagle
00181 
00182     switch( edge / 3 )
00183     {
00184         case 0:  // triangle at xi = 0
00185             opp = ( edge + 2 ) % 3;
00186 
00187             num_vtx = 5;
00188             // vertices in this xi = 0 triagnle
00189             vertex_indices_out[0] = 0;
00190             vertex_indices_out[1] = 1;
00191             vertex_indices_out[2] = 2;
00192             // adjacent vertices in xi = 1 triangle
00193             vertex_indices_out[3] = 3 + edge;
00194             vertex_indices_out[4] = 3 + ( edge + 1 ) % 3;
00195 
00196             // vertices in this xi = 0 triagnle
00197             d_coeff_d_xi_out[0][0] = -0.5;
00198             d_coeff_d_xi_out[0][1] = -1.0;
00199             d_coeff_d_xi_out[0][2] = -1.0;
00200             d_coeff_d_xi_out[1][0] = -0.5;
00201             d_coeff_d_xi_out[1][1] = 1.0;
00202             d_coeff_d_xi_out[1][2] = 0.0;
00203             d_coeff_d_xi_out[2][0] = -0.5;
00204             d_coeff_d_xi_out[2][1] = 0.0;
00205             d_coeff_d_xi_out[2][2] = 1.0;
00206             // clear dxi for vertex opposite edge in xi = 0 triangle
00207             d_coeff_d_xi_out[opp][0] = 0.0;
00208             // adjacent vertices in xi = 1 triangle
00209             d_coeff_d_xi_out[3][0] = 0.5;
00210             d_coeff_d_xi_out[3][1] = 0.0;
00211             d_coeff_d_xi_out[3][2] = 0.0;
00212             d_coeff_d_xi_out[4][0] = 0.5;
00213             d_coeff_d_xi_out[4][1] = 0.0;
00214             d_coeff_d_xi_out[4][2] = 0.0;
00215             break;
00216 
00217         case 1:  // lateral edges (not in either triangle)
00218             num_vtx               = 6;
00219             vertex_indices_out[0] = 0;
00220             vertex_indices_out[1] = 1;
00221             vertex_indices_out[2] = 2;
00222             vertex_indices_out[3] = 3;
00223             vertex_indices_out[4] = 4;
00224             vertex_indices_out[5] = 5;
00225 
00226             // set all deta & dzeta values, zero all dxi values
00227             d_coeff_d_xi_out[0][0] = 0.0;
00228             d_coeff_d_xi_out[0][1] = -0.5;
00229             d_coeff_d_xi_out[0][2] = -0.5;
00230             d_coeff_d_xi_out[1][0] = 0.0;
00231             d_coeff_d_xi_out[1][1] = 0.5;
00232             d_coeff_d_xi_out[1][2] = 0.0;
00233             d_coeff_d_xi_out[2][0] = 0.0;
00234             d_coeff_d_xi_out[2][1] = 0.0;
00235             d_coeff_d_xi_out[2][2] = 0.5;
00236             d_coeff_d_xi_out[3][0] = 0.0;
00237             d_coeff_d_xi_out[3][1] = -0.5;
00238             d_coeff_d_xi_out[3][2] = -0.5;
00239             d_coeff_d_xi_out[4][0] = 0.0;
00240             d_coeff_d_xi_out[4][1] = 0.5;
00241             d_coeff_d_xi_out[4][2] = 0.0;
00242             d_coeff_d_xi_out[5][0] = 0.0;
00243             d_coeff_d_xi_out[5][1] = 0.0;
00244             d_coeff_d_xi_out[5][2] = 0.5;
00245 
00246             // set dxi values for end points of edge
00247             d_coeff_d_xi_out[( edge - 3 )][0] = -1;
00248             d_coeff_d_xi_out[edge][0]         = 1;
00249             break;
00250 
00251         case 2:  // triangle at xi = 1
00252             opp = ( edge + 2 ) % 3;
00253 
00254             num_vtx = 5;
00255             // vertices in this xi = 1 triagnle
00256             vertex_indices_out[0] = 3;
00257             vertex_indices_out[1] = 4;
00258             vertex_indices_out[2] = 5;
00259             // adjacent vertices in xi = 1 triangle
00260             vertex_indices_out[3] = edge - 6;
00261             vertex_indices_out[4] = ( edge - 5 ) % 3;
00262 
00263             // vertices in this xi = 1 triagnle
00264             d_coeff_d_xi_out[0][0] = 0.5;
00265             d_coeff_d_xi_out[0][1] = -1.0;
00266             d_coeff_d_xi_out[0][2] = -1.0;
00267             d_coeff_d_xi_out[1][0] = 0.5;
00268             d_coeff_d_xi_out[1][1] = 1.0;
00269             d_coeff_d_xi_out[1][2] = 0.0;
00270             d_coeff_d_xi_out[2][0] = 0.5;
00271             d_coeff_d_xi_out[2][1] = 0.0;
00272             d_coeff_d_xi_out[2][2] = 1.0;
00273             // clear dxi for vertex opposite edge in xi = 1 triangle
00274             d_coeff_d_xi_out[opp][0] = 0.0;
00275             // adjacent vertices in xi = 0 triangle
00276             d_coeff_d_xi_out[3][0] = -0.5;
00277             d_coeff_d_xi_out[3][1] = 0.0;
00278             d_coeff_d_xi_out[3][2] = 0.0;
00279             d_coeff_d_xi_out[4][0] = -0.5;
00280             d_coeff_d_xi_out[4][1] = 0.0;
00281             d_coeff_d_xi_out[4][2] = 0.0;
00282             break;
00283     }
00284 }
00285 static void derivatives_at_mid_face( unsigned face, size_t* vertex_indices_out, MsqVector< 3 >* d_coeff_d_xi_out,
00286                                      size_t& num_vtx )
00287 {
00288     num_vtx               = 6;
00289     vertex_indices_out[0] = 0;
00290     vertex_indices_out[1] = 1;
00291     vertex_indices_out[2] = 2;
00292     vertex_indices_out[3] = 3;
00293     vertex_indices_out[4] = 4;
00294     vertex_indices_out[5] = 5;
00295 
00296     int opp;         // start vtx of edge opposite from quad face
00297     int tri_offset;  // offset in d_coeff_d_xi_out for triangle containing edge
00298 
00299     if( face < 3 )
00300     {  // quad face
00301         // set all values
00302         d_coeff_d_xi_out[0][0] = -0.5;
00303         d_coeff_d_xi_out[0][1] = -0.5;
00304         d_coeff_d_xi_out[0][2] = -0.5;
00305         d_coeff_d_xi_out[1][0] = -0.5;
00306         d_coeff_d_xi_out[1][1] = 0.5;
00307         d_coeff_d_xi_out[1][2] = 0.0;
00308         d_coeff_d_xi_out[2][0] = -0.5;
00309         d_coeff_d_xi_out[2][1] = 0.0;
00310         d_coeff_d_xi_out[2][2] = 0.5;
00311         d_coeff_d_xi_out[3][0] = 0.5;
00312         d_coeff_d_xi_out[3][1] = -0.5;
00313         d_coeff_d_xi_out[3][2] = -0.5;
00314         d_coeff_d_xi_out[4][0] = 0.5;
00315         d_coeff_d_xi_out[4][1] = 0.5;
00316         d_coeff_d_xi_out[4][2] = 0.0;
00317         d_coeff_d_xi_out[5][0] = 0.5;
00318         d_coeff_d_xi_out[5][1] = 0.0;
00319         d_coeff_d_xi_out[5][2] = 0.5;
00320         // clear dxi for ends of edge opposite from face
00321         opp                              = ( face + 2 ) % 3;
00322         d_coeff_d_xi_out[opp][0]         = 0.0;
00323         d_coeff_d_xi_out[( opp + 3 )][0] = 0.0;
00324     }
00325     else
00326     {  // triangular faces
00327         // set all xi values, zero all other values
00328         const double third     = 1. / 3;
00329         d_coeff_d_xi_out[0][0] = -third;
00330         d_coeff_d_xi_out[0][1] = 0;
00331         d_coeff_d_xi_out[0][2] = 0;
00332         d_coeff_d_xi_out[1][0] = -third;
00333         d_coeff_d_xi_out[1][1] = 0;
00334         d_coeff_d_xi_out[1][2] = 0;
00335         d_coeff_d_xi_out[2][0] = -third;
00336         d_coeff_d_xi_out[2][1] = 0;
00337         d_coeff_d_xi_out[2][2] = 0;
00338         d_coeff_d_xi_out[3][0] = third;
00339         d_coeff_d_xi_out[3][1] = 0;
00340         d_coeff_d_xi_out[3][2] = 0;
00341         d_coeff_d_xi_out[4][0] = third;
00342         d_coeff_d_xi_out[4][1] = 0;
00343         d_coeff_d_xi_out[4][2] = 0;
00344         d_coeff_d_xi_out[5][0] = third;
00345         d_coeff_d_xi_out[5][1] = 0;
00346         d_coeff_d_xi_out[5][2] = 0;
00347         // set deta and dzeta values for vertices in same triangle as edge
00348         tri_offset                          = 3 * ( face - 3 );  // either 0 or 3
00349         d_coeff_d_xi_out[tri_offset][1]     = -1.0;
00350         d_coeff_d_xi_out[tri_offset][2]     = -1.0;
00351         d_coeff_d_xi_out[tri_offset + 1][1] = 1.0;
00352         d_coeff_d_xi_out[tri_offset + 2][2] = 1.0;
00353     }
00354 }
00355 static void derivatives_at_mid_elem( size_t* vertex_indices_out, MsqVector< 3 >* d_coeff_d_xi_out, size_t& num_vtx )
00356 {
00357     const double third = 1. / 3;
00358 
00359     num_vtx = 6;
00360     ;
00361     vertex_indices_out[0] = 0;
00362     vertex_indices_out[1] = 1;
00363     vertex_indices_out[2] = 2;
00364     vertex_indices_out[3] = 3;
00365     vertex_indices_out[4] = 4;
00366     vertex_indices_out[5] = 5;
00367 
00368     d_coeff_d_xi_out[0][0] = -third;
00369     d_coeff_d_xi_out[0][1] = -0.5;
00370     d_coeff_d_xi_out[0][2] = -0.5;
00371     d_coeff_d_xi_out[1][0] = -third;
00372     d_coeff_d_xi_out[1][1] = 0.5;
00373     d_coeff_d_xi_out[1][2] = 0.0;
00374     d_coeff_d_xi_out[2][0] = -third;
00375     d_coeff_d_xi_out[2][1] = 0.0;
00376     d_coeff_d_xi_out[2][2] = 0.5;
00377     d_coeff_d_xi_out[3][0] = third;
00378     d_coeff_d_xi_out[3][1] = -0.5;
00379     d_coeff_d_xi_out[3][2] = -0.5;
00380     d_coeff_d_xi_out[4][0] = third;
00381     d_coeff_d_xi_out[4][1] = 0.5;
00382     d_coeff_d_xi_out[4][2] = 0.0;
00383     d_coeff_d_xi_out[5][0] = third;
00384     d_coeff_d_xi_out[5][1] = 0.0;
00385     d_coeff_d_xi_out[5][2] = 0.5;
00386 }
00387 
00388 void LinearPrism::derivatives( Sample loc, NodeSet nodeset, size_t* vertex_indices_out,
00389                                MsqVector< 3 >* d_coeff_d_xi_out, size_t& num_vtx, MsqError& err ) const
00390 {
00391     if( nodeset.have_any_mid_node() )
00392     {
00393         MSQ_SETERR( err )( nonlinear_error, MsqError::UNSUPPORTED_ELEMENT );
00394         return;
00395     }
00396 
00397     switch( loc.dimension )
00398     {
00399         case 0:
00400             derivatives_at_corner( loc.number, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00401             break;
00402         case 1:
00403             derivatives_at_mid_edge( loc.number, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00404             break;
00405         case 2:
00406             derivatives_at_mid_face( loc.number, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00407             break;
00408         case 3:
00409             derivatives_at_mid_elem( vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00410             break;
00411         default:
00412             MSQ_SETERR( err )( "Invalid/unsupported logical dimension", MsqError::INVALID_ARG );
00413     }
00414 }
00415 
00416 void LinearPrism::ideal( Sample, MsqMatrix< 3, 3 >& J, MsqError& ) const
00417 {
00418     const double a = 0.52455753171082409;  // 2^(-2/3) * 3^(-1/6)
00419     const double b = 0.90856029641606983;  // a * sqrt(3) = 1/2 cbrt(6)
00420 
00421     J( 0, 0 ) = 2 * a;
00422     J( 0, 1 ) = 0.0;
00423     J( 0, 2 ) = 0.0;
00424     J( 1, 0 ) = 0.0;
00425     J( 1, 1 ) = 2 * a;
00426     J( 1, 2 ) = a;
00427     J( 2, 0 ) = 0.0;
00428     J( 2, 1 ) = 0.0;
00429     J( 2, 2 ) = b;
00430 }
00431 
00432 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines