MOAB: Mesh Oriented datABase  (version 5.2.1)
QuadLagrangeShapeTest.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     retian 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 QuadLagrangeShapeTest.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "QuadLagrangeShape.hpp"
00034 #include "TopologyInfo.hpp"
00035 #include "MsqError.hpp"
00036 #include "IdealElements.hpp"
00037 #include "JacobianCalculator.hpp"
00038 
00039 #include "UnitUtil.hpp"
00040 
00041 #include <vector>
00042 #include <algorithm>
00043 #include <iostream>
00044 #include <sstream>
00045 
00046 using namespace MBMesquite;
00047 using namespace std;
00048 const double epsilon = 1e-6;
00049 #define ASSERT_VALUES_EQUAL( v1, v2, location, bits ) \
00050     ASSERT_MESSAGE( value_message( ( location ), ( bits ), ( v1 ), ( v2 ) ), ( fabs( ( v1 ) - ( v2 ) ) < epsilon ) )
00051 
00052 static inline CppUnit::Message value_message( unsigned location, NodeSet bits, double v1, double v2 )
00053 {
00054     CppUnit::Message m( "equality assertion failed" );
00055 
00056     std::ostringstream buffer1;
00057     buffer1 << "Expected : " << v1;
00058     m.addDetail( buffer1.str() );
00059 
00060     std::ostringstream buffer2;
00061     buffer2 << "Actual   : " << v2;
00062     m.addDetail( buffer2.str() );
00063 
00064     std::ostringstream buffer3;
00065     buffer3 << "Location : ";
00066     if( location < 4 )
00067         buffer3 << "Corner " << location;
00068     else if( location < 8 )
00069         buffer3 << "Edge " << location - 4;
00070     else if( location == 8 )
00071         buffer3 << "Mid-element";
00072     else
00073         buffer3 << "INVALID!!";
00074     m.addDetail( buffer3.str() );
00075 
00076     std::ostringstream buffer4;
00077     buffer4 << "Node Bits: " << bits;
00078     m.addDetail( buffer4.str() );
00079     return m;
00080 }
00081 
00082 class QuadLagrangeShapeTest : public CppUnit::TestFixture
00083 {
00084   private:
00085     CPPUNIT_TEST_SUITE( QuadLagrangeShapeTest );
00086 
00087     CPPUNIT_TEST( test_coeff_corners );
00088     CPPUNIT_TEST( test_coeff_edges );
00089     CPPUNIT_TEST( test_coeff_center );
00090 
00091     CPPUNIT_TEST( test_deriv_corners );
00092     CPPUNIT_TEST( test_deriv_edges );
00093     CPPUNIT_TEST( test_deriv_center );
00094 
00095     CPPUNIT_TEST( test_ideal_jacobian );
00096 
00097     CPPUNIT_TEST_SUITE_END();
00098 
00099     QuadLagrangeShape sf;
00100 
00101     void test_corner_coeff( int corner, NodeSet nodeset );
00102     void test_edge_coeff( int edge, NodeSet nodeset );
00103     void test_mid_coeff( NodeSet nodeset );
00104 
00105     void test_corner_derivs( int corner, NodeSet nodeset );
00106     void test_edge_derivs( int edge, NodeSet nodeset );
00107     void test_mid_derivs( NodeSet nodeset );
00108 
00109   public:
00110     void test_coeff_corners();
00111     void test_coeff_edges();
00112     void test_coeff_center();
00113 
00114     void test_deriv_corners();
00115     void test_deriv_edges();
00116     void test_deriv_center();
00117 
00118     void test_ideal_jacobian();
00119 };
00120 
00121 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( QuadLagrangeShapeTest, "QuadLagrangeShapeTest" );
00122 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( QuadLagrangeShapeTest, "Unit" );
00123 
00124 // Indices
00125 enum
00126 {
00127     XI  = 0,
00128     ETA = 1
00129 };
00130 // Parameter range
00131 // const double min_xi = -1;
00132 // const double max_xi =  1;
00133 // reparameterize in [0,1]
00134 const double min_xi = 0;
00135 const double max_xi = 1;
00136 const double mid_xi = 0.5 * ( min_xi + max_xi );
00137 
00138 // Some lagrange polynomial terms
00139 // static double l11( double xi ) { return 0.5 * (1 - xi); }
00140 // static double l12( double xi ) { return 0.5 * (1 + xi); }
00141 // static double l22( double xi ) { return 1 - xi*xi; }
00142 // reparameterize in [0,1]
00143 static double l11( double xi )
00144 {
00145     return 1 - xi;
00146 }
00147 static double l12( double xi )
00148 {
00149     return xi;
00150 }
00151 static double l22( double xi )
00152 {
00153     return 4 * xi * ( 1 - xi );
00154 }
00155 // First derivatives of lagrange polynomial terms
00156 // static double dl11( double    ) { return -0.5; }
00157 // static double dl12( double    ) { return  0.5; }
00158 // static double dl22( double xi ) { return -2 * xi; }
00159 // reparameterize in [0,1]
00160 static double dl11( double )
00161 {
00162     return -1;
00163 }
00164 static double dl12( double )
00165 {
00166     return 1;
00167 }
00168 static double dl22( double xi )
00169 {
00170     return 4 - 8 * xi;
00171 }
00172 
00173 // Mapping function weights.  See pg. 135 of Hughes and 'eval' method below
00174 // for an explanation of why we begin with linear element weight functions
00175 // for the corner nodes and a mixture of quadratic and linear weight
00176 // functions for mid-edge nodes.
00177 
00178 // Linear quad weights at corners
00179 static double N0( double xi, double eta )
00180 {
00181     return l11( xi ) * l11( eta );
00182 }
00183 static double N1( double xi, double eta )
00184 {
00185     return l12( xi ) * l11( eta );
00186 }
00187 static double N2( double xi, double eta )
00188 {
00189     return l12( xi ) * l12( eta );
00190 }
00191 static double N3( double xi, double eta )
00192 {
00193     return l11( xi ) * l12( eta );
00194 }
00195 // Mixture of linear and quadratic weight functions for mid-side nodes
00196 static double N4( double xi, double eta )
00197 {
00198     return l22( xi ) * l11( eta );
00199 }
00200 static double N5( double xi, double eta )
00201 {
00202     return l12( xi ) * l22( eta );
00203 }
00204 static double N6( double xi, double eta )
00205 {
00206     return l22( xi ) * l12( eta );
00207 }
00208 static double N7( double xi, double eta )
00209 {
00210     return l11( xi ) * l22( eta );
00211 }
00212 // Quadratic element weight function for mid-element node
00213 static double N8( double xi, double eta )
00214 {
00215     return l22( xi ) * l22( eta );
00216 }
00217 // Derivatives of above weight functions wrt xi
00218 static double dN0dxi( double xi, double eta )
00219 {
00220     return dl11( xi ) * l11( eta );
00221 }
00222 static double dN1dxi( double xi, double eta )
00223 {
00224     return dl12( xi ) * l11( eta );
00225 }
00226 static double dN2dxi( double xi, double eta )
00227 {
00228     return dl12( xi ) * l12( eta );
00229 }
00230 static double dN3dxi( double xi, double eta )
00231 {
00232     return dl11( xi ) * l12( eta );
00233 }
00234 static double dN4dxi( double xi, double eta )
00235 {
00236     return dl22( xi ) * l11( eta );
00237 }
00238 static double dN5dxi( double xi, double eta )
00239 {
00240     return dl12( xi ) * l22( eta );
00241 }
00242 static double dN6dxi( double xi, double eta )
00243 {
00244     return dl22( xi ) * l12( eta );
00245 }
00246 static double dN7dxi( double xi, double eta )
00247 {
00248     return dl11( xi ) * l22( eta );
00249 }
00250 static double dN8dxi( double xi, double eta )
00251 {
00252     return dl22( xi ) * l22( eta );
00253 }
00254 // Derivatives of above weight functions wrt eta
00255 static double dN0deta( double xi, double eta )
00256 {
00257     return l11( xi ) * dl11( eta );
00258 }
00259 static double dN1deta( double xi, double eta )
00260 {
00261     return l12( xi ) * dl11( eta );
00262 }
00263 static double dN2deta( double xi, double eta )
00264 {
00265     return l12( xi ) * dl12( eta );
00266 }
00267 static double dN3deta( double xi, double eta )
00268 {
00269     return l11( xi ) * dl12( eta );
00270 }
00271 static double dN4deta( double xi, double eta )
00272 {
00273     return l22( xi ) * dl11( eta );
00274 }
00275 static double dN5deta( double xi, double eta )
00276 {
00277     return l12( xi ) * dl22( eta );
00278 }
00279 static double dN6deta( double xi, double eta )
00280 {
00281     return l22( xi ) * dl12( eta );
00282 }
00283 static double dN7deta( double xi, double eta )
00284 {
00285     return l11( xi ) * dl22( eta );
00286 }
00287 static double dN8deta( double xi, double eta )
00288 {
00289     return l22( xi ) * dl22( eta );
00290 }
00291 
00292 typedef double ( *N_t )( double, double );
00293 static const N_t N_a[9]      = { &N0, &N1, &N2, &N3, &N4, &N5, &N6, &N7, &N8 };
00294 static const N_t dNdxi_a[9]  = { &dN0dxi, &dN1dxi, &dN2dxi, &dN3dxi, &dN4dxi, &dN5dxi, &dN6dxi, &dN7dxi, &dN8dxi };
00295 static const N_t dNdeta_a[9] = { &dN0deta, &dN1deta, &dN2deta, &dN3deta, &dN4deta,
00296                                  &dN5deta, &dN6deta, &dN7deta, &dN8deta };
00297 
00298 // Mapping function coefficient functions (i<-vertex number)
00299 static double N( unsigned i, double xi, double eta )
00300 {
00301     return N_a[i]( xi, eta );
00302 }
00303 
00304 // Partial derivatives of mapping function coefficients
00305 static double dNdxi( unsigned i, double xi, double eta )
00306 {
00307     return dNdxi_a[i]( xi, eta );
00308 }
00309 static double dNdeta( unsigned i, double xi, double eta )
00310 {
00311     return dNdeta_a[i]( xi, eta );
00312 }
00313 
00314 // Evaluate one of N, dNdxi or dNdeta for each vertex
00315 typedef double ( *f_t )( unsigned, double, double );
00316 static void eval( f_t function, NodeSet nodeset, double xi, double eta, double results[9] )
00317 {
00318     // Initial values are a) linear values for corners,
00319     // b) values for mid edges assuming no mid-face node,
00320     // and c) the mid-face value for a fully quadratic element.
00321     for( unsigned i = 0; i < 9; ++i )
00322         results[i] = function( i, xi, eta );
00323 
00324     // if center node is present, adjust mid-edge coefficients
00325     if( !nodeset.mid_face_node( 0 ) ) { results[8] = 0; }
00326     else
00327     {
00328         for( unsigned i = 0; i < 4; ++i )
00329             results[i] -= 0.25 * results[8];
00330         for( unsigned i = 4; i < 8; ++i )
00331             results[i] -= 0.5 * results[8];
00332     }
00333 
00334     // if mid-edge nodes are present, adjust values for adjacent corners
00335     for( unsigned i = 0; i < 4; ++i )
00336     {
00337         if( !nodeset.mid_edge_node( i ) ) { results[i + 4] = 0.0; }
00338         else
00339         {
00340             results[i] -= 0.5 * results[i + 4];              // 1st adjacent corner
00341             results[( i + 1 ) % 4] -= 0.5 * results[i + 4];  // 2nd adjacent corner
00342         }
00343     }
00344 }
00345 
00346 // Finally, what all the above stuff was building up to:
00347 // functions to query mapping function.
00348 
00349 static void get_coeffs( NodeSet nodebits, double xi, double eta, double coeffs_out[9] )
00350 {
00351     eval( &N, nodebits, xi, eta, coeffs_out );
00352 }
00353 
00354 static void get_partial_wrt_xi( NodeSet nodebits, double xi, double eta, double derivs_out[9] )
00355 {
00356     eval( &dNdxi, nodebits, xi, eta, derivs_out );
00357 }
00358 
00359 static void get_partial_wrt_eta( NodeSet nodebits, double xi, double eta, double derivs_out[9] )
00360 {
00361     eval( &dNdeta, nodebits, xi, eta, derivs_out );
00362 }
00363 
00364 // Pre-defined sample points (xi and eta values)
00365 static const double corners[4][2] = { { min_xi, min_xi }, { max_xi, min_xi }, { max_xi, max_xi }, { min_xi, max_xi } };
00366 static const double midedge[4][2] = { { mid_xi, min_xi }, { max_xi, mid_xi }, { mid_xi, max_xi }, { min_xi, mid_xi } };
00367 static const double midelem[2]    = { mid_xi, mid_xi };
00368 
00369 static void check_valid_indices( const size_t* vertices, size_t num_vtx, NodeSet nodeset )
00370 {
00371     // check valid size of list (at least three, at most all nodes)
00372     CPPUNIT_ASSERT( num_vtx <= 9 );
00373     CPPUNIT_ASSERT( num_vtx >= 3 );
00374     // make sure vertex indices are valid (in [0,8])
00375     size_t vertcopy[9];
00376     std::copy( vertices, vertices + num_vtx, vertcopy );
00377     std::sort( vertcopy, vertcopy + num_vtx );
00378     CPPUNIT_ASSERT( vertcopy[num_vtx - 1] <= 8 );  // max value less than 9
00379                                                    // make sure there are no duplicates in the list
00380     const size_t* iter = std::unique( vertcopy, vertcopy + num_vtx );
00381     CPPUNIT_ASSERT( iter == vertcopy + num_vtx );
00382 
00383     // make all vertices are present in element
00384     for( unsigned i = 0; i < num_vtx; ++i )
00385     {
00386         if( vertcopy[i] == 8 )
00387             CPPUNIT_ASSERT( nodeset.mid_face_node( 0 ) );
00388         else if( vertcopy[i] >= 4 )
00389             CPPUNIT_ASSERT( nodeset.mid_edge_node( vertcopy[i] - 4 ) );
00390     }
00391 }
00392 
00393 static void check_no_zeros( const MsqVector< 2 >* derivs, size_t num_vtx )
00394 {
00395     for( unsigned i = 0; i < num_vtx; ++i )
00396     {
00397         double dxi  = derivs[i][0];
00398         double deta = derivs[i][1];
00399         CPPUNIT_ASSERT( fabs( dxi ) > 1e-6 || fabs( deta ) > 1e-6 );
00400     }
00401 }
00402 
00403 static void compare_coefficients( const double* coeffs, const size_t* indices, size_t num_coeff,
00404                                   const double* expected_coeffs, unsigned loc, NodeSet bits )
00405 {
00406     // find the location in the returned list for each node
00407     size_t revidx[9];
00408     double test_vals[9];
00409     for( size_t i = 0; i < 9; ++i )
00410     {
00411         revidx[i]    = std::find( indices, indices + num_coeff, i ) - indices;
00412         test_vals[i] = ( revidx[i] == num_coeff ) ? 0.0 : coeffs[revidx[i]];
00413     }
00414 
00415     // Check that index list doesn't contain any nodes not actually
00416     // present in the element.
00417     CPPUNIT_ASSERT( bits.mid_edge_node( 0 ) || ( revidx[4] == num_coeff ) );
00418     CPPUNIT_ASSERT( bits.mid_edge_node( 1 ) || ( revidx[5] == num_coeff ) );
00419     CPPUNIT_ASSERT( bits.mid_edge_node( 2 ) || ( revidx[6] == num_coeff ) );
00420     CPPUNIT_ASSERT( bits.mid_edge_node( 3 ) || ( revidx[7] == num_coeff ) );
00421     CPPUNIT_ASSERT( bits.mid_face_node( 0 ) || ( revidx[8] == num_coeff ) );
00422 
00423     // compare expected and actual coefficient values
00424     ASSERT_VALUES_EQUAL( expected_coeffs[0], test_vals[0], loc, bits );
00425     ASSERT_VALUES_EQUAL( expected_coeffs[1], test_vals[1], loc, bits );
00426     ASSERT_VALUES_EQUAL( expected_coeffs[2], test_vals[2], loc, bits );
00427     ASSERT_VALUES_EQUAL( expected_coeffs[3], test_vals[3], loc, bits );
00428     ASSERT_VALUES_EQUAL( expected_coeffs[4], test_vals[4], loc, bits );
00429     ASSERT_VALUES_EQUAL( expected_coeffs[5], test_vals[5], loc, bits );
00430     ASSERT_VALUES_EQUAL( expected_coeffs[6], test_vals[6], loc, bits );
00431     ASSERT_VALUES_EQUAL( expected_coeffs[7], test_vals[7], loc, bits );
00432     ASSERT_VALUES_EQUAL( expected_coeffs[8], test_vals[8], loc, bits );
00433 }
00434 
00435 static void compare_derivatives( const size_t* vertices, size_t num_vtx, const MsqVector< 2 >* derivs,
00436                                  const double* expected_dxi, const double* expected_deta, unsigned loc, NodeSet bits )
00437 {
00438     check_valid_indices( vertices, num_vtx, bits );
00439     check_no_zeros( derivs, num_vtx );
00440 
00441     // Input has values in dxi & deta only for nodes in 'vertices'
00442     // Convert to values for every possible node, with zero's for
00443     // nodes that are not present.
00444     CPPUNIT_ASSERT( num_vtx <= 9 );
00445     double expanded_dxi[9], expanded_deta[9];
00446     std::fill( expanded_dxi, expanded_dxi + 9, 0.0 );
00447     std::fill( expanded_deta, expanded_deta + 9, 0.0 );
00448     for( unsigned i = 0; i < num_vtx; ++i )
00449     {
00450         CPPUNIT_ASSERT( vertices[i] <= 9 );
00451         expanded_dxi[vertices[i]]  = derivs[i][0];
00452         expanded_deta[vertices[i]] = derivs[i][1];
00453     }
00454 
00455     ASSERT_VALUES_EQUAL( expected_dxi[0], expanded_dxi[0], loc, bits );
00456     ASSERT_VALUES_EQUAL( expected_dxi[1], expanded_dxi[1], loc, bits );
00457     ASSERT_VALUES_EQUAL( expected_dxi[2], expanded_dxi[2], loc, bits );
00458     ASSERT_VALUES_EQUAL( expected_dxi[3], expanded_dxi[3], loc, bits );
00459     ASSERT_VALUES_EQUAL( expected_dxi[4], expanded_dxi[4], loc, bits );
00460     ASSERT_VALUES_EQUAL( expected_dxi[5], expanded_dxi[5], loc, bits );
00461     ASSERT_VALUES_EQUAL( expected_dxi[6], expanded_dxi[6], loc, bits );
00462     ASSERT_VALUES_EQUAL( expected_dxi[7], expanded_dxi[7], loc, bits );
00463     ASSERT_VALUES_EQUAL( expected_dxi[8], expanded_dxi[8], loc, bits );
00464 
00465     ASSERT_VALUES_EQUAL( expected_deta[0], expanded_deta[0], loc, bits );
00466     ASSERT_VALUES_EQUAL( expected_deta[1], expanded_deta[1], loc, bits );
00467     ASSERT_VALUES_EQUAL( expected_deta[2], expanded_deta[2], loc, bits );
00468     ASSERT_VALUES_EQUAL( expected_deta[3], expanded_deta[3], loc, bits );
00469     ASSERT_VALUES_EQUAL( expected_deta[4], expanded_deta[4], loc, bits );
00470     ASSERT_VALUES_EQUAL( expected_deta[5], expanded_deta[5], loc, bits );
00471     ASSERT_VALUES_EQUAL( expected_deta[6], expanded_deta[6], loc, bits );
00472     ASSERT_VALUES_EQUAL( expected_deta[7], expanded_deta[7], loc, bits );
00473     ASSERT_VALUES_EQUAL( expected_deta[8], expanded_deta[8], loc, bits );
00474 }
00475 
00476 void QuadLagrangeShapeTest::test_corner_coeff( int corner, NodeSet nodebits )
00477 {
00478     MsqPrintError err( std::cout );
00479 
00480     double expected[9];
00481     get_coeffs( nodebits, corners[corner][XI], corners[corner][ETA], expected );
00482 
00483     double coeff[100];
00484     size_t num_coeff = 11, indices[100];
00485     sf.coefficients( Sample( 0, corner ), nodebits, coeff, indices, num_coeff, err );
00486     CPPUNIT_ASSERT( !err );
00487 
00488     compare_coefficients( coeff, indices, num_coeff, expected, corner, nodebits );
00489 }
00490 
00491 void QuadLagrangeShapeTest::test_edge_coeff( int edge, NodeSet nodebits )
00492 {
00493     MsqPrintError err( std::cout );
00494 
00495     double expected[9];
00496     get_coeffs( nodebits, midedge[edge][XI], midedge[edge][ETA], expected );
00497 
00498     double coeff[100];
00499     size_t num_coeff = 11, indices[100];
00500     sf.coefficients( Sample( 1, edge ), nodebits, coeff, indices, num_coeff, err );
00501     CPPUNIT_ASSERT( !err );
00502 
00503     compare_coefficients( coeff, indices, num_coeff, expected, edge + 4, nodebits );
00504 }
00505 
00506 void QuadLagrangeShapeTest::test_mid_coeff( NodeSet nodebits )
00507 {
00508     MsqPrintError err( std::cout );
00509 
00510     double expected[9];
00511     get_coeffs( nodebits, midelem[XI], midelem[ETA], expected );
00512 
00513     double coeff[100];
00514     size_t num_coeff = 11, indices[100];
00515     sf.coefficients( Sample( 2, 0 ), nodebits, coeff, indices, num_coeff, err );
00516     CPPUNIT_ASSERT( !err );
00517 
00518     compare_coefficients( coeff, indices, num_coeff, expected, 8, nodebits );
00519 }
00520 
00521 void QuadLagrangeShapeTest::test_corner_derivs( int corner, NodeSet nodebits )
00522 {
00523     MsqPrintError err( std::cout );
00524 
00525     double expected_dxi[9], expected_deta[9];
00526     get_partial_wrt_xi( nodebits, corners[corner][XI], corners[corner][ETA], expected_dxi );
00527     get_partial_wrt_eta( nodebits, corners[corner][XI], corners[corner][ETA], expected_deta );
00528 
00529     size_t vertices[100], num_vtx = 23;
00530     MsqVector< 2 > derivs[100];
00531     sf.derivatives( Sample( 0, corner ), nodebits, vertices, derivs, num_vtx, err );
00532     CPPUNIT_ASSERT( !err );
00533 
00534     compare_derivatives( vertices, num_vtx, derivs, expected_dxi, expected_deta, corner, nodebits );
00535 }
00536 
00537 void QuadLagrangeShapeTest::test_edge_derivs( int edge, NodeSet nodebits )
00538 {
00539     MsqPrintError err( std::cout );
00540 
00541     double expected_dxi[9], expected_deta[9];
00542     get_partial_wrt_xi( nodebits, midedge[edge][XI], midedge[edge][ETA], expected_dxi );
00543     get_partial_wrt_eta( nodebits, midedge[edge][XI], midedge[edge][ETA], expected_deta );
00544 
00545     size_t vertices[100], num_vtx = 23;
00546     MsqVector< 2 > derivs[100];
00547     sf.derivatives( Sample( 1, edge ), nodebits, vertices, derivs, num_vtx, err );
00548     CPPUNIT_ASSERT( !err );
00549 
00550     compare_derivatives( vertices, num_vtx, derivs, expected_dxi, expected_deta, edge + 4, nodebits );
00551 }
00552 
00553 void QuadLagrangeShapeTest::test_mid_derivs( NodeSet nodebits )
00554 {
00555     MsqPrintError err( std::cout );
00556 
00557     double expected_dxi[9], expected_deta[9];
00558     get_partial_wrt_xi( nodebits, midelem[XI], midelem[ETA], expected_dxi );
00559     get_partial_wrt_eta( nodebits, midelem[XI], midelem[ETA], expected_deta );
00560 
00561     size_t vertices[100], num_vtx = 23;
00562     MsqVector< 2 > derivs[100];
00563     sf.derivatives( Sample( 2, 0 ), nodebits, vertices, derivs, num_vtx, err );
00564     CPPUNIT_ASSERT( !err );
00565 
00566     compare_derivatives( vertices, num_vtx, derivs, expected_dxi, expected_deta, 8, nodebits );
00567 }
00568 
00569 static NodeSet nodeset_from_bits( unsigned bits )
00570 {
00571     NodeSet result;
00572     for( unsigned i = 0; i < 4; ++i )
00573         if( bits & ( 1 << i ) ) result.set_mid_edge_node( i );
00574     if( bits & ( 1 << 4 ) ) result.set_mid_face_node( 0 );
00575     return result;
00576 }
00577 
00578 void QuadLagrangeShapeTest::test_coeff_corners()
00579 {
00580     // for every possible combination of higher-order nodes
00581     // (0x1F = 11111 : five possible higher-order nodes in quad)
00582     for( unsigned j = 0; j <= 0x1Fu; ++j )
00583         // for every corner
00584         for( unsigned i = 0; i < 4; ++i )
00585             test_corner_coeff( i, nodeset_from_bits( j ) );
00586 }
00587 
00588 void QuadLagrangeShapeTest::test_coeff_edges()
00589 {
00590     // for every possible combination of higher-order nodes
00591     // (0x1F = 11111 : five possible higher-order nodes in quad)
00592     for( unsigned j = 0; j <= 0x1Fu; ++j )
00593         // for every edge
00594         for( unsigned i = 0; i < 4; ++i )
00595             test_edge_coeff( i, nodeset_from_bits( j ) );
00596 }
00597 
00598 void QuadLagrangeShapeTest::test_coeff_center()
00599 {
00600     // for every possible combination of higher-order nodes
00601     // (0x1F = 11111 : five possible higher-order nodes in quad)
00602     for( unsigned j = 0; j <= 0x1Fu; ++j )
00603         test_mid_coeff( nodeset_from_bits( j ) );
00604 }
00605 
00606 void QuadLagrangeShapeTest::test_deriv_corners()
00607 {
00608     // for every possible combination of higher-order nodes
00609     // (0x1F = 11111 : five possible higher-order nodes in quad)
00610     for( unsigned j = 0; j <= 0x1Fu; ++j )
00611         // for every corner
00612         for( unsigned i = 0; i < 4; ++i )
00613             test_corner_derivs( i, nodeset_from_bits( j ) );
00614 }
00615 
00616 void QuadLagrangeShapeTest::test_deriv_edges()
00617 {
00618     // for every possible combination of higher-order nodes
00619     // (0x1F = 11111 : five possible higher-order nodes in quad)
00620     for( unsigned j = 0; j <= 0x1Fu; ++j )
00621         // for every edge
00622         for( unsigned i = 0; i < 4; ++i )
00623             test_edge_derivs( i, nodeset_from_bits( j ) );
00624 }
00625 
00626 void QuadLagrangeShapeTest::test_deriv_center()
00627 {
00628     // for every possible combination of higher-order nodes
00629     // (0x1F = 11111 : five possible higher-order nodes in quad)
00630     for( unsigned j = 0; j <= 0x1Fu; ++j )
00631         test_mid_derivs( nodeset_from_bits( j ) );
00632 }
00633 
00634 void QuadLagrangeShapeTest::test_ideal_jacobian()
00635 {
00636     MsqError err;
00637     MsqMatrix< 3, 2 > J_prime;
00638     sf.ideal( Sample( 2, 0 ), J_prime, err );
00639     ASSERT_NO_ERROR( err );
00640 
00641     // for this test that everything is in the xy-plane
00642     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, J_prime( 2, 0 ), 1e-12 );
00643     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, J_prime( 2, 1 ), 1e-12 );
00644     MsqMatrix< 2, 2 > J_act( J_prime.data() );
00645     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( J_act ), 1e-6 );
00646 
00647     const Vector3D* verts = unit_edge_element( QUADRILATERAL );
00648     CPPUNIT_ASSERT( verts );
00649 
00650     JacobianCalculator jc;
00651     jc.get_Jacobian_2D( &sf, NodeSet(), Sample( 2, 0 ), verts, 4, J_prime, err );
00652     ASSERT_NO_ERROR( err );
00653 
00654     // for this test that everything is in the xy-plane
00655     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, J_prime( 2, 0 ), 1e-12 );
00656     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, J_prime( 2, 1 ), 1e-12 );
00657     MsqMatrix< 2, 2 > J_exp( J_prime.data() );
00658     J_exp /= sqrt( det( J_exp ) );
00659 
00660     // Matrices should be a rotation of each other.
00661     // First, calculate tentative rotation matrix
00662     MsqMatrix< 2, 2 > R = inverse( J_exp ) * J_act;
00663     // next check that it is a rotation
00664     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( R ), 1e-6 );          // no scaling
00665     ASSERT_MATRICES_EQUAL( transpose( R ), inverse( R ), 1e-6 );  // orthogonal
00666 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines