MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 #ifndef MOAB_LINEAR_TET_HPP 00002 #define MOAB_LINEAR_TET_HPP 00003 00004 #include "moab/Matrix3.hpp" 00005 00006 namespace moab 00007 { 00008 namespace element_utility 00009 { 00010 00011 template < typename Entity_handle, typename Matrix > 00012 class Linear_tet_map 00013 { 00014 private: 00015 typedef Linear_tet_map< Entity_handle, Matrix > Self; 00016 00017 public: 00018 // Constructor 00019 Linear_tet_map() : Tinv(), eh() {} 00020 // Copy constructor 00021 Linear_tet_map( const Self& f ) : Tinv( f.Tinv ), eh( f.eh ) {} 00022 // Natural coordinates 00023 template < typename Moab, typename Points, typename Point > 00024 std::pair< bool, Point > operator()( const Moab& moab, const Entity_handle _eh, const Points& v, const Point& p, 00025 const double tol = 1e-6 ) 00026 { 00027 // Remove the warning about unused parameter 00028 if( NULL != &moab ) {} 00029 00030 set_tet( _eh, v ); 00031 // TODO: Make sure this is correct 00032 Point result = Tinv * p; 00033 return std::make_pair( is_contained( result, tol ), result ); 00034 } 00035 00036 private: 00037 template < typename Point > 00038 bool is_contained( const Point& result, const double tol = 1e-6 ) 00039 { 00040 double sum = 0.0; 00041 for( std::size_t i = 0; i < 3; ++i ) 00042 { 00043 sum += result[i]; 00044 if( result[i] < -tol ) { return false; } 00045 } 00046 return sum < 1.0 + tol; 00047 } 00048 template < typename Point, typename Field > 00049 double evaluate_scalar_field( const Point& p, const Field& field_values ) const 00050 { 00051 double f0 = field_values[0]; 00052 double f = f0; 00053 for( std::size_t i = 1; i < 5; ++i ) 00054 { 00055 f += ( field_values[i] - f0 ) * p[i - 1]; 00056 } 00057 return f; 00058 } 00059 template < typename Points, typename Field > 00060 double integrate_scalar_field( const Points& v, const Field& field_values ) const 00061 { 00062 double I( 0.0 ); 00063 for( unsigned int i = 0; i < 4; ++i ) 00064 { 00065 I += field_values[i]; 00066 } 00067 double det = 00068 Matrix( v[1][0] - v[0][0], v[2][0] - v[0][0], v[3][0] - v[0][0], v[1][1] - v[0][1], v[2][1] - v[0][1], 00069 v[3][1] - v[0][1], v[1][2] - v[0][2], v[2][2] - v[0][2], v[3][2] - v[0][2] ) 00070 .determinant(); 00071 I *= det / 24.0; 00072 return I; 00073 } 00074 00075 template < typename Points > 00076 void set_tet( const Entity_handle _eh, const Points& v ) 00077 { 00078 if( eh != _eh ) 00079 { 00080 eh = _eh; 00081 Tinv = moab::Matrix::inverse( Matrix( v[1][0] - v[0][0], v[2][0] - v[0][0], v[3][0] - v[0][0], 00082 v[1][1] - v[0][1], v[2][1] - v[0][1], v[3][1] - v[0][1], 00083 v[1][2] - v[0][2], v[2][2] - v[0][2], v[3][2] - v[0][2] ) ); 00084 } 00085 } 00086 00087 private: 00088 Matrix Tinv; 00089 Entity_handle eh; 00090 /* We don't need this! 00091 static const double reference_points[ 4][ 3] = { {0,0,0}, 00092 {1,0,0}, 00093 {0,1,0}, 00094 {0,0,1} };*/ 00095 00096 }; // Class Linear_tet_map 00097 00098 } // namespace element_utility 00099 00100 } // namespace moab 00101 #endif // MOAB_LINEAR_TET_HPP