MOAB: Mesh Oriented datABase  (version 5.4.0)
test_Matrix3.cpp
Go to the documentation of this file.
00001 #include <vector>
00002 #include <cassert>
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,
00025                           const moab::CartVect& B,
00026                           double eps,
00027                           const char* sA,
00028                           const char* sB,
00029                           int line,
00030                           const char* file )
00031 {
00032     check_equal( A.length(), B.length(), eps, sA, sB, line, file );
00033 
00034     double angle = find_angle( A, B );
00035 
00036     if( ( fabs( A[0] - B[0] ) <= eps || fabs( A[0] + B[0] ) <= eps ) &&
00037         ( fabs( A[1] - B[1] ) <= eps || fabs( A[1] + B[1] ) <= eps ) &&
00038         ( fabs( A[2] - B[2] ) <= eps || fabs( A[2] + B[2] ) <= eps ) &&
00039         ( angle <= eps || fabs( angle - 180.0 ) <= eps ) )
00040         return;
00041 
00042     std::cout << "Equality Test Failed: " << sA << " == " << sB << std::endl;
00043     std::cout << "  at line " << line << " of '" << file << "'" << std::endl;
00044 
00045     std::cout << "  Expected: ";
00046     std::cout << A << std::endl;
00047 
00048     std::cout << "  Actual:   ";
00049     std::cout << B << std::endl;
00050 
00051     std::cout << "Angle between vectors := " << angle << std::endl;
00052 
00053     flag_error();
00054 }
00055 
00056 void test_EigenDecomp();
00057 
00058 int main()
00059 {
00060 
00061     int result = 0;
00062 
00063     result += RUN_TEST( test_EigenDecomp );
00064 
00065     return result;
00066 }
00067 
00068 // test to ensure the Eigenvalues/vectors are calculated correctly and returned properly
00069 // from the Matrix3 class for a simple case
00070 void test_EigenDecomp()
00071 {
00072     // Create a matrix
00073     moab::Matrix3 mat;
00074 
00075     mat( 0 ) = 2;
00076     mat( 1 ) = -1;
00077     mat( 2 ) = 0;
00078     mat( 3 ) = -1;
00079     mat( 4 ) = 2;
00080     mat( 5 ) = -1;
00081     mat( 6 ) = 0;
00082     mat( 7 ) = -1;
00083     mat( 8 ) = 2;
00084 
00085     // now do the Eigen Decomposition of this Matrix
00086 
00087     CartVect lamda;
00088     Matrix3 vectors;
00089     moab::ErrorCode rval = mat.eigen_decomposition( lamda, vectors );CHECK_ERR( rval );
00090     for( int i = 0; i < 3; ++i )
00091         vectors.col( i ).normalize();
00092 
00093     // Hardcoded check values for the results
00094     double lamda_check[3];
00095     lamda_check[0] = 0.585786;
00096     lamda_check[1] = 2.0;
00097     lamda_check[2] = 3.41421;
00098 
00099     moab::CartVect vec0_check( 0.5, 0.707107, 0.5 );
00100     moab::CartVect vec1_check( 0.707107, 3.37748e-17, -0.707107 );
00101     moab::CartVect vec2_check( 0.5, -0.707107, 0.5 );
00102 
00103     // now verfy that the returns Eigenvalues and Eigenvectors are correct (within some tolerance)
00104     double tol = 1e-04;
00105 
00106     // use a deprecated constructor
00107     Matrix3 mat3( CartVect( 2, -1, 0 ), CartVect( -1, 2, -1 ), CartVect( 0, -1, 2 ), true );
00108     CHECK_REAL_EQUAL( mat( 1 ), mat3( 1 ), tol );
00109 
00110     vec0_check.normalize();
00111     vec1_check.normalize();
00112     vec2_check.normalize();
00113 
00114     // check that the correct Eigenvalues are returned correctly (in order)
00115     CHECK_REAL_EQUAL( lamda[0], lamda_check[0], tol );
00116     CHECK_REAL_EQUAL( lamda[1], lamda_check[1], tol );
00117     CHECK_REAL_EQUAL( lamda[2], lamda_check[2], tol );
00118 
00119     // check the Eigenvector values (order should correspond to the Eigenvalues)
00120     // first vector
00121     CHECK_EIGVECREAL_EQUAL( vectors.col( 0 ), vec0_check, tol );
00122 
00123     // sceond vector
00124     CHECK_EIGVECREAL_EQUAL( vectors.col( 1 ), vec1_check, tol );
00125 
00126     // third vector
00127     CHECK_EIGVECREAL_EQUAL( vectors.col( 2 ), vec2_check, tol );
00128 
00129     // another check to ensure the result is valid (AM-kM = 0)
00130     for( unsigned i = 0; i < 3; ++i )
00131     {
00132         moab::CartVect v = moab::Matrix::matrix_vector( mat, vectors.col( i ) ) - lamda[i] * vectors.col( i );
00133         CHECK_REAL_EQUAL( v.length(), 0, tol );
00134     }
00135 
00136     // for a real, symmetric matrix the Eigenvectors should be orthogonal
00137     CHECK_REAL_EQUAL( vectors.col( 0 ) % vectors.col( 1 ), 0, tol );
00138     CHECK_REAL_EQUAL( vectors.col( 0 ) % vectors.col( 2 ), 0, tol );
00139 
00140     return;
00141 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines