 |
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/SpectralHex.hpp"
2 #include "moab/Forward.hpp"
3
4 namespace moab
5 {
6
7 // SpectralHex
8
9 // filescope for static member data that is cached
10 int SpectralHex::_n;
11 real* SpectralHex::_z[3];
12 lagrange_data SpectralHex::_ld[3];
13 opt_data_3 SpectralHex::_data;
14 real* SpectralHex::_odwork;
15
16 bool SpectralHex::_init = false;
17
18 void SpectralHex::initFcn( const double* verts, const int nverts, double*& work )
19 {
20 if( _init && _n == order ) return;
21 if( _init && _n != order )
22 {
23 // TODO: free data cached
24 freedata();
25 }
26 // compute stuff that depends only on order
27 _init = true;
28 _n = order;
29 // triplicates! n is the same in all directions !!!
30 for( int d = 0; d < 3; d++ )
31 {
32 _z[d] = tmalloc( double, _n );
33 lobatto_nodes( _z[d], _n );
34 lagrange_setup( &_ld[d], _z[d], _n );
35 }
36 opt_alloc_3( &_data, _ld );
37
38 unsigned int nf = _n * _n, ne = _n, nw = 2 * _n * _n + 3 * _n;
39 _odwork = tmalloc( double, 6 * nf + 9 * ne + nw );
40 }
41 void SpectralHex::freedata()
42 {
43 for( int d = 0; d < 3; d++ )
44 {
45 free( _z[d] );
46 lagrange_free( &_ld[d] );
47 }
48 opt_free_3( &_data );
49 free( _odwork );
50 }
51
52 void SpectralHex::set_gl_points( double* x, double* y, double* z )
53 {
54 _xyz[0] = x;
55 _xyz[1] = y;
56 _xyz[2] = z;
57 }
58 CartVect SpectralHex::evaluate( const CartVect& params ) const
59 {
60 // piece that we shouldn't want to cache
61 int d = 0;
62 for( d = 0; d < 3; d++ )
63 {
64 lagrange_0( &_ld[d], params[d] );
65 }
66 CartVect result;
67 for( d = 0; d < 3; d++ )
68 {
69 result[d] = tensor_i3( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _ld[2].J, _ld[2].n,
70 _xyz[d], // this is the "field"
71 _odwork );
72 }
73 return result;
74 }
75 // replicate the functionality of hex_findpt
76 bool SpectralHex::evaluate_reverse( CartVect const& xyz,
77 CartVect& params,
78 double iter_tol,
79 const double inside_tol,
80 const CartVect& init ) const
81 {
82 params = init;
83
84 // find nearest point
85 double x_star[3];
86 xyz.get( x_star );
87
88 double r[3] = { 0, 0, 0 }; // initial guess for parametric coords
89 unsigned c = opt_no_constraints_3;
90 double dist = opt_findpt_3( &_data, (const double**)_xyz, x_star, r, &c );
91 // if it did not converge, get out with throw...
92 if( dist > 0.9e+30 ) return false;
93 // c tells us if we landed inside the element or exactly on a face, edge, or node
94 // also, dist shows the distance to the computed point.
95 // copy parametric coords back
96 params = r;
97
98 return is_inside( params, inside_tol );
99 }
100 Matrix3 SpectralHex::jacobian( const CartVect& params ) const
101 {
102 real x_i[3];
103 params.get( x_i );
104 // set the positions of GL nodes, before evaluations
105 _data.elx[0] = _xyz[0];
106 _data.elx[1] = _xyz[1];
107 _data.elx[2] = _xyz[2];
108 opt_vol_set_intp_3( &_data, x_i );
109 Matrix3 J( 0. );
110 // it is organized differently
111 J( 0, 0 ) = _data.jac[0]; // dx/dr
112 J( 0, 1 ) = _data.jac[1]; // dx/ds
113 J( 0, 2 ) = _data.jac[2]; // dx/dt
114 J( 1, 0 ) = _data.jac[3]; // dy/dr
115 J( 1, 1 ) = _data.jac[4]; // dy/ds
116 J( 1, 2 ) = _data.jac[5]; // dy/dt
117 J( 2, 0 ) = _data.jac[6]; // dz/dr
118 J( 2, 1 ) = _data.jac[7]; // dz/ds
119 J( 2, 2 ) = _data.jac[8]; // dz/dt
120 return J;
121 }
122 void SpectralHex::evaluate_vector( const CartVect& params, const double* field, int num_tuples, double* eval ) const
123 {
124 // piece that we shouldn't want to cache
125 int d;
126 for( d = 0; d < 3; d++ )
127 {
128 lagrange_0( &_ld[d], params[d] );
129 }
130
131 *eval = tensor_i3( _ld[0].J, _ld[0].n, _ld[1].J, _ld[1].n, _ld[2].J, _ld[2].n, field, _odwork );
132 }
133 void SpectralHex::integrate_vector( const double* field_values, int num_tuples, double* integral ) const
134 {
135 // set the position of GL points
136 // set the positions of GL nodes, before evaluations
137 _data.elx[0] = _xyz[0];
138 _data.elx[1] = _xyz[1];
139 _data.elx[2] = _xyz[2];
140 real params[3];
141 // triple loop; the most inner loop is in r direction, then s, then t
142 for( int l = 0; l < num_tuples; l++ )
143 integral[l] = 0.0;
144 // double volume = 0;
145 int index = 0; // used fr the inner loop
146 for( int k = 0; k < _n; k++ )
147 {
148 params[2] = _ld[2].z[k];
149 // double wk= _ld[2].w[k];
150 for( int j = 0; j < _n; j++ )
151 {
152 params[1] = _ld[1].z[j];
153 // double wj= _ld[1].w[j];
154 for( int i = 0; i < _n; i++ )
155 {
156 params[0] = _ld[0].z[i];
157 // double wi= _ld[0].w[i];
158 opt_vol_set_intp_3( &_data, params );
159 double wk = _ld[2].J[k];
160 double wj = _ld[1].J[j];
161 double wi = _ld[0].J[i];
162 Matrix3 J( 0. );
163 // it is organized differently
164 J( 0, 0 ) = _data.jac[0]; // dx/dr
165 J( 0, 1 ) = _data.jac[1]; // dx/ds
166 J( 0, 2 ) = _data.jac[2]; // dx/dt
167 J( 1, 0 ) = _data.jac[3]; // dy/dr
168 J( 1, 1 ) = _data.jac[4]; // dy/ds
169 J( 1, 2 ) = _data.jac[5]; // dy/dt
170 J( 2, 0 ) = _data.jac[6]; // dz/dr
171 J( 2, 1 ) = _data.jac[7]; // dz/ds
172 J( 2, 2 ) = _data.jac[8]; // dz/dt
173 double bm = wk * wj * wi * J.determinant();
174 for( int l = 0; l < num_tuples; l++ )
175 integral[l] += bm * field_values[num_tuples * index + l];
176 // volume +=bm;
177 }
178 }
179 }
180 // std::cout << "volume: " << volume << "\n";
181 }
182
183 int SpectralHex::insideFcn( const double* params, const int ndim, const double tol )
184 {
185 return EvalSet::inside( params, ndim, tol );
186 }
187
188 // SpectralHex
189
190 } // namespace moab
1.9.1.