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