MOAB: Mesh Oriented datABase  (version 5.3.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, 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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines