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