MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 #include <vector> 00002 #include <assert.h> 00003 #include "moab/Matrix3.hpp" 00004 #include "moab/Core.hpp" 00005 #include "moab/CartVect.hpp" 00006 #include "TestUtil.hpp" 00007 00008 using namespace moab; 00009 00010 #define ACOS( x ) acos( std::min( std::max( x, -1.0 ), 1.0 ) ) 00011 00012 double find_angle( const moab::CartVect& A, const moab::CartVect& B ) 00013 { 00014 const double lA = A.length(), lB = B.length(); 00015 assert( lA > 0.0 ); 00016 assert( lB > 0.0 ); 00017 const double dPI = 3.14159265; 00018 const double dot = ( A[0] * B[0] + A[1] * B[1] + A[2] * B[2] ); 00019 return ACOS( dot / ( lA * lB ) ) * 180.0 / dPI; 00020 } 00021 00022 #define CHECK_EIGVECREAL_EQUAL( EXP, ACT, EPS ) \ 00023 check_equal_eigvect( ( EXP ), ( ACT ), ( EPS ), #EXP, #ACT, __LINE__, __FILE__ ) 00024 void check_equal_eigvect( const moab::CartVect& A, const moab::CartVect& B, double eps, const char* sA, const char* sB, 00025 int line, const char* file ) 00026 { 00027 check_equal( A.length(), B.length(), eps, sA, sB, line, file ); 00028 00029 double angle = find_angle( A, B ); 00030 00031 if( ( fabs( A[0] - B[0] ) <= eps || fabs( A[0] + B[0] ) <= eps ) && 00032 ( fabs( A[1] - B[1] ) <= eps || fabs( A[1] + B[1] ) <= eps ) && 00033 ( fabs( A[2] - B[2] ) <= eps || fabs( A[2] + B[2] ) <= eps ) && 00034 ( angle <= eps || fabs( angle - 180.0 ) <= eps ) ) 00035 return; 00036 00037 std::cout << "Equality Test Failed: " << sA << " == " << sB << std::endl; 00038 std::cout << " at line " << line << " of '" << file << "'" << std::endl; 00039 00040 std::cout << " Expected: "; 00041 std::cout << A << std::endl; 00042 00043 std::cout << " Actual: "; 00044 std::cout << B << std::endl; 00045 00046 std::cout << "Angle between vectors := " << angle << std::endl; 00047 00048 flag_error(); 00049 } 00050 00051 void test_EigenDecomp(); 00052 00053 int main() 00054 { 00055 00056 int result = 0; 00057 00058 result += RUN_TEST( test_EigenDecomp ); 00059 00060 return result; 00061 } 00062 00063 // test to ensure the Eigenvalues/vectors are calculated correctly and returned properly 00064 // from the Matrix3 class for a simple case 00065 void test_EigenDecomp() 00066 { 00067 // Create a matrix 00068 moab::Matrix3 mat; 00069 00070 mat( 0 ) = 2; 00071 mat( 1 ) = -1; 00072 mat( 2 ) = 0; 00073 mat( 3 ) = -1; 00074 mat( 4 ) = 2; 00075 mat( 5 ) = -1; 00076 mat( 6 ) = 0; 00077 mat( 7 ) = -1; 00078 mat( 8 ) = 2; 00079 00080 // now do the Eigen Decomposition of this Matrix 00081 00082 CartVect lamda; 00083 Matrix3 vectors; 00084 moab::ErrorCode rval = mat.eigen_decomposition( lamda, vectors );CHECK_ERR( rval ); 00085 for( int i = 0; i < 3; ++i ) 00086 vectors.col( i ).normalize(); 00087 00088 // Hardcoded check values for the results 00089 double lamda_check[3]; 00090 lamda_check[0] = 0.585786; 00091 lamda_check[1] = 2.0; 00092 lamda_check[2] = 3.41421; 00093 00094 moab::CartVect vec0_check( 0.5, 0.707107, 0.5 ); 00095 moab::CartVect vec1_check( 0.707107, 3.37748e-17, -0.707107 ); 00096 moab::CartVect vec2_check( 0.5, -0.707107, 0.5 ); 00097 00098 // now verfy that the returns Eigenvalues and Eigenvectors are correct (within some tolerance) 00099 double tol = 1e-04; 00100 00101 // use a deprecated constructor 00102 Matrix3 mat3( CartVect( 2, -1, 0 ), CartVect( -1, 2, -1 ), CartVect( 0, -1, 2 ), true ); 00103 CHECK_REAL_EQUAL( mat( 1 ), mat3( 1 ), tol ); 00104 00105 vec0_check.normalize(); 00106 vec1_check.normalize(); 00107 vec2_check.normalize(); 00108 00109 // check that the correct Eigenvalues are returned correctly (in order) 00110 CHECK_REAL_EQUAL( lamda[0], lamda_check[0], tol ); 00111 CHECK_REAL_EQUAL( lamda[1], lamda_check[1], tol ); 00112 CHECK_REAL_EQUAL( lamda[2], lamda_check[2], tol ); 00113 00114 // check the Eigenvector values (order should correspond to the Eigenvalues) 00115 // first vector 00116 CHECK_EIGVECREAL_EQUAL( vectors.col( 0 ), vec0_check, tol ); 00117 00118 // sceond vector 00119 CHECK_EIGVECREAL_EQUAL( vectors.col( 1 ), vec1_check, tol ); 00120 00121 // third vector 00122 CHECK_EIGVECREAL_EQUAL( vectors.col( 2 ), vec2_check, tol ); 00123 00124 // another check to ensure the result is valid (AM-kM = 0) 00125 for( unsigned i = 0; i < 3; ++i ) 00126 { 00127 moab::CartVect v = moab::Matrix::matrix_vector( mat, vectors.col( i ) ) - lamda[i] * vectors.col( i ); 00128 CHECK_REAL_EQUAL( v.length(), 0, tol ); 00129 } 00130 00131 // for a real, symmetric matrix the Eigenvectors should be orthogonal 00132 CHECK_REAL_EQUAL( vectors.col( 0 ) % vectors.col( 1 ), 0, tol ); 00133 CHECK_REAL_EQUAL( vectors.col( 0 ) % vectors.col( 2 ), 0, tol ); 00134 00135 return; 00136 }