MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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 (2010) kraftche@cae.wisc.edu 00024 00025 ***************************************************************** */ 00026 00027 /** \file HexLagrangeShapeTest.cpp 00028 * \brief 00029 * \author Jason Kraftcheck 00030 */ 00031 00032 #include "Mesquite.hpp" 00033 #include "HexLagrangeShape.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 ) \ 00050 ASSERT_MESSAGE( value_message( ( location ), ( v1 ), ( v2 ) ), test_value( ( v1 ), ( v2 ) ) ) 00051 00052 static bool test_value( double v1, double v2 ) 00053 { 00054 return fabs( v1 - v2 ) < epsilon; 00055 } 00056 00057 static inline CppUnit::Message value_message( unsigned location, double v1, double v2 ) 00058 { 00059 CppUnit::Message m( "equality assertion failed" ); 00060 00061 std::ostringstream buffer1; 00062 buffer1 << "Expected : " << v1; 00063 m.addDetail( buffer1.str() ); 00064 00065 std::ostringstream buffer2; 00066 buffer2 << "Actual : " << v2; 00067 m.addDetail( buffer2.str() ); 00068 00069 std::ostringstream buffer3; 00070 buffer3 << "Location : "; 00071 if( location < 9 ) 00072 buffer3 << "Corner " << location; 00073 else if( location < 20 ) 00074 buffer3 << "Edge " << location - 8; 00075 else if( location < 26 ) 00076 buffer3 << "Face " << location - 20; 00077 else if( location == 26 ) 00078 buffer3 << "Mid-element"; 00079 else 00080 buffer3 << "INVALID!!"; 00081 m.addDetail( buffer3.str() ); 00082 return m; 00083 } 00084 00085 static bool test_value( MsqVector< 3 > v1, MsqVector< 3 > v2 ) 00086 { 00087 return length( v1 - v2 ) < epsilon; 00088 } 00089 00090 static inline CppUnit::Message value_message( unsigned location, MsqVector< 3 > v1, MsqVector< 3 > v2 ) 00091 { 00092 CppUnit::Message m( "equality assertion failed" ); 00093 00094 std::ostringstream buffer1; 00095 buffer1 << "Expected : [" << v1[0] << ", " << v1[1] << ", " << v1[2] << "]"; 00096 m.addDetail( buffer1.str() ); 00097 00098 std::ostringstream buffer2; 00099 buffer2 << "Actual : [" << v2[0] << ", " << v2[1] << ", " << v2[2] << "]"; 00100 m.addDetail( buffer2.str() ); 00101 00102 std::ostringstream buffer3; 00103 buffer3 << "Location : "; 00104 if( location < 9 ) 00105 buffer3 << "Corner " << location; 00106 else if( location < 20 ) 00107 buffer3 << "Edge " << location - 8; 00108 else if( location < 26 ) 00109 buffer3 << "Face " << location - 20; 00110 else if( location == 26 ) 00111 buffer3 << "Mid-element"; 00112 else 00113 buffer3 << "INVALID!!"; 00114 m.addDetail( buffer3.str() ); 00115 return m; 00116 } 00117 00118 class HexLagrangeShapeTest : public CppUnit::TestFixture 00119 { 00120 private: 00121 CPPUNIT_TEST_SUITE( HexLagrangeShapeTest ); 00122 00123 CPPUNIT_TEST( test_corners_coeff ); 00124 CPPUNIT_TEST( test_edges_coeff ); 00125 CPPUNIT_TEST( test_faces_coeff ); 00126 CPPUNIT_TEST( test_mid_coeff ); 00127 00128 CPPUNIT_TEST( test_corners_derivs ); 00129 CPPUNIT_TEST( test_edges_derivs ); 00130 CPPUNIT_TEST( test_faces_derivs ); 00131 CPPUNIT_TEST( test_mid_derivs ); 00132 00133 CPPUNIT_TEST( test_ideal_jacobian ); 00134 00135 CPPUNIT_TEST_SUITE_END(); 00136 00137 HexLagrangeShape sf; 00138 00139 public: 00140 void test_corner_coeff( int corner ); 00141 void test_edge_coeff( int edge ); 00142 void test_face_coeff( int face ); 00143 void test_mid_coeff(); 00144 00145 void test_corner_derivs( int corner ); 00146 void test_edge_derivs( int edge ); 00147 void test_face_derivs( int face ); 00148 void test_mid_derivs(); 00149 00150 NodeSet allNodes; 00151 void setUp() 00152 { 00153 allNodes.set_all_nodes( HEXAHEDRON ); 00154 } 00155 00156 void test_corners_coeff(); 00157 void test_edges_coeff(); 00158 void test_faces_coeff(); 00159 00160 void test_corners_derivs(); 00161 void test_edges_derivs(); 00162 void test_faces_derivs(); 00163 00164 void test_ideal_jacobian(); 00165 }; 00166 00167 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( HexLagrangeShapeTest, "HexLagrangeShapeTest" ); 00168 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( HexLagrangeShapeTest, "Unit" ); 00169 00170 // Indices 00171 enum 00172 { 00173 XI = 0, 00174 ETA = 1, 00175 ZETA = 2 00176 }; 00177 // Parameter range 00178 const double min_xi = 0; 00179 const double max_xi = 1; 00180 const double mid_xi = 0.5 * ( min_xi + max_xi ); 00181 00182 // Some lagrange polynomial terms 00183 static double l21( double xi ) 00184 { 00185 return ( 2 * xi - 1 ) * ( xi - 1 ); 00186 } 00187 static double l22( double xi ) 00188 { 00189 return 4 * xi * ( 1 - xi ); 00190 } 00191 static double l23( double xi ) 00192 { 00193 return xi * ( 2 * xi - 1 ); 00194 } 00195 // First derivatives of lagrange polynomial terms 00196 static double dl21( double xi ) 00197 { 00198 return 4 * xi - 3; 00199 } 00200 static double dl22( double xi ) 00201 { 00202 return 4 - 8 * xi; 00203 } 00204 static double dl23( double xi ) 00205 { 00206 return 4 * xi - 1; 00207 } 00208 // Indexible polynomail terms 00209 typedef double ( *f_t )( double ); 00210 f_t l[4] = { 0, &l21, &l22, &l23 }; 00211 f_t dl[4] = { 0, &dl21, &dl22, &dl23 }; 00212 const int C[][3] = { 00213 { 1, 1, 1 }, // 0 00214 { 3, 1, 1 }, // 1 00215 { 3, 3, 1 }, // 2 00216 { 1, 3, 1 }, // 3 00217 00218 { 1, 1, 3 }, // 4 00219 { 3, 1, 3 }, // 5 00220 { 3, 3, 3 }, // 6 00221 { 1, 3, 3 }, // 7 00222 00223 { 2, 1, 1 }, // 8 00224 { 3, 2, 1 }, // 9 00225 { 2, 3, 1 }, // 10 00226 { 1, 2, 1 }, // 11 00227 00228 { 1, 1, 2 }, // 12 00229 { 3, 1, 2 }, // 13 00230 { 3, 3, 2 }, // 14 00231 { 1, 3, 2 }, // 15 00232 00233 { 2, 1, 3 }, // 16 00234 { 3, 2, 3 }, // 17 00235 { 2, 3, 3 }, // 18 00236 { 1, 2, 3 }, // 19 00237 00238 { 2, 1, 2 }, // 20 00239 { 3, 2, 2 }, // 21 00240 { 2, 3, 2 }, // 22 00241 { 1, 2, 2 }, // 23 00242 00243 { 2, 2, 1 }, // 24 00244 { 2, 2, 3 }, // 25 00245 00246 { 2, 2, 2 } // 26 00247 }; 00248 00249 const double corners[8][3] = { { min_xi, min_xi, min_xi }, { max_xi, min_xi, min_xi }, { max_xi, max_xi, min_xi }, 00250 { min_xi, max_xi, min_xi }, { min_xi, min_xi, max_xi }, { max_xi, min_xi, max_xi }, 00251 { max_xi, max_xi, max_xi }, { min_xi, max_xi, max_xi } }; 00252 00253 static MsqVector< 3 > XI_corner( int i ) 00254 { 00255 return MsqVector< 3 >( corners[i] ); 00256 } 00257 00258 static MsqVector< 3 > XI_edge( int i ) 00259 { 00260 const unsigned* idx = TopologyInfo::edge_vertices( HEXAHEDRON, i ); 00261 return 0.5 * ( XI_corner( idx[0] ) + XI_corner( idx[1] ) ); 00262 } 00263 00264 static MsqVector< 3 > XI_face( int i ) 00265 { 00266 unsigned n; 00267 const unsigned* idx = TopologyInfo::face_vertices( HEXAHEDRON, i, n ); 00268 MsqVector< 3 > result = XI_corner( idx[0] ); 00269 for( unsigned j = 1; j < n; ++j ) 00270 result += XI_corner( idx[j] ); 00271 result *= 1.0 / n; 00272 return result; 00273 } 00274 00275 static MsqVector< 3 > XI_elem() 00276 { 00277 const double vals[] = { mid_xi, mid_xi, mid_xi }; 00278 return MsqVector< 3 >( vals ); 00279 } 00280 00281 static double N( int i, MsqVector< 3 > xi ) 00282 { 00283 const int* c = C[i]; 00284 return l[c[XI]]( xi[XI] ) * l[c[ETA]]( xi[ETA] ) * l[c[ZETA]]( xi[ZETA] ); 00285 } 00286 00287 static MsqVector< 3 > dN( int i, MsqVector< 3 > xi ) 00288 { 00289 MsqVector< 3 > result; 00290 const int* c = C[i]; 00291 result[XI] = dl[c[XI]]( xi[XI] ) * l[c[ETA]]( xi[ETA] ) * l[c[ZETA]]( xi[ZETA] ); 00292 result[ETA] = l[c[XI]]( xi[XI] ) * dl[c[ETA]]( xi[ETA] ) * l[c[ZETA]]( xi[ZETA] ); 00293 result[ZETA] = l[c[XI]]( xi[XI] ) * l[c[ETA]]( xi[ETA] ) * dl[c[ZETA]]( xi[ZETA] ); 00294 return result; 00295 } 00296 00297 static void get_coeffs( MsqVector< 3 > xi, double coeffs_out[27] ) 00298 { 00299 for( int i = 0; i < 27; ++i ) 00300 coeffs_out[i] = N( i, xi ); 00301 } 00302 00303 static void get_derivs( MsqVector< 3 > xi, MsqVector< 3 > derivs_out[27] ) 00304 { 00305 for( int i = 0; i < 27; ++i ) 00306 derivs_out[i] = dN( i, xi ); 00307 } 00308 00309 static void check_valid_indices( const size_t* vertices, size_t num_vtx ) 00310 { 00311 // check valid size of list (at least fout, at most all nodes) 00312 CPPUNIT_ASSERT( num_vtx <= 27 ); 00313 CPPUNIT_ASSERT( num_vtx >= 4 ); 00314 // make sure vertex indices are valid (in [0,26]) 00315 size_t vertcopy[27]; 00316 std::copy( vertices, vertices + num_vtx, vertcopy ); 00317 std::sort( vertcopy, vertcopy + num_vtx ); 00318 CPPUNIT_ASSERT( vertcopy[num_vtx - 1] <= 26 ); // max value less than 27 00319 // make sure there are no duplicates in the list 00320 const size_t* iter = std::unique( vertcopy, vertcopy + num_vtx ); 00321 CPPUNIT_ASSERT( iter == vertcopy + num_vtx ); 00322 } 00323 00324 static void check_no_zeros( const MsqVector< 3 >* derivs, size_t num_vtx ) 00325 { 00326 for( unsigned i = 0; i < num_vtx; ++i ) 00327 { 00328 CPPUNIT_ASSERT( length( derivs[i] ) > 1e-6 ); 00329 } 00330 } 00331 00332 static void compare_coefficients( const double* coeffs, const size_t* indices, size_t num_coeff, 00333 const double* expected_coeffs, unsigned loc ) 00334 { 00335 // find the location in the returned list for each node 00336 size_t revidx[27]; 00337 double test_vals[27]; 00338 for( size_t i = 0; i < 27; ++i ) 00339 { 00340 revidx[i] = std::find( indices, indices + num_coeff, i ) - indices; 00341 test_vals[i] = ( revidx[i] == num_coeff ) ? 0.0 : coeffs[revidx[i]]; 00342 } 00343 00344 // compare expected and actual coefficient values 00345 ASSERT_VALUES_EQUAL( expected_coeffs[0], test_vals[0], loc ); 00346 ASSERT_VALUES_EQUAL( expected_coeffs[1], test_vals[1], loc ); 00347 ASSERT_VALUES_EQUAL( expected_coeffs[2], test_vals[2], loc ); 00348 ASSERT_VALUES_EQUAL( expected_coeffs[3], test_vals[3], loc ); 00349 ASSERT_VALUES_EQUAL( expected_coeffs[4], test_vals[4], loc ); 00350 ASSERT_VALUES_EQUAL( expected_coeffs[5], test_vals[5], loc ); 00351 ASSERT_VALUES_EQUAL( expected_coeffs[6], test_vals[6], loc ); 00352 ASSERT_VALUES_EQUAL( expected_coeffs[7], test_vals[7], loc ); 00353 ASSERT_VALUES_EQUAL( expected_coeffs[8], test_vals[8], loc ); 00354 ASSERT_VALUES_EQUAL( expected_coeffs[9], test_vals[9], loc ); 00355 ASSERT_VALUES_EQUAL( expected_coeffs[10], test_vals[10], loc ); 00356 ASSERT_VALUES_EQUAL( expected_coeffs[11], test_vals[11], loc ); 00357 ASSERT_VALUES_EQUAL( expected_coeffs[12], test_vals[12], loc ); 00358 ASSERT_VALUES_EQUAL( expected_coeffs[13], test_vals[13], loc ); 00359 ASSERT_VALUES_EQUAL( expected_coeffs[14], test_vals[14], loc ); 00360 ASSERT_VALUES_EQUAL( expected_coeffs[15], test_vals[15], loc ); 00361 ASSERT_VALUES_EQUAL( expected_coeffs[16], test_vals[16], loc ); 00362 ASSERT_VALUES_EQUAL( expected_coeffs[17], test_vals[17], loc ); 00363 ASSERT_VALUES_EQUAL( expected_coeffs[18], test_vals[18], loc ); 00364 ASSERT_VALUES_EQUAL( expected_coeffs[19], test_vals[19], loc ); 00365 ASSERT_VALUES_EQUAL( expected_coeffs[20], test_vals[20], loc ); 00366 ASSERT_VALUES_EQUAL( expected_coeffs[21], test_vals[21], loc ); 00367 ASSERT_VALUES_EQUAL( expected_coeffs[22], test_vals[22], loc ); 00368 ASSERT_VALUES_EQUAL( expected_coeffs[23], test_vals[23], loc ); 00369 ASSERT_VALUES_EQUAL( expected_coeffs[24], test_vals[24], loc ); 00370 ASSERT_VALUES_EQUAL( expected_coeffs[25], test_vals[25], loc ); 00371 ASSERT_VALUES_EQUAL( expected_coeffs[26], test_vals[26], loc ); 00372 } 00373 00374 static void compare_derivatives( const size_t* vertices, size_t num_vtx, const MsqVector< 3 >* actual, 00375 const MsqVector< 3 >* expected, unsigned loc ) 00376 { 00377 check_valid_indices( vertices, num_vtx ); 00378 check_no_zeros( actual, num_vtx ); 00379 00380 // Input has values in dxi & deta & dzeta only for nodes in 'vertices' 00381 // Convert to values for every possible node, with zero's for 00382 // nodes that are not present. 00383 CPPUNIT_ASSERT( num_vtx <= 9 ); 00384 MsqVector< 3 > expanded[27]; 00385 std::fill( expanded, expanded + 27, MsqVector< 3 >( 0.0 ) ); 00386 for( unsigned i = 0; i < num_vtx; ++i ) 00387 { 00388 CPPUNIT_ASSERT( vertices[i] <= 26 ); 00389 expanded[vertices[i]] = actual[i]; 00390 } 00391 00392 ASSERT_VALUES_EQUAL( expected[0], expanded[0], loc ); 00393 ASSERT_VALUES_EQUAL( expected[1], expanded[1], loc ); 00394 ASSERT_VALUES_EQUAL( expected[2], expanded[2], loc ); 00395 ASSERT_VALUES_EQUAL( expected[3], expanded[3], loc ); 00396 ASSERT_VALUES_EQUAL( expected[4], expanded[4], loc ); 00397 ASSERT_VALUES_EQUAL( expected[5], expanded[5], loc ); 00398 ASSERT_VALUES_EQUAL( expected[6], expanded[6], loc ); 00399 ASSERT_VALUES_EQUAL( expected[7], expanded[7], loc ); 00400 ASSERT_VALUES_EQUAL( expected[8], expanded[8], loc ); 00401 ASSERT_VALUES_EQUAL( expected[9], expanded[9], loc ); 00402 ASSERT_VALUES_EQUAL( expected[10], expanded[10], loc ); 00403 ASSERT_VALUES_EQUAL( expected[11], expanded[11], loc ); 00404 ASSERT_VALUES_EQUAL( expected[12], expanded[12], loc ); 00405 ASSERT_VALUES_EQUAL( expected[13], expanded[13], loc ); 00406 ASSERT_VALUES_EQUAL( expected[14], expanded[14], loc ); 00407 ASSERT_VALUES_EQUAL( expected[15], expanded[15], loc ); 00408 ASSERT_VALUES_EQUAL( expected[16], expanded[16], loc ); 00409 ASSERT_VALUES_EQUAL( expected[17], expanded[17], loc ); 00410 ASSERT_VALUES_EQUAL( expected[18], expanded[18], loc ); 00411 ASSERT_VALUES_EQUAL( expected[19], expanded[19], loc ); 00412 ASSERT_VALUES_EQUAL( expected[20], expanded[20], loc ); 00413 ASSERT_VALUES_EQUAL( expected[21], expanded[21], loc ); 00414 ASSERT_VALUES_EQUAL( expected[22], expanded[22], loc ); 00415 ASSERT_VALUES_EQUAL( expected[23], expanded[23], loc ); 00416 ASSERT_VALUES_EQUAL( expected[24], expanded[24], loc ); 00417 ASSERT_VALUES_EQUAL( expected[25], expanded[25], loc ); 00418 ASSERT_VALUES_EQUAL( expected[26], expanded[26], loc ); 00419 } 00420 00421 void HexLagrangeShapeTest::test_corner_coeff( int corner ) 00422 { 00423 MsqPrintError err( std::cout ); 00424 00425 double expected[27]; 00426 get_coeffs( XI_corner( corner ), expected ); 00427 00428 double coeff[100]; 00429 size_t num_coeff = 11, indices[100]; 00430 sf.coefficients( Sample( 0, corner ), allNodes, coeff, indices, num_coeff, err ); 00431 CPPUNIT_ASSERT( !err ); 00432 00433 compare_coefficients( coeff, indices, num_coeff, expected, corner ); 00434 } 00435 00436 void HexLagrangeShapeTest::test_edge_coeff( int edge ) 00437 { 00438 MsqPrintError err( std::cout ); 00439 00440 double expected[27]; 00441 get_coeffs( XI_edge( edge ), expected ); 00442 00443 double coeff[100]; 00444 size_t num_coeff = 11, indices[100]; 00445 sf.coefficients( Sample( 1, edge ), allNodes, coeff, indices, num_coeff, err ); 00446 CPPUNIT_ASSERT( !err ); 00447 00448 compare_coefficients( coeff, indices, num_coeff, expected, edge + 8 ); 00449 } 00450 00451 void HexLagrangeShapeTest::test_face_coeff( int face ) 00452 { 00453 MsqPrintError err( std::cout ); 00454 00455 double expected[27]; 00456 get_coeffs( XI_face( face ), expected ); 00457 00458 double coeff[100]; 00459 size_t num_coeff = 11, indices[100]; 00460 sf.coefficients( Sample( 2, face ), allNodes, coeff, indices, num_coeff, err ); 00461 CPPUNIT_ASSERT( !err ); 00462 00463 compare_coefficients( coeff, indices, num_coeff, expected, face + 20 ); 00464 } 00465 00466 void HexLagrangeShapeTest::test_mid_coeff() 00467 { 00468 MsqPrintError err( std::cout ); 00469 00470 double expected[27]; 00471 get_coeffs( XI_elem(), expected ); 00472 00473 double coeff[100]; 00474 size_t num_coeff = 11, indices[100]; 00475 sf.coefficients( Sample( 3, 0 ), allNodes, coeff, indices, num_coeff, err ); 00476 CPPUNIT_ASSERT( !err ); 00477 00478 compare_coefficients( coeff, indices, num_coeff, expected, 26 ); 00479 } 00480 00481 void HexLagrangeShapeTest::test_corner_derivs( int corner ) 00482 { 00483 MsqPrintError err( std::cout ); 00484 00485 MsqVector< 3 > expected[27]; 00486 get_derivs( XI_corner( corner ), expected ); 00487 00488 size_t vertices[100], num_vtx = 23; 00489 MsqVector< 3 > derivs[100]; 00490 sf.derivatives( Sample( 0, corner ), allNodes, vertices, derivs, num_vtx, err ); 00491 CPPUNIT_ASSERT( !err ); 00492 00493 compare_derivatives( vertices, num_vtx, derivs, expected, corner ); 00494 } 00495 00496 void HexLagrangeShapeTest::test_edge_derivs( int edge ) 00497 { 00498 MsqPrintError err( std::cout ); 00499 00500 MsqVector< 3 > expected[27]; 00501 get_derivs( XI_edge( edge ), expected ); 00502 00503 size_t vertices[100], num_vtx = 23; 00504 MsqVector< 3 > derivs[100]; 00505 sf.derivatives( Sample( 1, edge ), allNodes, vertices, derivs, num_vtx, err ); 00506 CPPUNIT_ASSERT( !err ); 00507 00508 compare_derivatives( vertices, num_vtx, derivs, expected, edge + 8 ); 00509 } 00510 00511 void HexLagrangeShapeTest::test_face_derivs( int face ) 00512 { 00513 MsqPrintError err( std::cout ); 00514 00515 MsqVector< 3 > expected[27]; 00516 get_derivs( XI_face( face ), expected ); 00517 00518 size_t vertices[100], num_vtx = 23; 00519 MsqVector< 3 > derivs[100]; 00520 sf.derivatives( Sample( 2, face ), allNodes, vertices, derivs, num_vtx, err ); 00521 CPPUNIT_ASSERT( !err ); 00522 00523 compare_derivatives( vertices, num_vtx, derivs, expected, face + 20 ); 00524 } 00525 00526 void HexLagrangeShapeTest::test_mid_derivs() 00527 { 00528 MsqPrintError err( std::cout ); 00529 00530 MsqVector< 3 > expected[27]; 00531 get_derivs( XI_elem(), expected ); 00532 00533 size_t vertices[100], num_vtx = 23; 00534 MsqVector< 3 > derivs[100]; 00535 sf.derivatives( Sample( 3, 0 ), allNodes, vertices, derivs, num_vtx, err ); 00536 CPPUNIT_ASSERT( !err ); 00537 00538 compare_derivatives( vertices, num_vtx, derivs, expected, 26 ); 00539 } 00540 00541 void HexLagrangeShapeTest::test_corners_coeff() 00542 { 00543 for( unsigned i = 0; i < 8; ++i ) 00544 test_corner_coeff( i ); 00545 } 00546 00547 void HexLagrangeShapeTest::test_edges_coeff() 00548 { 00549 for( unsigned i = 0; i < 12; ++i ) 00550 test_edge_coeff( i ); 00551 } 00552 00553 void HexLagrangeShapeTest::test_faces_coeff() 00554 { 00555 for( unsigned i = 0; i < 6; ++i ) 00556 test_face_coeff( i ); 00557 } 00558 00559 void HexLagrangeShapeTest::test_corners_derivs() 00560 { 00561 for( unsigned i = 0; i < 8; ++i ) 00562 test_corner_derivs( i ); 00563 } 00564 00565 void HexLagrangeShapeTest::test_edges_derivs() 00566 { 00567 for( unsigned i = 0; i < 12; ++i ) 00568 test_edge_derivs( i ); 00569 } 00570 00571 void HexLagrangeShapeTest::test_faces_derivs() 00572 { 00573 for( unsigned i = 0; i < 6; ++i ) 00574 test_face_derivs( i ); 00575 } 00576 00577 void HexLagrangeShapeTest::test_ideal_jacobian() 00578 { 00579 MsqError err; 00580 MsqMatrix< 3, 3 > J_act, J_exp( 1.0 ); 00581 sf.ideal( Sample( 3, 0 ), J_act, err ); 00582 ASSERT_NO_ERROR( err ); 00583 00584 // Matrices should be a rotation of each other. 00585 // First, calculate tentative rotation matrix 00586 MsqMatrix< 3, 3 > R = inverse( J_exp ) * J_act; 00587 // next check that it is a rotation 00588 CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( R ), 1e-6 ); // no scaling 00589 ASSERT_MATRICES_EQUAL( transpose( R ), inverse( R ), 1e-6 ); // orthogonal 00590 }