linear_tet_map.hpp

Go to the documentation of this file.

#ifndef MOAB_LINEAR_TET_HPP
#define MOAB_LINEAR_TET_HPP
 
#include "moab/Matrix3.hpp"
 
namespace moab
{
namespace element_utility
{
 
 template < typename Entity_handle, typename Matrix >
 class Linear_tet_map
 {
 private:
 typedef Linear_tet_map< Entity_handle, Matrix > Self;
 
 public:
 // Constructor
 Linear_tet_map() : Tinv(), eh() {}
 // Copy constructor
 Linear_tet_map( const Self& f ) : Tinv( f.Tinv ), eh( f.eh ) {}
 // Natural coordinates
 template < typename Moab, typename Points, typename Point >
 std::pair< bool, Point > operator()( const Moab& moab,
 const Entity_handle _eh,
 const Points& v,
 const Point& p,
 const double tol = 1e-6 )
 {
 // Remove the warning about unused parameter
 if( NULL != &moab )
 {
 }
 
 set_tet( _eh, v );
 // TODO: Make sure this is correct
 Point result = Tinv * p;
 return std::make_pair( is_contained( result, tol ), result );
 }
 
 private:
 template < typename Point >
 bool is_contained( const Point& result, const double tol = 1e-6 )
 {
 double sum = 0.0;
 for( std::size_t i = 0; i < 3; ++i )
 {
 sum += result[i];
 if( result[i] < -tol )
 {
 return false;
 }
 }
 return sum < 1.0 + tol;
 }
 template < typename Point, typename Field >
 double evaluate_scalar_field( const Point& p, const Field& field_values ) const
 {
 double f0 = field_values[0];
 double f = f0;
 for( std::size_t i = 1; i < 5; ++i )
 {
 f += ( field_values[i] - f0 ) * p[i - 1];
 }
 return f;
 }
 template < typename Points, typename Field >
 double integrate_scalar_field( const Points& v, const Field& field_values ) const
 {
 double I( 0.0 );
 for( unsigned int i = 0; i < 4; ++i )
 {
 I += field_values[i];
 }
 double det =
 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],
 v[3][1] - v[0][1], v[1][2] - v[0][2], v[2][2] - v[0][2], v[3][2] - v[0][2] )
 .determinant();
 I *= det / 24.0;
 return I;
 }
 
 template < typename Points >
 void set_tet( const Entity_handle _eh, const Points& v )
 {
 if( eh != _eh )
 {
 eh = _eh;
 Tinv = moab::Matrix::inverse( 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], v[3][1] - v[0][1],
 v[1][2] - v[0][2], v[2][2] - v[0][2], v[3][2] - v[0][2] ) );
 }
 }
 
 private:
 Matrix Tinv;
 Entity_handle eh;
 /* We don't need this!
 static const double reference_points[ 4][ 3] = { {0,0,0},
 {1,0,0},
 {0,1,0},
 {0,0,1} };*/
 
 }; // Class Linear_tet_map
 
} // namespace element_utility
 
} // namespace moab
#endif // MOAB_LINEAR_TET_HPP