MOAB: Mesh Oriented datABase  (version 5.4.0)
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,
00267                                   const size_t* indices,
00268                                   const double* expected_coeffs,
00269                                   size_t num_coeff,
00270                                   unsigned loc,
00271                                   NodeSet bits )
00272 {
00273     // find the location in the returned list for each node
00274     size_t revidx[6];
00275     double test_vals[6];
00276     for( size_t i = 0; i < 6; ++i )
00277     {
00278         revidx[i]    = std::find( indices, indices + num_coeff, i ) - indices;
00279         test_vals[i] = ( revidx[i] == num_coeff ) ? 0.0 : coeffs[revidx[i]];
00280     }
00281 
00282     // Check that index list doesn't contain any nodes not actually
00283     // present in the element.
00284     CPPUNIT_ASSERT( bits.mid_edge_node( 0 ) || ( revidx[3] == num_coeff ) );
00285     CPPUNIT_ASSERT( bits.mid_edge_node( 1 ) || ( revidx[4] == num_coeff ) );
00286     CPPUNIT_ASSERT( bits.mid_edge_node( 2 ) || ( revidx[5] == num_coeff ) );
00287 
00288     // compare expected and actual coefficient values
00289     ASSERT_VALUES_EQUAL( expected_coeffs[0], test_vals[0], loc, bits );
00290     ASSERT_VALUES_EQUAL( expected_coeffs[1], test_vals[1], loc, bits );
00291     ASSERT_VALUES_EQUAL( expected_coeffs[2], test_vals[2], loc, bits );
00292     ASSERT_VALUES_EQUAL( expected_coeffs[3], test_vals[3], loc, bits );
00293     ASSERT_VALUES_EQUAL( expected_coeffs[4], test_vals[4], loc, bits );
00294     ASSERT_VALUES_EQUAL( expected_coeffs[5], test_vals[5], loc, bits );
00295 }
00296 
00297 static void compare_derivatives( const size_t* vertices,
00298                                  size_t num_vtx,
00299                                  const MsqVector< 2 >* derivs,
00300                                  const double* expected_derivs,
00301                                  unsigned loc,
00302                                  NodeSet bits )
00303 {
00304     check_valid_indices( vertices, num_vtx );
00305     check_no_zeros( derivs, num_vtx );
00306     double expanded_derivs[12] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
00307     for( unsigned i = 0; i < num_vtx; ++i )
00308     {
00309         expanded_derivs[2 * vertices[i]]     = derivs[i][0];
00310         expanded_derivs[2 * vertices[i] + 1] = derivs[i][1];
00311     }
00312 
00313     ASSERT_VALUES_EQUAL( expected_derivs[0], expanded_derivs[0], loc, bits );
00314     ASSERT_VALUES_EQUAL( expected_derivs[1], expanded_derivs[1], loc, bits );
00315     ASSERT_VALUES_EQUAL( expected_derivs[2], expanded_derivs[2], loc, bits );
00316     ASSERT_VALUES_EQUAL( expected_derivs[3], expanded_derivs[3], loc, bits );
00317     ASSERT_VALUES_EQUAL( expected_derivs[4], expanded_derivs[4], loc, bits );
00318     ASSERT_VALUES_EQUAL( expected_derivs[5], expanded_derivs[5], loc, bits );
00319     ASSERT_VALUES_EQUAL( expected_derivs[6], expanded_derivs[6], loc, bits );
00320     ASSERT_VALUES_EQUAL( expected_derivs[7], expanded_derivs[7], loc, bits );
00321     ASSERT_VALUES_EQUAL( expected_derivs[8], expanded_derivs[8], loc, bits );
00322     ASSERT_VALUES_EQUAL( expected_derivs[9], expanded_derivs[9], loc, bits );
00323     ASSERT_VALUES_EQUAL( expected_derivs[10], expanded_derivs[10], loc, bits );
00324     ASSERT_VALUES_EQUAL( expected_derivs[11], expanded_derivs[11], loc, bits );
00325 }
00326 
00327 void TriLagrangeShapeTest::test_corner_coeff( int corner, NodeSet nodebits )
00328 {
00329     MsqPrintError err( std::cout );
00330 
00331     double expected[6];
00332     get_coeff( nodebits, rs_corner[corner], expected );
00333 
00334     double coeff[27];
00335     size_t num_coeff = 17, indices[27];
00336     sf.coefficients( Sample( 0, corner ), nodebits, coeff, indices, num_coeff, err );
00337     CPPUNIT_ASSERT( !err );
00338 
00339     compare_coefficients( coeff, indices, expected, num_coeff, corner, nodebits );
00340 }
00341 
00342 void TriLagrangeShapeTest::test_edge_coeff( int edge, NodeSet nodebits )
00343 {
00344     MsqPrintError err( std::cout );
00345 
00346     double expected[6];
00347     get_coeff( nodebits, rs_edge[edge], expected );
00348 
00349     double coeff[27];
00350     size_t num_coeff = 17, indices[27];
00351     sf.coefficients( Sample( 1, edge ), nodebits, coeff, indices, num_coeff, err );
00352     CPPUNIT_ASSERT( !err );
00353 
00354     compare_coefficients( coeff, indices, expected, num_coeff, edge + 3, nodebits );
00355 }
00356 
00357 void TriLagrangeShapeTest::test_mid_coeff( NodeSet nodebits )
00358 {
00359     MsqPrintError err( std::cout );
00360 
00361     double expected[6];
00362     get_coeff( nodebits, rs_mid, expected );
00363 
00364     double coeff[27];
00365     size_t num_coeff = 17, indices[27];
00366     sf.coefficients( Sample( 2, 0 ), nodebits, coeff, indices, num_coeff, err );
00367     CPPUNIT_ASSERT( !err );
00368 
00369     compare_coefficients( coeff, indices, expected, num_coeff, 6, nodebits );
00370 }
00371 
00372 void TriLagrangeShapeTest::test_corner_derivs( int corner, NodeSet nodebits )
00373 {
00374     MsqPrintError err( std::cout );
00375 
00376     double expected[12];
00377     get_derivs( nodebits, rs_corner[corner], expected );
00378 
00379     size_t n = 19, vertices[100];
00380     MsqVector< 2 > derivs[100];
00381     sf.derivatives( Sample( 0, corner ), nodebits, vertices, derivs, n, err );
00382     CPPUNIT_ASSERT( !err );
00383 
00384     compare_derivatives( vertices, n, derivs, expected, corner, nodebits );
00385 }
00386 
00387 void TriLagrangeShapeTest::test_edge_derivs( int edge, NodeSet nodebits )
00388 {
00389     MsqPrintError err( std::cout );
00390 
00391     double expected[12];
00392     get_derivs( nodebits, rs_edge[edge], expected );
00393 
00394     size_t n = 19, vertices[100];
00395     MsqVector< 2 > derivs[100];
00396     sf.derivatives( Sample( 1, edge ), nodebits, vertices, derivs, n, err );
00397     CPPUNIT_ASSERT( !err );
00398 
00399     compare_derivatives( vertices, n, derivs, expected, edge + 3, nodebits );
00400 }
00401 
00402 void TriLagrangeShapeTest::test_mid_derivs( NodeSet nodebits )
00403 {
00404     MsqPrintError err( std::cout );
00405 
00406     double expected[12];
00407     get_derivs( nodebits, rs_mid, expected );
00408 
00409     size_t n = 19, vertices[100];
00410     MsqVector< 2 > derivs[100];
00411     sf.derivatives( Sample( 2, 0 ), nodebits, vertices, derivs, n, err );
00412     CPPUNIT_ASSERT( !err );
00413 
00414     compare_derivatives( vertices, n, derivs, expected, 6, nodebits );
00415 }
00416 
00417 void TriLagrangeShapeTest::test_coeff_corners()
00418 {
00419     NodeSet ns;
00420 
00421     ns.clear();
00422     test_corner_coeff( 0, ns );
00423     test_corner_coeff( 1, ns );
00424     test_corner_coeff( 2, ns );
00425 
00426     ns.set_mid_edge_node( 0 );
00427     test_corner_coeff( 0, ns );
00428     test_corner_coeff( 1, ns );
00429     test_corner_coeff( 2, ns );
00430 
00431     ns.clear();
00432     ns.set_mid_edge_node( 1 );
00433     test_corner_coeff( 0, ns );
00434     test_corner_coeff( 1, ns );
00435     test_corner_coeff( 2, ns );
00436 
00437     ns.set_mid_edge_node( 0 );
00438     test_corner_coeff( 0, ns );
00439     test_corner_coeff( 1, ns );
00440     test_corner_coeff( 2, ns );
00441 
00442     ns.clear();
00443     ns.set_mid_edge_node( 2 );
00444     test_corner_coeff( 0, ns );
00445     test_corner_coeff( 1, ns );
00446     test_corner_coeff( 2, ns );
00447 
00448     ns.set_mid_edge_node( 0 );
00449     test_corner_coeff( 0, ns );
00450     test_corner_coeff( 1, ns );
00451     test_corner_coeff( 2, ns );
00452 
00453     ns.set_mid_edge_node( 1 );
00454     test_corner_coeff( 0, ns );
00455     test_corner_coeff( 1, ns );
00456     test_corner_coeff( 2, ns );
00457 
00458     ns.clear_mid_edge_node( 0 );
00459     test_corner_coeff( 0, ns );
00460     test_corner_coeff( 1, ns );
00461     test_corner_coeff( 2, ns );
00462 }
00463 
00464 void TriLagrangeShapeTest::test_coeff_edges()
00465 {
00466     NodeSet ns;
00467 
00468     ns.clear();
00469     test_edge_coeff( 0, ns );
00470     test_edge_coeff( 1, ns );
00471     test_edge_coeff( 2, ns );
00472 
00473     ns.set_mid_edge_node( 0 );
00474     test_edge_coeff( 0, ns );
00475     test_edge_coeff( 1, ns );
00476     test_edge_coeff( 2, ns );
00477 
00478     ns.clear();
00479     ns.set_mid_edge_node( 1 );
00480     test_edge_coeff( 0, ns );
00481     test_edge_coeff( 1, ns );
00482     test_edge_coeff( 2, ns );
00483 
00484     ns.set_mid_edge_node( 0 );
00485     test_edge_coeff( 0, ns );
00486     test_edge_coeff( 1, ns );
00487     test_edge_coeff( 2, ns );
00488 
00489     ns.clear();
00490     ns.set_mid_edge_node( 2 );
00491     test_edge_coeff( 0, ns );
00492     test_edge_coeff( 1, ns );
00493     test_edge_coeff( 2, ns );
00494 
00495     ns.set_mid_edge_node( 0 );
00496     test_edge_coeff( 0, ns );
00497     test_edge_coeff( 1, ns );
00498     test_edge_coeff( 2, ns );
00499 
00500     ns.set_mid_edge_node( 1 );
00501     test_edge_coeff( 0, ns );
00502     test_edge_coeff( 1, ns );
00503     test_edge_coeff( 2, ns );
00504 
00505     ns.clear_mid_edge_node( 0 );
00506     test_edge_coeff( 0, ns );
00507     test_edge_coeff( 1, ns );
00508     test_edge_coeff( 2, ns );
00509 }
00510 
00511 void TriLagrangeShapeTest::test_coeff_center()
00512 {
00513     NodeSet ns;
00514 
00515     ns.clear();
00516     test_mid_coeff( ns );
00517     ns.set_mid_edge_node( 0 );
00518     test_mid_coeff( ns );
00519     ns.clear();
00520     ns.set_mid_edge_node( 1 );
00521     test_mid_coeff( ns );
00522     ns.set_mid_edge_node( 0 );
00523     test_mid_coeff( ns );
00524     ns.clear();
00525     ns.set_mid_edge_node( 2 );
00526     test_mid_coeff( ns );
00527     ns.set_mid_edge_node( 0 );
00528     test_mid_coeff( ns );
00529     ns.set_mid_edge_node( 1 );
00530     test_mid_coeff( ns );
00531     ns.clear_mid_edge_node( 0 );
00532     test_mid_coeff( ns );
00533 }
00534 
00535 void TriLagrangeShapeTest::test_deriv_corners()
00536 {
00537     NodeSet ns;
00538 
00539     ns.clear();
00540     test_corner_derivs( 0, ns );
00541     test_corner_derivs( 1, ns );
00542     test_corner_derivs( 2, ns );
00543 
00544     ns.set_mid_edge_node( 0 );
00545     test_corner_derivs( 0, ns );
00546     test_corner_derivs( 1, ns );
00547     test_corner_derivs( 2, ns );
00548 
00549     ns.clear();
00550     ns.set_mid_edge_node( 1 );
00551     test_corner_derivs( 0, ns );
00552     test_corner_derivs( 1, ns );
00553     test_corner_derivs( 2, ns );
00554 
00555     ns.set_mid_edge_node( 0 );
00556     test_corner_derivs( 0, ns );
00557     test_corner_derivs( 1, ns );
00558     test_corner_derivs( 2, ns );
00559 
00560     ns.clear();
00561     ns.set_mid_edge_node( 2 );
00562     test_corner_derivs( 0, ns );
00563     test_corner_derivs( 1, ns );
00564     test_corner_derivs( 2, ns );
00565 
00566     ns.set_mid_edge_node( 0 );
00567     test_corner_derivs( 0, ns );
00568     test_corner_derivs( 1, ns );
00569     test_corner_derivs( 2, ns );
00570 
00571     ns.set_mid_edge_node( 1 );
00572     test_corner_derivs( 0, ns );
00573     test_corner_derivs( 1, ns );
00574     test_corner_derivs( 2, ns );
00575 
00576     ns.clear_mid_edge_node( 0 );
00577     test_corner_derivs( 0, ns );
00578     test_corner_derivs( 1, ns );
00579     test_corner_derivs( 2, ns );
00580 }
00581 
00582 void TriLagrangeShapeTest::test_deriv_edges()
00583 {
00584     NodeSet ns;
00585 
00586     ns.clear();
00587     test_edge_derivs( 0, ns );
00588     test_edge_derivs( 1, ns );
00589     test_edge_derivs( 2, ns );
00590 
00591     ns.set_mid_edge_node( 0 );
00592     test_edge_derivs( 0, ns );
00593     test_edge_derivs( 1, ns );
00594     test_edge_derivs( 2, ns );
00595 
00596     ns.clear();
00597     ns.set_mid_edge_node( 1 );
00598     test_edge_derivs( 0, ns );
00599     test_edge_derivs( 1, ns );
00600     test_edge_derivs( 2, ns );
00601 
00602     ns.set_mid_edge_node( 0 );
00603     test_edge_derivs( 0, ns );
00604     test_edge_derivs( 1, ns );
00605     test_edge_derivs( 2, ns );
00606 
00607     ns.clear();
00608     ns.set_mid_edge_node( 2 );
00609     test_edge_derivs( 0, ns );
00610     test_edge_derivs( 1, ns );
00611     test_edge_derivs( 2, ns );
00612 
00613     ns.set_mid_edge_node( 0 );
00614     test_edge_derivs( 0, ns );
00615     test_edge_derivs( 1, ns );
00616     test_edge_derivs( 2, ns );
00617 
00618     ns.set_mid_edge_node( 1 );
00619     test_edge_derivs( 0, ns );
00620     test_edge_derivs( 1, ns );
00621     test_edge_derivs( 2, ns );
00622 
00623     ns.clear_mid_edge_node( 0 );
00624     test_edge_derivs( 0, ns );
00625     test_edge_derivs( 1, ns );
00626     test_edge_derivs( 2, ns );
00627 }
00628 
00629 void TriLagrangeShapeTest::test_deriv_center()
00630 {
00631     NodeSet ns;
00632 
00633     ns.clear();
00634     test_mid_derivs( ns );
00635     ns.set_mid_edge_node( 0 );
00636     test_mid_derivs( ns );
00637     ns.clear();
00638     ns.set_mid_edge_node( 1 );
00639     test_mid_derivs( ns );
00640     ns.set_mid_edge_node( 0 );
00641     test_mid_derivs( ns );
00642     ns.clear();
00643     ns.set_mid_edge_node( 2 );
00644     test_mid_derivs( ns );
00645     ns.set_mid_edge_node( 0 );
00646     test_mid_derivs( ns );
00647     ns.set_mid_edge_node( 1 );
00648     test_mid_derivs( ns );
00649     ns.clear_mid_edge_node( 0 );
00650     test_mid_derivs( ns );
00651 }
00652 
00653 void TriLagrangeShapeTest::test_ideal_jacobian()
00654 {
00655     MsqError err;
00656     MsqMatrix< 3, 2 > J_prime;
00657     sf.ideal( Sample( 2, 0 ), J_prime, err );
00658     ASSERT_NO_ERROR( err );
00659 
00660     // for this test that everything is in the xy-plane
00661     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, J_prime( 2, 0 ), 1e-12 );
00662     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, J_prime( 2, 1 ), 1e-12 );
00663     MsqMatrix< 2, 2 > J_act( J_prime.data() );
00664     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( J_act ), 1e-6 );
00665 
00666     const Vector3D* verts = unit_edge_element( TRIANGLE );
00667     CPPUNIT_ASSERT( verts );
00668 
00669     JacobianCalculator jc;
00670     jc.get_Jacobian_2D( &sf, NodeSet(), Sample( 2, 0 ), verts, 3, J_prime, err );
00671     ASSERT_NO_ERROR( err );
00672 
00673     // for this test that everything is in the xy-plane
00674     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, J_prime( 2, 0 ), 1e-12 );
00675     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, J_prime( 2, 1 ), 1e-12 );
00676     MsqMatrix< 2, 2 > J_exp( J_prime.data() );
00677     J_exp /= sqrt( det( J_exp ) );
00678 
00679     // Matrices should be a rotation of each other.
00680     // First, calculate tentative rotation matrix
00681     MsqMatrix< 2, 2 > R = inverse( J_exp ) * J_act;
00682     // next check that it is a rotation
00683     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( R ), 1e-6 );          // no scaling
00684     ASSERT_MATRICES_EQUAL( transpose( R ), inverse( R ), 1e-6 );  // orthogonal
00685 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines