MOAB: Mesh Oriented datABase  (version 5.4.0)
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 ) )
00326     {
00327         results[8] = 0;
00328     }
00329     else
00330     {
00331         for( unsigned i = 0; i < 4; ++i )
00332             results[i] -= 0.25 * results[8];
00333         for( unsigned i = 4; i < 8; ++i )
00334             results[i] -= 0.5 * results[8];
00335     }
00336 
00337     // if mid-edge nodes are present, adjust values for adjacent corners
00338     for( unsigned i = 0; i < 4; ++i )
00339     {
00340         if( !nodeset.mid_edge_node( i ) )
00341         {
00342             results[i + 4] = 0.0;
00343         }
00344         else
00345         {
00346             results[i] -= 0.5 * results[i + 4];              // 1st adjacent corner
00347             results[( i + 1 ) % 4] -= 0.5 * results[i + 4];  // 2nd adjacent corner
00348         }
00349     }
00350 }
00351 
00352 // Finally, what all the above stuff was building up to:
00353 // functions to query mapping function.
00354 
00355 static void get_coeffs( NodeSet nodebits, double xi, double eta, double coeffs_out[9] )
00356 {
00357     eval( &N, nodebits, xi, eta, coeffs_out );
00358 }
00359 
00360 static void get_partial_wrt_xi( NodeSet nodebits, double xi, double eta, double derivs_out[9] )
00361 {
00362     eval( &dNdxi, nodebits, xi, eta, derivs_out );
00363 }
00364 
00365 static void get_partial_wrt_eta( NodeSet nodebits, double xi, double eta, double derivs_out[9] )
00366 {
00367     eval( &dNdeta, nodebits, xi, eta, derivs_out );
00368 }
00369 
00370 // Pre-defined sample points (xi and eta values)
00371 static const double corners[4][2] = { { min_xi, min_xi }, { max_xi, min_xi }, { max_xi, max_xi }, { min_xi, max_xi } };
00372 static const double midedge[4][2] = { { mid_xi, min_xi }, { max_xi, mid_xi }, { mid_xi, max_xi }, { min_xi, mid_xi } };
00373 static const double midelem[2]    = { mid_xi, mid_xi };
00374 
00375 static void check_valid_indices( const size_t* vertices, size_t num_vtx, NodeSet nodeset )
00376 {
00377     // check valid size of list (at least three, at most all nodes)
00378     CPPUNIT_ASSERT( num_vtx <= 9 );
00379     CPPUNIT_ASSERT( num_vtx >= 3 );
00380     // make sure vertex indices are valid (in [0,8])
00381     size_t vertcopy[9];
00382     std::copy( vertices, vertices + num_vtx, vertcopy );
00383     std::sort( vertcopy, vertcopy + num_vtx );
00384     CPPUNIT_ASSERT( vertcopy[num_vtx - 1] <= 8 );  // max value less than 9
00385                                                    // make sure there are no duplicates in the list
00386     const size_t* iter = std::unique( vertcopy, vertcopy + num_vtx );
00387     CPPUNIT_ASSERT( iter == vertcopy + num_vtx );
00388 
00389     // make all vertices are present in element
00390     for( unsigned i = 0; i < num_vtx; ++i )
00391     {
00392         if( vertcopy[i] == 8 )
00393             CPPUNIT_ASSERT( nodeset.mid_face_node( 0 ) );
00394         else if( vertcopy[i] >= 4 )
00395             CPPUNIT_ASSERT( nodeset.mid_edge_node( vertcopy[i] - 4 ) );
00396     }
00397 }
00398 
00399 static void check_no_zeros( const MsqVector< 2 >* derivs, size_t num_vtx )
00400 {
00401     for( unsigned i = 0; i < num_vtx; ++i )
00402     {
00403         double dxi  = derivs[i][0];
00404         double deta = derivs[i][1];
00405         CPPUNIT_ASSERT( fabs( dxi ) > 1e-6 || fabs( deta ) > 1e-6 );
00406     }
00407 }
00408 
00409 static void compare_coefficients( const double* coeffs,
00410                                   const size_t* indices,
00411                                   size_t num_coeff,
00412                                   const double* expected_coeffs,
00413                                   unsigned loc,
00414                                   NodeSet bits )
00415 {
00416     // find the location in the returned list for each node
00417     size_t revidx[9];
00418     double test_vals[9];
00419     for( size_t i = 0; i < 9; ++i )
00420     {
00421         revidx[i]    = std::find( indices, indices + num_coeff, i ) - indices;
00422         test_vals[i] = ( revidx[i] == num_coeff ) ? 0.0 : coeffs[revidx[i]];
00423     }
00424 
00425     // Check that index list doesn't contain any nodes not actually
00426     // present in the element.
00427     CPPUNIT_ASSERT( bits.mid_edge_node( 0 ) || ( revidx[4] == num_coeff ) );
00428     CPPUNIT_ASSERT( bits.mid_edge_node( 1 ) || ( revidx[5] == num_coeff ) );
00429     CPPUNIT_ASSERT( bits.mid_edge_node( 2 ) || ( revidx[6] == num_coeff ) );
00430     CPPUNIT_ASSERT( bits.mid_edge_node( 3 ) || ( revidx[7] == num_coeff ) );
00431     CPPUNIT_ASSERT( bits.mid_face_node( 0 ) || ( revidx[8] == num_coeff ) );
00432 
00433     // compare expected and actual coefficient values
00434     ASSERT_VALUES_EQUAL( expected_coeffs[0], test_vals[0], loc, bits );
00435     ASSERT_VALUES_EQUAL( expected_coeffs[1], test_vals[1], loc, bits );
00436     ASSERT_VALUES_EQUAL( expected_coeffs[2], test_vals[2], loc, bits );
00437     ASSERT_VALUES_EQUAL( expected_coeffs[3], test_vals[3], loc, bits );
00438     ASSERT_VALUES_EQUAL( expected_coeffs[4], test_vals[4], loc, bits );
00439     ASSERT_VALUES_EQUAL( expected_coeffs[5], test_vals[5], loc, bits );
00440     ASSERT_VALUES_EQUAL( expected_coeffs[6], test_vals[6], loc, bits );
00441     ASSERT_VALUES_EQUAL( expected_coeffs[7], test_vals[7], loc, bits );
00442     ASSERT_VALUES_EQUAL( expected_coeffs[8], test_vals[8], loc, bits );
00443 }
00444 
00445 static void compare_derivatives( const size_t* vertices,
00446                                  size_t num_vtx,
00447                                  const MsqVector< 2 >* derivs,
00448                                  const double* expected_dxi,
00449                                  const double* expected_deta,
00450                                  unsigned loc,
00451                                  NodeSet bits )
00452 {
00453     check_valid_indices( vertices, num_vtx, bits );
00454     check_no_zeros( derivs, num_vtx );
00455 
00456     // Input has values in dxi & deta only for nodes in 'vertices'
00457     // Convert to values for every possible node, with zero's for
00458     // nodes that are not present.
00459     CPPUNIT_ASSERT( num_vtx <= 9 );
00460     double expanded_dxi[9], expanded_deta[9];
00461     std::fill( expanded_dxi, expanded_dxi + 9, 0.0 );
00462     std::fill( expanded_deta, expanded_deta + 9, 0.0 );
00463     for( unsigned i = 0; i < num_vtx; ++i )
00464     {
00465         CPPUNIT_ASSERT( vertices[i] <= 9 );
00466         expanded_dxi[vertices[i]]  = derivs[i][0];
00467         expanded_deta[vertices[i]] = derivs[i][1];
00468     }
00469 
00470     ASSERT_VALUES_EQUAL( expected_dxi[0], expanded_dxi[0], loc, bits );
00471     ASSERT_VALUES_EQUAL( expected_dxi[1], expanded_dxi[1], loc, bits );
00472     ASSERT_VALUES_EQUAL( expected_dxi[2], expanded_dxi[2], loc, bits );
00473     ASSERT_VALUES_EQUAL( expected_dxi[3], expanded_dxi[3], loc, bits );
00474     ASSERT_VALUES_EQUAL( expected_dxi[4], expanded_dxi[4], loc, bits );
00475     ASSERT_VALUES_EQUAL( expected_dxi[5], expanded_dxi[5], loc, bits );
00476     ASSERT_VALUES_EQUAL( expected_dxi[6], expanded_dxi[6], loc, bits );
00477     ASSERT_VALUES_EQUAL( expected_dxi[7], expanded_dxi[7], loc, bits );
00478     ASSERT_VALUES_EQUAL( expected_dxi[8], expanded_dxi[8], loc, bits );
00479 
00480     ASSERT_VALUES_EQUAL( expected_deta[0], expanded_deta[0], loc, bits );
00481     ASSERT_VALUES_EQUAL( expected_deta[1], expanded_deta[1], loc, bits );
00482     ASSERT_VALUES_EQUAL( expected_deta[2], expanded_deta[2], loc, bits );
00483     ASSERT_VALUES_EQUAL( expected_deta[3], expanded_deta[3], loc, bits );
00484     ASSERT_VALUES_EQUAL( expected_deta[4], expanded_deta[4], loc, bits );
00485     ASSERT_VALUES_EQUAL( expected_deta[5], expanded_deta[5], loc, bits );
00486     ASSERT_VALUES_EQUAL( expected_deta[6], expanded_deta[6], loc, bits );
00487     ASSERT_VALUES_EQUAL( expected_deta[7], expanded_deta[7], loc, bits );
00488     ASSERT_VALUES_EQUAL( expected_deta[8], expanded_deta[8], loc, bits );
00489 }
00490 
00491 void QuadLagrangeShapeTest::test_corner_coeff( int corner, NodeSet nodebits )
00492 {
00493     MsqPrintError err( std::cout );
00494 
00495     double expected[9];
00496     get_coeffs( nodebits, corners[corner][XI], corners[corner][ETA], expected );
00497 
00498     double coeff[100];
00499     size_t num_coeff = 11, indices[100];
00500     sf.coefficients( Sample( 0, corner ), nodebits, coeff, indices, num_coeff, err );
00501     CPPUNIT_ASSERT( !err );
00502 
00503     compare_coefficients( coeff, indices, num_coeff, expected, corner, nodebits );
00504 }
00505 
00506 void QuadLagrangeShapeTest::test_edge_coeff( int edge, NodeSet nodebits )
00507 {
00508     MsqPrintError err( std::cout );
00509 
00510     double expected[9];
00511     get_coeffs( nodebits, midedge[edge][XI], midedge[edge][ETA], expected );
00512 
00513     double coeff[100];
00514     size_t num_coeff = 11, indices[100];
00515     sf.coefficients( Sample( 1, edge ), nodebits, coeff, indices, num_coeff, err );
00516     CPPUNIT_ASSERT( !err );
00517 
00518     compare_coefficients( coeff, indices, num_coeff, expected, edge + 4, nodebits );
00519 }
00520 
00521 void QuadLagrangeShapeTest::test_mid_coeff( NodeSet nodebits )
00522 {
00523     MsqPrintError err( std::cout );
00524 
00525     double expected[9];
00526     get_coeffs( nodebits, midelem[XI], midelem[ETA], expected );
00527 
00528     double coeff[100];
00529     size_t num_coeff = 11, indices[100];
00530     sf.coefficients( Sample( 2, 0 ), nodebits, coeff, indices, num_coeff, err );
00531     CPPUNIT_ASSERT( !err );
00532 
00533     compare_coefficients( coeff, indices, num_coeff, expected, 8, nodebits );
00534 }
00535 
00536 void QuadLagrangeShapeTest::test_corner_derivs( int corner, NodeSet nodebits )
00537 {
00538     MsqPrintError err( std::cout );
00539 
00540     double expected_dxi[9], expected_deta[9];
00541     get_partial_wrt_xi( nodebits, corners[corner][XI], corners[corner][ETA], expected_dxi );
00542     get_partial_wrt_eta( nodebits, corners[corner][XI], corners[corner][ETA], expected_deta );
00543 
00544     size_t vertices[100], num_vtx = 23;
00545     MsqVector< 2 > derivs[100];
00546     sf.derivatives( Sample( 0, corner ), nodebits, vertices, derivs, num_vtx, err );
00547     CPPUNIT_ASSERT( !err );
00548 
00549     compare_derivatives( vertices, num_vtx, derivs, expected_dxi, expected_deta, corner, nodebits );
00550 }
00551 
00552 void QuadLagrangeShapeTest::test_edge_derivs( int edge, NodeSet nodebits )
00553 {
00554     MsqPrintError err( std::cout );
00555 
00556     double expected_dxi[9], expected_deta[9];
00557     get_partial_wrt_xi( nodebits, midedge[edge][XI], midedge[edge][ETA], expected_dxi );
00558     get_partial_wrt_eta( nodebits, midedge[edge][XI], midedge[edge][ETA], expected_deta );
00559 
00560     size_t vertices[100], num_vtx = 23;
00561     MsqVector< 2 > derivs[100];
00562     sf.derivatives( Sample( 1, edge ), nodebits, vertices, derivs, num_vtx, err );
00563     CPPUNIT_ASSERT( !err );
00564 
00565     compare_derivatives( vertices, num_vtx, derivs, expected_dxi, expected_deta, edge + 4, nodebits );
00566 }
00567 
00568 void QuadLagrangeShapeTest::test_mid_derivs( NodeSet nodebits )
00569 {
00570     MsqPrintError err( std::cout );
00571 
00572     double expected_dxi[9], expected_deta[9];
00573     get_partial_wrt_xi( nodebits, midelem[XI], midelem[ETA], expected_dxi );
00574     get_partial_wrt_eta( nodebits, midelem[XI], midelem[ETA], expected_deta );
00575 
00576     size_t vertices[100], num_vtx = 23;
00577     MsqVector< 2 > derivs[100];
00578     sf.derivatives( Sample( 2, 0 ), nodebits, vertices, derivs, num_vtx, err );
00579     CPPUNIT_ASSERT( !err );
00580 
00581     compare_derivatives( vertices, num_vtx, derivs, expected_dxi, expected_deta, 8, nodebits );
00582 }
00583 
00584 static NodeSet nodeset_from_bits( unsigned bits )
00585 {
00586     NodeSet result;
00587     for( unsigned i = 0; i < 4; ++i )
00588         if( bits & ( 1 << i ) ) result.set_mid_edge_node( i );
00589     if( bits & ( 1 << 4 ) ) result.set_mid_face_node( 0 );
00590     return result;
00591 }
00592 
00593 void QuadLagrangeShapeTest::test_coeff_corners()
00594 {
00595     // for every possible combination of higher-order nodes
00596     // (0x1F = 11111 : five possible higher-order nodes in quad)
00597     for( unsigned j = 0; j <= 0x1Fu; ++j )
00598         // for every corner
00599         for( unsigned i = 0; i < 4; ++i )
00600             test_corner_coeff( i, nodeset_from_bits( j ) );
00601 }
00602 
00603 void QuadLagrangeShapeTest::test_coeff_edges()
00604 {
00605     // for every possible combination of higher-order nodes
00606     // (0x1F = 11111 : five possible higher-order nodes in quad)
00607     for( unsigned j = 0; j <= 0x1Fu; ++j )
00608         // for every edge
00609         for( unsigned i = 0; i < 4; ++i )
00610             test_edge_coeff( i, nodeset_from_bits( j ) );
00611 }
00612 
00613 void QuadLagrangeShapeTest::test_coeff_center()
00614 {
00615     // for every possible combination of higher-order nodes
00616     // (0x1F = 11111 : five possible higher-order nodes in quad)
00617     for( unsigned j = 0; j <= 0x1Fu; ++j )
00618         test_mid_coeff( nodeset_from_bits( j ) );
00619 }
00620 
00621 void QuadLagrangeShapeTest::test_deriv_corners()
00622 {
00623     // for every possible combination of higher-order nodes
00624     // (0x1F = 11111 : five possible higher-order nodes in quad)
00625     for( unsigned j = 0; j <= 0x1Fu; ++j )
00626         // for every corner
00627         for( unsigned i = 0; i < 4; ++i )
00628             test_corner_derivs( i, nodeset_from_bits( j ) );
00629 }
00630 
00631 void QuadLagrangeShapeTest::test_deriv_edges()
00632 {
00633     // for every possible combination of higher-order nodes
00634     // (0x1F = 11111 : five possible higher-order nodes in quad)
00635     for( unsigned j = 0; j <= 0x1Fu; ++j )
00636         // for every edge
00637         for( unsigned i = 0; i < 4; ++i )
00638             test_edge_derivs( i, nodeset_from_bits( j ) );
00639 }
00640 
00641 void QuadLagrangeShapeTest::test_deriv_center()
00642 {
00643     // for every possible combination of higher-order nodes
00644     // (0x1F = 11111 : five possible higher-order nodes in quad)
00645     for( unsigned j = 0; j <= 0x1Fu; ++j )
00646         test_mid_derivs( nodeset_from_bits( j ) );
00647 }
00648 
00649 void QuadLagrangeShapeTest::test_ideal_jacobian()
00650 {
00651     MsqError err;
00652     MsqMatrix< 3, 2 > J_prime;
00653     sf.ideal( Sample( 2, 0 ), J_prime, err );
00654     ASSERT_NO_ERROR( err );
00655 
00656     // for this test that everything is in the xy-plane
00657     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, J_prime( 2, 0 ), 1e-12 );
00658     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, J_prime( 2, 1 ), 1e-12 );
00659     MsqMatrix< 2, 2 > J_act( J_prime.data() );
00660     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( J_act ), 1e-6 );
00661 
00662     const Vector3D* verts = unit_edge_element( QUADRILATERAL );
00663     CPPUNIT_ASSERT( verts );
00664 
00665     JacobianCalculator jc;
00666     jc.get_Jacobian_2D( &sf, NodeSet(), Sample( 2, 0 ), verts, 4, J_prime, err );
00667     ASSERT_NO_ERROR( err );
00668 
00669     // for this test that everything is in the xy-plane
00670     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, J_prime( 2, 0 ), 1e-12 );
00671     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, J_prime( 2, 1 ), 1e-12 );
00672     MsqMatrix< 2, 2 > J_exp( J_prime.data() );
00673     J_exp /= sqrt( det( J_exp ) );
00674 
00675     // Matrices should be a rotation of each other.
00676     // First, calculate tentative rotation matrix
00677     MsqMatrix< 2, 2 > R = inverse( J_exp ) * J_act;
00678     // next check that it is a rotation
00679     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( R ), 1e-6 );          // no scaling
00680     ASSERT_MATRICES_EQUAL( transpose( R ), inverse( R ), 1e-6 );  // orthogonal
00681 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines