MOAB: Mesh Oriented datABase  (version 5.4.1)
LinearQuadrilateral.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) [email protected]
00024 
00025   ***************************************************************** */
00026 /** \file LinearQuadrilateral.cpp
00027  *  \author Jason Kraftcheck
00028  */
00029 
00030 #include "Mesquite.hpp"
00031 #include "MsqError.hpp"
00032 #include "LinearQuadrilateral.hpp"
00033 
00034 namespace MBMesquite
00035 {
00036 
00037 static const char* nonlinear_error = "Attempt to use LinearQuadrilateral mapping function for a nonlinear element\n";
00038 
00039 EntityTopology LinearQuadrilateral::element_topology() const
00040 {
00041     return QUADRILATERAL;
00042 }
00043 
00044 int LinearQuadrilateral::num_nodes() const
00045 {
00046     return 4;
00047 }
00048 
00049 void LinearQuadrilateral::coefficients( Sample location,
00050                                         NodeSet nodeset,
00051                                         double* coeff_out,
00052                                         size_t* indices_out,
00053                                         size_t& num_coeff,
00054                                         MsqError& err ) const
00055 {
00056     if( nodeset.have_any_mid_node() )
00057     {
00058         MSQ_SETERR( err )( nonlinear_error, MsqError::UNSUPPORTED_ELEMENT );
00059         return;
00060     }
00061 
00062     coefficients( location, nodeset, coeff_out, indices_out, num_coeff );
00063 }
00064 
00065 void LinearQuadrilateral::coefficients( Sample location,
00066                                         NodeSet,
00067                                         double* coeff_out,
00068                                         size_t* indices_out,
00069                                         size_t& num_coeff )
00070 {
00071     switch( location.dimension )
00072     {
00073         case 0:
00074             num_coeff      = 1;
00075             indices_out[0] = location.number;
00076             coeff_out[0]   = 1.0;
00077             break;
00078         case 1:
00079             num_coeff      = 2;
00080             indices_out[0] = location.number;
00081             indices_out[1] = ( location.number + 1 ) % 4;
00082             coeff_out[0]   = 0.5;
00083             coeff_out[1]   = 0.5;
00084             break;
00085         case 2:
00086             num_coeff      = 4;
00087             indices_out[0] = 0;
00088             indices_out[1] = 1;
00089             indices_out[2] = 2;
00090             indices_out[3] = 3;
00091             coeff_out[0]   = 0.25;
00092             coeff_out[1]   = 0.25;
00093             coeff_out[2]   = 0.25;
00094             coeff_out[3]   = 0.25;
00095             break;
00096         default:
00097             num_coeff = 0;
00098     }
00099 }
00100 
00101 const unsigned xi    = 0;
00102 const unsigned eta   = 1;
00103 const int sign[2][4] = { { -1, 1, 1, -1 },    // xi
00104                          { -1, -1, 1, 1 } };  // eta
00105 
00106 static void derivatives_at_corner( unsigned corner, size_t* vertex_indices, MsqVector< 2 >* derivs, size_t& num_vtx )
00107 {
00108     const unsigned adj_in_xi  = ( 5 - corner ) % 4;
00109     const unsigned adj_in_eta = 3 - corner;
00110 
00111     num_vtx           = 3;
00112     vertex_indices[0] = corner;
00113     vertex_indices[1] = adj_in_xi;
00114     vertex_indices[2] = adj_in_eta;
00115 
00116     derivs[0][0] = sign[xi][corner];
00117     derivs[0][1] = sign[eta][corner];
00118     derivs[1][0] = sign[xi][adj_in_xi];
00119     derivs[1][1] = 0.0;
00120     derivs[2][0] = 0.0;
00121     derivs[2][1] = sign[eta][adj_in_eta];
00122 }
00123 
00124 static void derivatives_at_mid_edge( unsigned edge, size_t* vertices, MsqVector< 2 >* derivs, size_t& num_vtx )
00125 {
00126     const unsigned start_vtx  = edge;
00127     const unsigned end_vtx    = ( edge + 1 ) % 4;
00128     const unsigned othr1_vtx  = ( edge + 2 ) % 4;
00129     const unsigned othr2_vtx  = ( edge + 3 ) % 4;
00130     const unsigned direction  = edge % 2;
00131     const unsigned orthogonal = 1 - direction;
00132 
00133     num_vtx     = 4;
00134     vertices[0] = 0;
00135     vertices[1] = 1;
00136     vertices[2] = 2;
00137     vertices[3] = 3;
00138 
00139     derivs[start_vtx][direction] = sign[direction][start_vtx];
00140     derivs[end_vtx][direction]   = sign[direction][end_vtx];
00141     derivs[othr1_vtx][direction] = 0.0;
00142     derivs[othr2_vtx][direction] = 0.0;
00143 
00144     derivs[0][orthogonal] = 0.5 * sign[orthogonal][0];
00145     derivs[1][orthogonal] = 0.5 * sign[orthogonal][1];
00146     derivs[2][orthogonal] = 0.5 * sign[orthogonal][2];
00147     derivs[3][orthogonal] = 0.5 * sign[orthogonal][3];
00148 }
00149 
00150 static void derivatives_at_mid_elem( size_t* vertices, MsqVector< 2 >* derivs, size_t& num_vtx )
00151 {
00152     num_vtx      = 4;
00153     vertices[0]  = 0;
00154     derivs[0][0] = -0.5;
00155     derivs[0][1] = -0.5;
00156     vertices[1]  = 1;
00157     derivs[1][0] = 0.5;
00158     derivs[1][1] = -0.5;
00159     vertices[2]  = 2;
00160     derivs[2][0] = 0.5;
00161     derivs[2][1] = 0.5;
00162     vertices[3]  = 3;
00163     derivs[3][0] = -0.5;
00164     derivs[3][1] = 0.5;
00165 }
00166 
00167 void LinearQuadrilateral::derivatives( Sample loc,
00168                                        NodeSet nodeset,
00169                                        size_t* vertex_indices_out,
00170                                        MsqVector< 2 >* d_coeff_d_xi_out,
00171                                        size_t& num_vtx,
00172                                        MsqError& err ) const
00173 {
00174     if( nodeset.have_any_mid_node() )
00175     {
00176         MSQ_SETERR( err )( nonlinear_error, MsqError::UNSUPPORTED_ELEMENT );
00177         return;
00178     }
00179 
00180     derivatives( loc, nodeset, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00181 }
00182 
00183 void LinearQuadrilateral::derivatives( Sample loc,
00184                                        NodeSet,
00185                                        size_t* vertex_indices_out,
00186                                        MsqVector< 2 >* d_coeff_d_xi_out,
00187                                        size_t& num_vtx )
00188 {
00189     switch( loc.dimension )
00190     {
00191         case 0:
00192             derivatives_at_corner( loc.number, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00193             break;
00194         case 1:
00195             derivatives_at_mid_edge( loc.number, vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00196             break;
00197         case 2:
00198             derivatives_at_mid_elem( vertex_indices_out, d_coeff_d_xi_out, num_vtx );
00199             break;
00200         default:
00201             num_vtx = 0;
00202     }
00203 }
00204 
00205 void LinearQuadrilateral::ideal( Sample, MsqMatrix< 3, 2 >& J, MsqError& ) const
00206 {
00207     J( 0, 0 ) = J( 1, 1 ) = 1.0;
00208     J( 0, 1 ) = J( 1, 0 ) = 0.0;
00209     J( 2, 0 ) = J( 2, 1 ) = 0.0;
00210 }
00211 
00212 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines