MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 #include "moab/LocalDiscretization/SpectralHex.hpp" 00002 #include "moab/Forward.hpp" 00003 00004 namespace moab 00005 { 00006 00007 // SpectralHex 00008 00009 // filescope for static member data that is cached 00010 int SpectralHex::_n; 00011 real* SpectralHex::_z[3]; 00012 lagrange_data SpectralHex::_ld[3]; 00013 opt_data_3 SpectralHex::_data; 00014 real* SpectralHex::_odwork; 00015 00016 bool SpectralHex::_init = false; 00017 00018 void SpectralHex::initFcn( const double* verts, const int nverts, double*& work ) 00019 { 00020 if( _init && _n == order ) return; 00021 if( _init && _n != order ) 00022 { 00023 // TODO: free data cached 00024 freedata(); 00025 } 00026 // compute stuff that depends only on order 00027 _init = true; 00028 _n = order; 00029 // triplicates! n is the same in all directions !!! 00030 for( int d = 0; d < 3; d++ ) 00031 { 00032 _z[d] = tmalloc( double, _n ); 00033 lobatto_nodes( _z[d], _n ); 00034 lagrange_setup( &_ld[d], _z[d], _n ); 00035 } 00036 opt_alloc_3( &_data, _ld ); 00037 00038 unsigned int nf = _n * _n, ne = _n, nw = 2 * _n * _n + 3 * _n; 00039 _odwork = tmalloc( double, 6 * nf + 9 * ne + nw ); 00040 } 00041 void SpectralHex::freedata() 00042 { 00043 for( int d = 0; d < 3; d++ ) 00044 { 00045 free( _z[d] ); 00046 lagrange_free( &_ld[d] ); 00047 } 00048 opt_free_3( &_data ); 00049 free( _odwork ); 00050 } 00051 00052 void SpectralHex::set_gl_points( double* x, double* y, double* z ) 00053 { 00054 _xyz[0] = x; 00055 _xyz[1] = y; 00056 _xyz[2] = z; 00057 } 00058 CartVect SpectralHex::evaluate( const CartVect& params ) const 00059 { 00060 // piece that we shouldn't want to cache 00061 int d = 0; 00062 for( d = 0; d < 3; d++ ) 00063 { 00064 lagrange_0( &_ld[d], params[d] ); 00065 } 00066 CartVect result; 00067 for( d = 0; d < 3; d++ ) 00068 { 00069 result[d] = tensor_i3( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _ld[2].J, _ld[2].n, 00070 _xyz[d], // this is the "field" 00071 _odwork ); 00072 } 00073 return result; 00074 } 00075 // replicate the functionality of hex_findpt 00076 bool SpectralHex::evaluate_reverse( CartVect const& xyz, CartVect& params, double iter_tol, const double inside_tol, 00077 const CartVect& init ) const 00078 { 00079 params = init; 00080 00081 // find nearest point 00082 double x_star[3]; 00083 xyz.get( x_star ); 00084 00085 double r[3] = { 0, 0, 0 }; // initial guess for parametric coords 00086 unsigned c = opt_no_constraints_3; 00087 double dist = opt_findpt_3( &_data, (const double**)_xyz, x_star, r, &c ); 00088 // if it did not converge, get out with throw... 00089 if( dist > 0.9e+30 ) return false; 00090 // c tells us if we landed inside the element or exactly on a face, edge, or node 00091 // also, dist shows the distance to the computed point. 00092 // copy parametric coords back 00093 params = r; 00094 00095 return is_inside( params, inside_tol ); 00096 } 00097 Matrix3 SpectralHex::jacobian( const CartVect& params ) const 00098 { 00099 real x_i[3]; 00100 params.get( x_i ); 00101 // set the positions of GL nodes, before evaluations 00102 _data.elx[0] = _xyz[0]; 00103 _data.elx[1] = _xyz[1]; 00104 _data.elx[2] = _xyz[2]; 00105 opt_vol_set_intp_3( &_data, x_i ); 00106 Matrix3 J( 0. ); 00107 // it is organized differently 00108 J( 0, 0 ) = _data.jac[0]; // dx/dr 00109 J( 0, 1 ) = _data.jac[1]; // dx/ds 00110 J( 0, 2 ) = _data.jac[2]; // dx/dt 00111 J( 1, 0 ) = _data.jac[3]; // dy/dr 00112 J( 1, 1 ) = _data.jac[4]; // dy/ds 00113 J( 1, 2 ) = _data.jac[5]; // dy/dt 00114 J( 2, 0 ) = _data.jac[6]; // dz/dr 00115 J( 2, 1 ) = _data.jac[7]; // dz/ds 00116 J( 2, 2 ) = _data.jac[8]; // dz/dt 00117 return J; 00118 } 00119 void SpectralHex::evaluate_vector( const CartVect& params, const double* field, int num_tuples, double* eval ) const 00120 { 00121 // piece that we shouldn't want to cache 00122 int d; 00123 for( d = 0; d < 3; d++ ) 00124 { 00125 lagrange_0( &_ld[d], params[d] ); 00126 } 00127 00128 *eval = tensor_i3( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _ld[2].J, _ld[2].n, field, _odwork ); 00129 } 00130 void SpectralHex::integrate_vector( const double* field_values, int num_tuples, double* integral ) const 00131 { 00132 // set the position of GL points 00133 // set the positions of GL nodes, before evaluations 00134 _data.elx[0] = _xyz[0]; 00135 _data.elx[1] = _xyz[1]; 00136 _data.elx[2] = _xyz[2]; 00137 real params[3]; 00138 // triple loop; the most inner loop is in r direction, then s, then t 00139 for( int l = 0; l < num_tuples; l++ ) 00140 integral[l] = 0.0; 00141 // double volume = 0; 00142 int index = 0; // used fr the inner loop 00143 for( int k = 0; k < _n; k++ ) 00144 { 00145 params[2] = _ld[2].z[k]; 00146 // double wk= _ld[2].w[k]; 00147 for( int j = 0; j < _n; j++ ) 00148 { 00149 params[1] = _ld[1].z[j]; 00150 // double wj= _ld[1].w[j]; 00151 for( int i = 0; i < _n; i++ ) 00152 { 00153 params[0] = _ld[0].z[i]; 00154 // double wi= _ld[0].w[i]; 00155 opt_vol_set_intp_3( &_data, params ); 00156 double wk = _ld[2].J[k]; 00157 double wj = _ld[1].J[j]; 00158 double wi = _ld[0].J[i]; 00159 Matrix3 J( 0. ); 00160 // it is organized differently 00161 J( 0, 0 ) = _data.jac[0]; // dx/dr 00162 J( 0, 1 ) = _data.jac[1]; // dx/ds 00163 J( 0, 2 ) = _data.jac[2]; // dx/dt 00164 J( 1, 0 ) = _data.jac[3]; // dy/dr 00165 J( 1, 1 ) = _data.jac[4]; // dy/ds 00166 J( 1, 2 ) = _data.jac[5]; // dy/dt 00167 J( 2, 0 ) = _data.jac[6]; // dz/dr 00168 J( 2, 1 ) = _data.jac[7]; // dz/ds 00169 J( 2, 2 ) = _data.jac[8]; // dz/dt 00170 double bm = wk * wj * wi * J.determinant(); 00171 for( int l = 0; l < num_tuples; l++ ) 00172 integral[l] += bm * field_values[num_tuples * index + l]; 00173 // volume +=bm; 00174 } 00175 } 00176 } 00177 // std::cout << "volume: " << volume << "\n"; 00178 } 00179 00180 int SpectralHex::insideFcn( const double* params, const int ndim, const double tol ) 00181 { 00182 return EvalSet::inside( params, ndim, tol ); 00183 } 00184 00185 // SpectralHex 00186 00187 } // namespace moab