![]() |
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
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] |
void moab::ElemUtil::hex_eval | ( | realType * | field, |
int | n, | ||
CartVect | rstCartVec, | ||
double & | value | ||
) |
Definition at line 246 of file ElemUtil.cpp.
References moab::CartVect::get(), lagrange_0(), lagrange_free(), lagrange_setup(), lobatto_nodes(), tensor_i3(), and tmalloc.
{
int d;
realType rst[3];
rstCartVec.get( rst );
// can cache stuff below
lagrange_data ld[3];
realType* z[3];
for( d = 0; d < 3; ++d )
{
z[d] = tmalloc( realType, n );
lobatto_nodes( z[d], n );
lagrange_setup( &ld[d], z[d], n );
}
// cut and paste -- see findpt.c
const unsigned nf = n * n, ne = n, nw = 2 * n * n + 3 * n;
realType* od_work = tmalloc( realType, 6 * nf + 9 * ne + nw );
// piece that we shouldn't want to cache
for( d = 0; d < 3; d++ )
{
lagrange_0( &ld[d], rst[d] );
}
value = tensor_i3( ld[0].J, ld[0].n, ld[1].J, ld[1].n, ld[2].J, ld[2].n, field, od_work );
// all this could be cached
for( d = 0; d < 3; d++ )
{
free( z[d] );
lagrange_free( &ld[d] );
}
free( od_work );
}
void moab::ElemUtil::hex_findpt | ( | realType * | xm[3], |
int | n, | ||
CartVect | xyz, | ||
CartVect & | rst, | ||
double & | dist | ||
) |
Definition at line 201 of file ElemUtil.cpp.
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().
{
// compute stuff that only depends on the order -- could be cached
realType* z[3];
lagrange_data ld[3];
opt_data_3 data;
// triplicates
for( int d = 0; d < 3; d++ )
{
z[d] = tmalloc( realType, n );
lobatto_nodes( z[d], n );
lagrange_setup( &ld[d], z[d], n );
}
opt_alloc_3( &data, ld );
// find nearest point
realType x_star[3];
xyz.get( x_star );
realType r[3] = { 0, 0, 0 }; // initial guess for parametric coords
unsigned c = opt_no_constraints_3;
dist = opt_findpt_3( &data, (const realType**)xm, x_star, r, &c );
// c tells us if we landed inside the element or exactly on a face, edge, or node
// copy parametric coords back
rst = r;
// Clean-up (move to destructor if we decide to cache)
opt_free_3( &data );
for( int d = 0; d < 3; ++d )
lagrange_free( &ld[d] );
for( int d = 0; d < 3; ++d )
free( z[d] );
}
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.
References debug, moab::Matrix3::determinant(), moab::ElemUtil::LinearHexMap::evaluate_scalar_field(), gauss_1, gauss_2, gauss_3, gauss_4, gauss_5, and moab::ElemUtil::LinearHexMap::jacobian().
{
// Create the LinearHexMap object using the hex_corners array of CartVects
LinearHexMap hex( hex_corners );
// Use the correct table of points and locations based on the num_pts parameter
const double( *g_pts )[2] = 0;
switch( num_pts )
{
case 1:
g_pts = gauss_1;
break;
case 2:
g_pts = gauss_2;
break;
case 3:
g_pts = gauss_3;
break;
case 4:
g_pts = gauss_4;
break;
case 5:
g_pts = gauss_5;
break;
default: // value out of range
return false;
}
// Test code - print Gaussian Quadrature data
if( debug )
{
for( int r = 0; r < num_pts; r++ )
for( int c = 0; c < 2; c++ )
std::cout << "g_pts[" << r << "][" << c << "]=" << g_pts[r][c] << std::endl;
}
// End Test code
double soln = 0.0;
for( int i = 0; i < num_pts; i++ )
{ // Loop for xi direction
double w_i = g_pts[i][0];
double xi_i = g_pts[i][1];
for( int j = 0; j < num_pts; j++ )
{ // Loop for eta direction
double w_j = g_pts[j][0];
double eta_j = g_pts[j][1];
for( int k = 0; k < num_pts; k++ )
{ // Loop for zeta direction
double w_k = g_pts[k][0];
double zeta_k = g_pts[k][1];
// Calculate the "realType" space point given the "normal" point
CartVect normal_pt( xi_i, eta_j, zeta_k );
// Calculate the value of F(x(xi,eta,zeta),y(xi,eta,zeta),z(xi,eta,zeta)
double field = hex.evaluate_scalar_field( normal_pt, corner_fields );
// Calculate the Jacobian for this "normal" point and its determinant
Matrix3 J = hex.jacobian( normal_pt );
double det = J.determinant();
// Calculate integral and update the solution
soln = soln + ( w_i * w_j * w_k * field * det );
}
}
}
// Set the output parameter
field_val = soln;
return true;
}
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.
References moab::ElemUtil::VolMap::solve_inverse().
Referenced by point_in_trilinear_hex().
{
return LinearHexMap( corner_coords ).solve_inverse( x, xi, tol );
}
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.
References moab::CartVect::get(), hex_findpt(), and MOAB_POLY_EPS.
{
const int ndim = 3;
const int nverts = 8;
const int vertMap[nverts] = { 0, 1, 3, 2, 4, 5, 7, 6 }; // Map from nat to lex ordering
const int n = 2; // linear
realType coords[ndim * nverts]; // buffer
realType* xm[ndim];
for( int i = 0; i < ndim; i++ )
xm[i] = coords + i * nverts;
// stuff hex into coords
for( int i = 0; i < nverts; i++ )
{
realType vcoord[ndim];
hex[i].get( vcoord );
for( int d = 0; d < ndim; d++ )
coords[d * nverts + vertMap[i]] = vcoord[d];
}
double dist = 0.0;
ElemUtil::hex_findpt( xm, n, xyz, ncoords, dist );
if( 3 * MOAB_POLY_EPS < dist )
{
// outside element, set extremal values to something outside range
for( int j = 0; j < 3; j++ )
{
if( ncoords[j] < ( -1.0 - etol ) || ncoords[j] > ( 1.0 + etol ) ) ncoords[j] *= 10;
}
}
}
bool moab::ElemUtil::point_in_trilinear_hex | ( | const CartVect * | hex, |
const CartVect & | xyz, | ||
double | etol | ||
) |
Definition at line 162 of file ElemUtil.cpp.
References nat_coords_trilinear_hex().
Referenced by main(), and point_in_trilinear_hex().
{
CartVect xi;
return nat_coords_trilinear_hex( hex, xyz, xi, etol ) && ( fabs( xi[0] - 1.0 ) < etol ) &&
( fabs( xi[1] - 1.0 ) < etol ) && ( fabs( xi[2] - 1.0 ) < etol );
}
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.
References box_min(), dim, and point_in_trilinear_hex().
{
// all values scaled by 2 (eliminates 3 flops)
const CartVect mid = box_max + box_min;
const CartVect dim = box_max - box_min;
const CartVect pt = 2 * xyz - mid;
return ( fabs( pt[0] - dim[0] ) < etol ) && ( fabs( pt[1] - dim[1] ) < etol ) &&
( fabs( pt[2] - dim[2] ) < etol ) && point_in_trilinear_hex( hex, xyz, etol );
}
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().