 |
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
|
Classes | |
| class | VolMap |
| Class representing a 3-D mapping function (e.g. shape function for volume element) More... | |
| class | LinearHexMap |
| Shape function for trilinear hexahedron. More... | |
Functions | |
| bool | nat_coords_trilinear_hex (const CartVect *corner_coords, const CartVect &x, CartVect &xi, double tol) |
| void | nat_coords_trilinear_hex2 (const CartVect hex[8], const CartVect &xyz, CartVect &ncoords, double etol) |
| bool | point_in_trilinear_hex (const CartVect *hex, const CartVect &xyz, double etol) |
| bool | point_in_trilinear_hex (const CartVect *hex, const CartVect &xyz, const CartVect &box_min, const CartVect &box_max, double etol) |
| void | hex_findpt (realType *xm[3], int n, CartVect xyz, CartVect &rst, double &dist) |
| void | hex_eval (realType *field, int n, CartVect rstCartVec, double &value) |
| bool | integrate_trilinear_hex (const CartVect *hex_corners, double *corner_fields, double &field_val, int num_pts) |
| void | nat_coords_trilinear_hex2 (const CartVect *hex_corners, const CartVect &x, CartVect &xi, double til) |
Variables | |
| bool | debug = false |
| const double | gauss_1 [1][2] = { { 2.0, 0.0 } } |
| const double | gauss_2 [2][2] = { { 1.0, -0.5773502691 }, { 1.0, 0.5773502691 } } |
| const double | gauss_3 [3][2] |
| const double | gauss_4 [4][2] |
| const double | gauss_5 [5][2] |
Definition at line 246 of file ElemUtil.cpp.
247 {
248 int d;
249 realType rst[3];
250 rstCartVec.get( rst );
251
252 // can cache stuff below
253 lagrange_data ld[3];
254 realType* z[3];
255 for( d = 0; d < 3; ++d )
256 {
257 z[d] = tmalloc( realType, n );
258 lobatto_nodes( z[d], n );
259 lagrange_setup( &ld[d], z[d], n );
260 }
261
262 // cut and paste -- see findpt.c
263 const unsigned nf = n * n, ne = n, nw = 2 * n * n + 3 * n;
264 realType* od_work = tmalloc( realType, 6 * nf + 9 * ne + nw );
265
266 // piece that we shouldn't want to cache
267 for( d = 0; d < 3; d++ )
268 {
269 lagrange_0( &ld[d], rst[d] );
270 }
271
272 value = tensor_i3( ld[0].J, ld[0].n, ld[1].J, ld[1].n, ld[2].J, ld[2].n, field, od_work );
273
274 // all this could be cached
275 for( d = 0; d < 3; d++ )
276 {
277 free( z[d] );
278 lagrange_free( &ld[d] );
279 }
280 free( od_work );
281 }
References moab::CartVect::get(), lagrange_0(), lagrange_free(), lagrange_setup(), lobatto_nodes(), tensor_i3(), and tmalloc.
| void moab::ElemUtil::hex_findpt | ( | realType * | xm[3], |
| int | n, | ||
| CartVect | xyz, | ||
| CartVect & | rst, | ||
| double & | dist | ||
| ) |
Definition at line 201 of file ElemUtil.cpp.
202 {
203
204 // compute stuff that only depends on the order -- could be cached
205 realType* z[3];
206 lagrange_data ld[3];
207 opt_data_3 data;
208
209 // triplicates
210 for( int d = 0; d < 3; d++ )
211 {
212 z[d] = tmalloc( realType, n );
213 lobatto_nodes( z[d], n );
214 lagrange_setup( &ld[d], z[d], n );
215 }
216
217 opt_alloc_3( &data, ld );
218
219 // find nearest point
220 realType x_star[3];
221 xyz.get( x_star );
222
223 realType r[3] = { 0, 0, 0 }; // initial guess for parametric coords
224 unsigned c = opt_no_constraints_3;
225 dist = opt_findpt_3( &data, (const realType**)xm, x_star, r, &c );
226 // c tells us if we landed inside the element or exactly on a face, edge, or node
227
228 // copy parametric coords back
229 rst = r;
230
231 // Clean-up (move to destructor if we decide to cache)
232 opt_free_3( &data );
233 for( int d = 0; d < 3; ++d )
234 lagrange_free( &ld[d] );
235 for( int d = 0; d < 3; ++d )
236 free( z[d] );
237 }
References moab::CartVect::get(), lagrange_free(), lagrange_setup(), lobatto_nodes(), opt_alloc_3(), opt_findpt_3(), opt_free_3(), opt_no_constraints_3, and tmalloc.
Referenced by nat_coords_trilinear_hex2().
| bool moab::ElemUtil::integrate_trilinear_hex | ( | const CartVect * | hex_corners, |
| double * | corner_fields, | ||
| double & | field_val, | ||
| int | num_pts | ||
| ) |
Definition at line 312 of file ElemUtil.cpp.
313 {
314 // Create the LinearHexMap object using the hex_corners array of CartVects
315 LinearHexMap hex( hex_corners );
316
317 // Use the correct table of points and locations based on the num_pts parameter
318 const double( *g_pts )[2] = 0;
319 switch( num_pts )
320 {
321 case 1:
322 g_pts = gauss_1;
323 break;
324
325 case 2:
326 g_pts = gauss_2;
327 break;
328
329 case 3:
330 g_pts = gauss_3;
331 break;
332
333 case 4:
334 g_pts = gauss_4;
335 break;
336
337 case 5:
338 g_pts = gauss_5;
339 break;
340
341 default: // value out of range
342 return false;
343 }
344
345 // Test code - print Gaussian Quadrature data
346 if( debug )
347 {
348 for( int r = 0; r < num_pts; r++ )
349 for( int c = 0; c < 2; c++ )
350 std::cout << "g_pts[" << r << "][" << c << "]=" << g_pts[r][c] << std::endl;
351 }
352 // End Test code
353
354 double soln = 0.0;
355
356 for( int i = 0; i < num_pts; i++ )
357 { // Loop for xi direction
358 double w_i = g_pts[i][0];
359 double xi_i = g_pts[i][1];
360 for( int j = 0; j < num_pts; j++ )
361 { // Loop for eta direction
362 double w_j = g_pts[j][0];
363 double eta_j = g_pts[j][1];
364 for( int k = 0; k < num_pts; k++ )
365 { // Loop for zeta direction
366 double w_k = g_pts[k][0];
367 double zeta_k = g_pts[k][1];
368
369 // Calculate the "realType" space point given the "normal" point
370 CartVect normal_pt( xi_i, eta_j, zeta_k );
371
372 // Calculate the value of F(x(xi,eta,zeta),y(xi,eta,zeta),z(xi,eta,zeta)
373 double field = hex.evaluate_scalar_field( normal_pt, corner_fields );
374
375 // Calculate the Jacobian for this "normal" point and its determinant
376 Matrix3 J = hex.jacobian( normal_pt );
377 double det = J.determinant();
378
379 // Calculate integral and update the solution
380 soln = soln + ( w_i * w_j * w_k * field * det );
381 }
382 }
383 }
384
385 // Set the output parameter
386 field_val = soln;
387
388 return true;
389 }
References debug, moab::Matrix3::determinant(), moab::ElemUtil::LinearHexMap::evaluate_scalar_field(), gauss_1, gauss_2, gauss_3, gauss_4, gauss_5, hex_corners, and moab::ElemUtil::LinearHexMap::jacobian().
| bool moab::ElemUtil::nat_coords_trilinear_hex | ( | const CartVect * | corner_coords, |
| const CartVect & | x, | ||
| CartVect & | xi, | ||
| double | tol | ||
| ) |
Definition at line 119 of file ElemUtil.cpp.
120 {
121 return LinearHexMap( corner_coords ).solve_inverse( x, xi, tol );
122 }
References moab::ElemUtil::VolMap::solve_inverse().
Referenced by point_in_trilinear_hex().
| void moab::ElemUtil::nat_coords_trilinear_hex2 | ( | const CartVect * | hex_corners, |
| const CartVect & | x, | ||
| CartVect & | xi, | ||
| double | til | ||
| ) |
| void moab::ElemUtil::nat_coords_trilinear_hex2 | ( | const CartVect | hex[8], |
| const CartVect & | xyz, | ||
| CartVect & | ncoords, | ||
| double | etol | ||
| ) |
Definition at line 127 of file ElemUtil.cpp.
129 {
130 const int ndim = 3;
131 const int nverts = 8;
132 const int vertMap[nverts] = { 0, 1, 3, 2, 4, 5, 7, 6 }; // Map from nat to lex ordering
133
134 const int n = 2; // linear
135 realType coords[ndim * nverts]; // buffer
136
137 realType* xm[ndim];
138 for( int i = 0; i < ndim; i++ )
139 xm[i] = coords + i * nverts;
140
141 // stuff hex into coords
142 for( int i = 0; i < nverts; i++ )
143 {
144 realType vcoord[ndim];
145 hex[i].get( vcoord );
146
147 for( int d = 0; d < ndim; d++ )
148 coords[d * nverts + vertMap[i]] = vcoord[d];
149 }
150
151 double dist = 0.0;
152 ElemUtil::hex_findpt( xm, n, xyz, ncoords, dist );
153 if( 3 * MOAB_POLY_EPS < dist )
154 {
155 // outside element, set extremal values to something outside range
156 for( int j = 0; j < 3; j++ )
157 {
158 if( ncoords[j] < ( -1.0 - etol ) || ncoords[j] > ( 1.0 + etol ) ) ncoords[j] *= 10;
159 }
160 }
161 }
References moab::CartVect::get(), hex_findpt(), and MOAB_POLY_EPS.
| bool moab::ElemUtil::point_in_trilinear_hex | ( | const CartVect * | hex, |
| const CartVect & | xyz, | ||
| const CartVect & | box_min, | ||
| const CartVect & | box_max, | ||
| double | etol | ||
| ) |
Definition at line 169 of file ElemUtil.cpp.
174 {
175 // all values scaled by 2 (eliminates 3 flops)
176 const CartVect mid = box_max + box_min;
177 const CartVect dim = box_max - box_min;
178 const CartVect pt = 2 * xyz - mid;
179 return ( fabs( pt[0] - dim[0] ) < etol ) && ( fabs( pt[1] - dim[1] ) < etol ) &&
180 ( fabs( pt[2] - dim[2] ) < etol ) && point_in_trilinear_hex( hex, xyz, etol );
181 }
References box_max(), box_min(), dim, and point_in_trilinear_hex().
| bool moab::ElemUtil::point_in_trilinear_hex | ( | const CartVect * | hex, |
| const CartVect & | xyz, | ||
| double | etol | ||
| ) |
Definition at line 162 of file ElemUtil.cpp.
163 {
164 CartVect xi;
165 return nat_coords_trilinear_hex( hex, xyz, xi, etol ) && ( fabs( xi[0] - 1.0 ) < etol ) &&
166 ( fabs( xi[1] - 1.0 ) < etol ) && ( fabs( xi[2] - 1.0 ) < etol );
167 }
References nat_coords_trilinear_hex().
Referenced by main(), and point_in_trilinear_hex().
| bool moab::ElemUtil::debug = false |
Definition at line 13 of file ElemUtil.cpp.
Referenced by integrate_trilinear_hex().
| const double moab::ElemUtil::gauss_1[1][2] = { { 2.0, 0.0 } } |
Definition at line 290 of file ElemUtil.cpp.
Referenced by integrate_trilinear_hex().
| const double moab::ElemUtil::gauss_2[2][2] = { { 1.0, -0.5773502691 }, { 1.0, 0.5773502691 } } |
Definition at line 292 of file ElemUtil.cpp.
Referenced by integrate_trilinear_hex().
| const double moab::ElemUtil::gauss_3[3][2] |
= { { 0.5555555555, -0.7745966692 }, { 0.8888888888, 0.0 }, { 0.5555555555, 0.7745966692 } }
Definition at line 294 of file ElemUtil.cpp.
Referenced by integrate_trilinear_hex().
| const double moab::ElemUtil::gauss_4[4][2] |
= { { 0.3478548451, -0.8611363116 }, { 0.6521451549, -0.3399810436 }, { 0.6521451549, 0.3399810436 }, { 0.3478548451, 0.8611363116 } }
Definition at line 298 of file ElemUtil.cpp.
Referenced by integrate_trilinear_hex().
| const double moab::ElemUtil::gauss_5[5][2] |
= { { 0.2369268851, -0.9061798459 }, { 0.4786286705, -0.5384693101 }, { 0.5688888889, 0.0 }, { 0.4786286705, 0.5384693101 }, { 0.2369268851, 0.9061798459 } }
Definition at line 303 of file ElemUtil.cpp.
Referenced by integrate_trilinear_hex().
1.9.1.