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