MOAB: Mesh Oriented datABase  (version 5.4.1)
Classes | Enumerations | Functions | Variables
moab::GeomUtil 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...

Enumerations

enum  intersection_type {
  NONE = 0, INTERIOR, NODE0, NODE1,
  NODE2, EDGE0, EDGE1, EDGE2
}
 Plücker test for intersection between a ray and a triangle. More...

Functions

static void min_max_3 (double a, double b, double c, double &min, double &max)
static double dot_abs (const CartVect &u, const CartVect &v)
bool segment_box_intersect (CartVect box_min, CartVect box_max, const CartVect &seg_pt, const CartVect &seg_unit_dir, double &seg_start, double &seg_end)
bool first (const CartVect &a, const CartVect &b)
double plucker_edge_test (const CartVect &vertexa, const CartVect &vertexb, const CartVect &ray, const CartVect &ray_normal)
bool plucker_ray_tri_intersect (const CartVect vertices[3], const CartVect &origin, const CartVect &direction, double &dist_out, const double *nonneg_ray_len, const double *neg_ray_len, const int *orientation, intersection_type *type)
bool ray_tri_intersect (const CartVect vertices[3], const CartVect &ray_point, const CartVect &ray_unit_direction, double &t_out, const double *ray_length=0)
 Test for intersection between a ray and a triangle.
bool ray_box_intersect (const CartVect &box_min, const CartVect &box_max, const CartVect &ray_pt, const CartVect &ray_dir, double &t_enter, double &t_exit)
bool box_plane_overlap (const CartVect &plane_normal, double plane_coeff, CartVect box_min_corner, CartVect box_max_corner)
 Test if plane intersects axis-aligned box.
bool box_tri_overlap (const CartVect triangle_corners[3], const CartVect &box_center, const CartVect &box_half_dims)
 Test if triangle intersects axis-aligned box.
bool box_tri_overlap (const CartVect triangle_corners[3], const CartVect &box_min_corner, const CartVect &box_max_corner, double tolerance)
 Test if triangle intersects axis-aligned box.
bool box_elem_overlap (const CartVect *elem_corners, EntityType elem_type, const CartVect &box_center, const CartVect &box_half_dims, int nodecount=0)
 Test if the specified element intersects an axis-aligned box.
static CartVect quad_norm (const CartVect &v1, const CartVect &v2, const CartVect &v3, const CartVect &v4)
static CartVect tri_norm (const CartVect &v1, const CartVect &v2, const CartVect &v3)
bool box_linear_elem_overlap (const CartVect *elem_corners, EntityType elem_type, const CartVect &box_center, const CartVect &box_half_dims)
 Test if the specified element intersects an axis-aligned box.
bool box_linear_elem_overlap (const CartVect *elem_corners, EntityType elem_type, const CartVect &box_half_dims)
 Test if the specified element intersects an axis-aligned box.
bool box_hex_overlap (const CartVect *elem_corners, const CartVect &center, const CartVect &dims)
static bool box_tet_overlap_edge (const CartVect &dims, const CartVect &edge, const CartVect &ve, const CartVect &v1, const CartVect &v2)
bool box_tet_overlap (const CartVect *corners_in, const CartVect &center, const CartVect &dims)
void closest_location_on_tri (const CartVect &location, const CartVect *vertices, CartVect &closest_out)
 find closest location on triangle
void closest_location_on_tri (const CartVect &location, const CartVect *vertices, double tolerance, CartVect &closest_out, int &closest_topo)
 find closest topological location on triangle
void closest_location_on_polygon (const CartVect &location, const CartVect *vertices, int num_vertices, CartVect &closest_out)
 find closest location on polygon
void closest_location_on_box (const CartVect &min, const CartVect &max, const CartVect &point, CartVect &closest)
bool box_point_overlap (const CartVect &box_min_corner, const CartVect &box_max_corner, const CartVect &point, double tolerance)
bool boxes_overlap (const CartVect &box_min1, const CartVect &box_max1, const CartVect &box_min2, const CartVect &box_max2, double tolerance)
bool bounding_boxes_overlap (const CartVect *list1, int num1, const CartVect *list2, int num2, double tolerance)
bool bounding_boxes_overlap_2d (const double *list1, int num1, const double *list2, int num2, double tolerance)
bool nat_coords_trilinear_hex (const CartVect *corner_coords, const CartVect &x, CartVect &xi, double tol)
bool point_in_trilinear_hex (const CartVect *hex, const CartVect &xyz, double etol)
bool box_hex_overlap (const CartVect hexv[8], const CartVect &box_center, const CartVect &box_dims)
bool box_tet_overlap (const CartVect tet_corners[4], const CartVect &box_center, const CartVect &box_dims)

Variables

const intersection_type type_list [] = { INTERIOR, EDGE0, EDGE1, NODE1, EDGE2, NODE0, NODE2 }

Enumeration Type Documentation

enum moab::GeomUtil::intersection_type

Plücker test for intersection between a ray and a triangle.

Parameters:
verticesNodes of the triangle.
ray_pointThe start point of the ray.
ray_unit_directionThe direction of the ray. Must be a unit vector.
t_outOutput: The distance along the ray from ray_point in the direction of ray_unit_direction at which the ray intersected the triangle.
nonneg_ray_lengthOptional: If non-null, a maximum length for the ray, converting this function to a segment-tri-intersect test.
neg_ray_lengthOptional: If non-null, a maximum length for the ray behind the origin, converting this function to a segment-tri-intersect test.
orientationOptional: Reject intersections without the specified orientation of ray with respect to triangle normal vector. Indicate desired orientation by passing 1 (forward), -1 (reverse), or 0 (no preference).
int_typeOptional Output: The type of intersection; used to identify edge/node intersections.
Returns:
true if intersection, false otherwise.
Enumerator:
NONE 
INTERIOR 
NODE0 
NODE1 
NODE2 
EDGE0 
EDGE1 
EDGE2 

Definition at line 105 of file GeomUtil.hpp.

    {
        NONE = 0,
        INTERIOR,
        NODE0,
        NODE1,
        NODE2,
        EDGE0,
        EDGE1,
        EDGE2
    };

Function Documentation

bool moab::GeomUtil::bounding_boxes_overlap ( const CartVect *  list1,
int  num1,
const CartVect *  list2,
int  num2,
double  tolerance 
)

Definition at line 1351 of file GeomUtil.cpp.

References boxes_overlap().

Referenced by moab::Intx2MeshInPlane::computeIntersectionBetweenTgtAndSrc(), moab::IntxRllCssphere::computeIntersectionBetweenTgtAndSrc(), and moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc().

    {
        assert( num1 >= 1 && num2 >= 1 );
        CartVect box_min1 = list1[0], box_max1 = list1[0];
        CartVect box_min2 = list2[0], box_max2 = list2[0];
        for( int i = 1; i < num1; i++ )
        {
            for( int k = 0; k < 3; k++ )
            {
                double val = list1[i][k];
                if( box_min1[k] > val ) box_min1[k] = val;
                if( box_max1[k] < val ) box_max1[k] = val;
            }
        }
        for( int i = 1; i < num2; i++ )
        {
            for( int k = 0; k < 3; k++ )
            {
                double val = list2[i][k];
                if( box_min2[k] > val ) box_min2[k] = val;
                if( box_max2[k] < val ) box_max2[k] = val;
            }
        }

        return boxes_overlap( box_min1, box_max1, box_min2, box_max2, tolerance );
    }
bool moab::GeomUtil::bounding_boxes_overlap_2d ( const double *  list1,
int  num1,
const double *  list2,
int  num2,
double  tolerance 
)

Definition at line 1379 of file GeomUtil.cpp.

Referenced by moab::Intx2MeshOnSphere::computeIntersectionBetweenTgtAndSrc().

    {
        /*
         * box1:
         *         (bmax11, bmax12)
         *      |-------|
         *      |       |
         *      |-------|
         *  (bmin11, bmin12)

         *         box2:
         *                (bmax21, bmax22)
         *             |----------|
         *             |          |
         *             |----------|
         *     (bmin21, bmin22)
         */
        double bmin11, bmax11, bmin12, bmax12;
        bmin11 = bmax11 = list1[0];
        bmin12 = bmax12 = list1[1];

        double bmin21, bmax21, bmin22, bmax22;
        bmin21 = bmax21 = list2[0];
        bmin22 = bmax22 = list2[1];

        for( int i = 1; i < num1; i++ )
        {
            if( bmin11 > list1[2 * i] ) bmin11 = list1[2 * i];
            if( bmax11 < list1[2 * i] ) bmax11 = list1[2 * i];
            if( bmin12 > list1[2 * i + 1] ) bmin12 = list1[2 * i + 1];
            if( bmax12 < list1[2 * i + 1] ) bmax12 = list1[2 * i + 1];
        }
        for( int i = 1; i < num2; i++ )
        {
            if( bmin21 > list2[2 * i] ) bmin21 = list2[2 * i];
            if( bmax21 < list2[2 * i] ) bmax21 = list2[2 * i];
            if( bmin22 > list2[2 * i + 1] ) bmin22 = list2[2 * i + 1];
            if( bmax22 < list2[2 * i + 1] ) bmax22 = list2[2 * i + 1];
        }

        if( ( bmax11 < bmin21 + tolerance ) || ( bmax21 < bmin11 + tolerance ) ) return false;

        if( ( bmax12 < bmin22 + tolerance ) || ( bmax22 < bmin12 + tolerance ) ) return false;

        return true;
    }
bool moab::GeomUtil::box_elem_overlap ( const CartVect *  elem_corners,
EntityType  elem_type,
const CartVect &  box_center,
const CartVect &  box_half_dims,
int  nodecount = 0 
)

Test if the specified element intersects an axis-aligned box.

Test if element intersects axis-aligned box. Use element-specific optimization if available, otherwise call box_general_elem_overlap.

Parameters:
elem_cornersThe coordinates of the element vertices
elem_typeThe toplogy of the element.
box_centerThe center of the axis-aligned box
box_half_dimsHalf of the width of the box in each axial direction.

Definition at line 519 of file GeomUtil.cpp.

References box_hex_overlap(), box_linear_elem_overlap(), box_tet_overlap(), box_tri_overlap(), MBHEX, MBPOLYGON, MBPOLYHEDRON, MBTET, and MBTRI.

Referenced by moab::AdaptiveKDTree::intersect_children_with_elems().

    {

        switch( elem_type )
        {
            case MBTRI:
                return box_tri_overlap( elem_corners, center, dims );
            case MBTET:
                return box_tet_overlap( elem_corners, center, dims );
            case MBHEX:
                return box_hex_overlap( elem_corners, center, dims );
            case MBPOLYGON: {
                CartVect vt[3];
                vt[0] = elem_corners[0];
                vt[1] = elem_corners[1];
                for( int j = 2; j < nodecount; j++ )
                {
                    vt[2] = elem_corners[j];
                    if( box_tri_overlap( vt, center, dims ) ) return true;
                }
                // none of the triangles overlap, then we return false
                return false;
            }
            case MBPOLYHEDRON:
                assert( false );
                return false;
            default:
                return box_linear_elem_overlap( elem_corners, elem_type, center, dims );
        }
    }
bool moab::GeomUtil::box_hex_overlap ( const CartVect  hexv[8],
const CartVect &  box_center,
const CartVect &  box_dims 
)
bool moab::GeomUtil::box_hex_overlap ( const CartVect *  elem_corners,
const CartVect &  center,
const CartVect &  dims 
)

Definition at line 744 of file GeomUtil.cpp.

References center(), moab::cross(), moab::dot(), dot_abs(), and quad_norm().

Referenced by box_elem_overlap(), and test_box_hex_overlap().

    {
        // Do Separating Axis Theorem:
        // If the element and the box overlap, then the 1D projections
        // onto at least one of the axes in the following three sets
        // must overlap (assuming convex polyhedral element).
        // 1) The normals of the faces of the box (the principal axes)
        // 2) The crossproduct of each element edge with each box edge
        //    (crossproduct of each edge with each principal axis)
        // 3) The normals of the faces of the element

        unsigned i, e, f;  // loop counters
        double dot, cross[2], tmp;
        CartVect norm;
        const CartVect corners[8] = { elem_corners[0] - center, elem_corners[1] - center, elem_corners[2] - center,
                                      elem_corners[3] - center, elem_corners[4] - center, elem_corners[5] - center,
                                      elem_corners[6] - center, elem_corners[7] - center };

        // test box face normals (principal axes)
        int not_less[3]    = { 8, 8, 8 };
        int not_greater[3] = { 8, 8, 8 };
        int not_inside;
        for( i = 0; i < 8; ++i )
        {  // for each element corner
            not_inside = 3;

            if( corners[i][0] < -dims[0] )
                --not_less[0];
            else if( corners[i][0] > dims[0] )
                --not_greater[0];
            else
                --not_inside;

            if( corners[i][1] < -dims[1] )
                --not_less[1];
            else if( corners[i][1] > dims[1] )
                --not_greater[1];
            else
                --not_inside;

            if( corners[i][2] < -dims[2] )
                --not_less[2];
            else if( corners[i][2] > dims[2] )
                --not_greater[2];
            else
                --not_inside;

            if( !not_inside ) return true;
        }
        // If all points less than min_x of box, then
        // not_less[0] == 0, and therefore
        // the following product is zero.
        if( not_greater[0] * not_greater[1] * not_greater[2] * not_less[0] * not_less[1] * not_less[2] == 0 )
            return false;

        // Test edge-edge crossproducts
        const unsigned edges[12][2] = { { 0, 1 }, { 0, 4 }, { 0, 3 }, { 2, 3 }, { 2, 1 }, { 2, 6 },
                                        { 5, 6 }, { 5, 1 }, { 5, 4 }, { 7, 4 }, { 7, 3 }, { 7, 6 } };

        // Edge directions for box are principal axis, so
        // for each element edge, check along the cross-product
        // of that edge with each of the tree principal axes.
        for( e = 0; e < 12; ++e )
        {  // for each element edge
            // get which element vertices bound the edge
            const CartVect& v0 = corners[edges[e][0]];
            const CartVect& v1 = corners[edges[e][1]];

            // X-Axis

            // calculate crossproduct: axis x (v1 - v0),
            // where v1 and v0 are edge vertices.
            cross[0] = v0[2] - v1[2];
            cross[1] = v1[1] - v0[1];
            // skip if parallel
            if( ( cross[0] * cross[0] + cross[1] * cross[1] ) >= std::numeric_limits< double >::epsilon() )
            {
                dot         = fabs( cross[0] * dims[1] ) + fabs( cross[1] * dims[2] );
                not_less[0] = not_greater[0] = 7;
                for( i = ( edges[e][0] + 1 ) % 8; i != edges[e][0]; i = ( i + 1 ) % 8 )
                {  // for each element corner
                    tmp = cross[0] * corners[i][1] + cross[1] * corners[i][2];
                    not_less[0] -= ( tmp < -dot );
                    not_greater[0] -= ( tmp > dot );
                }

                if( not_less[0] * not_greater[0] == 0 ) return false;
            }

            // Y-Axis

            // calculate crossproduct: axis x (v1 - v0),
            // where v1 and v0 are edge vertices.
            cross[0] = v0[0] - v1[0];
            cross[1] = v1[2] - v0[2];
            // skip if parallel
            if( ( cross[0] * cross[0] + cross[1] * cross[1] ) >= std::numeric_limits< double >::epsilon() )
            {
                dot         = fabs( cross[0] * dims[2] ) + fabs( cross[1] * dims[0] );
                not_less[0] = not_greater[0] = 7;
                for( i = ( edges[e][0] + 1 ) % 8; i != edges[e][0]; i = ( i + 1 ) % 8 )
                {  // for each element corner
                    tmp = cross[0] * corners[i][2] + cross[1] * corners[i][0];
                    not_less[0] -= ( tmp < -dot );
                    not_greater[0] -= ( tmp > dot );
                }

                if( not_less[0] * not_greater[0] == 0 ) return false;
            }

            // Z-Axis

            // calculate crossproduct: axis x (v1 - v0),
            // where v1 and v0 are edge vertices.
            cross[0] = v0[1] - v1[1];
            cross[1] = v1[0] - v0[0];
            // skip if parallel
            if( ( cross[0] * cross[0] + cross[1] * cross[1] ) >= std::numeric_limits< double >::epsilon() )
            {
                dot         = fabs( cross[0] * dims[0] ) + fabs( cross[1] * dims[1] );
                not_less[0] = not_greater[0] = 7;
                for( i = ( edges[e][0] + 1 ) % 8; i != edges[e][0]; i = ( i + 1 ) % 8 )
                {  // for each element corner
                    tmp = cross[0] * corners[i][0] + cross[1] * corners[i][1];
                    not_less[0] -= ( tmp < -dot );
                    not_greater[0] -= ( tmp > dot );
                }

                if( not_less[0] * not_greater[0] == 0 ) return false;
            }
        }

        // test element face normals
        const unsigned faces[6][4] = { { 0, 1, 2, 3 }, { 0, 1, 5, 4 }, { 1, 2, 6, 5 },
                                       { 2, 6, 7, 3 }, { 3, 7, 4, 0 }, { 7, 4, 5, 6 } };
        for( f = 0; f < 6; ++f )
        {
            norm = quad_norm( corners[faces[f][0]], corners[faces[f][1]], corners[faces[f][2]], corners[faces[f][3]] );

            dot = dot_abs( norm, dims );

            // for each element vertex
            not_less[0] = not_greater[0] = 8;
            for( i = 0; i < 8; ++i )
            {
                tmp = norm % corners[i];
                not_less[0] -= ( tmp < -dot );
                not_greater[0] -= ( tmp > dot );
            }

            if( not_less[0] * not_greater[0] == 0 ) return false;
        }

        // Overlap on all tested axes.
        return true;
    }
bool moab::GeomUtil::box_linear_elem_overlap ( const CartVect *  elem_corners,
EntityType  elem_type,
const CartVect &  box_center,
const CartVect &  box_half_dims 
)

Test if the specified element intersects an axis-aligned box.

Uses MBCN and separating axis theorem for general algorithm that works for all fixed-size elements (not poly*).

Parameters:
elem_cornersThe coordinates of the element vertices
elem_typeThe toplogy of the element.
box_centerThe center of the axis-aligned box
box_half_dimsHalf of the width of the box in each axial direction.

Definition at line 564 of file GeomUtil.cpp.

References moab::CN::VerticesPerEntity().

Referenced by box_elem_overlap(), and LinearElemOverlapTest::operator()().

    {
        CartVect corners[8];
        const unsigned num_corner = CN::VerticesPerEntity( type );
        assert( num_corner <= sizeof( corners ) / sizeof( corners[0] ) );
        for( unsigned i = 0; i < num_corner; ++i )
            corners[i] = elem_corners[i] - box_center;
        return box_linear_elem_overlap( corners, type, box_halfdims );
    }
bool moab::GeomUtil::box_linear_elem_overlap ( const CartVect *  elem_corners,
EntityType  elem_type,
const CartVect &  box_half_dims 
)

Test if the specified element intersects an axis-aligned box.

Uses MBCN and separating axis theorem for general algorithm that works for all fixed-size elements (not poly*). Box and element vertices must be translated such that box center is at origin.

Parameters:
elem_cornersThe coordinates of the element vertices, in local coordinate system of box.
elem_typeThe toplogy of the element.
box_half_dimsHalf of the width of the box in each axial direction.

Definition at line 577 of file GeomUtil.cpp.

References moab::cross(), moab::dot(), dot_abs(), MBQUAD, MBTRI, moab::CN::NumSubEntities(), quad_norm(), moab::CN::SubEntityType(), moab::CN::SubEntityVertexIndices(), tri_norm(), and moab::CN::VerticesPerEntity().

    {
        // Do Separating Axis Theorem:
        // If the element and the box overlap, then the 1D projections
        // onto at least one of the axes in the following three sets
        // must overlap (assuming convex polyhedral element).
        // 1) The normals of the faces of the box (the principal axes)
        // 2) The crossproduct of each element edge with each box edge
        //    (crossproduct of each edge with each principal axis)
        // 3) The normals of the faces of the element

        int e, f;  // loop counters
        int i;
        double dot, cross[2], tmp;
        CartVect norm;
        int indices[4];  // element edge/face vertex indices

        // test box face normals (principal axes)
        const int num_corner = CN::VerticesPerEntity( type );
        int not_less[3]      = { num_corner, num_corner, num_corner };
        int not_greater[3]   = { num_corner, num_corner, num_corner };
        int not_inside;
        for( i = 0; i < num_corner; ++i )
        {  // for each element corner
            not_inside = 3;

            if( elem_corners[i][0] < -dims[0] )
                --not_less[0];
            else if( elem_corners[i][0] > dims[0] )
                --not_greater[0];
            else
                --not_inside;

            if( elem_corners[i][1] < -dims[1] )
                --not_less[1];
            else if( elem_corners[i][1] > dims[1] )
                --not_greater[1];
            else
                --not_inside;

            if( elem_corners[i][2] < -dims[2] )
                --not_less[2];
            else if( elem_corners[i][2] > dims[2] )
                --not_greater[2];
            else
                --not_inside;

            if( !not_inside ) return true;
        }
        // If all points less than min_x of box, then
        // not_less[0] == 0, and therefore
        // the following product is zero.
        if( not_greater[0] * not_greater[1] * not_greater[2] * not_less[0] * not_less[1] * not_less[2] == 0 )
            return false;

        // Test edge-edge crossproducts

        // Edge directions for box are principal axis, so
        // for each element edge, check along the cross-product
        // of that edge with each of the tree principal axes.
        const int num_edge = CN::NumSubEntities( type, 1 );
        for( e = 0; e < num_edge; ++e )
        {  // for each element edge
            // get which element vertices bound the edge
            CN::SubEntityVertexIndices( type, 1, e, indices );

            // X-Axis

            // calculate crossproduct: axis x (v1 - v0),
            // where v1 and v0 are edge vertices.
            cross[0] = elem_corners[indices[0]][2] - elem_corners[indices[1]][2];
            cross[1] = elem_corners[indices[1]][1] - elem_corners[indices[0]][1];
            // skip if parallel
            if( ( cross[0] * cross[0] + cross[1] * cross[1] ) >= std::numeric_limits< double >::epsilon() )
            {
                dot         = fabs( cross[0] * dims[1] ) + fabs( cross[1] * dims[2] );
                not_less[0] = not_greater[0] = num_corner - 1;
                for( i = ( indices[0] + 1 ) % num_corner; i != indices[0]; i = ( i + 1 ) % num_corner )
                {  // for each element corner
                    tmp = cross[0] * elem_corners[i][1] + cross[1] * elem_corners[i][2];
                    not_less[0] -= ( tmp < -dot );
                    not_greater[0] -= ( tmp > dot );
                }

                if( not_less[0] * not_greater[0] == 0 ) return false;
            }

            // Y-Axis

            // calculate crossproduct: axis x (v1 - v0),
            // where v1 and v0 are edge vertices.
            cross[0] = elem_corners[indices[0]][0] - elem_corners[indices[1]][0];
            cross[1] = elem_corners[indices[1]][2] - elem_corners[indices[0]][2];
            // skip if parallel
            if( ( cross[0] * cross[0] + cross[1] * cross[1] ) >= std::numeric_limits< double >::epsilon() )
            {
                dot         = fabs( cross[0] * dims[2] ) + fabs( cross[1] * dims[0] );
                not_less[0] = not_greater[0] = num_corner - 1;
                for( i = ( indices[0] + 1 ) % num_corner; i != indices[0]; i = ( i + 1 ) % num_corner )
                {  // for each element corner
                    tmp = cross[0] * elem_corners[i][2] + cross[1] * elem_corners[i][0];
                    not_less[0] -= ( tmp < -dot );
                    not_greater[0] -= ( tmp > dot );
                }

                if( not_less[0] * not_greater[0] == 0 ) return false;
            }

            // Z-Axis

            // calculate crossproduct: axis x (v1 - v0),
            // where v1 and v0 are edge vertices.
            cross[0] = elem_corners[indices[0]][1] - elem_corners[indices[1]][1];
            cross[1] = elem_corners[indices[1]][0] - elem_corners[indices[0]][0];
            // skip if parallel
            if( ( cross[0] * cross[0] + cross[1] * cross[1] ) >= std::numeric_limits< double >::epsilon() )
            {
                dot         = fabs( cross[0] * dims[0] ) + fabs( cross[1] * dims[1] );
                not_less[0] = not_greater[0] = num_corner - 1;
                for( i = ( indices[0] + 1 ) % num_corner; i != indices[0]; i = ( i + 1 ) % num_corner )
                {  // for each element corner
                    tmp = cross[0] * elem_corners[i][0] + cross[1] * elem_corners[i][1];
                    not_less[0] -= ( tmp < -dot );
                    not_greater[0] -= ( tmp > dot );
                }

                if( not_less[0] * not_greater[0] == 0 ) return false;
            }
        }

        // test element face normals
        const int num_face = CN::NumSubEntities( type, 2 );
        for( f = 0; f < num_face; ++f )
        {
            CN::SubEntityVertexIndices( type, 2, f, indices );
            switch( CN::SubEntityType( type, 2, f ) )
            {
                case MBTRI:
                    norm = tri_norm( elem_corners[indices[0]], elem_corners[indices[1]], elem_corners[indices[2]] );
                    break;
                case MBQUAD:
                    norm = quad_norm( elem_corners[indices[0]], elem_corners[indices[1]], elem_corners[indices[2]],
                                      elem_corners[indices[3]] );
                    break;
                default:
                    assert( false );
                    continue;
            }

            dot = dot_abs( norm, dims );

            // for each element vertex
            not_less[0] = not_greater[0] = num_corner;
            for( i = 0; i < num_corner; ++i )
            {
                tmp = norm % elem_corners[i];
                not_less[0] -= ( tmp < -dot );
                not_greater[0] -= ( tmp > dot );
            }

            if( not_less[0] * not_greater[0] == 0 ) return false;
        }

        // Overlap on all tested axes.
        return true;
    }
bool moab::GeomUtil::box_plane_overlap ( const CartVect &  plane_normal,
double  plane_coeff,
CartVect  box_min_corner,
CartVect  box_max_corner 
)

Test if plane intersects axis-aligned box.

Test for intersection between an unbounded plane and an axis-aligned box.

Parameters:
plane_normalVector in plane normal direction (need *not* be a unit vector). The N in the plane equation: N . X + D = 0
plane_coeffThe scalar 'D' term in the plane equation: N . X + D = 0
box_min_cornerThe smallest coordinates of the box along each axis. The corner of the box for which all three coordinate values are smaller than those of any other corner. The X, Y, Z values for the planes normal to those axes and bounding the box on the -X, -Y, and -Z sides respectively.
box_max_cornerThe largest coordinates of the box along each axis. The corner of the box for which all three coordinate values are larger than those of any other corner. The X, Y, Z values for the planes normal to those axes and bounding the box on the +X, +Y, and +Z sides respectively.
Returns:
true if overlap, false otherwise.

Definition at line 399 of file GeomUtil.cpp.

References swap().

Referenced by box_tri_overlap(), and test_box_plane_norm().

    {
        if( normal[0] < 0.0 ) std::swap( min[0], max[0] );
        if( normal[1] < 0.0 ) std::swap( min[1], max[1] );
        if( normal[2] < 0.0 ) std::swap( min[2], max[2] );

        return ( normal % min <= -d ) && ( normal % max >= -d );
    }
bool moab::GeomUtil::box_point_overlap ( const CartVect &  box_min_corner,
const CartVect &  box_max_corner,
const CartVect &  point,
double  tolerance 
)

Definition at line 1322 of file GeomUtil.cpp.

References closest_location_on_box(), and moab::tolerance.

Referenced by tag_vertices().

    {
        CartVect closest;
        closest_location_on_box( box_min_corner, box_max_corner, point, closest );
        closest -= point;
        return closest % closest < tolerance * tolerance;
    }
bool moab::GeomUtil::box_tet_overlap ( const CartVect  tet_corners[4],
const CartVect &  box_center,
const CartVect &  box_dims 
)
bool moab::GeomUtil::box_tet_overlap ( const CartVect *  corners_in,
const CartVect &  center,
const CartVect &  dims 
)

Definition at line 945 of file GeomUtil.cpp.

References box_tet_overlap_edge(), center(), moab::dot(), dot_abs(), and swap().

Referenced by box_elem_overlap(), and test_box_tet_overlap().

    {
        // Do Separating Axis Theorem:
        // If the element and the box overlap, then the 1D projections
        // onto at least one of the axes in the following three sets
        // must overlap (assuming convex polyhedral element).
        // 1) The normals of the faces of the box (the principal axes)
        // 2) The crossproduct of each element edge with each box edge
        //    (crossproduct of each edge with each principal axis)
        // 3) The normals of the faces of the element

        // Translate problem such that box center is at origin.
        const CartVect corners[4] = { corners_in[0] - center, corners_in[1] - center, corners_in[2] - center,
                                      corners_in[3] - center };

        // 0) Check if any vertex is within the box
        if( fabs( corners[0][0] ) <= dims[0] && fabs( corners[0][1] ) <= dims[1] && fabs( corners[0][2] ) <= dims[2] )
            return true;
        if( fabs( corners[1][0] ) <= dims[0] && fabs( corners[1][1] ) <= dims[1] && fabs( corners[1][2] ) <= dims[2] )
            return true;
        if( fabs( corners[2][0] ) <= dims[0] && fabs( corners[2][1] ) <= dims[1] && fabs( corners[2][2] ) <= dims[2] )
            return true;
        if( fabs( corners[3][0] ) <= dims[0] && fabs( corners[3][1] ) <= dims[1] && fabs( corners[3][2] ) <= dims[2] )
            return true;

        // 1) Check for overlap on each principal axis (box face normal)
        // X
        if( corners[0][0] < -dims[0] && corners[1][0] < -dims[0] && corners[2][0] < -dims[0] &&
            corners[3][0] < -dims[0] )
            return false;
        if( corners[0][0] > dims[0] && corners[1][0] > dims[0] && corners[2][0] > dims[0] && corners[3][0] > dims[0] )
            return false;
        // Y
        if( corners[0][1] < -dims[1] && corners[1][1] < -dims[1] && corners[2][1] < -dims[1] &&
            corners[3][1] < -dims[1] )
            return false;
        if( corners[0][1] > dims[1] && corners[1][1] > dims[1] && corners[2][1] > dims[1] && corners[3][1] > dims[1] )
            return false;
        // Z
        if( corners[0][2] < -dims[2] && corners[1][2] < -dims[2] && corners[2][2] < -dims[2] &&
            corners[3][2] < -dims[2] )
            return false;
        if( corners[0][2] > dims[2] && corners[1][2] > dims[2] && corners[2][2] > dims[2] && corners[3][2] > dims[2] )
            return false;

        // 3) test element face normals
        CartVect norm;
        double dot, dot1, dot2;

        const CartVect v01 = corners[1] - corners[0];
        const CartVect v02 = corners[2] - corners[0];
        norm               = v01 * v02;
        dot                = dot_abs( norm, dims );
        dot1               = norm % corners[0];
        dot2               = norm % corners[3];
        if( dot1 > dot2 ) std::swap( dot1, dot2 );
        if( dot2 < -dot || dot1 > dot ) return false;

        const CartVect v03 = corners[3] - corners[0];
        norm               = v03 * v01;
        dot                = dot_abs( norm, dims );
        dot1               = norm % corners[0];
        dot2               = norm % corners[2];
        if( dot1 > dot2 ) std::swap( dot1, dot2 );
        if( dot2 < -dot || dot1 > dot ) return false;

        norm = v02 * v03;
        dot  = dot_abs( norm, dims );
        dot1 = norm % corners[0];
        dot2 = norm % corners[1];
        if( dot1 > dot2 ) std::swap( dot1, dot2 );
        if( dot2 < -dot || dot1 > dot ) return false;

        const CartVect v12 = corners[2] - corners[1];
        const CartVect v13 = corners[3] - corners[1];
        norm               = v13 * v12;
        dot                = dot_abs( norm, dims );
        dot1               = norm % corners[0];
        dot2               = norm % corners[1];
        if( dot1 > dot2 ) std::swap( dot1, dot2 );
        if( dot2 < -dot || dot1 > dot ) return false;

        // 2) test edge-edge cross products

        const CartVect v23 = corners[3] - corners[2];
        return box_tet_overlap_edge( dims, v01, corners[0], corners[2], corners[3] ) &&
               box_tet_overlap_edge( dims, v02, corners[0], corners[1], corners[3] ) &&
               box_tet_overlap_edge( dims, v03, corners[0], corners[1], corners[2] ) &&
               box_tet_overlap_edge( dims, v12, corners[1], corners[0], corners[3] ) &&
               box_tet_overlap_edge( dims, v13, corners[1], corners[0], corners[2] ) &&
               box_tet_overlap_edge( dims, v23, corners[2], corners[0], corners[1] );
    }
static bool moab::GeomUtil::box_tet_overlap_edge ( const CartVect &  dims,
const CartVect &  edge,
const CartVect &  ve,
const CartVect &  v1,
const CartVect &  v2 
) [inline, static]

Definition at line 901 of file GeomUtil.cpp.

References moab::dot(), and min_max_3().

Referenced by box_tet_overlap().

    {
        double dot, dot1, dot2, dot3, min, max;

        // edge x X
        if( fabs( edge[1] * edge[2] ) > std::numeric_limits< double >::epsilon() )
        {
            dot  = fabs( edge[2] ) * dims[1] + fabs( edge[1] ) * dims[2];
            dot1 = edge[2] * ve[1] - edge[1] * ve[2];
            dot2 = edge[2] * v1[1] - edge[1] * v1[2];
            dot3 = edge[2] * v2[1] - edge[1] * v2[2];
            min_max_3( dot1, dot2, dot3, min, max );
            if( max < -dot || min > dot ) return false;
        }

        // edge x Y
        if( fabs( edge[1] * edge[2] ) > std::numeric_limits< double >::epsilon() )
        {
            dot  = fabs( edge[2] ) * dims[0] + fabs( edge[0] ) * dims[2];
            dot1 = -edge[2] * ve[0] + edge[0] * ve[2];
            dot2 = -edge[2] * v1[0] + edge[0] * v1[2];
            dot3 = -edge[2] * v2[0] + edge[0] * v2[2];
            min_max_3( dot1, dot2, dot3, min, max );
            if( max < -dot || min > dot ) return false;
        }

        // edge x Z
        if( fabs( edge[1] * edge[2] ) > std::numeric_limits< double >::epsilon() )
        {
            dot  = fabs( edge[1] ) * dims[0] + fabs( edge[0] ) * dims[1];
            dot1 = edge[1] * ve[0] - edge[0] * ve[1];
            dot2 = edge[1] * v1[0] - edge[0] * v1[1];
            dot3 = edge[1] * v2[0] - edge[0] * v2[1];
            min_max_3( dot1, dot2, dot3, min, max );
            if( max < -dot || min > dot ) return false;
        }

        return true;
    }
