MOAB: Mesh Oriented datABase  (version 5.2.1)
linear_tet_map.hpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines