moab::ElemUtil Namespace Reference

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]

Function Documentation

hex_eval()

void moab::ElemUtil::hex_eval	(	realType *	field,
		int	n,
		CartVect	rstCartVec,
		double &	value
	)

Definition at line 246 of file ElemUtil.cpp.

 {
 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 );
 }

References moab::CartVect::get(), lagrange_0(), lagrange_free(), lagrange_setup(), lobatto_nodes(), tensor_i3(), and tmalloc.

hex_findpt()

void moab::ElemUtil::hex_findpt	(	realType *	xm[3],
		int	n,
		CartVect	xyz,
		CartVect &	rst,
		double &	dist
	)

Definition at line 201 of file ElemUtil.cpp.

 {
 
 // 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] );
 }

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().

integrate_trilinear_hex()

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.

 {
 // 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;
 }

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().

nat_coords_trilinear_hex()

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.

 {
 return LinearHexMap( corner_coords ).solve_inverse( x, xi, tol );
 }

References moab::ElemUtil::VolMap::solve_inverse().

Referenced by point_in_trilinear_hex().

nat_coords_trilinear_hex2() [1/2]

void moab::ElemUtil::nat_coords_trilinear_hex2	(	const CartVect *	hex_corners,
		const CartVect &	x,
		CartVect &	xi,
		double	til
	)

nat_coords_trilinear_hex2() [2/2]

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.

 {
 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;
 }
 }
 }

References moab::CartVect::get(), hex_findpt(), and MOAB_POLY_EPS.

point_in_trilinear_hex() [1/2]

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.

 {
 // 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 );
 }

References box_max(), box_min(), dim, and point_in_trilinear_hex().

point_in_trilinear_hex() [2/2]

bool moab::ElemUtil::point_in_trilinear_hex	(	const CartVect *	hex,
		const CartVect &	xyz,
		double	etol
	)

Definition at line 162 of file ElemUtil.cpp.

 {
 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 );
 }

References nat_coords_trilinear_hex().

Referenced by main(), and point_in_trilinear_hex().

Variable Documentation

debug

bool moab::ElemUtil::debug = false

Definition at line 13 of file ElemUtil.cpp.

Referenced by integrate_trilinear_hex().

gauss_1

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().

gauss_2

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().

gauss_3

const double moab::ElemUtil::gauss_3[3][2]

Initial value:

= { { 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().

gauss_4

const double moab::ElemUtil::gauss_4[4][2]

Initial value:

= { { 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().

gauss_5

const double moab::ElemUtil::gauss_5[5][2]

Initial value:

= { { 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().