MOAB: Mesh Oriented datABase  (version 5.3.1)
TriLagrangeShapeTest.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 TriLagrangeShapeTest.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "TriLagrangeShape.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 < 3 )
00067         buffer3 << "Corner " << location;
00068     else if( location < 6 )
00069         buffer3 << "Edge " << location - 3;
00070     else if( location == 6 )
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 TriLagrangeShapeTest : public CppUnit::TestFixture
00083 {
00084   private:
00085     CPPUNIT_TEST_SUITE( TriLagrangeShapeTest );
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     TriLagrangeShape 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( TriLagrangeShapeTest, "TriLagrangeShapeTest" );
00122 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( TriLagrangeShapeTest, "Unit" );
00123 
00124 double N0( double r, double s )
00125 {
00126     double t = 1 - r - s;
00127     return t * ( 2 * t - 1 );
00128 }
00129 double N1( double r, double )
00130 {
00131     return r * ( 2 * r - 1 );
00132 }
00133 double N2( double, double s )
00134 {
00135     return s * ( 2 * s - 1 );
00136 }
00137 double N3( double r, double s )
00138 {
00139     double t = 1 - r - s;
00140     return 4 * r * t;
00141 }
00142 double N4( double r, double s )
00143 {
00144     return 4 * r * s;
00145 }
00146 double N5( double r, double s )
00147 {
00148     double t = 1 - r - s;
00149     return 4 * s * t;
00150 }
00151 double dN0dr( double r, double s )
00152 {
00153     return 4 * r + 4 * s - 3;
00154 }
00155 double dN0ds( double r, double s )
00156 {
00157     return 4 * r + 4 * s - 3;
00158 }
00159 double dN1dr( double r, double )
00160 {
00161     return 4 * r - 1;
00162 }
00163 double dN1ds( double, double )
00164 {
00165     return 0;
00166 }
00167 double dN2dr( double, double )
00168 {
00169     return 0;
00170 }
00171 double dN2ds( double, double s )
00172 {
00173     return 4 * s - 1;
00174 }
00175 double dN3dr( double r, double s )
00176 {
00177     return 4 * ( 1 - s - 2 * r );
00178 }
00179 double dN3ds( double r, double )
00180 {
00181     return -4 * r;
00182 }
00183 double dN4dr( double, double s )
00184 {
00185     return 4 * s;
00186 }
00187 double dN4ds( double r, double )
00188 {
00189     return 4 * r;
00190 }
00191 double dN5dr( double, double s )
00192 {
00193     return -4 * s;
00194 }
00195 double dN5ds( double r, double s )
00196 {
00197     return 4 * ( 1 - r - 2 * s );
00198 }
00199 
00200 typedef double ( *N_t )( double, double );
00201 const N_t N[]    = { &N0, &N1, &N2, &N3, &N4, &N5 };
00202 const N_t dNdr[] = { &dN0dr, &dN1dr, &dN2dr, &dN3dr, &dN4dr, &dN5dr };
00203 const N_t dNds[] = { &dN0ds, &dN1ds, &dN2ds, &dN3ds, &dN4ds, &dN5ds };
00204 
00205 const double rs_corner[][2] = { { 0, 0 }, { 1, 0 }, { 0, 1 } };
00206 const double rs_edge[][2]   = { { 0.5, 0.0 }, { 0.5, 0.5 }, { 0.0, 0.5 } };
00207 const double rs_mid[2]      = { 1.0 / 3.0, 1.0 / 3.0 };
00208 
00209 static void get_coeff( NodeSet nodeset, const double* rs, double* coeffs )
00210 {
00211     for( int i = 0; i < 6; ++i )
00212         coeffs[i] = ( *N[i] )( rs[0], rs[1] );
00213     for( int i = 0; i < 3; ++i )
00214         if( !nodeset.mid_edge_node( i ) )
00215         {
00216             coeffs[i] += 0.5 * coeffs[i + 3];
00217             coeffs[( i + 1 ) % 3] += 0.5 * coeffs[i + 3];
00218             coeffs[i + 3] = 0;
00219         }
00220 }
00221 
00222 static void get_derivs( NodeSet nodeset, const double* rs, double* derivs )
00223 {
00224     for( int i = 0; i < 6; ++i )
00225     {
00226         derivs[2 * i]     = ( *dNdr[i] )( rs[0], rs[1] );
00227         derivs[2 * i + 1] = ( *dNds[i] )( rs[0], rs[1] );
00228     }
00229     for( int i = 0; i < 3; ++i )
00230         if( !nodeset.mid_edge_node( i ) )
00231         {
00232             derivs[2 * i] += 0.5 * derivs[2 * i + 6];
00233             derivs[2 * i + 1] += 0.5 * derivs[2 * i + 7];
00234             int j = ( i + 1 ) % 3;
00235             derivs[2 * j] += 0.5 * derivs[2 * i + 6];
00236             derivs[2 * j + 1] += 0.5 * derivs[2 * i + 7];
00237             derivs[2 * i + 6] = 0;
00238             derivs[2 * i + 7] = 0;
00239         }
00240 }
00241 
00242 static void check_valid_indices( const size_t* vtx_in, size_t num_vtx )
00243 {
00244     CPPUNIT_ASSERT( num_vtx < 7 );
00245     CPPUNIT_ASSERT( num_vtx > 2 );
00246     size_t vertices[6];
00247     std::copy( vtx_in, vtx_in + num_vtx, vertices );
00248     std::sort( vertices, vertices + num_vtx );
00249     for( unsigned i = 1; i < num_vtx; ++i )
00250     {
00251         CPPUNIT_ASSERT( vertices[i] != vertices[i - 1] );
00252         CPPUNIT_ASSERT( vertices[i] < 6 );
00253     }
00254 }
00255 
00256 static void check_no_zeros( const MsqVector< 2 >* derivs, size_t num_vtx )
00257 {
00258     for( unsigned i = 0; i < num_vtx; ++i )
00259     {
00260         double dr = derivs[i][0];
00261         double ds = derivs[i][1];
00262         CPPUNIT_ASSERT( ( fabs( dr ) > 1e-6 ) || ( fabs( ds ) > 1e-6 ) );
00263     }
00264 }
00265 
00266 static void compare_coefficients( const double* coeffs, const size_t* indices, const double* expected_coeffs,
00267                                   size_t num_coeff, unsigned loc, NodeSet bits )
00268 {
00269     // find the location in the returned list for each node
00270     size_t revidx[6];
00271     double test_vals[6];
00272     for( size_t i = 0; i < 6; ++i )
00273     {
00274         revidx[i]    = std::find( indices, indices + num_coeff, i ) - indices;
00275         test_vals[i] = ( revidx[i] == num_coeff ) ? 0.0 : coeffs[revidx[i]];
00276     }
00277 
00278     // Check that index list doesn't contain any nodes not actually
00279     // present in the element.
00280     CPPUNIT_ASSERT( bits.mid_edge_node( 0 ) || ( revidx[3] == num_coeff ) );
00281     CPPUNIT_ASSERT( bits.mid_edge_node( 1 ) || ( revidx[4] == num_coeff ) );
00282     CPPUNIT_ASSERT( bits.mid_edge_node( 2 ) || ( revidx[5] == num_coeff ) );
00283 
00284     // compare expected and actual coefficient values
00285     ASSERT_VALUES_EQUAL( expected_coeffs[0], test_vals[0], loc, bits );
00286     ASSERT_VALUES_EQUAL( expected_coeffs[1], test_vals[1], loc, bits );
00287     ASSERT_VALUES_EQUAL( expected_coeffs[2], test_vals[2], loc, bits );
00288     ASSERT_VALUES_EQUAL( expected_coeffs[3], test_vals[3], loc, bits );
00289     ASSERT_VALUES_EQUAL( expected_coeffs[4], test_vals[4], loc, bits );
00290     ASSERT_VALUES_EQUAL( expected_coeffs[5], test_vals[5], loc, bits );
00291 }
00292 
00293 static void compare_derivatives( const size_t* vertices, size_t num_vtx, const MsqVector< 2 >* derivs,
00294                                  const double* expected_derivs, unsigned loc, NodeSet bits )
00295 {
00296     check_valid_indices( vertices, num_vtx );
00297     check_no_zeros( derivs, num_vtx );
00298     double expanded_derivs[12] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
00299     for( unsigned i = 0; i < num_vtx; ++i )
00300     {
00301         expanded_derivs[2 * vertices[i]]     = derivs[i][0];
00302         expanded_derivs[2 * vertices[i] + 1] = derivs[i][1];
00303     }
00304 
00305     ASSERT_VALUES_EQUAL( expected_derivs[0], expanded_derivs[0], loc, bits );
00306     ASSERT_VALUES_EQUAL( expected_derivs[1], expanded_derivs[1], loc, bits );
00307     ASSERT_VALUES_EQUAL( expected_derivs[2], expanded_derivs[2], loc, bits );
00308     ASSERT_VALUES_EQUAL( expected_derivs[3], expanded_derivs[3], loc, bits );
00309     ASSERT_VALUES_EQUAL( expected_derivs[4], expanded_derivs[4], loc, bits );
00310     ASSERT_VALUES_EQUAL( expected_derivs[5], expanded_derivs[5], loc, bits );
00311     ASSERT_VALUES_EQUAL( expected_derivs[6], expanded_derivs[6], loc, bits );
00312     ASSERT_VALUES_EQUAL( expected_derivs[7], expanded_derivs[7], loc, bits );
00313     ASSERT_VALUES_EQUAL( expected_derivs[8], expanded_derivs[8], loc, bits );
00314     ASSERT_VALUES_EQUAL( expected_derivs[9], expanded_derivs[9], loc, bits );
00315     ASSERT_VALUES_EQUAL( expected_derivs[10], expanded_derivs[10], loc, bits );
00316     ASSERT_VALUES_EQUAL( expected_derivs[11], expanded_derivs[11], loc, bits );
00317 }
00318 
00319 void TriLagrangeShapeTest::test_corner_coeff( int corner, NodeSet nodebits )
00320 {
00321     MsqPrintError err( std::cout );
00322 
00323     double expected[6];
00324     get_coeff( nodebits, rs_corner[corner], expected );
00325 
00326     double coeff[27];
00327     size_t num_coeff = 17, indices[27];
00328     sf.coefficients( Sample( 0, corner ), nodebits, coeff, indices, num_coeff, err );
00329     CPPUNIT_ASSERT( !err );
00330 
00331     compare_coefficients( coeff, indices, expected, num_coeff, corner, nodebits );
00332 }
00333 
00334 void TriLagrangeShapeTest::test_edge_coeff( int edge, NodeSet nodebits )
00335 {
00336     MsqPrintError err( std::cout );
00337 
00338     double expected[6];
00339     get_coeff( nodebits, rs_edge[edge], expected );
00340 
00341     double coeff[27];
00342     size_t num_coeff = 17, indices[27];
00343     sf.coefficients( Sample( 1, edge ), nodebits, coeff, indices, num_coeff, err );
00344     CPPUNIT_ASSERT( !err );
00345 
00346     compare_coefficients( coeff, indices, expected, num_coeff, edge + 3, nodebits );
00347 }
00348 
00349 void TriLagrangeShapeTest::test_mid_coeff( NodeSet nodebits )
00350 {
00351     MsqPrintError err( std::cout );
00352 
00353     double expected[6];
00354     get_coeff( nodebits, rs_mid, expected );
00355 
00356     double coeff[27];
00357     size_t num_coeff = 17, indices[27];
00358     sf.coefficients( Sample( 2, 0 ), nodebits, coeff, indices, num_coeff, err );
00359     CPPUNIT_ASSERT( !err );
00360 
00361     compare_coefficients( coeff, indices, expected, num_coeff, 6, nodebits );
00362 }
00363 
00364 void TriLagrangeShapeTest::test_corner_derivs( int corner, NodeSet nodebits )
00365 {
00366     MsqPrintError err( std::cout );
00367 
00368     double expected[12];
00369     get_derivs( nodebits, rs_corner[corner], expected );
00370 
00371     size_t n = 19, vertices[100];
00372     MsqVector< 2 > derivs[100];
00373     sf.derivatives( Sample( 0, corner ), nodebits, vertices, derivs, n, err );
00374     CPPUNIT_ASSERT( !err );
00375 
00376     compare_derivatives( vertices, n, derivs, expected, corner, nodebits );
00377 }
00378 
00379 void TriLagrangeShapeTest::test_edge_derivs( int edge, NodeSet nodebits )
00380 {
00381     MsqPrintError err( std::cout );
00382 
00383     double expected[12];
00384     get_derivs( nodebits, rs_edge[edge], expected );
00385 
00386     size_t n = 19, vertices[100];
00387     MsqVector< 2 > derivs[100];
00388     sf.derivatives( Sample( 1, edge ), nodebits, vertices, derivs, n, err );
00389     CPPUNIT_ASSERT( !err );
00390 
00391     compare_derivatives( vertices, n, derivs, expected, edge + 3, nodebits );
00392 }
00393 
00394 void TriLagrangeShapeTest::test_mid_derivs( NodeSet nodebits )
00395 {
00396     MsqPrintError err( std::cout );
00397 
00398     double expected[12];
00399     get_derivs( nodebits, rs_mid, expected );
00400 
00401     size_t n = 19, vertices[100];
00402     MsqVector< 2 > derivs[100];
00403     sf.derivatives( Sample( 2, 0 ), nodebits, vertices, derivs, n, err );
00404     CPPUNIT_ASSERT( !err );
00405 
00406     compare_derivatives( vertices, n, derivs, expected, 6, nodebits );
00407 }
00408 
00409 void TriLagrangeShapeTest::test_coeff_corners()
00410 {
00411     NodeSet ns;
00412 
00413     ns.clear();
00414     test_corner_coeff( 0, ns );
00415     test_corner_coeff( 1, ns );
00416     test_corner_coeff( 2, ns );
00417 
00418     ns.set_mid_edge_node( 0 );
00419     test_corner_coeff( 0, ns );
00420     test_corner_coeff( 1, ns );
00421     test_corner_coeff( 2, ns );
00422 
00423     ns.clear();
00424     ns.set_mid_edge_node( 1 );
00425     test_corner_coeff( 0, ns );
00426     test_corner_coeff( 1, ns );
00427     test_corner_coeff( 2, ns );
00428 
00429     ns.set_mid_edge_node( 0 );
00430     test_corner_coeff( 0, ns );
00431     test_corner_coeff( 1, ns );
00432     test_corner_coeff( 2, ns );
00433 
00434     ns.clear();
00435     ns.set_mid_edge_node( 2 );
00436     test_corner_coeff( 0, ns );
00437     test_corner_coeff( 1, ns );
00438     test_corner_coeff( 2, ns );
00439 
00440     ns.set_mid_edge_node( 0 );
00441     test_corner_coeff( 0, ns );
00442     test_corner_coeff( 1, ns );
00443     test_corner_coeff( 2, ns );
00444 
00445     ns.set_mid_edge_node( 1 );
00446     test_corner_coeff( 0, ns );
00447     test_corner_coeff( 1, ns );
00448     test_corner_coeff( 2, ns );
00449 
00450     ns.clear_mid_edge_node( 0 );
00451     test_corner_coeff( 0, ns );
00452     test_corner_coeff( 1, ns );
00453     test_corner_coeff( 2, ns );
00454 }
00455 
00456 void TriLagrangeShapeTest::test_coeff_edges()
00457 {
00458     NodeSet ns;
00459 
00460     ns.clear();
00461     test_edge_coeff( 0, ns );
00462     test_edge_coeff( 1, ns );
00463     test_edge_coeff( 2, ns );
00464 
00465     ns.set_mid_edge_node( 0 );
00466     test_edge_coeff( 0, ns );
00467     test_edge_coeff( 1, ns );
00468     test_edge_coeff( 2, ns );
00469 
00470     ns.clear();
00471     ns.set_mid_edge_node( 1 );
00472     test_edge_coeff( 0, ns );
00473     test_edge_coeff( 1, ns );
00474     test_edge_coeff( 2, ns );
00475 
00476     ns.set_mid_edge_node( 0 );
00477     test_edge_coeff( 0, ns );
00478     test_edge_coeff( 1, ns );
00479     test_edge_coeff( 2, ns );
00480 
00481     ns.clear();
00482     ns.set_mid_edge_node( 2 );
00483     test_edge_coeff( 0, ns );
00484     test_edge_coeff( 1, ns );
00485     test_edge_coeff( 2, ns );
00486 
00487     ns.set_mid_edge_node( 0 );
00488     test_edge_coeff( 0, ns );
00489     test_edge_coeff( 1, ns );
00490     test_edge_coeff( 2, ns );
00491 
00492     ns.set_mid_edge_node( 1 );
00493     test_edge_coeff( 0, ns );
00494     test_edge_coeff( 1, ns );
00495     test_edge_coeff( 2, ns );
00496 
00497     ns.clear_mid_edge_node( 0 );
00498     test_edge_coeff( 0, ns );
00499     test_edge_coeff( 1, ns );
00500     test_edge_coeff( 2, ns );
00501 }
00502 
00503 void TriLagrangeShapeTest::test_coeff_center()
00504 {
00505     NodeSet ns;
00506 
00507     ns.clear();
00508     test_mid_coeff( ns );
00509     ns.set_mid_edge_node( 0 );
00510     test_mid_coeff( ns );
00511     ns.clear();
00512     ns.set_mid_edge_node( 1 );
00513     test_mid_coeff( ns );
00514     ns.set_mid_edge_node( 0 );
00515     test_mid_coeff( ns );
00516     ns.clear();
00517     ns.set_mid_edge_node( 2 );
00518     test_mid_coeff( ns );
00519     ns.set_mid_edge_node( 0 );
00520     test_mid_coeff( ns );
00521     ns.set_mid_edge_node( 1 );
00522     test_mid_coeff( ns );
00523     ns.clear_mid_edge_node( 0 );
00524     test_mid_coeff( ns );
00525 }
00526 
00527 void TriLagrangeShapeTest::test_deriv_corners()
00528 {
00529     NodeSet ns;
00530 
00531     ns.clear();
00532     test_corner_derivs( 0, ns );
00533     test_corner_derivs( 1, ns );
00534     test_corner_derivs( 2, ns );
00535 
00536     ns.set_mid_edge_node( 0 );
00537     test_corner_derivs( 0, ns );
00538     test_corner_derivs( 1, ns );
00539     test_corner_derivs( 2, ns );
00540 
00541     ns.clear();
00542     ns.set_mid_edge_node( 1 );
00543     test_corner_derivs( 0, ns );
00544     test_corner_derivs( 1, ns );
00545     test_corner_derivs( 2, ns );
00546 
00547     ns.set_mid_edge_node( 0 );
00548     test_corner_derivs( 0, ns );
00549     test_corner_derivs( 1, ns );
00550     test_corner_derivs( 2, ns );
00551 
00552     ns.clear();
00553     ns.set_mid_edge_node( 2 );
00554     test_corner_derivs( 0, ns );
00555     test_corner_derivs( 1, ns );
00556     test_corner_derivs( 2, ns );
00557 
00558     ns.set_mid_edge_node( 0 );
00559     test_corner_derivs( 0, ns );
00560     test_corner_derivs( 1, ns );
00561     test_corner_derivs( 2, ns );
00562 
00563     ns.set_mid_edge_node( 1 );
00564     test_corner_derivs( 0, ns );
00565     test_corner_derivs( 1, ns );
00566     test_corner_derivs( 2, ns );
00567 
00568     ns.clear_mid_edge_node( 0 );
00569     test_corner_derivs( 0, ns );
00570     test_corner_derivs( 1, ns );
00571     test_corner_derivs( 2, ns );
00572 }
00573 
00574 void TriLagrangeShapeTest::test_deriv_edges()
00575 {
00576     NodeSet ns;
00577 
00578     ns.clear();
00579     test_edge_derivs( 0, ns );
00580     test_edge_derivs( 1, ns );
00581     test_edge_derivs( 2, ns );
00582 
00583     ns.set_mid_edge_node( 0 );
00584     test_edge_derivs( 0, ns );
00585     test_edge_derivs( 1, ns );
00586     test_edge_derivs( 2, ns );
00587 
00588     ns.clear();
00589     ns.set_mid_edge_node( 1 );
00590     test_edge_derivs( 0, ns );
00591     test_edge_derivs( 1, ns );
00592     test_edge_derivs( 2, ns );
00593 
00594     ns.set_mid_edge_node( 0 );
00595     test_edge_derivs( 0, ns );
00596     test_edge_derivs( 1, ns );
00597     test_edge_derivs( 2, ns );
00598 
00599     ns.clear();
00600     ns.set_mid_edge_node( 2 );
00601     test_edge_derivs( 0, ns );
00602     test_edge_derivs( 1, ns );
00603     test_edge_derivs( 2, ns );
00604 
00605     ns.set_mid_edge_node( 0 );
00606     test_edge_derivs( 0, ns );
00607     test_edge_derivs( 1, ns );
00608     test_edge_derivs( 2, ns );
00609 
00610     ns.set_mid_edge_node( 1 );
00611     test_edge_derivs( 0, ns );
00612     test_edge_derivs( 1, ns );
00613     test_edge_derivs( 2, ns );
00614 
00615     ns.clear_mid_edge_node( 0 );
00616     test_edge_derivs( 0, ns );
00617     test_edge_derivs( 1, ns );
00618     test_edge_derivs( 2, ns );
00619 }
00620 
00621 void TriLagrangeShapeTest::test_deriv_center()
00622 {
00623     NodeSet ns;
00624 
00625     ns.clear();
00626     test_mid_derivs( ns );
00627     ns.set_mid_edge_node( 0 );
00628     test_mid_derivs( ns );
00629     ns.clear();
00630     ns.set_mid_edge_node( 1 );
00631     test_mid_derivs( ns );
00632     ns.set_mid_edge_node( 0 );
00633     test_mid_derivs( ns );
00634     ns.clear();
00635     ns.set_mid_edge_node( 2 );
00636     test_mid_derivs( ns );
00637     ns.set_mid_edge_node( 0 );
00638     test_mid_derivs( ns );
00639     ns.set_mid_edge_node( 1 );
00640     test_mid_derivs( ns );
00641     ns.clear_mid_edge_node( 0 );
00642     test_mid_derivs( ns );
00643 }
00644 
00645 void TriLagrangeShapeTest::test_ideal_jacobian()
00646 {
00647     MsqError err;
00648     MsqMatrix< 3, 2 > J_prime;
00649     sf.ideal( Sample( 2, 0 ), J_prime, err );
00650     ASSERT_NO_ERROR( err );
00651 
00652     // for this test that everything is in the xy-plane
00653     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, J_prime( 2, 0 ), 1e-12 );
00654     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, J_prime( 2, 1 ), 1e-12 );
00655     MsqMatrix< 2, 2 > J_act( J_prime.data() );
00656     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( J_act ), 1e-6 );
00657 
00658     const Vector3D* verts = unit_edge_element( TRIANGLE );
00659     CPPUNIT_ASSERT( verts );
00660 
00661     JacobianCalculator jc;
00662     jc.get_Jacobian_2D( &sf, NodeSet(), Sample( 2, 0 ), verts, 3, J_prime, err );
00663     ASSERT_NO_ERROR( err );
00664 
00665     // for this test that everything is in the xy-plane
00666     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, J_prime( 2, 0 ), 1e-12 );
00667     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, J_prime( 2, 1 ), 1e-12 );
00668     MsqMatrix< 2, 2 > J_exp( J_prime.data() );
00669     J_exp /= sqrt( det( J_exp ) );
00670 
00671     // Matrices should be a rotation of each other.
00672     // First, calculate tentative rotation matrix
00673     MsqMatrix< 2, 2 > R = inverse( J_exp ) * J_act;
00674     // next check that it is a rotation
00675     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( R ), 1e-6 );          // no scaling
00676     ASSERT_MATRICES_EQUAL( transpose( R ), inverse( R ), 1e-6 );  // orthogonal
00677 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines