Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
quadratic_hex_map.hpp
Go to the documentation of this file.
1 #ifndef MOAB_QUADRATIC_HEX_HPP
2 #define MOAB_QUADRATIC_HEX_HPP
3 
4 #include "moab/Matrix3.hpp"
5 #include "moab/CartVect.hpp"
6 #include <sstream>
7 #include <iomanip>
8 #include <iostream>
9 
10 namespace moab
11 {
12 
13 namespace element_utility
14 {
15 
16  namespace
17  {
18 
19  double SH( const int i, const double xi )
20  {
21  switch( i )
22  {
23  case -1:
24  return ( xi * xi - xi ) / 2;
25  case 0:
26  return 1 - xi * xi;
27  case 1:
28  return ( xi * xi + xi ) / 2;
29  default:
30  return 0.;
31  }
32  }
33  double DSH( const int i, const double xi )
34  {
35  switch( i )
36  {
37  case -1:
38  return xi - 0.5;
39  case 0:
40  return -2 * xi;
41  case 1:
42  return xi + 0.5;
43  default:
44  return 0.;
45  }
46  }
47 
48  } // namespace
49 
50  template < typename _Matrix >
51  class Quadratic_hex_map
52  {
53  public:
54  typedef _Matrix Matrix;
55 
56  private:
57  typedef Quadratic_hex_map< Matrix > Self;
58 
59  public:
60  // Constructor
61  Quadratic_hex_map() {}
62  // Copy constructor
63  Quadratic_hex_map( const Self& f ) {}
64 
65  public:
66  // Natural coordinates
67  template < typename Moab, typename Entity_handle, typename Points, typename Point >
68  std::pair< bool, Point > operator()( const Moab& /* moab */,
69  const Entity_handle& /* h */,
70  const Points& v,
71  const Point& p,
72  const double tol = 1.e-6 ) const
73  {
74  Point result( 3, 0.0 );
75  bool point_found = solve_inverse( p, result, v, tol ) && is_contained( result, tol );
76  return std::make_pair( point_found, result );
77  }
78 
79  private:
80  // This is a hack to avoid a .cpp file and C++11
81  // reference_points(i,j) will be a 1 or -1;
82  // This should unroll..
83  inline double reference_points( const std::size_t& i, const std::size_t& j ) const
84  {
85  const double rpts[27][3] = { { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 }, // reference_points nodes: 0-7
86  { -1, 1, -1 }, // mid-edge nodes: 8-19
87  { -1, -1, 1 }, // center-face nodes 20-25 center node 26
88  { 1, -1, 1 }, //
89  { 1, 1, 1 }, { -1, 1, 1 }, // 4 ----- 19 ----- 7
90  { 0, -1, -1 }, // . | . |
91  { 1, 0, -1 }, // 16 25 18 |
92  { 0, 1, -1 }, // . | . |
93  { -1, 0, -1 }, // 5 ----- 17 ----- 6 |
94  { -1, -1, 0 }, // | 12 | 23 15
95  { 1, -1, 0 }, // | | |
96  { 1, 1, 0 }, // | 20 | 26 | 22 |
97  { -1, 1, 0 }, // | | |
98  { 0, -1, 1 }, // 13 21 | 14 |
99  { 1, 0, 1 }, // | 0 ----- 11 ----- 3
100  { 0, 1, 1 }, // | . | .
101  { -1, 0, 1 }, // | 8 24 | 10
102  { 0, -1, 0 }, // | . | .
103  { 1, 0, 0 }, // 1 ----- 9 ----- 2
104  { 0, 1, 0 }, //
105  { -1, 0, 0 }, { 0, 0, -1 }, { 0, 0, 1 }, { 0, 0, 0 } };
106  return rpts[i][j];
107  }
108 
109  template < typename Point >
110  bool is_contained( const Point& p, const double tol ) const
111  {
112  // just look at the box+tol here
113  return ( p[0] >= -1. - tol ) && ( p[0] <= 1. + tol ) && ( p[1] >= -1. - tol ) && ( p[1] <= 1. + tol ) &&
114  ( p[2] >= -1. - tol ) && ( p[2] <= 1. + tol );
115  }
116 
117  template < typename Point, typename Points >
118  bool solve_inverse( const Point& x, Point& xi, const Points& points, const double tol = 1.e-6 ) const
119  {
120  const double error_tol_sqr = tol * tol;
121  Point delta( 3, 0.0 );
122  xi = delta;
123  evaluate( xi, points, delta );
124  vec_subtract( delta, x );
125  std::size_t num_iterations = 0;
126 #ifdef QUADRATIC_HEX_DEBUG
127  std::stringstream ss;
128  ss << "Point: ";
129  ss << x[0] << ", " << x[1] << ", " << x[2] << std::endl;
130  ss << "Hex: ";
131  for( int i = 0; i < 8; ++i )
132  {
133  ss << points[i][0] << ", " << points[i][1] << ", " << points[i][2] << std::endl;
134  }
135  ss << std::endl;
136 #endif
137  while( normsq( delta ) > error_tol_sqr )
138  {
139 #ifdef QUADRATIC_HEX_DEBUG
140  ss << "Iter #: " << num_iterations << " Err: " << sqrt( normsq( delta ) ) << " Iterate: ";
141  ss << xi[0] << ", " << xi[1] << ", " << xi[2] << std::endl;
142 #endif
143  if( ++num_iterations >= 5 )
144  {
145  return false;
146  }
147  Matrix J;
148  jacobian( xi, points, J );
149  double det = moab::Matrix::determinant3( J );
150  if( fabs( det ) < 1.e-10 )
151  {
152 #ifdef QUADRATIC_HEX_DEBUG
153  std::cerr << ss.str();
154 #endif
155 #ifndef QUADRATIC_HEX_DEBUG
156  std::cerr << x[0] << ", " << x[1] << ", " << x[2] << std::endl;
157 #endif
158  std::cerr << "inverse solve failure: det: " << det << std::endl;
159  exit( -1 );
160  }
161  vec_subtract( xi, moab::Matrix::inverse( J, 1.0 / det ) * delta );
162  evaluate( xi, points, delta );
163  vec_subtract( delta, x );
164  }
165  return true;
166  }
167 
168  template < typename Point, typename Points >
169  Point& evaluate( const Point& p, const Points& points, Point& f ) const
170  {
171  typedef typename Points::value_type Vector;
172  Vector result;
173  for( int i = 0; i < 3; ++i )
174  {
175  result[i] = 0;
176  }
177  for( unsigned i = 0; i < 27; ++i )
178  {
179  const double sh = SH( reference_points( i, 0 ), p[0] ) * SH( reference_points( i, 1 ), p[1] ) *
180  SH( reference_points( i, 2 ), p[2] );
181  result += sh * points[i];
182  }
183  for( int i = 0; i < 3; ++i )
184  {
185  f[i] = result[i];
186  }
187  return f;
188  }
189  template < typename Point, typename Field >
190  double evaluate_scalar_field( const Point& p, const Field& field ) const
191  {
192  double x = 0.0;
193  for( int i = 0; i < 27; i++ )
194  {
195  const double sh = SH( reference_points( i, 0 ), p[0] ) * SH( reference_points( i, 1 ), p[1] ) *
196  SH( reference_points( i, 2 ), p[2] );
197  x += sh * field[i];
198  }
199  return x;
200  }
201  template < typename Field, typename Points >
202  double integrate_scalar_field( const Points& p, const Field& field_values ) const
203  {
204  // TODO: gaussian integration , probably 2x2x2
205  return 0.;
206  }
207 
208  template < typename Point, typename Points >
209  Matrix& jacobian( const Point& p, const Points& /* points */, Matrix& J ) const
210  {
211  J = Matrix( 0.0 );
212  for( int i = 0; i < 27; i++ )
213  {
214  const double sh[3] = { SH( reference_points( i, 0 ), p[0] ), SH( reference_points( i, 1 ), p[1] ),
215  SH( reference_points( i, 2 ), p[2] ) };
216  const double dsh[3] = { DSH( reference_points( i, 0 ), p[0] ), DSH( reference_points( i, 1 ), p[1] ),
217  DSH( reference_points( i, 2 ), p[2] ) };
218  for( int j = 0; j < 3; j++ )
219  {
220  // dxj/dr first column
221  J( j, 0 ) += dsh[0] * sh[1] * sh[2] * reference_points( i, j );
222  J( j, 1 ) += sh[0] * dsh[1] * sh[2] * reference_points( i, j ); // dxj/ds
223  J( j, 2 ) += sh[0] * sh[1] * dsh[2] * reference_points( i, j ); // dxj/dt
224  }
225  }
226  return J;
227  }
228 
229  private:
230  }; // Class Quadratic_hex_map
231 
232 } // namespace element_utility
233 
234 } // namespace moab
235 #endif // MOAB_QUADRATIC_HEX_nPP