 |
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
|
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
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
00007 #include
00008 #include
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.
Dark theme by Tilen Majerle. All rights reserved.