MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 #include "moab/LocalDiscretization/SpectralQuad.hpp" 00002 #include "moab/Forward.hpp" 00003 00004 namespace moab 00005 { 00006 00007 // filescope for static member data that is cached 00008 int SpectralQuad::_n; 00009 real* SpectralQuad::_z[2]; 00010 lagrange_data SpectralQuad::_ld[2]; 00011 opt_data_2 SpectralQuad::_data; 00012 real* SpectralQuad::_odwork; 00013 real* SpectralQuad::_glpoints; 00014 bool SpectralQuad::_init = false; 00015 00016 SpectralQuad::SpectralQuad() : Map( 0 ) {} 00017 // the preferred constructor takes pointers to GL blocked positions 00018 SpectralQuad::SpectralQuad( int order, double* x, double* y, double* z ) : Map( 0 ) 00019 { 00020 Init( order ); 00021 _xyz[0] = x; 00022 _xyz[1] = y; 00023 _xyz[2] = z; 00024 } 00025 SpectralQuad::SpectralQuad( int order ) : Map( 4 ) 00026 { 00027 Init( order ); 00028 _xyz[0] = _xyz[1] = _xyz[2] = NULL; 00029 } 00030 SpectralQuad::~SpectralQuad() 00031 { 00032 if( _init ) freedata(); 00033 _init = false; 00034 } 00035 void SpectralQuad::Init( int order ) 00036 { 00037 if( _init && _n == order ) return; 00038 if( _init && _n != order ) 00039 { 00040 // TODO: free data cached 00041 freedata(); 00042 } 00043 // compute stuff that depends only on order 00044 _init = true; 00045 _n = order; 00046 // duplicates! n is the same in all directions !!! 00047 for( int d = 0; d < 2; d++ ) 00048 { 00049 _z[d] = tmalloc( double, _n ); 00050 lobatto_nodes( _z[d], _n ); 00051 lagrange_setup( &_ld[d], _z[d], _n ); 00052 } 00053 opt_alloc_2( &_data, _ld ); 00054 00055 unsigned int nf = _n * _n, ne = _n, nw = 2 * _n * _n + 3 * _n; 00056 _odwork = tmalloc( double, 6 * nf + 9 * ne + nw ); 00057 _glpoints = tmalloc( double, 3 * nf ); 00058 } 00059 00060 void SpectralQuad::freedata() 00061 { 00062 for( int d = 0; d < 2; d++ ) 00063 { 00064 free( _z[d] ); 00065 lagrange_free( &_ld[d] ); 00066 } 00067 opt_free_2( &_data ); 00068 free( _odwork ); 00069 free( _glpoints ); 00070 } 00071 00072 void SpectralQuad::set_gl_points( double* x, double* y, double* z ) 00073 { 00074 _xyz[0] = x; 00075 _xyz[1] = y; 00076 _xyz[2] = z; 00077 } 00078 CartVect SpectralQuad::evalFcn( const double* params, const double* field, const int ndim, const int num_tuples, 00079 double* work, double* result ) 00080 { 00081 // piece that we shouldn't want to cache 00082 int d = 0; 00083 for( d = 0; d < 2; d++ ) 00084 { 00085 lagrange_0( &_ld[d], params[d] ); 00086 } 00087 CartVect result; 00088 for( d = 0; d < 3; d++ ) 00089 { 00090 result[d] = tensor_i2( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _xyz[d], _odwork ); 00091 } 00092 return result; 00093 } 00094 // replicate the functionality of hex_findpt 00095 bool SpectralQuad::reverseEvalFcn( const double* posn, const double* verts, const int nverts, const int ndim, 00096 const double iter_tol, const double inside_tol, double* work, double* params, 00097 int* is_inside ) 00098 { 00099 params = init; 00100 00101 // find nearest point 00102 double x_star[3]; 00103 xyz.get( x_star ); 00104 00105 double r[2] = { 0, 0 }; // initial guess for parametric coords 00106 unsigned c = opt_no_constraints_3; 00107 double dist = opt_findpt_2( &_data, (const double**)_xyz, x_star, r, &c ); 00108 // if it did not converge, get out with throw... 00109 if( dist > 0.9e+30 ) 00110 { 00111 std::vector< CartVect > dummy; 00112 throw Map::EvaluationError( CartVect( x_star ), dummy ); 00113 } 00114 // c tells us if we landed inside the element or exactly on a face, edge, or node 00115 // also, dist shows the distance to the computed point. 00116 // copy parametric coords back 00117 params = r; 00118 00119 return insideFcn( params, 2, inside_tol ); 00120 } 00121 00122 Matrix3 SpectralQuad::jacobian( const double* params, const double* verts, const int nverts, const int ndim, 00123 double* work, double* result ) 00124 { 00125 // not implemented 00126 Matrix3 J( 0. ); 00127 return J; 00128 } 00129 00130 void SpectralQuad::evaluate_vector( const CartVect& params, const double* field, int num_tuples, double* eval ) const 00131 { 00132 // piece that we shouldn't want to cache 00133 int d; 00134 for( d = 0; d < 2; d++ ) 00135 { 00136 lagrange_0( &_ld[d], params[d] ); 00137 } 00138 00139 *eval = tensor_i2( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, field, _odwork ); 00140 } 00141 void SpectralQuad::integrate_vector( const double* field, const double* verts, const int nverts, const int ndim, 00142 const int num_tuples, double* work, double* result ) 00143 { 00144 // not implemented 00145 } 00146 00147 int SpectralQuad::insideFcn( const double* params, const int ndim, const double tol ) 00148 { 00149 return EvalSet::inside( params, ndim, tol ); 00150 } 00151 00152 // something we don't do for spectral hex; we do it here because 00153 // we do not store the position of gl points in a tag yet 00154 void SpectralQuad::compute_gl_positions() 00155 { 00156 // will need to use shape functions on a simple linear quad to compute gl points 00157 // so we know the position of gl points in parametric space, we will just compute those 00158 // from the 3d vertex position (corner nodes of the quad), using simple mapping 00159 assert( this->vertex.size() == 4 ); 00160 static double corner_params[4][2] = { { -1., -1. }, { 1., -1. }, { 1., 1. }, { -1., 1. } }; 00161 // we will use the cached lobatto nodes in parametric space _z[d] (the same in both directions) 00162 int indexGL = 0; 00163 int n2 = _n * _n; 00164 for( int i = 0; i < _n; i++ ) 00165 { 00166 double csi = _z[0][i]; 00167 for( int j = 0; j < _n; j++ ) 00168 { 00169 double eta = _z[1][j]; // we could really use the same _z[0] array of lobatto nodes 00170 CartVect pos( 0.0 ); 00171 for( int k = 0; k < 4; k++ ) 00172 { 00173 const double N_k = ( 1 + csi * corner_params[k][0] ) * ( 1 + eta * corner_params[k][1] ); 00174 pos += N_k * vertex[k]; 00175 } 00176 pos *= 0.25; // these are x, y, z of gl points; reorder them 00177 _glpoints[indexGL] = pos[0]; // x 00178 _glpoints[indexGL + n2] = pos[1]; // y 00179 _glpoints[indexGL + 2 * n2] = pos[2]; // z 00180 indexGL++; 00181 } 00182 } 00183 // now, we can set the _xyz pointers to internal memory allocated to these! 00184 _xyz[0] = &( _glpoints[0] ); 00185 _xyz[1] = &( _glpoints[n2] ); 00186 _xyz[2] = &( _glpoints[2 * n2] ); 00187 } 00188 void SpectralQuad::get_gl_points( double*& x, double*& y, double*& z, int& size ) 00189 { 00190 x = (double*)_xyz[0]; 00191 y = (double*)_xyz[1]; 00192 z = (double*)_xyz[2]; 00193 size = _n * _n; 00194 } 00195 00196 } // namespace moab