linear_tet_map.hpp

Go to the documentation of this file.

    1 #ifndef MOAB_LINEAR_TET_HPP
    2 #define MOAB_LINEAR_TET_HPP
    3  
    4 #include "moab/Matrix3.hpp"
    5  
    6 namespace moab
    7 {
    8 namespace element_utility
    9 {
   10  
   11     template < typename Entity_handle, typename Matrix >
   12     class Linear_tet_map
   13     {
   14       private:
   15         typedef Linear_tet_map< Entity_handle, Matrix > Self;
   16  
   17       public:
   18         // Constructor
   19         Linear_tet_map() : Tinv(), eh() {}
   20         // Copy constructor
   21         Linear_tet_map( const Self& f ) : Tinv( f.Tinv ), eh( f.eh ) {}
   22         // Natural coordinates
   23         template < typename Moab, typename Points, typename Point >
   24         std::pair< bool, Point > operator()( const Moab& moab,
   25                                              const Entity_handle _eh,
   26                                              const Points& v,
   27                                              const Point& p,
   28                                              const double tol = 1e-6 )
   29         {
   30             // Remove the warning about unused parameter
   31             if( NULL != &moab )
   32             {
   33             }
   34  
   35             set_tet( _eh, v );
   36             // TODO: Make sure this is correct
   37             Point result = Tinv * p;
   38             return std::make_pair( is_contained( result, tol ), result );
   39         }
   40  
   41       private:
   42         template < typename Point >
   43         bool is_contained( const Point& result, const double tol = 1e-6 )
   44         {
   45             double sum = 0.0;
   46             for( std::size_t i = 0; i < 3; ++i )
   47             {
   48                 sum += result[i];
   49                 if( result[i] < -tol )
   50                 {
   51                     return false;
   52                 }
   53             }
   54             return sum < 1.0 + tol;
   55         }
   56         template < typename Point, typename Field >
   57         double evaluate_scalar_field( const Point& p, const Field& field_values ) const
   58         {
   59             double f0 = field_values[0];
   60             double f  = f0;
   61             for( std::size_t i = 1; i < 5; ++i )
   62             {
   63                 f += ( field_values[i] - f0 ) * p[i - 1];
   64             }
   65             return f;
   66         }
   67         template < typename Points, typename Field >
   68         double integrate_scalar_field( const Points& v, const Field& field_values ) const
   69         {
   70             double I( 0.0 );
   71             for( unsigned int i = 0; i < 4; ++i )
   72             {
   73                 I += field_values[i];
   74             }
   75             double det =
   76                 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],
   77                         v[3][1] - v[0][1], v[1][2] - v[0][2], v[2][2] - v[0][2], v[3][2] - v[0][2] )
   78                     .determinant();
   79             I *= det / 24.0;
   80             return I;
   81         }
   82  
   83         template < typename Points >
   84         void set_tet( const Entity_handle _eh, const Points& v )
   85         {
   86             if( eh != _eh )
   87             {
   88                 eh   = _eh;
   89                 Tinv = moab::Matrix::inverse( Matrix( v[1][0] - v[0][0], v[2][0] - v[0][0], v[3][0] - v[0][0],
   90                                                       v[1][1] - v[0][1], v[2][1] - v[0][1], v[3][1] - v[0][1],
   91                                                       v[1][2] - v[0][2], v[2][2] - v[0][2], v[3][2] - v[0][2] ) );
   92             }
   93         }
   94  
   95       private:
   96         Matrix Tinv;
   97         Entity_handle eh;
   98         /* We don't need this!
   99         static const double reference_points[ 4][ 3] = { {0,0,0},
  100                                      {1,0,0},
  101                                      {0,1,0},
  102                                      {0,0,1} };*/
  103  
  104     };  // Class Linear_tet_map
  105  
  106 }  // namespace element_utility
  107  
  108 }  // namespace moab
  109 #endif  // MOAB_LINEAR_TET_HPP