Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
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 >
13  {
14  private:
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