bool moab::GeomUtil::box_tri_overlap ( const CartVect  triangle_corners[3],
const CartVect &  box_center,
const CartVect &  box_half_dims 
)

Test if triangle intersects axis-aligned box.

Test if a triangle intersects an axis-aligned box.

Parameters:
triangle_cornersThe corners of the triangle.
box_centerThe center of the box.
box_hanf_dimsThe distance along each axis, respectively, from the box_center to the boundary of the box.
Returns:
true if overlap, false otherwise.

Definition at line 425 of file GeomUtil.cpp.

References box_plane_overlap(), and CHECK_RANGE.

Referenced by box_elem_overlap(), box_tri_overlap(), general_box_tri_overlap_test(), test_box_tri_overlap(), and test_valid_tree().

    {
        // translate everything such that box is centered at origin
        const CartVect v0( vertices[0] - box_center );
        const CartVect v1( vertices[1] - box_center );
        const CartVect v2( vertices[2] - box_center );

        // do case 1) tests
        if( v0[0] > box_dims[0] && v1[0] > box_dims[0] && v2[0] > box_dims[0] ) return false;
        if( v0[1] > box_dims[1] && v1[1] > box_dims[1] && v2[1] > box_dims[1] ) return false;
        if( v0[2] > box_dims[2] && v1[2] > box_dims[2] && v2[2] > box_dims[2] ) return false;
        if( v0[0] < -box_dims[0] && v1[0] < -box_dims[0] && v2[0] < -box_dims[0] ) return false;
        if( v0[1] < -box_dims[1] && v1[1] < -box_dims[1] && v2[1] < -box_dims[1] ) return false;
        if( v0[2] < -box_dims[2] && v1[2] < -box_dims[2] && v2[2] < -box_dims[2] ) return false;

        // compute triangle edge vectors
        const CartVect e0( vertices[1] - vertices[0] );
        const CartVect e1( vertices[2] - vertices[1] );
        const CartVect e2( vertices[0] - vertices[2] );

        // do case 3) tests
        double fex, fey, fez, p0, p1, p2, rad;
        fex = fabs( e0[0] );
        fey = fabs( e0[1] );
        fez = fabs( e0[2] );

        p0  = e0[2] * v0[1] - e0[1] * v0[2];
        p2  = e0[2] * v2[1] - e0[1] * v2[2];
        rad = fez * box_dims[1] + fey * box_dims[2];
        CHECK_RANGE( p0, p2, rad );

        p0  = -e0[2] * v0[0] + e0[0] * v0[2];
        p2  = -e0[2] * v2[0] + e0[0] * v2[2];
        rad = fez * box_dims[0] + fex * box_dims[2];
        CHECK_RANGE( p0, p2, rad );

        p1  = e0[1] * v1[0] - e0[0] * v1[1];
        p2  = e0[1] * v2[0] - e0[0] * v2[1];
        rad = fey * box_dims[0] + fex * box_dims[1];
        CHECK_RANGE( p1, p2, rad );

        fex = fabs( e1[0] );
        fey = fabs( e1[1] );
        fez = fabs( e1[2] );

        p0  = e1[2] * v0[1] - e1[1] * v0[2];
        p2  = e1[2] * v2[1] - e1[1] * v2[2];
        rad = fez * box_dims[1] + fey * box_dims[2];
        CHECK_RANGE( p0, p2, rad );

        p0  = -e1[2] * v0[0] + e1[0] * v0[2];
        p2  = -e1[2] * v2[0] + e1[0] * v2[2];
        rad = fez * box_dims[0] + fex * box_dims[2];
        CHECK_RANGE( p0, p2, rad );

        p0  = e1[1] * v0[0] - e1[0] * v0[1];
        p1  = e1[1] * v1[0] - e1[0] * v1[1];
        rad = fey * box_dims[0] + fex * box_dims[1];
        CHECK_RANGE( p0, p1, rad );

        fex = fabs( e2[0] );
        fey = fabs( e2[1] );
        fez = fabs( e2[2] );

        p0  = e2[2] * v0[1] - e2[1] * v0[2];
        p1  = e2[2] * v1[1] - e2[1] * v1[2];
        rad = fez * box_dims[1] + fey * box_dims[2];
        CHECK_RANGE( p0, p1, rad );

        p0  = -e2[2] * v0[0] + e2[0] * v0[2];
        p1  = -e2[2] * v1[0] + e2[0] * v1[2];
        rad = fez * box_dims[0] + fex * box_dims[2];
        CHECK_RANGE( p0, p1, rad );

        p1  = e2[1] * v1[0] - e2[0] * v1[1];
        p2  = e2[1] * v2[0] - e2[0] * v2[1];
        rad = fey * box_dims[0] + fex * box_dims[1];
        CHECK_RANGE( p1, p2, rad );

        // do case 2) test
        CartVect n = e0 * e1;
        return box_plane_overlap( n, -( n % v0 ), -box_dims, box_dims );
    }
bool moab::GeomUtil::box_tri_overlap ( const CartVect  triangle_corners[3],
const CartVect &  box_min_corner,
const CartVect &  box_max_corner,
double  tolerance 
)

Test if triangle intersects axis-aligned box.

Test if a triangle intersects an axis-aligned box.

Parameters:
triangle_cornersThe corners of the triangle.
box_min_cornerThe smallest coordinates of the box along each axis. The corner of the box for which all three coordinate values are smaller than those of any other corner. The X, Y, Z values for the planes normal to those axes and bounding the box on the -X, -Y, and -Z sides respectively.
box_max_cornerThe largest coordinates of the box along each axis. The corner of the box for which all three coordinate values are larger than those of any other corner. The X, Y, Z values for the planes normal to those axes and bounding the box on the +X, +Y, and +Z sides respectively.
toleranceThe tolerance used in the intersection test. The box size is increased by this amount before the intersection test.
Returns:
true if overlap, false otherwise.

Definition at line 509 of file GeomUtil.cpp.

References box_tri_overlap().

    {
        const CartVect box_center = 0.5 * ( box_max_corner + box_min_corner );
        const CartVect box_hf_dim = 0.5 * ( box_max_corner - box_min_corner );
        return box_tri_overlap( triangle_corners, box_center, box_hf_dim + CartVect( tolerance ) );
    }
bool moab::GeomUtil::boxes_overlap ( const CartVect &  box_min1,
const CartVect &  box_max1,
const CartVect &  box_min2,
const CartVect &  box_max2,
double  tolerance 
)

Definition at line 1333 of file GeomUtil.cpp.

Referenced by bounding_boxes_overlap().

    {

        for( int k = 0; k < 3; k++ )
        {
            double b1min = box_min1[k], b1max = box_max1[k];
            double b2min = box_min2[k], b2max = box_max2[k];
            if( b1min - tolerance > b2max ) return false;
            if( b2min - tolerance > b1max ) return false;
        }
        return true;
    }
void moab::GeomUtil::closest_location_on_box ( const CartVect &  min,
const CartVect &  max,
const CartVect &  point,
CartVect &  closest 
)

Definition at line 1315 of file GeomUtil.cpp.

Referenced by box_point_overlap(), and test_closest_location_on_box().

    {
        closest[0] = point[0] < min[0] ? min[0] : point[0] > max[0] ? max[0] : point[0];
        closest[1] = point[1] < min[1] ? min[1] : point[1] > max[1] ? max[1] : point[1];
        closest[2] = point[2] < min[2] ? min[2] : point[2] > max[2] ? max[2] : point[2];
    }
void moab::GeomUtil::closest_location_on_polygon ( const CartVect &  location,
const CartVect *  vertices,
int  num_vertices,
CartVect &  closest_out 
)

find closest location on polygon

Find closest location on polygon

Parameters:
locationInput position to evaluate from
verticesArray of corner vertex coordinates.
num_verticesLength of 'vertices' array.
closest_outResult position

Definition at line 1244 of file GeomUtil.cpp.

References t.

Referenced by moab::OrientedBoxTreeTool::closest_to_location(), and test_closest_location_on_polygon().

    {
        const int n = num_vertices;
        CartVect d, p, v;
        double shortest_sqr, dist_sqr, t_closest, t;
        int i, e;

        // Find closest edge of polygon.
        e         = n - 1;
        v         = vertices[0] - vertices[e];
        t_closest = ( v % ( location - vertices[e] ) ) / ( v % v );
        if( t_closest < 0.0 )
            d = location - vertices[e];
        else if( t_closest > 1.0 )
            d = location - vertices[0];
        else
            d = location - vertices[e] - t_closest * v;
        shortest_sqr = d % d;
        for( i = 0; i < n - 1; ++i )
        {
            v = vertices[i + 1] - vertices[i];
            t = ( v % ( location - vertices[i] ) ) / ( v % v );
            if( t < 0.0 )
                d = location - vertices[i];
            else if( t > 1.0 )
                d = location - vertices[i + 1];
            else
                d = location - vertices[i] - t * v;
            dist_sqr = d % d;
            if( dist_sqr < shortest_sqr )
            {
                e            = i;
                shortest_sqr = dist_sqr;
                t_closest    = t;
            }
        }

        // If we are beyond the bounds of the edge, then
        // the point is outside and closest to a vertex
        if( t_closest <= 0.0 )
        {
            closest_out = vertices[e];
            return;
        }
        else if( t_closest >= 1.0 )
        {
            closest_out = vertices[( e + 1 ) % n];
            return;
        }

        // Now check which side of the edge we are one
        const CartVect v0   = vertices[e] - vertices[( e + n - 1 ) % n];
        const CartVect v1   = vertices[( e + 1 ) % n] - vertices[e];
        const CartVect v2   = vertices[( e + 2 ) % n] - vertices[( e + 1 ) % n];
        const CartVect norm = ( 1.0 - t_closest ) * ( v0 * v1 ) + t_closest * ( v1 * v2 );
        // if on outside of edge, result is closest point on edge
        if( ( norm % ( ( vertices[e] - location ) * v1 ) ) <= 0.0 )
        {
            closest_out = vertices[e] + t_closest * v1;
            return;
        }

        // Inside.  Project to plane defined by point and normal at
        // closest edge
        const double D = -( norm % ( vertices[e] + t_closest * v1 ) );
        closest_out    = ( location - ( norm % location + D ) * norm ) / ( norm % norm );
    }
void moab::GeomUtil::closest_location_on_tri ( const CartVect &  location,
const CartVect *  vertices,
CartVect &  closest_out 
)

find closest location on triangle

Find closest location on linear triangle.

Parameters:
locationInput position to evaluate from
verticesArray of three corner vertex coordinates.
closest_outResult position

Definition at line 1060 of file GeomUtil.cpp.

References t.

Referenced by check_point_in_triangles(), closest_location_on_tri(), closest_point_in_triangles(), moab::OrientedBoxTreeTool::closest_to_location(), moab::closest_to_triangles(), do_closest_point_test(), moab::AdaptiveKDTree::sphere_intersect_triangles(), moab::OrientedBoxTreeTool::sphere_intersect_triangles(), and test_closest_location_on_tri().

    {                                                    // ops      comparisons
        const CartVect sv( vertices[1] - vertices[0] );  // +3 = 3
        const CartVect tv( vertices[2] - vertices[0] );  // +3 = 6
        const CartVect pv( vertices[0] - location );     // +3 = 9
        const double ss  = sv % sv;                      // +5 = 14
        const double st  = sv % tv;                      // +5 = 19
        const double tt  = tv % tv;                      // +5 = 24
        const double sp  = sv % pv;                      // +5 = 29
        const double tp  = tv % pv;                      // +5 = 34
        const double det = ss * tt - st * st;            // +3 = 37
        double s         = st * tp - tt * sp;            // +3 = 40
        double t         = st * sp - ss * tp;            // +3 = 43
        if( s + t < det )
        {  // +1 = 44, +1 = 1
            if( s < 0 )
            {  //          +1 = 2
                if( t < 0 )
                {  //          +1 = 3
                    // region 4
                    if( sp < 0 )
                    {                                   //          +1 = 4
                        if( -sp > ss )                  //          +1 = 5
                            closest_out = vertices[1];  //      44       5
                        else
                            closest_out = vertices[0] - ( sp / ss ) * sv;  // +7 = 51,      5
                    }
                    else if( tp < 0 )
                    {                                   //          +1 = 5
                        if( -tp > tt )                  //          +1 = 6
                            closest_out = vertices[2];  //      44,      6
                        else
                            closest_out = vertices[0] - ( tp / tt ) * tv;  // +7 = 51,      6
                    }
                    else
                    {
                        closest_out = vertices[0];  //      44,      5
                    }
                }
                else
                {
                    // region 3
                    if( tp >= 0 )                   //          +1 = 4
                        closest_out = vertices[0];  //      44,      4
                    else if( -tp >= tt )            //          +1 = 5
                        closest_out = vertices[2];  //      44,      5
                    else
                        closest_out = vertices[0] - ( tp / tt ) * tv;  // +7 = 51,      5
                }
            }
            else if( t < 0 )
            {  //          +1 = 3
                // region 5;
                if( sp >= 0.0 )                 //          +1 = 4
                    closest_out = vertices[0];  //      44,      4
                else if( -sp >= ss )            //          +1 = 5
                    closest_out = vertices[1];  //      44       5
                else
                    closest_out = vertices[0] - ( sp / ss ) * sv;  // +7 = 51,      5
            }
            else
            {
                // region 0
                const double inv_det = 1.0 / det;             // +1 = 45
                s *= inv_det;                                 // +1 = 46
                t *= inv_det;                                 // +1 = 47
                closest_out = vertices[0] + s * sv + t * tv;  //+12 = 59,      3
            }
        }
        else
        {
            if( s < 0 )
            {  //          +1 = 2
                // region 2
                s = st + sp;  // +1 = 45
                t = tt + tp;  // +1 = 46
                if( t > s )
                {                                         //          +1 = 3
                    const double num = t - s;             // +1 = 47
                    const double den = ss - 2 * st + tt;  // +3 = 50
                    if( num > den )                       //          +1 = 4
                        closest_out = vertices[1];        //      50,      4
                    else
                    {
                        s           = num / den;                          // +1 = 51
                        t           = 1 - s;                              // +1 = 52
                        closest_out = s * vertices[1] + t * vertices[2];  // +9 = 61,      4
                    }
                }
                else if( t <= 0 )               //          +1 = 4
                    closest_out = vertices[2];  //      46,      4
                else if( tp >= 0 )              //          +1 = 5
                    closest_out = vertices[0];  //      46,      5
                else
                    closest_out = vertices[0] - ( tp / tt ) * tv;  // +7 = 53,      5
            }
            else if( t < 0 )
            {  //          +1 = 3
                // region 6
                t = st + tp;  // +1 = 45
                s = ss + sp;  // +1 = 46
                if( s > t )
                {                                         //          +1 = 4
                    const double num = t - s;             // +1 = 47
                    const double den = tt - 2 * st + ss;  // +3 = 50
                    if( num > den )                       //          +1 = 5
                        closest_out = vertices[2];        //      50,      5
                    else
                    {
                        t           = num / den;                          // +1 = 51
                        s           = 1 - t;                              // +1 = 52
                        closest_out = s * vertices[1] + t * vertices[2];  // +9 = 61,      5
                    }
                }
                else if( s <= 0 )               //          +1 = 5
                    closest_out = vertices[1];  //      46,      5
                else if( sp >= 0 )              //          +1 = 6
                    closest_out = vertices[0];  //      46,      6
                else
                    closest_out = vertices[0] - ( sp / ss ) * sv;  // +7 = 53,      6
            }
            else
            {
                // region 1
                const double num = tt + tp - st - sp;  // +3 = 47
                if( num <= 0 )
                {                               //          +1 = 4
                    closest_out = vertices[2];  //      47,      4
                }
                else
                {
                    const double den = ss - 2 * st + tt;  // +3 = 50
                    if( num >= den )                      //          +1 = 5
                        closest_out = vertices[1];        //      50,      5
                    else
                    {
                        s           = num / den;                          // +1 = 51
                        t           = 1 - s;                              // +1 = 52
                        closest_out = s * vertices[1] + t * vertices[2];  // +9 = 61,      5
                    }
                }
            }
        }
    }
void moab::GeomUtil::closest_location_on_tri ( const CartVect &  location,
const CartVect *  vertices,
double  tolerance,
CartVect &  closest_out,
int &  closest_topo 
)

find closest topological location on triangle

Find closest location on linear triangle.

Parameters:
locationInput position to evaluate from
verticesArray of three corner vertex coordinates.
toleranceTolerance to use when comparing to corners and edges
closest_outResult position
closest_topoClosest topological entity 0-2 : vertex index 3-5 : edge beginning at closest_topo - 3 6 : triangle interior

Definition at line 1205 of file GeomUtil.cpp.

References closest_location_on_tri(), t, and moab::tolerance.

    {
        const double tsqr = tolerance * tolerance;
        int i;
        CartVect pv[3], ev, ep;
        double t;

        closest_location_on_tri( location, vertices, closest_out );

        for( i = 0; i < 3; ++i )
        {
            pv[i] = vertices[i] - closest_out;
            if( ( pv[i] % pv[i] ) <= tsqr )
            {
                closest_topo = i;
                return;
            }
        }

        for( i = 0; i < 3; ++i )
        {
            ev = vertices[( i + 1 ) % 3] - vertices[i];
            t  = ( ev % pv[i] ) / ( ev % ev );
            ep = closest_out - ( vertices[i] + t * ev );
            if( ( ep % ep ) <= tsqr )
            {
                closest_topo = i + 3;
                return;
            }
        }

        closest_topo = 6;
    }
static double moab::GeomUtil::dot_abs ( const CartVect &  u,
const CartVect &  v 
) [inline, static]

Definition at line 65 of file GeomUtil.cpp.

Referenced by box_hex_overlap(), box_linear_elem_overlap(), and box_tet_overlap().

    {
        return fabs( u[0] * v[0] ) + fabs( u[1] * v[1] ) + fabs( u[2] * v[2] );
    }
bool moab::GeomUtil::first ( const CartVect &  a,
const CartVect &  b 
) [inline]

Definition at line 114 of file GeomUtil.cpp.

Referenced by MetisPartitioner::assemble_taggedsets_graph(), check_order_by_sets_and_adj(), moab::check_range(), moab::TypeSequenceManager::check_valid_handles(), moab::OrientedBox::covariance_data_from_tris(), moab::ParallelComm::create_interface_sets(), moab::ReadHDF5::delete_non_side_elements(), dump_sets(), moab::TypeSequenceManager::erase(), filter_options(), filter_options1(), find_coincident_elements(), moab::MeshSet::FIRST_OF_DIM(), moab::get_adjacencies_union(), moab::Core::get_connectivity(), moab::Core::get_coords(), moab::Core::get_entities_by_dimension(), moab::RangeSetIterator::get_next_by_dimension(), moab::Core::get_number_entities_by_dimension(), moab::AEntityFactory::get_zero_to_n_elements(), moab::WriteHDF5::initialize_mesh(), moab::Range::insert(), moab::intersect(), moab::WriteGMV::local_write_mesh(), moab::MergeMesh::merge_using_integer_tag(), moab::Range::num_of_dimension(), plucker_edge_test(), moab::DualTool::print_cell(), ProgOptions::printUsage(), moab::ReadHDF5::read_set_data(), moab::ParallelComm::resolve_shared_ents(), moab::BVHTree::set_interval(), moab::Bvh_tree< _Entity_handles, _Box, _Moab, _Parametrizer >::set_interval(), moab::ParCommGraph::settle_comm_by_ids(), moab::FBEngine::split_surface(), moab::Range::subset_by_dimension(), tag_info_test(), moab::ParallelMergeMesh::TagSharedElements(), test_coords_connect_iterate(), test_delete_tag_data(), test_is_free_handle(), test_iterates(), v_quad_distortion(), v_quad_oddy(), v_tri_scaled_jacobian(), and moab::ReadVtk::vtk_read_polygons().

    {
        if( a[0] < b[0] )
        {
            return true;
        }
        else if( a[0] == b[0] )
        {
            if( a[1] < b[1] )
            {
                return true;
            }
            else if( a[1] == b[1] )
            {
                if( a[2] < b[2] )
                {
                    return true;
                }
                else
                {
                    return false;
                }
            }
            else
            {
                return false;
            }
        }
        else
        {
            return false;
        }
    }
static void moab::GeomUtil::min_max_3 ( double  a,
double  b,
double  c,
double &  min,
double &  max 
) [inline, static]

Definition at line 38 of file GeomUtil.cpp.

Referenced by box_tet_overlap_edge().

    {
        if( a < b )
        {
            if( a < c )
            {
                min = a;
                max = b > c ? b : c;
            }
            else
            {
                min = c;
                max = b;
            }
        }
        else if( b < c )
        {
            min = b;
            max = a > c ? a : c;
        }
        else
        {
            min = c;
            max = a;
        }
    }
bool moab::GeomUtil::nat_coords_trilinear_hex ( const CartVect *  corner_coords,
const CartVect &  x,
CartVect &  xi,
double  tol 
)

Definition at line 1516 of file GeomUtil.cpp.

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

Referenced by point_in_trilinear_hex().

    {
        return LinearHexMap( corner_coords ).solve_inverse( x, xi, tol );
    }
double moab::GeomUtil::plucker_edge_test ( const CartVect &  vertexa,
const CartVect &  vertexb,
const CartVect &  ray,
const CartVect &  ray_normal 
)

Definition at line 148 of file GeomUtil.cpp.

References first().

Referenced by plucker_ray_tri_intersect().

    {

        double pip;
        const double near_zero = 10 * std::numeric_limits< double >::epsilon();

        if( first( vertexa, vertexb ) )
        {
            const CartVect edge        = vertexb - vertexa;
            const CartVect edge_normal = edge * vertexa;
            pip                        = ray % edge_normal + ray_normal % edge;
        }
        else
        {
            const CartVect edge        = vertexa - vertexb;
            const CartVect edge_normal = edge * vertexb;
            pip                        = ray % edge_normal + ray_normal % edge;
            pip                        = -pip;
        }

        if( near_zero > fabs( pip ) ) pip = 0.0;

        return pip;
    }
bool moab::GeomUtil::plucker_ray_tri_intersect ( const CartVect  vertices[3],
const CartVect &  origin,
const CartVect &  direction,
double &  dist_out,
const double *  nonneg_ray_len,
const double *  neg_ray_len,
const int *  orientation,
intersection_type *  type 
)

Definition at line 193 of file GeomUtil.cpp.

References EXIT_EARLY, origin, plucker_edge_test(), and type_list.

Referenced by moab::RayIntersectSets::leaf(), moab::OrientedBoxTreeTool::ray_intersect_triangles(), and test_plucker_ray_tri_intersect().

    {

        const CartVect raya = direction;
        const CartVect rayb = direction * origin;

        // Determine the value of the first Plucker coordinate from edge 0
        double plucker_coord0 = plucker_edge_test( vertices[0], vertices[1], raya, rayb );

        // If orientation is set, confirm that sign of plucker_coordinate indicate
        // correct orientation of intersection
        if( orientation && ( *orientation ) * plucker_coord0 > 0 )
        {
            EXIT_EARLY
        }

        // Determine the value of the second Plucker coordinate from edge 1
        double plucker_coord1 = plucker_edge_test( vertices[1], vertices[2], raya, rayb );

        // If orientation is set, confirm that sign of plucker_coordinate indicate
        // correct orientation of intersection
        if( orientation )
        {
            if( ( *orientation ) * plucker_coord1 > 0 )
            {
                EXIT_EARLY
            }
            // If the orientation is not specified, all plucker_coords must be the same sign or
            // zero.
        }
        else if( ( 0.0 < plucker_coord0 && 0.0 > plucker_coord1 ) || ( 0.0 > plucker_coord0 && 0.0 < plucker_coord1 ) )
        {
            EXIT_EARLY
        }

        // Determine the value of the second Plucker coordinate from edge 2
        double plucker_coord2 = plucker_edge_test( vertices[2], vertices[0], raya, rayb );

        // If orientation is set, confirm that sign of plucker_coordinate indicate
        // correct orientation of intersection
        if( orientation )
        {
            if( ( *orientation ) * plucker_coord2 > 0 )
            {
                EXIT_EARLY
            }
            // If the orientation is not specified, all plucker_coords must be the same sign or
            // zero.
        }
        else if( ( 0.0 < plucker_coord1 && 0.0 > plucker_coord2 ) || ( 0.0 > plucker_coord1 && 0.0 < plucker_coord2 ) ||
                 ( 0.0 < plucker_coord0 && 0.0 > plucker_coord2 ) || ( 0.0 > plucker_coord0 && 0.0 < plucker_coord2 ) )
        {
            EXIT_EARLY
        }

        // check for coplanar case to avoid dividing by zero
        if( 0.0 == plucker_coord0 && 0.0 == plucker_coord1 && 0.0 == plucker_coord2 )
        {
            EXIT_EARLY
        }

        // get the distance to intersection
        const double inverse_sum = 1.0 / ( plucker_coord0 + plucker_coord1 + plucker_coord2 );
        assert( 0.0 != inverse_sum );
        const CartVect intersection( plucker_coord0 * inverse_sum * vertices[2] +
                                     plucker_coord1 * inverse_sum * vertices[0] +
                                     plucker_coord2 * inverse_sum * vertices[1] );

        // To minimize numerical error, get index of largest magnitude direction.
        int idx            = 0;
        double max_abs_dir = 0;
        for( unsigned int i = 0; i < 3; ++i )
        {
            if( fabs( direction[i] ) > max_abs_dir )
            {
                idx         = i;
                max_abs_dir = fabs( direction[i] );
            }
        }
        const double dist = ( intersection[idx] - origin[idx] ) / direction[idx];

        // is the intersection within distance limits?
        if( ( nonneg_ray_len && *nonneg_ray_len < dist ) ||  // intersection is beyond positive limit
            ( neg_ray_len && *neg_ray_len >= dist ) ||       // intersection is behind negative limit
            ( !neg_ray_len && 0 > dist ) )
        {  // Unless a neg_ray_len is used, don't return negative distances
            EXIT_EARLY
        }

        dist_out = dist;

        if( type )
            *type = type_list[( ( 0.0 == plucker_coord2 ) << 2 ) + ( ( 0.0 == plucker_coord1 ) << 1 ) +
                              ( 0.0 == plucker_coord0 )];

        return true;
    }
bool moab::GeomUtil::point_in_trilinear_hex ( const CartVect *  hex,
const CartVect &  xyz,
double  etol 
)

Definition at line 1521 of file GeomUtil.cpp.

References nat_coords_trilinear_hex().

Referenced by do_kdtree_test(), do_linear_test(), and hex_containing_point().

    {
        CartVect xi;
        return nat_coords_trilinear_hex( hex, xyz, xi, etol ) && fabs( xi[0] ) - 1 < etol && fabs( xi[1] ) - 1 < etol &&
               fabs( xi[2] ) - 1 < etol;
    }
static CartVect moab::GeomUtil::quad_norm ( const CartVect &  v1,
const CartVect &  v2,
const CartVect &  v3,
const CartVect &  v4 
) [inline, static]

Definition at line 554 of file GeomUtil.cpp.

Referenced by box_hex_overlap(), and box_linear_elem_overlap().

    {
        return ( -v1 + v2 + v3 - v4 ) * ( -v1 - v2 + v3 + v4 );
    }
bool moab::GeomUtil::ray_box_intersect ( const CartVect &  box_min,
const CartVect &  box_max,
const CartVect &  ray_pt,
const CartVect &  ray_dir,
double &  t_enter,
double &  t_exit 
)

Find range of overlap between ray and axis-aligned box.

Parameters:
box_minBox corner with minimum coordinate values
box_maxBox corner with minimum coordinate values
ray_ptCoordinates of start point of ray
ray_dirDirectionion vector for ray such that the ray is defined by r(t) = ray_point + t * ray_vect for t > 0.
t_enterOutput: if return value is true, this value is the parameter location along the ray at which the ray entered the leaf. If return value is false, then this value is undefined.
t_exitOutput: if return value is true, this value is the parameter location along the ray at which the ray exited the leaf. If return value is false, then this value is undefined.
Returns:
true if ray intersects leaf, false otherwise.

Definition at line 347 of file GeomUtil.cpp.

Referenced by moab::AdaptiveKDTreeIter::intersect_ray().

    {
        const double epsilon = 1e-12;
        double t1, t2;

        // Use 'slabs' method from 13.6.1 of Akenine-Moller
        t_enter = 0.0;
        t_exit  = std::numeric_limits< double >::infinity();

        // Intersect with each pair of axis-aligned planes bounding
        // opposite faces of the leaf box
        bool ray_is_valid = false;  // is ray direction vector zero?
        for( int axis = 0; axis < 3; ++axis )
        {
            if( fabs( ray_dir[axis] ) < epsilon )
            {  // ray parallel to planes
                if( ray_pt[axis] >= box_min[axis] && ray_pt[axis] <= box_max[axis] )
                    continue;
                else
                    return false;
            }

            // find t values at which ray intersects each plane
            ray_is_valid = true;
            t1           = ( box_min[axis] - ray_pt[axis] ) / ray_dir[axis];
            t2           = ( box_max[axis] - ray_pt[axis] ) / ray_dir[axis];

            // t_enter = max( t_enter_x, t_enter_y, t_enter_z )
            // t_exit  = min( t_exit_x, t_exit_y, t_exit_z )
            //   where
            // t_enter_x = min( t1_x, t2_x );
            // t_exit_x  = max( t1_x, t2_x )
            if( t1 < t2 )
            {
                if( t_enter < t1 ) t_enter = t1;
                if( t_exit > t2 ) t_exit = t2;
            }
            else
            {
                if( t_enter < t2 ) t_enter = t2;
                if( t_exit > t1 ) t_exit = t1;
            }
        }

        return ray_is_valid && ( t_enter <= t_exit );
    }
bool moab::GeomUtil::ray_tri_intersect ( const CartVect  vertices[3],
const CartVect &  ray_point,
const CartVect &  ray_unit_direction,
double &  t_out,
const double *  ray_length = 0 
)

Test for intersection between a ray and a triangle.

Parameters:
ray_pointThe start point of the ray.
ray_unit_direcitonThe direction of the ray. Must be a unit vector.
t_outOutput: The distance along the ray from ray_point in the direction of ray_unit_direction at which the ray itersected the triangle.
ray_lengthOptional: If non-null, a pointer a maximum length for the ray, converting this function to a segment-tri- intersect test.
Returns:
true if intersection, false otherwise.

Definition at line 299 of file GeomUtil.cpp.

References t.

Referenced by moab::AdaptiveKDTree::ray_intersect_triangles(), and test_ray_tri_intersect().

    {
        const CartVect p0 = vertices[0] - vertices[1];  // abc
        const CartVect p1 = vertices[0] - vertices[2];  // def
                                                        // ghi<-v
        const CartVect p   = vertices[0] - b;           // jkl
        const CartVect c   = p1 * v;                    // eiMinushf,gfMinusdi,dhMinuseg
        const double mP    = p0 % c;
        const double betaP = p % c;
        if( mP > 0 )
        {
            if( betaP < 0 ) return false;
        }
        else if( mP < 0 )
        {
            if( betaP > 0 ) return false;
        }
        else
        {
            return false;
        }

        const CartVect d = p0 * p;  // jcMinusal,blMinuskc,akMinusjb
        double gammaP    = v % d;
        if( mP > 0 )
        {
            if( gammaP < 0 || betaP + gammaP > mP ) return false;
        }
        else if( betaP + gammaP < mP || gammaP > 0 )
            return false;

        const double tP    = p1 % d;
        const double m     = 1.0 / mP;
        const double beta  = betaP * m;
        const double gamma = gammaP * m;
        const double t     = -tP * m;
        if( ray_length && t > *ray_length ) return false;

        if( beta < 0 || gamma < 0 || beta + gamma > 1 || t < 0.0 ) return false;

        t_out = t;
        return true;
    }
bool moab::GeomUtil::segment_box_intersect ( CartVect  box_min,
CartVect  box_max,
const CartVect &  seg_pt,
const CartVect &  seg_unit_dir,
double &  seg_start,
double &  seg_end 
)

Given a line segment and an axis-aligned box, return the sub-segment of the line segment that itersects the box.

Can be used to intersect ray with box by passing seg_end as HUGE_VAL or std::numeric_limits<double>::maximum().

Parameters:
box_minMinimum corner of axis-aligned box
box_maxMaximum corner of axis-aligned box
seg_ptA point in the line containing the segement
seg_unit_dirA unit vector in the direction of the line containing the semgent.
seg_startThe distance from seg_pt in the direction of seg_unit_dir at which the segment begins. As input, the start of the original segment, as output, the start of the sub-segment intersecting the box. Note: seg_start must be less than seg_end
seg_endThe distance from seg_pt in the direction of seg_unit_dir at which the segment ends. As input, the end of the original segment, as output, the end of the sub-segment intersecting the box. Note: seg_start must be less than seg_end
Returns:
true if line semgent intersects box, false otherwise.

Definition at line 70 of file GeomUtil.cpp.

References moab::Util::is_finite().

Referenced by moab::AdaptiveKDTree::ray_intersect_triangles(), and test_segment_box_intersect().

    {
        // translate so that seg_pt is at origin
        box_min -= seg_pt;
        box_max -= seg_pt;

        for( unsigned i = 0; i < 3; ++i )
        {  // X, Y, and Z slabs

            // intersect line with slab planes
            const double t_min = box_min[i] / seg_unit_dir[i];
            const double t_max = box_max[i] / seg_unit_dir[i];

            // check if line is parallel to planes
            if( !Util::is_finite( t_min ) )
            {
                if( box_min[i] > 0.0 || box_max[i] < 0.0 ) return false;
                continue;
            }

            if( seg_unit_dir[i] < 0 )
            {
                if( t_min < seg_end ) seg_end = t_min;
                if( t_max > seg_start ) seg_start = t_max;
            }
            else
            {  // seg_unit_dir[i] > 0
                if( t_min > seg_start ) seg_start = t_min;
                if( t_max < seg_end ) seg_end = t_max;
            }
        }

        return seg_start <= seg_end;
    }
static CartVect moab::GeomUtil::tri_norm ( const CartVect &  v1,
const CartVect &  v2,
const CartVect &  v3 
) [inline, static]

Definition at line 559 of file GeomUtil.cpp.

Referenced by box_linear_elem_overlap().

    {
        return ( v2 - v1 ) * ( v3 - v1 );
    }

Variable Documentation

const intersection_type moab::GeomUtil::type_list[] = { INTERIOR, EDGE0, EDGE1, NODE1, EDGE2, NODE0, NODE2 }

Definition at line 125 of file GeomUtil.hpp.

Referenced by plucker_ray_tri_intersect().

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines