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 (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 }