MOAB: Mesh Oriented datABase  (version 5.4.1)
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

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

Variable Documentation

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

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

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

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines