 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 #include "moab/LocalDiscretization/SpectralQuad.hpp"
2 #include "moab/Forward.hpp"
3
4 namespace moab
5 {
6
7 // filescope for static member data that is cached
8 int SpectralQuad::_n;
9 real* SpectralQuad::_z[2];
10 lagrange_data SpectralQuad::_ld[2];
11 opt_data_2 SpectralQuad::_data;
12 real* SpectralQuad::_odwork;
13 real* SpectralQuad::_glpoints;
14 bool SpectralQuad::_init = false;
15
16 SpectralQuad::SpectralQuad() : Map( 0 ) {}
17 // the preferred constructor takes pointers to GL blocked positions
18 SpectralQuad::SpectralQuad( int order, double* x, double* y, double* z ) : Map( 0 )
19 {
20 Init( order );
21 _xyz[0] = x;
22 _xyz[1] = y;
23 _xyz[2] = z;
24 }
25 SpectralQuad::SpectralQuad( int order ) : Map( 4 )
26 {
27 Init( order );
28 _xyz[0] = _xyz[1] = _xyz[2] = NULL;
29 }
30 SpectralQuad::~SpectralQuad()
31 {
32 if( _init ) freedata();
33 _init = false;
34 }
35 void SpectralQuad::Init( int order )
36 {
37 if( _init && _n == order ) return;
38 if( _init && _n != order )
39 {
40 // TODO: free data cached
41 freedata();
42 }
43 // compute stuff that depends only on order
44 _init = true;
45 _n = order;
46 // duplicates! n is the same in all directions !!!
47 for( int d = 0; d < 2; d++ )
48 {
49 _z[d] = tmalloc( double, _n );
50 lobatto_nodes( _z[d], _n );
51 lagrange_setup( &_ld[d], _z[d], _n );
52 }
53 opt_alloc_2( &_data, _ld );
54
55 unsigned int nf = _n * _n, ne = _n, nw = 2 * _n * _n + 3 * _n;
56 _odwork = tmalloc( double, 6 * nf + 9 * ne + nw );
57 _glpoints = tmalloc( double, 3 * nf );
58 }
59
60 void SpectralQuad::freedata()
61 {
62 for( int d = 0; d < 2; d++ )
63 {
64 free( _z[d] );
65 lagrange_free( &_ld[d] );
66 }
67 opt_free_2( &_data );
68 free( _odwork );
69 free( _glpoints );
70 }
71
72 void SpectralQuad::set_gl_points( double* x, double* y, double* z )
73 {
74 _xyz[0] = x;
75 _xyz[1] = y;
76 _xyz[2] = z;
77 }
78 CartVect SpectralQuad::evalFcn( const double* params,
79 const double* field,
80 const int ndim,
81 const int num_tuples,
82 double* work,
83 double* result )
84 {
85 // piece that we shouldn't want to cache
86 int d = 0;
87 for( d = 0; d < 2; d++ )
88 {
89 lagrange_0( &_ld[d], params[d] );
90 }
91 CartVect result;
92 for( d = 0; d < 3; d++ )
93 {
94 result[d] = tensor_i2( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _xyz[d], _odwork );
95 }
96 return result;
97 }
98 // replicate the functionality of hex_findpt
99 bool SpectralQuad::reverseEvalFcn( const double* posn,
100 const double* verts,
101 const int nverts,
102 const int ndim,
103 const double iter_tol,
104 const double inside_tol,
105 double* work,
106 double* params,
107 int* is_inside )
108 {
109 params = init;
110
111 // find nearest point
112 double x_star[3];
113 xyz.get( x_star );
114
115 double r[2] = { 0, 0 }; // initial guess for parametric coords
116 unsigned c = opt_no_constraints_3;
117 double dist = opt_findpt_2( &_data, (const double**)_xyz, x_star, r, &c );
118 // if it did not converge, get out with throw...
119 if( dist > 0.9e+30 )
120 {
121 std::vector< CartVect > dummy;
122 throw Map::EvaluationError( CartVect( x_star ), dummy );
123 }
124 // c tells us if we landed inside the element or exactly on a face, edge, or node
125 // also, dist shows the distance to the computed point.
126 // copy parametric coords back
127 params = r;
128
129 return insideFcn( params, 2, inside_tol );
130 }
131
132 Matrix3 SpectralQuad::jacobian( const double* params,
133 const double* verts,
134 const int nverts,
135 const int ndim,
136 double* work,
137 double* result )
138 {
139 // not implemented
140 Matrix3 J( 0. );
141 return J;
142 }
143
144 void SpectralQuad::evaluate_vector( const CartVect& params, const double* field, int num_tuples, double* eval ) const
145 {
146 // piece that we shouldn't want to cache
147 int d;
148 for( d = 0; d < 2; d++ )
149 {
150 lagrange_0( &_ld[d], params[d] );
151 }
152
153 *eval = tensor_i2( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, field, _odwork );
154 }
155 void SpectralQuad::integrate_vector( const double* field,
156 const double* verts,
157 const int nverts,
158 const int ndim,
159 const int num_tuples,
160 double* work,
161 double* result )
162 {
163 // not implemented
164 }
165
166 int SpectralQuad::insideFcn( const double* params, const int ndim, const double tol )
167 {
168 return EvalSet::inside( params, ndim, tol );
169 }
170
171 // something we don't do for spectral hex; we do it here because
172 // we do not store the position of gl points in a tag yet
173 void SpectralQuad::compute_gl_positions()
174 {
175 // will need to use shape functions on a simple linear quad to compute gl points
176 // so we know the position of gl points in parametric space, we will just compute those
177 // from the 3d vertex position (corner nodes of the quad), using simple mapping
178 assert( this->vertex.size() == 4 );
179 static double corner_params[4][2] = { { -1., -1. }, { 1., -1. }, { 1., 1. }, { -1., 1. } };
180 // we will use the cached lobatto nodes in parametric space _z[d] (the same in both directions)
181 int indexGL = 0;
182 int n2 = _n * _n;
183 for( int i = 0; i < _n; i++ )
184 {
185 double csi = _z[0][i];
186 for( int j = 0; j < _n; j++ )
187 {
188 double eta = _z[1][j]; // we could really use the same _z[0] array of lobatto nodes
189 CartVect pos( 0.0 );
190 for( int k = 0; k < 4; k++ )
191 {
192 const double N_k = ( 1 + csi * corner_params[k][0] ) * ( 1 + eta * corner_params[k][1] );
193 pos += N_k * vertex[k];
194 }
195 pos *= 0.25; // these are x, y, z of gl points; reorder them
196 _glpoints[indexGL] = pos[0]; // x
197 _glpoints[indexGL + n2] = pos[1]; // y
198 _glpoints[indexGL + 2 * n2] = pos[2]; // z
199 indexGL++;
200 }
201 }
202 // now, we can set the _xyz pointers to internal memory allocated to these!
203 _xyz[0] = &( _glpoints[0] );
204 _xyz[1] = &( _glpoints[n2] );
205 _xyz[2] = &( _glpoints[2 * n2] );
206 }
207 void SpectralQuad::get_gl_points( double*& x, double*& y, double*& z, int& size )
208 {
209 x = (double*)_xyz[0];
210 y = (double*)_xyz[1];
211 z = (double*)_xyz[2];
212 size = _n * _n;
213 }
214
215 } // namespace moab
1.9.1.