|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
00001 #ifndef MOAB_QUADRATIC_HEX_HPP 00002 #define MOAB_QUADRATIC_HEX_HPP 00003 00004 #include "moab/Matrix3.hpp" 00005 #include "moab/CartVect.hpp" 00006 #include <sstream> 00007 #include <iomanip> 00008 #include <iostream> 00009 00010 namespace moab 00011 { 00012 00013 namespace element_utility 00014 { 00015 00016 namespace 00017 { 00018 00019 double SH( const int i, const double xi ) 00020 { 00021 switch( i ) 00022 { 00023 case -1: 00024 return ( xi * xi - xi ) / 2; 00025 case 0: 00026 return 1 - xi * xi; 00027 case 1: 00028 return ( xi * xi + xi ) / 2; 00029 default: 00030 return 0.; 00031 } 00032 } 00033 double DSH( const int i, const double xi ) 00034 { 00035 switch( i ) 00036 { 00037 case -1: 00038 return xi - 0.5; 00039 case 0: 00040 return -2 * xi; 00041 case 1: 00042 return xi + 0.5; 00043 default: 00044 return 0.; 00045 } 00046 } 00047 00048 } // namespace 00049 00050 template < typename _Matrix > 00051 class Quadratic_hex_map 00052 { 00053 public: 00054 typedef _Matrix Matrix; 00055 00056 private: 00057 typedef Quadratic_hex_map< Matrix > Self; 00058 00059 public: 00060 // Constructor 00061 Quadratic_hex_map() {} 00062 // Copy constructor 00063 Quadratic_hex_map( const Self& f ) {} 00064 00065 public: 00066 // Natural coordinates 00067 template < typename Moab, typename Entity_handle, typename Points, typename Point > 00068 std::pair< bool, Point > operator()( const Moab& /* moab */, 00069 const Entity_handle& /* h */, 00070 const Points& v, 00071 const Point& p, 00072 const double tol = 1.e-6 ) const 00073 { 00074 Point result( 3, 0.0 ); 00075 bool point_found = solve_inverse( p, result, v, tol ) && is_contained( result, tol ); 00076 return std::make_pair( point_found, result ); 00077 } 00078 00079 private: 00080 // This is a hack to avoid a .cpp file and C++11 00081 // reference_points(i,j) will be a 1 or -1; 00082 // This should unroll.. 00083 inline double reference_points( const std::size_t& i, const std::size_t& j ) const 00084 { 00085 const double rpts[27][3] = { { -1, -1, -1 }, { 1, -1, -1 }, { 1, 1, -1 }, // reference_points nodes: 0-7 00086 { -1, 1, -1 }, // mid-edge nodes: 8-19 00087 { -1, -1, 1 }, // center-face nodes 20-25 center node 26 00088 { 1, -1, 1 }, // 00089 { 1, 1, 1 }, { -1, 1, 1 }, // 4 ----- 19 ----- 7 00090 { 0, -1, -1 }, // . | . | 00091 { 1, 0, -1 }, // 16 25 18 | 00092 { 0, 1, -1 }, // . | . | 00093 { -1, 0, -1 }, // 5 ----- 17 ----- 6 | 00094 { -1, -1, 0 }, // | 12 | 23 15 00095 { 1, -1, 0 }, // | | | 00096 { 1, 1, 0 }, // | 20 | 26 | 22 | 00097 { -1, 1, 0 }, // | | | 00098 { 0, -1, 1 }, // 13 21 | 14 | 00099 { 1, 0, 1 }, // | 0 ----- 11 ----- 3 00100 { 0, 1, 1 }, // | . | . 00101 { -1, 0, 1 }, // | 8 24 | 10 00102 { 0, -1, 0 }, // | . | . 00103 { 1, 0, 0 }, // 1 ----- 9 ----- 2 00104 { 0, 1, 0 }, // 00105 { -1, 0, 0 }, { 0, 0, -1 }, { 0, 0, 1 }, { 0, 0, 0 } }; 00106 return rpts[i][j]; 00107 } 00108 00109 template < typename Point > 00110 bool is_contained( const Point& p, const double tol ) const 00111 { 00112 // just look at the box+tol here 00113 return ( p[0] >= -1. - tol ) && ( p[0] <= 1. + tol ) && ( p[1] >= -1. - tol ) && ( p[1] <= 1. + tol ) && 00114 ( p[2] >= -1. - tol ) && ( p[2] <= 1. + tol ); 00115 } 00116 00117 template < typename Point, typename Points > 00118 bool solve_inverse( const Point& x, Point& xi, const Points& points, const double tol = 1.e-6 ) const 00119 { 00120 const double error_tol_sqr = tol * tol; 00121 Point delta( 3, 0.0 ); 00122 xi = delta; 00123 evaluate( xi, points, delta ); 00124 vec_subtract( delta, x ); 00125 std::size_t num_iterations = 0; 00126 #ifdef QUADRATIC_HEX_DEBUG 00127 std::stringstream ss; 00128 ss << "Point: "; 00129 ss << x[0] << ", " << x[1] << ", " << x[2] << std::endl; 00130 ss << "Hex: "; 00131 for( int i = 0; i < 8; ++i ) 00132 { 00133 ss << points[i][0] << ", " << points[i][1] << ", " << points[i][2] << std::endl; 00134 } 00135 ss << std::endl; 00136 #endif 00137 while( normsq( delta ) > error_tol_sqr ) 00138 { 00139 #ifdef QUADRATIC_HEX_DEBUG 00140 ss << "Iter #: " << num_iterations << " Err: " << sqrt( normsq( delta ) ) << " Iterate: "; 00141 ss << xi[0] << ", " << xi[1] << ", " << xi[2] << std::endl; 00142 #endif 00143 if( ++num_iterations >= 5 ) 00144 { 00145 return false; 00146 } 00147 Matrix J; 00148 jacobian( xi, points, J ); 00149 double det = moab::Matrix::determinant3( J ); 00150 if( fabs( det ) < 1.e-10 ) 00151 { 00152 #ifdef QUADRATIC_HEX_DEBUG 00153 std::cerr << ss.str(); 00154 #endif 00155 #ifndef QUADRATIC_HEX_DEBUG 00156 std::cerr << x[0] << ", " << x[1] << ", " << x[2] << std::endl; 00157 #endif 00158 std::cerr << "inverse solve failure: det: " << det << std::endl; 00159 exit( -1 ); 00160 } 00161 vec_subtract( xi, moab::Matrix::inverse( J, 1.0 / det ) * delta ); 00162 evaluate( xi, points, delta ); 00163 vec_subtract( delta, x ); 00164 } 00165 return true; 00166 } 00167 00168 template < typename Point, typename Points > 00169 Point& evaluate( const Point& p, const Points& points, Point& f ) const 00170 { 00171 typedef typename Points::value_type Vector; 00172 Vector result; 00173 for( int i = 0; i < 3; ++i ) 00174 { 00175 result[i] = 0; 00176 } 00177 for( unsigned i = 0; i < 27; ++i ) 00178 { 00179 const double sh = SH( reference_points( i, 0 ), p[0] ) * SH( reference_points( i, 1 ), p[1] ) * 00180 SH( reference_points( i, 2 ), p[2] ); 00181 result += sh * points[i]; 00182 } 00183 for( int i = 0; i < 3; ++i ) 00184 { 00185 f[i] = result[i]; 00186 } 00187 return f; 00188 } 00189 template < typename Point, typename Field > 00190 double evaluate_scalar_field( const Point& p, const Field& field ) const 00191 { 00192 double x = 0.0; 00193 for( int i = 0; i < 27; i++ ) 00194 { 00195 const double sh = SH( reference_points( i, 0 ), p[0] ) * SH( reference_points( i, 1 ), p[1] ) * 00196 SH( reference_points( i, 2 ), p[2] ); 00197 x += sh * field[i]; 00198 } 00199 return x; 00200 } 00201 template < typename Field, typename Points > 00202 double integrate_scalar_field( const Points& p, const Field& field_values ) const 00203 { 00204 // TODO: gaussian integration , probably 2x2x2 00205 return 0.; 00206 } 00207 00208 template < typename Point, typename Points > 00209 Matrix& jacobian( const Point& p, const Points& /* points */, Matrix& J ) const 00210 { 00211 J = Matrix( 0.0 ); 00212 for( int i = 0; i < 27; i++ ) 00213 { 00214 const double sh[3] = { SH( reference_points( i, 0 ), p[0] ), SH( reference_points( i, 1 ), p[1] ), 00215 SH( reference_points( i, 2 ), p[2] ) }; 00216 const double dsh[3] = { DSH( reference_points( i, 0 ), p[0] ), DSH( reference_points( i, 1 ), p[1] ), 00217 DSH( reference_points( i, 2 ), p[2] ) }; 00218 for( int j = 0; j < 3; j++ ) 00219 { 00220 // dxj/dr first column 00221 J( j, 0 ) += dsh[0] * sh[1] * sh[2] * reference_points( i, j ); 00222 J( j, 1 ) += sh[0] * dsh[1] * sh[2] * reference_points( i, j ); // dxj/ds 00223 J( j, 2 ) += sh[0] * sh[1] * dsh[2] * reference_points( i, j ); // dxj/dt 00224 } 00225 } 00226 return J; 00227 } 00228 00229 private: 00230 }; // Class Quadratic_hex_map 00231 00232 } // namespace element_utility 00233 00234 } // namespace moab 00235 #endif // MOAB_QUADRATIC_HEX_nPP
1.7.6.1 
