moab
|
00001 #ifndef MOAB_SPECTRAL_HEX_HPP 00002 #define MOAB_SPECTRAL_HEX_HPP 00003 00004 #include "moab/Matrix3.hpp" 00005 #include "moab/CartVect.hpp" 00006 #include "moab/FindPtFuncs.h" 00007 #include <sstream> 00008 #include <iomanip> 00009 #include <iostream> 00010 00011 namespace moab { 00012 00013 namespace element_utility { 00014 00015 namespace {} //non-exported functionality 00016 00017 template< typename _Matrix> 00018 class Spectral_hex_map { 00019 public: 00020 typedef _Matrix Matrix; 00021 private: 00022 typedef Spectral_hex_map< Matrix> Self; 00023 public: 00024 //Constructor 00025 Spectral_hex_map() {}; 00026 Spectral_hex_map( int order){ initialize_spectral_hex( order); } 00027 //Copy constructor 00028 Spectral_hex_map( const Self & f ) {} 00029 private: 00030 void initialize_spectral_hex( int order){ 00031 if (_init && _n==order){ return; } 00032 if( _init && _n != order){ free_data();} 00033 _init = true; 00034 _n = order; 00035 for( int d = 0; d < 3; d++){ 00036 lobatto_nodes(_z[ d], _n); 00037 lagrange_setup(&_ld[ d], _z[ d], _n); 00038 } 00039 opt_alloc_3(&_data, _ld); 00040 std::size_t nf = _n*_n, ne = _n, nw = 2*_n*_n + 3*_n; 00041 _odwork = tmalloc(real, 6*nf + 9*ne + nw); 00042 } 00043 00044 void free_data(){ 00045 for(int d=0; d<3; d++){ 00046 free(_z[d]); 00047 lagrange_free(&_ld[d]); 00048 } 00049 opt_free_3(&_data); 00050 free(_odwork); 00051 } 00052 00053 public: 00054 //Natural coordinates 00055 template< typename Moab, typename Entity_handle, 00056 typename Points, typename Point> 00057 std::pair< bool, Point> operator()( const Moab & /* moab */, 00058 const Entity_handle & /* h */, 00059 const Points & v, 00060 const Point & p, 00061 const double tol = 1.e-6) { 00062 Point result(3, 0.0); 00063 /* 00064 moab.tag_get_by_ptr(_xm1Tag, &eh, 1,(const void **) &_xyz[ 0] ); 00065 moab.tag_get_by_ptr(_ym1Tag, &eh, 1,(const void **) &_xyz[ 1] ); 00066 moab.tag_get_by_ptr(_zm1Tag, &eh, 1,(const void **) &_xyz[ 2] ); 00067 */ 00068 bool point_found = solve_inverse( p, result, v, tol) && 00069 is_contained( result, tol); 00070 return std::make_pair( point_found, result); 00071 } 00072 00073 private: 00074 void set_gl_points( double * x, double * y, double * z){ 00075 _xyz[ 0] = x; _xyz[ 1] = y; _xyz[ 2] = z; 00076 } 00077 template< typename Point> 00078 bool is_contained( const Point & p, const double tol) const{ 00079 //just look at the box+tol here 00080 return ( p[0]>=-1.-tol) && (p[0]<=1.+tol) && 00081 ( p[1]>=-1.-tol) && (p[1]<=1.+tol) && 00082 ( p[2]>=-1.-tol) && (p[2]<=1.+tol); 00083 } 00084 00085 template< typename Point, typename Points> 00086 bool solve_inverse( const Point & x, 00087 Point & xi, 00088 const Points & points, 00089 const double tol=1.e-6) { 00090 const double error_tol_sqr = tol*tol; 00091 Point delta(3,0.0); 00092 xi = delta; 00093 evaluate( xi, points, delta); 00094 vec_subtract( delta, x); 00095 std::size_t num_iterations=0; 00096 #ifdef SPECTRAL_HEX_DEBUG 00097 std::stringstream ss; 00098 ss << "Point: "; 00099 ss << x[ 0 ] << ", " << x[ 1] 00100 << ", " << x [ 2] << std::endl; 00101 ss << "Hex: "; 00102 for(int i = 0; i < 8; ++i){ 00103 ss << points[ i][ 0] << ", " << points[ i][ 1] << ", " 00104 << points[ i][ 2] << std::endl; 00105 } 00106 ss << std::endl; 00107 #endif 00108 while ( normsq( delta) > error_tol_sqr) { 00109 #ifdef SPECTRAL_HEX_DEBUG 00110 ss << "Iter #: " << num_iterations 00111 << " Err: " << sqrt( normsq( delta)) << " Iterate: "; 00112 ss << xi[ 0 ] << ", " << xi[ 1] 00113 << ", " << xi[ 2] << std::endl; 00114 #endif 00115 if( ++num_iterations >= 5){ return false; } 00116 Matrix J; 00117 jacobian( xi, points, J); 00118 double det = moab::Matrix::determinant3( J); 00119 if (fabs(det) < 1.e-10){ 00120 #ifdef SPECTRAL_HEX_DEBUG 00121 std::cerr << ss.str(); 00122 #endif 00123 #ifndef SPECTRAL_HEX_DEBUG 00124 std::cerr << x[ 0 ] << ", " << x[ 1] 00125 << ", " << x [ 2] << std::endl; 00126 #endif 00127 std::cerr << "inverse solve failure: det: " << det << std::endl; 00128 exit( -1); 00129 } 00130 vec_subtract( xi, moab::Matrix::inverse(J, 1.0/det) * delta); 00131 evaluate( xi, points, delta); 00132 vec_subtract( delta, x); 00133 } 00134 return true; 00135 } 00136 00137 template< typename Point, typename Points> 00138 Point& evaluate( const Point & p, const Points & /* points */, Point & f) { 00139 for (int d = 0; d < 3; ++d) { lagrange_0(&_ld[ d], p[ 0]); } 00140 for (int d = 0; d < 3; ++d) { 00141 f[ d] = tensor_i3( _ld[ 0].J, _ld[ 0].n, 00142 _ld[1].J, _ld[1].n, 00143 _ld[2].J, _ld[2].n, 00144 _xyz[ d], 00145 _odwork); 00146 } 00147 return f; 00148 } 00149 00150 template< typename Point, typename Field> 00151 double evaluate_scalar_field(const Point & p, const Field & field) const { 00152 int d; 00153 for(d=0; d<3; d++){ lagrange_0(&_ld[d], p[d]); } 00154 return tensor_i3( _ld[0].J,_ld[0].n, 00155 _ld[1].J,_ld[1].n, 00156 _ld[2].J,_ld[2].n, 00157 field, _odwork); 00158 } 00159 template< typename Points, typename Field> 00160 double integrate_scalar_field(const Points & p, 00161 const Field & field) const { 00162 // set the position of GL points 00163 // set the positions of GL nodes, before evaluations 00164 _data.elx[0]=_xyz[0]; 00165 _data.elx[1]=_xyz[1]; 00166 _data.elx[2]=_xyz[2]; 00167 double xi[3]; 00168 //triple loop; the most inner loop is in r direction, then s, then t 00169 double integral = 0.; 00170 //double volume = 0; 00171 int index=0; // used fr the inner loop 00172 for (int k=0; k<_n; k++ ) { 00173 xi[2]=_ld[2].z[k]; 00174 //double wk= _ld[2].w[k]; 00175 for (int j=0; j<_n; j++) { 00176 xi[1]=_ld[1].z[j]; 00177 //double wj= _ld[1].w[j]; 00178 for (int i=0; i<_n; i++) { 00179 xi[0]=_ld[0].z[i]; 00180 //double wi= _ld[0].w[i]; 00181 opt_vol_set_intp_3((opt_data_3 *)&_data,xi); // cast away const-ness 00182 double wk= _ld[2].J[k]; 00183 double wj= _ld[1].J[j]; 00184 double wi= _ld[0].J[i]; 00185 Matrix3 J(0.); 00186 for (int n = 0; n < 8; n++) 00187 J(n/3, n%3) = _data.jac[n]; 00188 double bm = wk*wj*wi* J.determinant(); 00189 integral+= bm*field[index++]; 00190 //volume +=bm; 00191 } 00192 } 00193 } 00194 //std::cout << "volume: " << volume << "\n"; 00195 return integral; 00196 } 00197 00198 template< typename Point, typename Points> 00199 Matrix& jacobian( const Point & /* p */, const Points & /* points */, Matrix & J) { 00200 real x[ 3]; 00201 for (int i = 0; i < 3; ++i) { _data.elx[ i] = _xyz[ i]; } 00202 opt_vol_set_intp_3(& _data, x); 00203 for (int i = 0; i < 9; ++i) { J(i%3, i/3) = _data.jac[ i]; } 00204 return J; 00205 } 00206 00207 private: 00208 bool _init; 00209 int _n; 00210 real * _z[ 3]; 00211 lagrange_data _ld[ 3]; 00212 opt_data_3 _data; 00213 real * _odwork; 00214 real * _xyz[ 3]; 00215 }; //Class Spectral_hex_map 00216 00217 }// namespace element_utility 00218 00219 }// namespace moab 00220 #endif //MOAB_SPECTRAL_HEX_nPP