Loading [MathJax]/extensions/tex2jax.js
Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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