MOAB: Mesh Oriented datABase  (version 5.2.1)
TetLagrangeShapeTest.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 TetLagrangeShapeTest.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "TetLagrangeShape.hpp"
00034 #include "TopologyInfo.hpp"
00035 #include "MsqError.hpp"
00036 #include "IdealElements.hpp"
00037 #include "JacobianCalculator.hpp"
00038 #include <math.h>
00039 
00040 #include "UnitUtil.hpp"
00041 
00042 #include <vector>
00043 #include <algorithm>
00044 #include <iostream>
00045 #include <sstream>
00046 
00047 using namespace MBMesquite;
00048 using namespace std;
00049 
00050 const double epsilon = 1e-6;
00051 #define ASSERT_VALUES_EQUAL( v1, v2, location, bits ) \
00052     ASSERT_MESSAGE( value_message( ( location ), ( bits ), ( v1 ), ( v2 ) ), ( fabs( ( v1 ) - ( v2 ) ) < epsilon ) )
00053 
00054 static inline CppUnit::Message value_message( unsigned location, NodeSet bits, double v1, double v2 )
00055 {
00056     CppUnit::Message m( "equality assertion failed" );
00057 
00058     std::ostringstream buffer1;
00059     buffer1 << "Expected : " << v1;
00060     m.addDetail( buffer1.str() );
00061 
00062     std::ostringstream buffer2;
00063     buffer2 << "Actual   : " << v2;
00064     m.addDetail( buffer2.str() );
00065 
00066     std::ostringstream buffer3;
00067     buffer3 << "Location : ";
00068     if( location < 4 )
00069         buffer3 << "Corner " << location;
00070     else if( location < 10 )
00071         buffer3 << "Edge " << location - 4;
00072     else if( location < 14 )
00073         buffer3 << "Face " << location - 10;
00074     else if( location == 14 )
00075         buffer3 << "Mid-element";
00076     else
00077         buffer3 << "INVALID!!";
00078     m.addDetail( buffer3.str() );
00079 
00080     std::ostringstream buffer4;
00081     buffer4 << "Node Bits: " << bits;
00082     m.addDetail( buffer4.str() );
00083     return m;
00084 }
00085 
00086 class TetLagrangeShapeTest : public CppUnit::TestFixture
00087 {
00088   private:
00089     CPPUNIT_TEST_SUITE( TetLagrangeShapeTest );
00090 
00091     CPPUNIT_TEST( test_coeff_corners );
00092     CPPUNIT_TEST( test_coeff_edges );
00093     CPPUNIT_TEST( test_coeff_faces );
00094     CPPUNIT_TEST( test_coeff_center );
00095 
00096     CPPUNIT_TEST( test_deriv_corners );
00097     CPPUNIT_TEST( test_deriv_edges );
00098     CPPUNIT_TEST( test_deriv_faces );
00099     CPPUNIT_TEST( test_deriv_center );
00100 
00101     CPPUNIT_TEST( test_mid_elem_node_coeff );
00102     CPPUNIT_TEST( test_mid_elem_node_deriv );
00103     CPPUNIT_TEST( test_mid_face_node_coeff );
00104     CPPUNIT_TEST( test_mid_face_node_deriv );
00105 
00106     CPPUNIT_TEST( test_ideal_jacobian );
00107 
00108     CPPUNIT_TEST_SUITE_END();
00109 
00110     TetLagrangeShape sf;
00111 
00112     void test_corner_coeff( int corner, NodeSet nodeset );
00113     void test_edge_coeff( int edge, NodeSet nodeset );
00114     void test_face_coeff( int face, NodeSet nodeset );
00115     void test_mid_coeff( NodeSet nodebits );
00116 
00117     void test_corner_derivs( int corner, NodeSet nodeset );
00118     void test_edge_derivs( int edge, NodeSet nodeset );
00119     void test_face_derivs( int face, NodeSet nodeset );
00120     void test_mid_derivs( NodeSet nodeset );
00121 
00122     void test_invalid_nodebits_coeff( NodeSet nodeset );
00123     void test_invalid_nodebits_deriv( NodeSet nodeset );
00124 
00125   public:
00126     void test_coeff_corners();
00127     void test_coeff_edges();
00128     void test_coeff_faces();
00129     void test_coeff_center();
00130 
00131     void test_deriv_corners();
00132     void test_deriv_edges();
00133     void test_deriv_faces();
00134     void test_deriv_center();
00135 
00136     void test_mid_elem_node_coeff();
00137     void test_mid_elem_node_deriv();
00138     void test_mid_face_node_coeff();
00139     void test_mid_face_node_deriv();
00140 
00141     void test_ideal_jacobian();
00142 };
00143 
00144 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( TetLagrangeShapeTest, "TetLagrangeShapeTest" );
00145 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( TetLagrangeShapeTest, "Unit" );
00146 
00147 static double N0( double r, double s, double t )
00148 {
00149     double u = 1 - r - s - t;
00150     return u * ( 2 * u - 1 );
00151 }
00152 static double N1( double r, double, double )
00153 {
00154     return r * ( 2 * r - 1 );
00155 }
00156 static double N2( double, double s, double )
00157 {
00158     return s * ( 2 * s - 1 );
00159 }
00160 static double N3( double, double, double t )
00161 {
00162     return t * ( 2 * t - 1 );
00163 }
00164 static double N4( double r, double s, double t )
00165 {
00166     double u = 1 - r - s - t;
00167     return 4 * r * u;
00168 }
00169 static double N5( double r, double s, double )
00170 {
00171     return 4 * r * s;
00172 }
00173 static double N6( double r, double s, double t )
00174 {
00175     double u = 1 - r - s - t;
00176     return 4 * s * u;
00177 }
00178 static double N7( double r, double s, double t )
00179 {
00180     double u = 1 - r - s - t;
00181     return 4 * u * t;
00182 }
00183 static double N8( double r, double, double t )
00184 {
00185     return 4 * r * t;
00186 }
00187 static double N9( double, double s, double t )
00188 {
00189     return 4 * s * t;
00190 }
00191 
00192 static double dN0dr( double r, double s, double t )
00193 {
00194     double u = 1 - r - s - t;
00195     return 1 - 4 * u;
00196 }
00197 static double dN0ds( double r, double s, double t )
00198 {
00199     double u = 1 - r - s - t;
00200     return 1 - 4 * u;
00201 }
00202 static double dN0dt( double r, double s, double t )
00203 {
00204     double u = 1 - r - s - t;
00205     return 1 - 4 * u;
00206 }
00207 
00208 static double dN1dr( double r, double, double )
00209 {
00210     return 4 * r - 1;
00211 }
00212 static double dN1ds( double, double, double )
00213 {
00214     return 0;
00215 }
00216 static double dN1dt( double, double, double )
00217 {
00218     return 0;
00219 }
00220 
00221 static double dN2dr( double, double, double )
00222 {
00223     return 0;
00224 }
00225 static double dN2ds( double, double s, double )
00226 {
00227     return 4 * s - 1;
00228 }
00229 static double dN2dt( double, double, double )
00230 {
00231     return 0;
00232 }
00233 
00234 static double dN3dr( double, double, double )
00235 {
00236     return 0;
00237 }
00238 static double dN3ds( double, double, double )
00239 {
00240     return 0;
00241 }
00242 static double dN3dt( double, double, double t )
00243 {
00244     return 4 * t - 1;
00245 }
00246 
00247 static double dN4dr( double r, double s, double t )
00248 {
00249     double u = 1 - r - s - t;
00250     return 4 * ( u - r );
00251 }
00252 static double dN4ds( double r, double, double )
00253 {
00254     return -4 * r;
00255 }
00256 static double dN4dt( double r, double, double )
00257 {
00258     return -4 * r;
00259 }
00260 
00261 static double dN5dr( double, double s, double )
00262 {
00263     return 4 * s;
00264 }
00265 static double dN5ds( double r, double, double )
00266 {
00267     return 4 * r;
00268 }
00269 static double dN5dt( double, double, double )
00270 {
00271     return 0;
00272 }
00273 
00274 static double dN6dr( double, double s, double )
00275 {
00276     return -4 * s;
00277 }
00278 static double dN6ds( double r, double s, double t )
00279 {
00280     double u = 1 - r - s - t;
00281     return 4 * ( u - s );
00282 }
00283 static double dN6dt( double, double s, double )
00284 {
00285     return -4 * s;
00286 }
00287 
00288 static double dN7dr( double, double, double t )
00289 {
00290     return -4 * t;
00291 }
00292 static double dN7ds( double, double, double t )
00293 {
00294     return -4 * t;
00295 }
00296 static double dN7dt( double r, double s, double t )
00297 {
00298     double u = 1 - r - s - t;
00299     return 4 * ( u - t );
00300 }
00301 
00302 static double dN8dr( double, double, double t )
00303 {
00304     return 4 * t;
00305 }
00306 static double dN8ds( double, double, double )
00307 {
00308     return 0;
00309 }
00310 static double dN8dt( double r, double, double )
00311 {
00312     return 4 * r;
00313 }
00314 
00315 static double dN9dr( double, double, double )
00316 {
00317     return 0;
00318 }
00319 static double dN9ds( double, double, double t )
00320 {
00321     return 4 * t;
00322 }
00323 static double dN9dt( double, double s, double )
00324 {
00325     return 4 * s;
00326 }
00327 
00328 typedef double ( *N_t )( double, double, double );
00329 static const N_t N[]    = { &N0, &N1, &N2, &N3, &N4, &N5, &N6, &N7, &N8, &N9 };
00330 static const N_t dNdr[] = { &dN0dr, &dN1dr, &dN2dr, &dN3dr, &dN4dr, &dN5dr, &dN6dr, &dN7dr, &dN8dr, &dN9dr };
00331 static const N_t dNds[] = { &dN0ds, &dN1ds, &dN2ds, &dN3ds, &dN4ds, &dN5ds, &dN6ds, &dN7ds, &dN8ds, &dN9ds };
00332 static const N_t dNdt[] = { &dN0dt, &dN1dt, &dN2dt, &dN3dt, &dN4dt, &dN5dt, &dN6dt, &dN7dt, &dN8dt, &dN9dt };
00333 
00334 static const double rst_corner[][3] = { { 0, 0, 0 }, { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } };
00335 static const double rst_edge[][3]   = { { 0.5, 0.0, 0.0 }, { 0.5, 0.5, 0.0 }, { 0.0, 0.5, 0.0 },
00336                                       { 0.0, 0.0, 0.5 }, { 0.5, 0.0, 0.5 }, { 0.0, 0.5, 0.5 } };
00337 static const double rst_face[][3]   = { { 1. / 3, 0.00, 1. / 3 },
00338                                       { 1. / 3, 1. / 3, 1. / 3 },
00339                                       { 0.00, 1. / 3, 1. / 3 },
00340                                       { 1. / 3, 1. / 3, 0.00 } };
00341 static const double rst_mid[3]      = { 0.25, 0.25, 0.25 };
00342 
00343 static unsigned edges[][2] = { { 0, 1 }, { 1, 2 }, { 2, 0 }, { 0, 3 }, { 1, 3 }, { 2, 3 } };
00344 
00345 static void get_coeff( NodeSet nodeset, const double* rst, double* coeffs )
00346 {
00347     for( int i = 0; i < 10; ++i )
00348         coeffs[i] = ( *N[i] )( rst[0], rst[1], rst[2] );
00349     for( int i = 0; i < 6; ++i )
00350         if( !nodeset.mid_edge_node( i ) )
00351         {
00352             coeffs[edges[i][0]] += 0.5 * coeffs[i + 4];
00353             coeffs[edges[i][1]] += 0.5 * coeffs[i + 4];
00354             coeffs[i + 4] = 0;
00355         }
00356 }
00357 
00358 static void get_derivs( NodeSet nodeset, const double* rst, MsqVector< 3 >* derivs )
00359 {
00360     for( int i = 0; i < 10; ++i )
00361     {
00362         derivs[i][0] = ( *dNdr[i] )( rst[0], rst[1], rst[2] );
00363         derivs[i][1] = ( *dNds[i] )( rst[0], rst[1], rst[2] );
00364         derivs[i][2] = ( *dNdt[i] )( rst[0], rst[1], rst[2] );
00365     }
00366     for( int i = 0; i < 6; ++i )
00367         if( !nodeset.mid_edge_node( i ) )
00368         {
00369             int j = edges[i][0];
00370             derivs[j][0] += 0.5 * derivs[i + 4][0];
00371             derivs[j][1] += 0.5 * derivs[i + 4][1];
00372             derivs[j][2] += 0.5 * derivs[i + 4][2];
00373             j = edges[i][1];
00374             derivs[j][0] += 0.5 * derivs[i + 4][0];
00375             derivs[j][1] += 0.5 * derivs[i + 4][1];
00376             derivs[j][2] += 0.5 * derivs[i + 4][2];
00377             derivs[i + 4][0] = 0.0;
00378             derivs[i + 4][1] = 0.0;
00379             derivs[i + 4][2] = 0.0;
00380         }
00381 }
00382 
00383 static void check_valid_indices( const size_t* vertices, size_t num_vtx )
00384 {
00385     CPPUNIT_ASSERT( num_vtx < 11 );
00386     CPPUNIT_ASSERT( num_vtx > 3 );
00387     size_t vtxcopy[10];
00388     std::copy( vertices, vertices + num_vtx, vtxcopy );
00389     std::sort( vtxcopy, vtxcopy + num_vtx );
00390     for( unsigned i = 1; i < num_vtx; ++i )
00391     {
00392         CPPUNIT_ASSERT( vtxcopy[i] != vtxcopy[i - 1] );
00393         CPPUNIT_ASSERT( vtxcopy[i] < 10 );
00394     }
00395 }
00396 
00397 static void check_no_zeros( const MsqVector< 3 >* derivs, size_t num_vtx )
00398 {
00399     for( unsigned i = 0; i < num_vtx; ++i )
00400     {
00401         double dr = derivs[i][0];
00402         double ds = derivs[i][1];
00403         double dt = derivs[i][2];
00404         CPPUNIT_ASSERT( ( fabs( dr ) > 1e-6 ) || ( fabs( ds ) > 1e-6 ) || ( fabs( dt ) > 1e-6 ) );
00405     }
00406 }
00407 
00408 static void compare_coefficients( const double* coeffs, const size_t* indices, const double* expected_coeffs,
00409                                   size_t num_coeff, unsigned loc, NodeSet bits )
00410 {
00411     // find the location in the returned list for each node
00412     size_t revidx[10];
00413     double test_vals[10];
00414     for( size_t i = 0; i < 10; ++i )
00415     {
00416         revidx[i]    = std::find( indices, indices + num_coeff, i ) - indices;
00417         test_vals[i] = ( revidx[i] == num_coeff ) ? 0.0 : coeffs[revidx[i]];
00418     }
00419 
00420     // Check that index list doesn't contain any nodes not actually
00421     // present in the element.
00422     CPPUNIT_ASSERT( bits.mid_edge_node( 0 ) || ( revidx[4] == num_coeff ) );
00423     CPPUNIT_ASSERT( bits.mid_edge_node( 1 ) || ( revidx[5] == num_coeff ) );
00424     CPPUNIT_ASSERT( bits.mid_edge_node( 2 ) || ( revidx[6] == num_coeff ) );
00425     CPPUNIT_ASSERT( bits.mid_edge_node( 3 ) || ( revidx[7] == num_coeff ) );
00426     CPPUNIT_ASSERT( bits.mid_edge_node( 4 ) || ( revidx[8] == num_coeff ) );
00427     CPPUNIT_ASSERT( bits.mid_edge_node( 5 ) || ( revidx[9] == num_coeff ) );
00428 
00429     // compare expected and actual coefficient values
00430     ASSERT_VALUES_EQUAL( expected_coeffs[0], test_vals[0], loc, bits );
00431     ASSERT_VALUES_EQUAL( expected_coeffs[1], test_vals[1], loc, bits );
00432     ASSERT_VALUES_EQUAL( expected_coeffs[2], test_vals[2], loc, bits );
00433     ASSERT_VALUES_EQUAL( expected_coeffs[3], test_vals[3], loc, bits );
00434     ASSERT_VALUES_EQUAL( expected_coeffs[4], test_vals[4], loc, bits );
00435     ASSERT_VALUES_EQUAL( expected_coeffs[5], test_vals[5], loc, bits );
00436     ASSERT_VALUES_EQUAL( expected_coeffs[6], test_vals[6], loc, bits );
00437     ASSERT_VALUES_EQUAL( expected_coeffs[7], test_vals[7], loc, bits );
00438     ASSERT_VALUES_EQUAL( expected_coeffs[8], test_vals[8], loc, bits );
00439     ASSERT_VALUES_EQUAL( expected_coeffs[9], test_vals[9], loc, bits );
00440 }
00441 
00442 static void compare_derivatives( const size_t* vertices, size_t num_vtx, const MsqVector< 3 >* derivs,
00443                                  const MsqVector< 3 >* expected_derivs, unsigned loc, NodeSet bits )
00444 {
00445     check_valid_indices( vertices, num_vtx );
00446     check_no_zeros( derivs, num_vtx );
00447     MsqVector< 3 > expanded_derivs[30];
00448     memset( expanded_derivs, 0, sizeof( expanded_derivs ) );
00449     for( unsigned i = 0; i < num_vtx; ++i )
00450         expanded_derivs[vertices[i]] = derivs[i];
00451 
00452     ASSERT_VALUES_EQUAL( expected_derivs[0][0], expanded_derivs[0][0], loc, bits );
00453     ASSERT_VALUES_EQUAL( expected_derivs[0][1], expanded_derivs[0][1], loc, bits );
00454     ASSERT_VALUES_EQUAL( expected_derivs[0][2], expanded_derivs[0][2], loc, bits );
00455 
00456     ASSERT_VALUES_EQUAL( expected_derivs[1][0], expanded_derivs[1][0], loc, bits );
00457     ASSERT_VALUES_EQUAL( expected_derivs[1][1], expanded_derivs[1][1], loc, bits );
00458     ASSERT_VALUES_EQUAL( expected_derivs[1][2], expanded_derivs[1][2], loc, bits );
00459 
00460     ASSERT_VALUES_EQUAL( expected_derivs[2][0], expanded_derivs[2][0], loc, bits );
00461     ASSERT_VALUES_EQUAL( expected_derivs[2][1], expanded_derivs[2][1], loc, bits );
00462     ASSERT_VALUES_EQUAL( expected_derivs[2][2], expanded_derivs[2][2], loc, bits );
00463 
00464     ASSERT_VALUES_EQUAL( expected_derivs[3][0], expanded_derivs[3][0], loc, bits );
00465     ASSERT_VALUES_EQUAL( expected_derivs[3][1], expanded_derivs[3][1], loc, bits );
00466     ASSERT_VALUES_EQUAL( expected_derivs[3][2], expanded_derivs[3][2], loc, bits );
00467 
00468     ASSERT_VALUES_EQUAL( expected_derivs[4][0], expanded_derivs[4][0], loc, bits );
00469     ASSERT_VALUES_EQUAL( expected_derivs[4][1], expanded_derivs[4][1], loc, bits );
00470     ASSERT_VALUES_EQUAL( expected_derivs[4][2], expanded_derivs[4][2], loc, bits );
00471 
00472     ASSERT_VALUES_EQUAL( expected_derivs[5][0], expanded_derivs[5][0], loc, bits );
00473     ASSERT_VALUES_EQUAL( expected_derivs[5][1], expanded_derivs[5][1], loc, bits );
00474     ASSERT_VALUES_EQUAL( expected_derivs[5][2], expanded_derivs[5][2], loc, bits );
00475 
00476     ASSERT_VALUES_EQUAL( expected_derivs[6][0], expanded_derivs[6][0], loc, bits );
00477     ASSERT_VALUES_EQUAL( expected_derivs[6][1], expanded_derivs[6][1], loc, bits );
00478     ASSERT_VALUES_EQUAL( expected_derivs[6][2], expanded_derivs[6][2], loc, bits );
00479 
00480     ASSERT_VALUES_EQUAL( expected_derivs[7][0], expanded_derivs[7][0], loc, bits );
00481     ASSERT_VALUES_EQUAL( expected_derivs[7][1], expanded_derivs[7][1], loc, bits );
00482     ASSERT_VALUES_EQUAL( expected_derivs[7][2], expanded_derivs[7][2], loc, bits );
00483 
00484     ASSERT_VALUES_EQUAL( expected_derivs[8][0], expanded_derivs[8][0], loc, bits );
00485     ASSERT_VALUES_EQUAL( expected_derivs[8][1], expanded_derivs[8][1], loc, bits );
00486     ASSERT_VALUES_EQUAL( expected_derivs[8][2], expanded_derivs[8][2], loc, bits );
00487 
00488     ASSERT_VALUES_EQUAL( expected_derivs[9][0], expanded_derivs[9][0], loc, bits );
00489     ASSERT_VALUES_EQUAL( expected_derivs[9][1], expanded_derivs[9][1], loc, bits );
00490     ASSERT_VALUES_EQUAL( expected_derivs[9][2], expanded_derivs[9][2], loc, bits );
00491 }
00492 
00493 void TetLagrangeShapeTest::test_corner_coeff( int corner, NodeSet nodebits )
00494 {
00495     MsqPrintError err( std::cout );
00496 
00497     double expected[10];
00498     get_coeff( nodebits, rst_corner[corner], expected );
00499 
00500     double coeff[100];
00501     size_t n = 29, indices[100];
00502     sf.coefficients( Sample( 0, corner ), nodebits, coeff, indices, n, err );
00503     CPPUNIT_ASSERT( !err );
00504 
00505     compare_coefficients( coeff, indices, expected, n, corner, nodebits );
00506 }
00507 
00508 void TetLagrangeShapeTest::test_edge_coeff( int edge, NodeSet nodebits )
00509 {
00510     MsqPrintError err( std::cout );
00511 
00512     double expected[10];
00513     get_coeff( nodebits, rst_edge[edge], expected );
00514 
00515     double coeff[100];
00516     size_t n = 29, indices[100];
00517     sf.coefficients( Sample( 1, edge ), nodebits, coeff, indices, n, err );
00518     CPPUNIT_ASSERT( !err );
00519 
00520     compare_coefficients( coeff, indices, expected, n, edge + 4, nodebits );
00521 }
00522 
00523 void TetLagrangeShapeTest::test_face_coeff( int face, NodeSet nodebits )
00524 {
00525     MsqPrintError err( std::cout );
00526 
00527     double expected[10];
00528     get_coeff( nodebits, rst_face[face], expected );
00529 
00530     double coeff[100];
00531     size_t n = 29, indices[100];
00532     sf.coefficients( Sample( 2, face ), nodebits, coeff, indices, n, err );
00533     CPPUNIT_ASSERT( !err );
00534 
00535     compare_coefficients( coeff, indices, expected, n, face + 10, nodebits );
00536 }
00537 
00538 void TetLagrangeShapeTest::test_mid_coeff( NodeSet nodebits )
00539 {
00540     MsqPrintError err( std::cout );
00541 
00542     double expected[10];
00543     get_coeff( nodebits, rst_mid, expected );
00544 
00545     double coeff[100];
00546     size_t n = 29, indices[100];
00547     sf.coefficients( Sample( 3, 0 ), nodebits, coeff, indices, n, err );
00548     CPPUNIT_ASSERT( !err );
00549 
00550     compare_coefficients( coeff, indices, expected, n, 14, nodebits );
00551 }
00552 
00553 void TetLagrangeShapeTest::test_corner_derivs( int corner, NodeSet nodebits )
00554 {
00555     MsqPrintError err( std::cout );
00556 
00557     MsqVector< 3 > expected[10];
00558     get_derivs( nodebits, rst_corner[corner], expected );
00559 
00560     MsqVector< 3 > derivs[100];
00561     size_t vertices[100], n = 29;
00562     sf.derivatives( Sample( 0, corner ), nodebits, vertices, derivs, n, err );
00563     CPPUNIT_ASSERT( !err );
00564 
00565     compare_derivatives( vertices, n, derivs, expected, corner, nodebits );
00566 }
00567 
00568 void TetLagrangeShapeTest::test_edge_derivs( int edge, NodeSet nodebits )
00569 {
00570     MsqPrintError err( std::cout );
00571 
00572     MsqVector< 3 > expected[10];
00573     get_derivs( nodebits, rst_edge[edge], expected );
00574 
00575     MsqVector< 3 > derivs[100];
00576     size_t vertices[100], n = 29;
00577     sf.derivatives( Sample( 1, edge ), nodebits, vertices, derivs, n, err );
00578     CPPUNIT_ASSERT( !err );
00579 
00580     compare_derivatives( vertices, n, derivs, expected, edge + 4, nodebits );
00581 }
00582 
00583 void TetLagrangeShapeTest::test_face_derivs( int face, NodeSet nodebits )
00584 {
00585     MsqPrintError err( std::cout );
00586 
00587     MsqVector< 3 > expected[10];
00588     get_derivs( nodebits, rst_face[face], expected );
00589 
00590     MsqVector< 3 > derivs[100];
00591     size_t vertices[100], n = 29;
00592     sf.derivatives( Sample( 2, face ), nodebits, vertices, derivs, n, err );
00593     CPPUNIT_ASSERT( !err );
00594 
00595     compare_derivatives( vertices, n, derivs, expected, face + 10, nodebits );
00596 }
00597 
00598 void TetLagrangeShapeTest::test_mid_derivs( NodeSet nodebits )
00599 {
00600     MsqPrintError err( std::cout );
00601 
00602     MsqVector< 3 > expected[10];
00603     get_derivs( nodebits, rst_mid, expected );
00604 
00605     MsqVector< 3 > derivs[100];
00606     size_t vertices[100], n = 29;
00607     sf.derivatives( Sample( 3, 0 ), nodebits, vertices, derivs, n, err );
00608     CPPUNIT_ASSERT( !err );
00609 
00610     compare_derivatives( vertices, n, derivs, expected, 14, nodebits );
00611 }
00612 
00613 static NodeSet nodeset_from_bits( unsigned bits )
00614 {
00615     NodeSet result;
00616     for( unsigned i = 0; i < 6; ++i )
00617         if( bits & ( 1 << i ) ) result.set_mid_edge_node( i );
00618     for( unsigned i = 6; i < 10; ++i )
00619         if( bits & ( 1 << i ) ) result.set_mid_face_node( i - 6 );
00620     if( bits & ( 1 << 10 ) ) result.set_mid_region_node();
00621     return result;
00622 }
00623 
00624 void TetLagrangeShapeTest::test_coeff_corners()
00625 {
00626     for( unsigned i = 0; i < 0x40; ++i )
00627     {
00628         test_corner_coeff( 0, nodeset_from_bits( i ) );
00629         test_corner_coeff( 1, nodeset_from_bits( i ) );
00630         test_corner_coeff( 2, nodeset_from_bits( i ) );
00631         test_corner_coeff( 3, nodeset_from_bits( i ) );
00632     }
00633 }
00634 
00635 void TetLagrangeShapeTest::test_coeff_edges()
00636 {
00637     for( unsigned i = 0; i < 0x40; ++i )
00638     {
00639         test_edge_coeff( 0, nodeset_from_bits( i ) );
00640         test_edge_coeff( 1, nodeset_from_bits( i ) );
00641         test_edge_coeff( 2, nodeset_from_bits( i ) );
00642         test_edge_coeff( 3, nodeset_from_bits( i ) );
00643         test_edge_coeff( 4, nodeset_from_bits( i ) );
00644         test_edge_coeff( 5, nodeset_from_bits( i ) );
00645     }
00646 }
00647 
00648 void TetLagrangeShapeTest::test_coeff_faces()
00649 {
00650     for( unsigned i = 0; i < 0x40; ++i )
00651     {
00652         test_face_coeff( 0, nodeset_from_bits( i ) );
00653         test_face_coeff( 1, nodeset_from_bits( i ) );
00654         test_face_coeff( 2, nodeset_from_bits( i ) );
00655         test_face_coeff( 3, nodeset_from_bits( i ) );
00656     }
00657 }
00658 
00659 void TetLagrangeShapeTest::test_coeff_center()
00660 {
00661     for( unsigned i = 0; i < 0x40; ++i )
00662     {
00663         test_mid_coeff( nodeset_from_bits( i ) );
00664     }
00665 }
00666 
00667 void TetLagrangeShapeTest::test_deriv_corners()
00668 {
00669     for( unsigned i = 0; i < 0x40; ++i )
00670     {
00671         test_corner_derivs( 0, nodeset_from_bits( i ) );
00672         test_corner_derivs( 1, nodeset_from_bits( i ) );
00673         test_corner_derivs( 2, nodeset_from_bits( i ) );
00674         test_corner_derivs( 3, nodeset_from_bits( i ) );
00675     }
00676 }
00677 
00678 void TetLagrangeShapeTest::test_deriv_edges()
00679 {
00680     for( unsigned i = 0; i < 0x40; ++i )
00681     {
00682         test_edge_derivs( 0, nodeset_from_bits( i ) );
00683         test_edge_derivs( 1, nodeset_from_bits( i ) );
00684         test_edge_derivs( 2, nodeset_from_bits( i ) );
00685         test_edge_derivs( 3, nodeset_from_bits( i ) );
00686         test_edge_derivs( 4, nodeset_from_bits( i ) );
00687         test_edge_derivs( 5, nodeset_from_bits( i ) );
00688     }
00689 }
00690 
00691 void TetLagrangeShapeTest::test_deriv_faces()
00692 {
00693     for( unsigned i = 0; i < 0x40; ++i )
00694     {
00695         test_face_derivs( 0, nodeset_from_bits( i ) );
00696         test_face_derivs( 1, nodeset_from_bits( i ) );
00697         test_face_derivs( 2, nodeset_from_bits( i ) );
00698         test_face_derivs( 3, nodeset_from_bits( i ) );
00699     }
00700 }
00701 
00702 void TetLagrangeShapeTest::test_deriv_center()
00703 {
00704     for( unsigned i = 0; i < 0x40; ++i )
00705     {
00706         test_mid_derivs( nodeset_from_bits( i ) );
00707     }
00708 }
00709 
00710 void TetLagrangeShapeTest::test_invalid_nodebits_coeff( NodeSet bits )
00711 {
00712     MsqError err;
00713     double coeff[100];
00714     size_t n, indices[100];
00715 
00716     sf.coefficients( Sample( 0, 0 ), bits, coeff, indices, n, err );
00717     CPPUNIT_ASSERT( err );
00718     sf.coefficients( Sample( 0, 1 ), bits, coeff, indices, n, err );
00719     CPPUNIT_ASSERT( err );
00720     sf.coefficients( Sample( 0, 2 ), bits, coeff, indices, n, err );
00721     CPPUNIT_ASSERT( err );
00722     sf.coefficients( Sample( 0, 3 ), bits, coeff, indices, n, err );
00723     CPPUNIT_ASSERT( err );
00724 
00725     sf.coefficients( Sample( 1, 0 ), bits, coeff, indices, n, err );
00726     CPPUNIT_ASSERT( err );
00727     sf.coefficients( Sample( 1, 1 ), bits, coeff, indices, n, err );
00728     CPPUNIT_ASSERT( err );
00729     sf.coefficients( Sample( 1, 2 ), bits, coeff, indices, n, err );
00730     CPPUNIT_ASSERT( err );
00731     sf.coefficients( Sample( 1, 3 ), bits, coeff, indices, n, err );
00732     CPPUNIT_ASSERT( err );
00733     sf.coefficients( Sample( 1, 4 ), bits, coeff, indices, n, err );
00734     CPPUNIT_ASSERT( err );
00735     sf.coefficients( Sample( 1, 5 ), bits, coeff, indices, n, err );
00736     CPPUNIT_ASSERT( err );
00737 
00738     sf.coefficients( Sample( 2, 0 ), bits, coeff, indices, n, err );
00739     CPPUNIT_ASSERT( err );
00740     sf.coefficients( Sample( 2, 1 ), bits, coeff, indices, n, err );
00741     CPPUNIT_ASSERT( err );
00742     sf.coefficients( Sample( 2, 2 ), bits, coeff, indices, n, err );
00743     CPPUNIT_ASSERT( err );
00744     sf.coefficients( Sample( 2, 3 ), bits, coeff, indices, n, err );
00745     CPPUNIT_ASSERT( err );
00746 
00747     sf.coefficients( Sample( 3, 0 ), bits, coeff, indices, n, err );
00748     CPPUNIT_ASSERT( err );
00749 }
00750 
00751 void TetLagrangeShapeTest::test_invalid_nodebits_deriv( NodeSet bits )
00752 {
00753     MsqError err;
00754     size_t verts[100], n;
00755     MsqVector< 3 > derivs[100];
00756 
00757     sf.derivatives( Sample( 0, 0 ), bits, verts, derivs, n, err );
00758     CPPUNIT_ASSERT( err );
00759     sf.derivatives( Sample( 0, 1 ), bits, verts, derivs, n, err );
00760     CPPUNIT_ASSERT( err );
00761     sf.derivatives( Sample( 0, 2 ), bits, verts, derivs, n, err );
00762     CPPUNIT_ASSERT( err );
00763     sf.derivatives( Sample( 0, 3 ), bits, verts, derivs, n, err );
00764     CPPUNIT_ASSERT( err );
00765 
00766     sf.derivatives( Sample( 1, 0 ), bits, verts, derivs, n, err );
00767     CPPUNIT_ASSERT( err );
00768     sf.derivatives( Sample( 1, 1 ), bits, verts, derivs, n, err );
00769     CPPUNIT_ASSERT( err );
00770     sf.derivatives( Sample( 1, 2 ), bits, verts, derivs, n, err );
00771     CPPUNIT_ASSERT( err );
00772     sf.derivatives( Sample( 1, 3 ), bits, verts, derivs, n, err );
00773     CPPUNIT_ASSERT( err );
00774     sf.derivatives( Sample( 1, 4 ), bits, verts, derivs, n, err );
00775     CPPUNIT_ASSERT( err );
00776     sf.derivatives( Sample( 1, 5 ), bits, verts, derivs, n, err );
00777     CPPUNIT_ASSERT( err );
00778 
00779     sf.derivatives( Sample( 2, 0 ), bits, verts, derivs, n, err );
00780     CPPUNIT_ASSERT( err );
00781     sf.derivatives( Sample( 2, 1 ), bits, verts, derivs, n, err );
00782     CPPUNIT_ASSERT( err );
00783     sf.derivatives( Sample( 2, 2 ), bits, verts, derivs, n, err );
00784     CPPUNIT_ASSERT( err );
00785     sf.derivatives( Sample( 2, 3 ), bits, verts, derivs, n, err );
00786     CPPUNIT_ASSERT( err );
00787 
00788     sf.derivatives( Sample( 3, 0 ), bits, verts, derivs, n, err );
00789     CPPUNIT_ASSERT( err );
00790 }
00791 
00792 void TetLagrangeShapeTest::test_mid_elem_node_coeff()
00793 {
00794     NodeSet nodeset;
00795     nodeset.set_mid_region_node();
00796     test_invalid_nodebits_coeff( nodeset );
00797 }
00798 
00799 void TetLagrangeShapeTest::test_mid_elem_node_deriv()
00800 {
00801     NodeSet nodeset;
00802     nodeset.set_mid_region_node();
00803     test_invalid_nodebits_deriv( nodeset );
00804 }
00805 
00806 void TetLagrangeShapeTest::test_mid_face_node_coeff()
00807 {
00808     NodeSet nodeset1, nodeset2;
00809     nodeset1.set_mid_face_node( 0 );
00810     test_invalid_nodebits_coeff( nodeset1 );
00811     nodeset2.set_mid_face_node( 2 );
00812     test_invalid_nodebits_coeff( nodeset2 );
00813 }
00814 
00815 void TetLagrangeShapeTest::test_mid_face_node_deriv()
00816 {
00817     NodeSet nodeset1, nodeset2;
00818     nodeset1.set_mid_face_node( 0 );
00819     test_invalid_nodebits_deriv( nodeset1 );
00820     nodeset2.set_mid_face_node( 2 );
00821     test_invalid_nodebits_deriv( nodeset2 );
00822 }
00823 
00824 void TetLagrangeShapeTest::test_ideal_jacobian()
00825 {
00826     MsqError err;
00827     MsqMatrix< 3, 3 > J_act, J_exp;
00828     sf.ideal( Sample( 3, 0 ), J_act, err );
00829     ASSERT_NO_ERROR( err );
00830     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( J_act ), 1e-6 );
00831 
00832     const Vector3D* verts = unit_edge_element( TETRAHEDRON );
00833     CPPUNIT_ASSERT( verts );
00834 
00835     JacobianCalculator jc;
00836     jc.get_Jacobian_3D( &sf, NodeSet(), Sample( 2, 0 ), verts, 4, J_exp, err );
00837     ASSERT_NO_ERROR( err );
00838     J_exp /= MBMesquite::cbrt( det( J_exp ) );
00839 
00840     // Matrices should be a rotation of each other.
00841     // First, calculate tentative rotation matrix
00842     MsqMatrix< 3, 3 > R = inverse( J_exp ) * J_act;
00843     // next check that it is a rotation
00844     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( R ), 1e-6 );          // no scaling
00845     ASSERT_MATRICES_EQUAL( transpose( R ), inverse( R ), 1e-6 );  // orthogonal
00846 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines