MOAB: Mesh Oriented datABase  (version 5.4.0)
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 <cmath>
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,
00409                                   const size_t* indices,
00410                                   const double* expected_coeffs,
00411                                   size_t num_coeff,
00412                                   unsigned loc,
00413                                   NodeSet bits )
00414 {
00415     // find the location in the returned list for each node
00416     size_t revidx[10];
00417     double test_vals[10];
00418     for( size_t i = 0; i < 10; ++i )
00419     {
00420         revidx[i]    = std::find( indices, indices + num_coeff, i ) - indices;
00421         test_vals[i] = ( revidx[i] == num_coeff ) ? 0.0 : coeffs[revidx[i]];
00422     }
00423 
00424     // Check that index list doesn't contain any nodes not actually
00425     // present in the element.
00426     CPPUNIT_ASSERT( bits.mid_edge_node( 0 ) || ( revidx[4] == num_coeff ) );
00427     CPPUNIT_ASSERT( bits.mid_edge_node( 1 ) || ( revidx[5] == num_coeff ) );
00428     CPPUNIT_ASSERT( bits.mid_edge_node( 2 ) || ( revidx[6] == num_coeff ) );
00429     CPPUNIT_ASSERT( bits.mid_edge_node( 3 ) || ( revidx[7] == num_coeff ) );
00430     CPPUNIT_ASSERT( bits.mid_edge_node( 4 ) || ( revidx[8] == num_coeff ) );
00431     CPPUNIT_ASSERT( bits.mid_edge_node( 5 ) || ( revidx[9] == 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     ASSERT_VALUES_EQUAL( expected_coeffs[9], test_vals[9], loc, bits );
00444 }
00445 
00446 static void compare_derivatives( const size_t* vertices,
00447                                  size_t num_vtx,
00448                                  const MsqVector< 3 >* derivs,
00449                                  const MsqVector< 3 >* expected_derivs,
00450                                  unsigned loc,
00451                                  NodeSet bits )
00452 {
00453     check_valid_indices( vertices, num_vtx );
00454     check_no_zeros( derivs, num_vtx );
00455     MsqVector< 3 > expanded_derivs[30];
00456     memset( expanded_derivs, 0, sizeof( expanded_derivs ) );
00457     for( unsigned i = 0; i < num_vtx; ++i )
00458         expanded_derivs[vertices[i]] = derivs[i];
00459 
00460     ASSERT_VALUES_EQUAL( expected_derivs[0][0], expanded_derivs[0][0], loc, bits );
00461     ASSERT_VALUES_EQUAL( expected_derivs[0][1], expanded_derivs[0][1], loc, bits );
00462     ASSERT_VALUES_EQUAL( expected_derivs[0][2], expanded_derivs[0][2], loc, bits );
00463 
00464     ASSERT_VALUES_EQUAL( expected_derivs[1][0], expanded_derivs[1][0], loc, bits );
00465     ASSERT_VALUES_EQUAL( expected_derivs[1][1], expanded_derivs[1][1], loc, bits );
00466     ASSERT_VALUES_EQUAL( expected_derivs[1][2], expanded_derivs[1][2], loc, bits );
00467 
00468     ASSERT_VALUES_EQUAL( expected_derivs[2][0], expanded_derivs[2][0], loc, bits );
00469     ASSERT_VALUES_EQUAL( expected_derivs[2][1], expanded_derivs[2][1], loc, bits );
00470     ASSERT_VALUES_EQUAL( expected_derivs[2][2], expanded_derivs[2][2], loc, bits );
00471 
00472     ASSERT_VALUES_EQUAL( expected_derivs[3][0], expanded_derivs[3][0], loc, bits );
00473     ASSERT_VALUES_EQUAL( expected_derivs[3][1], expanded_derivs[3][1], loc, bits );
00474     ASSERT_VALUES_EQUAL( expected_derivs[3][2], expanded_derivs[3][2], loc, bits );
00475 
00476     ASSERT_VALUES_EQUAL( expected_derivs[4][0], expanded_derivs[4][0], loc, bits );
00477     ASSERT_VALUES_EQUAL( expected_derivs[4][1], expanded_derivs[4][1], loc, bits );
00478     ASSERT_VALUES_EQUAL( expected_derivs[4][2], expanded_derivs[4][2], loc, bits );
00479 
00480     ASSERT_VALUES_EQUAL( expected_derivs[5][0], expanded_derivs[5][0], loc, bits );
00481     ASSERT_VALUES_EQUAL( expected_derivs[5][1], expanded_derivs[5][1], loc, bits );
00482     ASSERT_VALUES_EQUAL( expected_derivs[5][2], expanded_derivs[5][2], loc, bits );
00483 
00484     ASSERT_VALUES_EQUAL( expected_derivs[6][0], expanded_derivs[6][0], loc, bits );
00485     ASSERT_VALUES_EQUAL( expected_derivs[6][1], expanded_derivs[6][1], loc, bits );
00486     ASSERT_VALUES_EQUAL( expected_derivs[6][2], expanded_derivs[6][2], loc, bits );
00487 
00488     ASSERT_VALUES_EQUAL( expected_derivs[7][0], expanded_derivs[7][0], loc, bits );
00489     ASSERT_VALUES_EQUAL( expected_derivs[7][1], expanded_derivs[7][1], loc, bits );
00490     ASSERT_VALUES_EQUAL( expected_derivs[7][2], expanded_derivs[7][2], loc, bits );
00491 
00492     ASSERT_VALUES_EQUAL( expected_derivs[8][0], expanded_derivs[8][0], loc, bits );
00493     ASSERT_VALUES_EQUAL( expected_derivs[8][1], expanded_derivs[8][1], loc, bits );
00494     ASSERT_VALUES_EQUAL( expected_derivs[8][2], expanded_derivs[8][2], loc, bits );
00495 
00496     ASSERT_VALUES_EQUAL( expected_derivs[9][0], expanded_derivs[9][0], loc, bits );
00497     ASSERT_VALUES_EQUAL( expected_derivs[9][1], expanded_derivs[9][1], loc, bits );
00498     ASSERT_VALUES_EQUAL( expected_derivs[9][2], expanded_derivs[9][2], loc, bits );
00499 }
00500 
00501 void TetLagrangeShapeTest::test_corner_coeff( int corner, NodeSet nodebits )
00502 {
00503     MsqPrintError err( std::cout );
00504 
00505     double expected[10];
00506     get_coeff( nodebits, rst_corner[corner], expected );
00507 
00508     double coeff[100];
00509     size_t n = 29, indices[100];
00510     sf.coefficients( Sample( 0, corner ), nodebits, coeff, indices, n, err );
00511     CPPUNIT_ASSERT( !err );
00512 
00513     compare_coefficients( coeff, indices, expected, n, corner, nodebits );
00514 }
00515 
00516 void TetLagrangeShapeTest::test_edge_coeff( int edge, NodeSet nodebits )
00517 {
00518     MsqPrintError err( std::cout );
00519 
00520     double expected[10];
00521     get_coeff( nodebits, rst_edge[edge], expected );
00522 
00523     double coeff[100];
00524     size_t n = 29, indices[100];
00525     sf.coefficients( Sample( 1, edge ), nodebits, coeff, indices, n, err );
00526     CPPUNIT_ASSERT( !err );
00527 
00528     compare_coefficients( coeff, indices, expected, n, edge + 4, nodebits );
00529 }
00530 
00531 void TetLagrangeShapeTest::test_face_coeff( int face, NodeSet nodebits )
00532 {
00533     MsqPrintError err( std::cout );
00534 
00535     double expected[10];
00536     get_coeff( nodebits, rst_face[face], expected );
00537 
00538     double coeff[100];
00539     size_t n = 29, indices[100];
00540     sf.coefficients( Sample( 2, face ), nodebits, coeff, indices, n, err );
00541     CPPUNIT_ASSERT( !err );
00542 
00543     compare_coefficients( coeff, indices, expected, n, face + 10, nodebits );
00544 }
00545 
00546 void TetLagrangeShapeTest::test_mid_coeff( NodeSet nodebits )
00547 {
00548     MsqPrintError err( std::cout );
00549 
00550     double expected[10];
00551     get_coeff( nodebits, rst_mid, expected );
00552 
00553     double coeff[100];
00554     size_t n = 29, indices[100];
00555     sf.coefficients( Sample( 3, 0 ), nodebits, coeff, indices, n, err );
00556     CPPUNIT_ASSERT( !err );
00557 
00558     compare_coefficients( coeff, indices, expected, n, 14, nodebits );
00559 }
00560 
00561 void TetLagrangeShapeTest::test_corner_derivs( int corner, NodeSet nodebits )
00562 {
00563     MsqPrintError err( std::cout );
00564 
00565     MsqVector< 3 > expected[10];
00566     get_derivs( nodebits, rst_corner[corner], expected );
00567 
00568     MsqVector< 3 > derivs[100];
00569     size_t vertices[100], n = 29;
00570     sf.derivatives( Sample( 0, corner ), nodebits, vertices, derivs, n, err );
00571     CPPUNIT_ASSERT( !err );
00572 
00573     compare_derivatives( vertices, n, derivs, expected, corner, nodebits );
00574 }
00575 
00576 void TetLagrangeShapeTest::test_edge_derivs( int edge, NodeSet nodebits )
00577 {
00578     MsqPrintError err( std::cout );
00579 
00580     MsqVector< 3 > expected[10];
00581     get_derivs( nodebits, rst_edge[edge], expected );
00582 
00583     MsqVector< 3 > derivs[100];
00584     size_t vertices[100], n = 29;
00585     sf.derivatives( Sample( 1, edge ), nodebits, vertices, derivs, n, err );
00586     CPPUNIT_ASSERT( !err );
00587 
00588     compare_derivatives( vertices, n, derivs, expected, edge + 4, nodebits );
00589 }
00590 
00591 void TetLagrangeShapeTest::test_face_derivs( int face, NodeSet nodebits )
00592 {
00593     MsqPrintError err( std::cout );
00594 
00595     MsqVector< 3 > expected[10];
00596     get_derivs( nodebits, rst_face[face], expected );
00597 
00598     MsqVector< 3 > derivs[100];
00599     size_t vertices[100], n = 29;
00600     sf.derivatives( Sample( 2, face ), nodebits, vertices, derivs, n, err );
00601     CPPUNIT_ASSERT( !err );
00602 
00603     compare_derivatives( vertices, n, derivs, expected, face + 10, nodebits );
00604 }
00605 
00606 void TetLagrangeShapeTest::test_mid_derivs( NodeSet nodebits )
00607 {
00608     MsqPrintError err( std::cout );
00609 
00610     MsqVector< 3 > expected[10];
00611     get_derivs( nodebits, rst_mid, expected );
00612 
00613     MsqVector< 3 > derivs[100];
00614     size_t vertices[100], n = 29;
00615     sf.derivatives( Sample( 3, 0 ), nodebits, vertices, derivs, n, err );
00616     CPPUNIT_ASSERT( !err );
00617 
00618     compare_derivatives( vertices, n, derivs, expected, 14, nodebits );
00619 }
00620 
00621 static NodeSet nodeset_from_bits( unsigned bits )
00622 {
00623     NodeSet result;
00624     for( unsigned i = 0; i < 6; ++i )
00625         if( bits & ( 1 << i ) ) result.set_mid_edge_node( i );
00626     for( unsigned i = 6; i < 10; ++i )
00627         if( bits & ( 1 << i ) ) result.set_mid_face_node( i - 6 );
00628     if( bits & ( 1 << 10 ) ) result.set_mid_region_node();
00629     return result;
00630 }
00631 
00632 void TetLagrangeShapeTest::test_coeff_corners()
00633 {
00634     for( unsigned i = 0; i < 0x40; ++i )
00635     {
00636         test_corner_coeff( 0, nodeset_from_bits( i ) );
00637         test_corner_coeff( 1, nodeset_from_bits( i ) );
00638         test_corner_coeff( 2, nodeset_from_bits( i ) );
00639         test_corner_coeff( 3, nodeset_from_bits( i ) );
00640     }
00641 }
00642 
00643 void TetLagrangeShapeTest::test_coeff_edges()
00644 {
00645     for( unsigned i = 0; i < 0x40; ++i )
00646     {
00647         test_edge_coeff( 0, nodeset_from_bits( i ) );
00648         test_edge_coeff( 1, nodeset_from_bits( i ) );
00649         test_edge_coeff( 2, nodeset_from_bits( i ) );
00650         test_edge_coeff( 3, nodeset_from_bits( i ) );
00651         test_edge_coeff( 4, nodeset_from_bits( i ) );
00652         test_edge_coeff( 5, nodeset_from_bits( i ) );
00653     }
00654 }
00655 
00656 void TetLagrangeShapeTest::test_coeff_faces()
00657 {
00658     for( unsigned i = 0; i < 0x40; ++i )
00659     {
00660         test_face_coeff( 0, nodeset_from_bits( i ) );
00661         test_face_coeff( 1, nodeset_from_bits( i ) );
00662         test_face_coeff( 2, nodeset_from_bits( i ) );
00663         test_face_coeff( 3, nodeset_from_bits( i ) );
00664     }
00665 }
00666 
00667 void TetLagrangeShapeTest::test_coeff_center()
00668 {
00669     for( unsigned i = 0; i < 0x40; ++i )
00670     {
00671         test_mid_coeff( nodeset_from_bits( i ) );
00672     }
00673 }
00674 
00675 void TetLagrangeShapeTest::test_deriv_corners()
00676 {
00677     for( unsigned i = 0; i < 0x40; ++i )
00678     {
00679         test_corner_derivs( 0, nodeset_from_bits( i ) );
00680         test_corner_derivs( 1, nodeset_from_bits( i ) );
00681         test_corner_derivs( 2, nodeset_from_bits( i ) );
00682         test_corner_derivs( 3, nodeset_from_bits( i ) );
00683     }
00684 }
00685 
00686 void TetLagrangeShapeTest::test_deriv_edges()
00687 {
00688     for( unsigned i = 0; i < 0x40; ++i )
00689     {
00690         test_edge_derivs( 0, nodeset_from_bits( i ) );
00691         test_edge_derivs( 1, nodeset_from_bits( i ) );
00692         test_edge_derivs( 2, nodeset_from_bits( i ) );
00693         test_edge_derivs( 3, nodeset_from_bits( i ) );
00694         test_edge_derivs( 4, nodeset_from_bits( i ) );
00695         test_edge_derivs( 5, nodeset_from_bits( i ) );
00696     }
00697 }
00698 
00699 void TetLagrangeShapeTest::test_deriv_faces()
00700 {
00701     for( unsigned i = 0; i < 0x40; ++i )
00702     {
00703         test_face_derivs( 0, nodeset_from_bits( i ) );
00704         test_face_derivs( 1, nodeset_from_bits( i ) );
00705         test_face_derivs( 2, nodeset_from_bits( i ) );
00706         test_face_derivs( 3, nodeset_from_bits( i ) );
00707     }
00708 }
00709 
00710 void TetLagrangeShapeTest::test_deriv_center()
00711 {
00712     for( unsigned i = 0; i < 0x40; ++i )
00713     {
00714         test_mid_derivs( nodeset_from_bits( i ) );
00715     }
00716 }
00717 
00718 void TetLagrangeShapeTest::test_invalid_nodebits_coeff( NodeSet bits )
00719 {
00720     MsqError err;
00721     double coeff[100];
00722     size_t n, indices[100];
00723 
00724     sf.coefficients( Sample( 0, 0 ), bits, coeff, indices, n, err );
00725     CPPUNIT_ASSERT( err );
00726     sf.coefficients( Sample( 0, 1 ), bits, coeff, indices, n, err );
00727     CPPUNIT_ASSERT( err );
00728     sf.coefficients( Sample( 0, 2 ), bits, coeff, indices, n, err );
00729     CPPUNIT_ASSERT( err );
00730     sf.coefficients( Sample( 0, 3 ), bits, coeff, indices, n, err );
00731     CPPUNIT_ASSERT( err );
00732 
00733     sf.coefficients( Sample( 1, 0 ), bits, coeff, indices, n, err );
00734     CPPUNIT_ASSERT( err );
00735     sf.coefficients( Sample( 1, 1 ), bits, coeff, indices, n, err );
00736     CPPUNIT_ASSERT( err );
00737     sf.coefficients( Sample( 1, 2 ), bits, coeff, indices, n, err );
00738     CPPUNIT_ASSERT( err );
00739     sf.coefficients( Sample( 1, 3 ), bits, coeff, indices, n, err );
00740     CPPUNIT_ASSERT( err );
00741     sf.coefficients( Sample( 1, 4 ), bits, coeff, indices, n, err );
00742     CPPUNIT_ASSERT( err );
00743     sf.coefficients( Sample( 1, 5 ), bits, coeff, indices, n, err );
00744     CPPUNIT_ASSERT( err );
00745 
00746     sf.coefficients( Sample( 2, 0 ), bits, coeff, indices, n, err );
00747     CPPUNIT_ASSERT( err );
00748     sf.coefficients( Sample( 2, 1 ), bits, coeff, indices, n, err );
00749     CPPUNIT_ASSERT( err );
00750     sf.coefficients( Sample( 2, 2 ), bits, coeff, indices, n, err );
00751     CPPUNIT_ASSERT( err );
00752     sf.coefficients( Sample( 2, 3 ), bits, coeff, indices, n, err );
00753     CPPUNIT_ASSERT( err );
00754 
00755     sf.coefficients( Sample( 3, 0 ), bits, coeff, indices, n, err );
00756     CPPUNIT_ASSERT( err );
00757 }
00758 
00759 void TetLagrangeShapeTest::test_invalid_nodebits_deriv( NodeSet bits )
00760 {
00761     MsqError err;
00762     size_t verts[100], n;
00763     MsqVector< 3 > derivs[100];
00764 
00765     sf.derivatives( Sample( 0, 0 ), bits, verts, derivs, n, err );
00766     CPPUNIT_ASSERT( err );
00767     sf.derivatives( Sample( 0, 1 ), bits, verts, derivs, n, err );
00768     CPPUNIT_ASSERT( err );
00769     sf.derivatives( Sample( 0, 2 ), bits, verts, derivs, n, err );
00770     CPPUNIT_ASSERT( err );
00771     sf.derivatives( Sample( 0, 3 ), bits, verts, derivs, n, err );
00772     CPPUNIT_ASSERT( err );
00773 
00774     sf.derivatives( Sample( 1, 0 ), bits, verts, derivs, n, err );
00775     CPPUNIT_ASSERT( err );
00776     sf.derivatives( Sample( 1, 1 ), bits, verts, derivs, n, err );
00777     CPPUNIT_ASSERT( err );
00778     sf.derivatives( Sample( 1, 2 ), bits, verts, derivs, n, err );
00779     CPPUNIT_ASSERT( err );
00780     sf.derivatives( Sample( 1, 3 ), bits, verts, derivs, n, err );
00781     CPPUNIT_ASSERT( err );
00782     sf.derivatives( Sample( 1, 4 ), bits, verts, derivs, n, err );
00783     CPPUNIT_ASSERT( err );
00784     sf.derivatives( Sample( 1, 5 ), bits, verts, derivs, n, err );
00785     CPPUNIT_ASSERT( err );
00786 
00787     sf.derivatives( Sample( 2, 0 ), bits, verts, derivs, n, err );
00788     CPPUNIT_ASSERT( err );
00789     sf.derivatives( Sample( 2, 1 ), bits, verts, derivs, n, err );
00790     CPPUNIT_ASSERT( err );
00791     sf.derivatives( Sample( 2, 2 ), bits, verts, derivs, n, err );
00792     CPPUNIT_ASSERT( err );
00793     sf.derivatives( Sample( 2, 3 ), bits, verts, derivs, n, err );
00794     CPPUNIT_ASSERT( err );
00795 
00796     sf.derivatives( Sample( 3, 0 ), bits, verts, derivs, n, err );
00797     CPPUNIT_ASSERT( err );
00798 }
00799 
00800 void TetLagrangeShapeTest::test_mid_elem_node_coeff()
00801 {
00802     NodeSet nodeset;
00803     nodeset.set_mid_region_node();
00804     test_invalid_nodebits_coeff( nodeset );
00805 }
00806 
00807 void TetLagrangeShapeTest::test_mid_elem_node_deriv()
00808 {
00809     NodeSet nodeset;
00810     nodeset.set_mid_region_node();
00811     test_invalid_nodebits_deriv( nodeset );
00812 }
00813 
00814 void TetLagrangeShapeTest::test_mid_face_node_coeff()
00815 {
00816     NodeSet nodeset1, nodeset2;
00817     nodeset1.set_mid_face_node( 0 );
00818     test_invalid_nodebits_coeff( nodeset1 );
00819     nodeset2.set_mid_face_node( 2 );
00820     test_invalid_nodebits_coeff( nodeset2 );
00821 }
00822 
00823 void TetLagrangeShapeTest::test_mid_face_node_deriv()
00824 {
00825     NodeSet nodeset1, nodeset2;
00826     nodeset1.set_mid_face_node( 0 );
00827     test_invalid_nodebits_deriv( nodeset1 );
00828     nodeset2.set_mid_face_node( 2 );
00829     test_invalid_nodebits_deriv( nodeset2 );
00830 }
00831 
00832 void TetLagrangeShapeTest::test_ideal_jacobian()
00833 {
00834     MsqError err;
00835     MsqMatrix< 3, 3 > J_act, J_exp;
00836     sf.ideal( Sample( 3, 0 ), J_act, err );
00837     ASSERT_NO_ERROR( err );
00838     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( J_act ), 1e-6 );
00839 
00840     const Vector3D* verts = unit_edge_element( TETRAHEDRON );
00841     CPPUNIT_ASSERT( verts );
00842 
00843     JacobianCalculator jc;
00844     jc.get_Jacobian_3D( &sf, NodeSet(), Sample( 2, 0 ), verts, 4, J_exp, err );
00845     ASSERT_NO_ERROR( err );
00846     J_exp /= MBMesquite::cbrt( det( J_exp ) );
00847 
00848     // Matrices should be a rotation of each other.
00849     // First, calculate tentative rotation matrix
00850     MsqMatrix< 3, 3 > R = inverse( J_exp ) * J_act;
00851     // next check that it is a rotation
00852     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( R ), 1e-6 );          // no scaling
00853     ASSERT_MATRICES_EQUAL( transpose( R ), inverse( R ), 1e-6 );  // orthogonal
00854 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines