Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
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,
00025                                              const Entity_handle _eh,
00026                                              const Points& v,
00027                                              const Point& p,
00028                                              const double tol = 1e-6 )
00029         {
00030             // Remove the warning about unused parameter
00031             if( NULL != &moab )
00032             {
00033             }
00034 
00035             set_tet( _eh, v );
00036             // TODO: Make sure this is correct
00037             Point result = Tinv * p;
00038             return std::make_pair( is_contained( result, tol ), result );
00039         }
00040 
00041       private:
00042         template < typename Point >
00043         bool is_contained( const Point& result, const double tol = 1e-6 )
00044         {
00045             double sum = 0.0;
00046             for( std::size_t i = 0; i < 3; ++i )
00047             {
00048                 sum += result[i];
00049                 if( result[i] < -tol )
00050                 {
00051                     return false;
00052                 }
00053             }
00054             return sum < 1.0 + tol;
00055         }
00056         template < typename Point, typename Field >
00057         double evaluate_scalar_field( const Point& p, const Field& field_values ) const
00058         {
00059             double f0 = field_values[0];
00060             double f  = f0;
00061             for( std::size_t i = 1; i < 5; ++i )
00062             {
00063                 f += ( field_values[i] - f0 ) * p[i - 1];
00064             }
00065             return f;
00066         }
00067         template < typename Points, typename Field >
00068         double integrate_scalar_field( const Points& v, const Field& field_values ) const
00069         {
00070             double I( 0.0 );
00071             for( unsigned int i = 0; i < 4; ++i )
00072             {
00073                 I += field_values[i];
00074             }
00075             double det =
00076                 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],
00077                         v[3][1] - v[0][1], v[1][2] - v[0][2], v[2][2] - v[0][2], v[3][2] - v[0][2] )
00078                     .determinant();
00079             I *= det / 24.0;
00080             return I;
00081         }
00082 
00083         template < typename Points >
00084         void set_tet( const Entity_handle _eh, const Points& v )
00085         {
00086             if( eh != _eh )
00087             {
00088                 eh   = _eh;
00089                 Tinv = moab::Matrix::inverse( Matrix( v[1][0] - v[0][0], v[2][0] - v[0][0], v[3][0] - v[0][0],
00090                                                       v[1][1] - v[0][1], v[2][1] - v[0][1], v[3][1] - v[0][1],
00091                                                       v[1][2] - v[0][2], v[2][2] - v[0][2], v[3][2] - v[0][2] ) );
00092             }
00093         }
00094 
00095       private:
00096         Matrix Tinv;
00097         Entity_handle eh;
00098         /* We don't need this!
00099         static const double reference_points[ 4][ 3] = { {0,0,0},
00100                                      {1,0,0},
00101                                      {0,1,0},
00102                                      {0,0,1} };*/
00103 
00104     };  // Class Linear_tet_map
00105 
00106 }  // namespace element_utility
00107 
00108 }  // namespace moab
00109 #endif  // MOAB_LINEAR_TET_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines