1 /* 2 * MOAB, a Mesh-Oriented datABase, is a software component for creating, 3 * storing and accessing finite element mesh data. 4 * 5 * Copyright 2004 Sandia Corporation. Under the terms of Contract 6 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government 7 * retains certain rights in this software. 8 * 9 * This library is free software; you can redistribute it and/or 10 * modify it under the terms of the GNU Lesser General Public 11 * License as published by the Free Software Foundation; either 12 * version 2.1 of the License, or (at your option) any later version. 13 * 14 */ 15 16 /**\file Geometry.hpp 17 *\author Jason Kraftcheck (kraftche@cae.wisc.edu) 18 *\date 2006-07-27 19 */ 20 21 #ifndef MB_GEOM_UTIL_HPP 22 #define MB_GEOM_UTIL_HPP 23 24 #include "moab/CartVect.hpp" 25 #include <cmath> 26 27 namespace moab { 28 29 /** \class GeomUtil 30 * \brief Functions for computational geometry on triangles, rays, and boxes 31 */ 32 namespace GeomUtil { 33 34 /** Given a line segment and an axis-aligned box, 35 * return the sub-segment of the line segment that 36 * itersects the box. 37 * 38 * Can be used to intersect ray with box by passing seg_end 39 * as HUGE_VAL or std::numeric_limits<double>::maximum(). 40 * 41 *\param box_min Minimum corner of axis-aligned box 42 *\param box_max Maximum corner of axis-aligned box 43 *\param seg_pt A point in the line containing the segement 44 *\param seg_unit_dir A unit vector in the direction of the line 45 * containing the semgent. 46 *\param seg_start The distance from seg_pt in the direction of 47 * seg_unit_dir at which the segment begins. 48 * As input, the start of the original segment, as output, the 49 * start of the sub-segment intersecting the box. 50 * Note: seg_start must be less than seg_end 51 *\param seg_end The distance from seg_pt in the direction of 52 * seg_unit_dir at which the segment ends. 53 * As input, the end of the original segment, as output, the 54 * end of the sub-segment intersecting the box. 55 * Note: seg_start must be less than seg_end 56 *\return true if line semgent intersects box, false otherwise. 57 */ 58 bool segment_box_intersect( CartVect box_min, 59 CartVect box_max, 60 const CartVect& seg_pt, 61 const CartVect& seg_unit_dir, 62 double& seg_start, double& seg_end ); 63 64 /**\brief Test for intersection between a ray and a triangle. 65 *\param ray_point The start point of the ray. 66 *\param ray_unit_direciton The direction of the ray. Must be a unit vector. 67 *\param t_out Output: The distance along the ray from ray_point in the 68 * direction of ray_unit_direction at which the ray 69 * itersected the triangle. 70 *\param ray_length Optional: If non-null, a pointer a maximum length 71 * for the ray, converting this function to a segment-tri- 72 * intersect test. 73 *\return true if intersection, false otherwise. 74 */ 75 bool ray_tri_intersect( const CartVect vertices[3], 76 const CartVect& ray_point, 77 const CartVect& ray_unit_direction, 78 double& t_out, 79 const double* ray_length = 0 ); 80 81 82 /**\brief Plücker test for intersection between a ray and a triangle. 83 *\param vertices Nodes of the triangle. 84 *\param ray_point The start point of the ray. 85 *\param ray_unit_direction The direction of the ray. Must be a unit vector. 86 *\param t_out Output: The distance along the ray from ray_point in the 87 * direction of ray_unit_direction at which the ray 88 * intersected the triangle. 89 *\param nonneg_ray_length Optional: If non-null, a maximum length for the ray, 90 * converting this function to a segment-tri-intersect 91 * test. 92 *\param neg_ray_length Optional: If non-null, a maximum length for the ray 93 * behind the origin, converting this function to a 94 * segment-tri-intersect test. 95 *\param orientation Optional: Reject intersections without the specified 96 * orientation of ray with respect to triangle normal 97 * vector. Indicate desired orientation by passing 98 * 1 (forward), -1 (reverse), or 0 (no preference). 99 *\param int_type Optional Output: The type of intersection; used to 100 * identify edge/node intersections. 101 *\return true if intersection, false otherwise. 102 */ 103 enum intersection_type {NONE=0, INTERIOR, NODE0, NODE1, NODE2, EDGE0, EDGE1, EDGE2}; 104 /* intersection type is determined by which of the intermediate values are == 0. There 105 are three such values that can therefore be encoded in a 3 bit integer. 106 0 = none are == 0 -> interior type 107 1 = pip0 == 0 -> EDGE0 108 2 = pip1 == 1 -> EDGE1 109 4 = pip2 == 2 -> EDGE2 110 5 = pip2 = pip0 == 0 -> NOEE0 111 3 = pip0 = pip1 == 0 -> NODE1 112 6 = pip1 = pip2 == 0 -> NODE2 */ 113 const intersection_type type_list[] = {INTERIOR, EDGE0, EDGE1, NODE1, EDGE2, NODE0, NODE2}; 114 115 bool plucker_ray_tri_intersect( const CartVect vertices[3], 116 const CartVect& ray_point, 117 const CartVect& ray_unit_direction, 118 double& dist_out, 119 const double* nonneg_ray_length = 0, 120 const double* neg_ray_length = 0, 121 const int* orientation = 0, 122 intersection_type* int_type = 0); 123 double plucker_edge_test( const CartVect& vertexa, 124 const CartVect& vertexb, 125 const CartVect& ray, 126 const CartVect& ray_normal); 127 128 //! Find range of overlap between ray and axis-aligned box. 129 //! 130 //!\param box_min Box corner with minimum coordinate values 131 //!\param box_max Box corner with minimum coordinate values 132 //!\param ray_pt Coordinates of start point of ray 133 //!\param ray_dir Directionion vector for ray such that 134 //! the ray is defined by r(t) = ray_point + t * ray_vect 135 //! for t > 0. 136 //!\param t_enter Output: if return value is true, this value 137 //! is the parameter location along the ray at which 138 //! the ray entered the leaf. If return value is false, 139 //! then this value is undefined. 140 //!\param t_exit Output: if return value is true, this value 141 //! is the parameter location along the ray at which 142 //! the ray exited the leaf. If return value is false, 143 //! then this value is undefined. 144 //!\return true if ray intersects leaf, false otherwise. 145 bool ray_box_intersect( const CartVect& box_min, 146 const CartVect& box_max, 147 const CartVect& ray_pt, 148 const CartVect& ray_dir, 149 double& t_enter, double& t_exit ); 150 151 /**\brief Test if plane intersects axis-aligned box 152 * 153 * Test for intersection between an unbounded plane and 154 * an axis-aligned box. 155 *\param plane_normal Vector in plane normal direction (need *not* 156 * be a unit vector). The N in 157 * the plane equation: N . X + D = 0 158 *\param plane_coeff The scalar 'D' term in the plane equation: 159 * N . X + D = 0 160 *\param box_min_corner The smallest coordinates of the box along each 161 * axis. The corner of the box for which all three 162 * coordinate values are smaller than those of any 163 * other corner. The X, Y, Z values for the planes 164 * normal to those axes and bounding the box on the 165 * -X, -Y, and -Z sides respectively. 166 *\param box_max_corner The largest coordinates of the box along each 167 * axis. The corner of the box for which all three 168 * coordinate values are larger than those of any 169 * other corner. The X, Y, Z values for the planes 170 * normal to those axes and bounding the box on the 171 * +X, +Y, and +Z sides respectively. 172 *\return true if overlap, false otherwise. 173 */ 174 bool box_plane_overlap( const CartVect& plane_normal, 175 double plane_coeff, 176 CartVect box_min_corner, 177 CartVect box_max_corner ); 178 179 /**\brief Test if triangle intersects axis-aligned box 180 * 181 * Test if a triangle intersects an axis-aligned box. 182 *\param triangle_corners The corners of the triangle. 183 *\param box_min_corner The smallest coordinates of the box along each 184 * axis. The corner of the box for which all three 185 * coordinate values are smaller than those of any 186 * other corner. The X, Y, Z values for the planes 187 * normal to those axes and bounding the box on the 188 * -X, -Y, and -Z sides respectively. 189 *\param box_max_corner The largest coordinates of the box along each 190 * axis. The corner of the box for which all three 191 * coordinate values are larger than those of any 192 * other corner. The X, Y, Z values for the planes 193 * normal to those axes and bounding the box on the 194 * +X, +Y, and +Z sides respectively. 195 *\param tolerance The tolerance used in the intersection test. The box 196 * size is increased by this amount before the intersection 197 * test. 198 *\return true if overlap, false otherwise. 199 */ 200 bool box_tri_overlap( const CartVect triangle_corners[3], 201 const CartVect& box_min_corner, 202 const CartVect& box_max_corner, 203 double tolerance ); 204 205 /**\brief Test if triangle intersects axis-aligned box 206 * 207 * Test if a triangle intersects an axis-aligned box. 208 *\param triangle_corners The corners of the triangle. 209 *\param box_center The center of the box. 210 *\param box_hanf_dims The distance along each axis, respectively, from the 211 * box_center to the boundary of the box. 212 *\return true if overlap, false otherwise. 213 */ 214 bool box_tri_overlap( const CartVect triangle_corners[3], 215 const CartVect& box_center, 216 const CartVect& box_half_dims ); 217 218 bool box_point_overlap( const CartVect& box_min_corner, 219 const CartVect& box_max_corner, 220 const CartVect& point, 221 double tolerance ); 222 223 /**\brief Test if the specified element intersects an axis-aligned box. 224 * 225 * Test if element intersects axis-aligned box. Use element-specific 226 * optimization if available, otherwise call box_general_elem_overlap. 227 * 228 *\param elem_corners The coordinates of the element vertices 229 *\param elem_type The toplogy of the element. 230 *\param box_center The center of the axis-aligned box 231 *\param box_half_dims Half of the width of the box in each axial 232 * direction. 233 */ 234 bool box_elem_overlap( const CartVect *elem_corners, 235 EntityType elem_type, 236 const CartVect& box_center, 237 const CartVect& box_half_dims ); 238 239 /**\brief Test if the specified element intersects an axis-aligned box. 240 * 241 * Uses MBCN and separating axis theorem for general algorithm that 242 * works for all fixed-size elements (not poly*). 243 * 244 *\param elem_corners The coordinates of the element vertices 245 *\param elem_type The toplogy of the element. 246 *\param box_center The center of the axis-aligned box 247 *\param box_half_dims Half of the width of the box in each axial 248 * direction. 249 */ 250 bool box_linear_elem_overlap( const CartVect *elem_corners, 251 EntityType elem_type, 252 const CartVect& box_center, 253 const CartVect& box_half_dims ); 254 255 /**\brief Test if the specified element intersects an axis-aligned box. 256 * 257 * Uses MBCN and separating axis theorem for general algorithm that 258 * works for all fixed-size elements (not poly*). Box and element 259 * vertices must be translated such that box center is at origin. 260 * 261 *\param elem_corners The coordinates of the element vertices, in 262 * local coordinate system of box. 263 *\param elem_type The toplogy of the element. 264 *\param box_half_dims Half of the width of the box in each axial 265 * direction. 266 */ 267 bool box_linear_elem_overlap( const CartVect *elem_corners, 268 EntityType elem_type, 269 const CartVect& box_half_dims ); 270 271 void closest_location_on_box( const CartVect& box_min_corner, 272 const CartVect& box_max_corner, 273 const CartVect& point, 274 CartVect& closest ); 275 276 /**\brief find closest location on triangle 277 * 278 * Find closest location on linear triangle. 279 *\param location Input position to evaluate from 280 *\param vertices Array of three corner vertex coordinates. 281 *\param closest_out Result position 282 */ 283 void closest_location_on_tri( const CartVect& location, 284 const CartVect* vertices, 285 CartVect& closest_out ); 286 287 /**\brief find closest location on polygon 288 * 289 * Find closest location on polygon 290 *\param location Input position to evaluate from 291 *\param vertices Array of corner vertex coordinates. 292 *\param num_vertices Length of 'vertices' array. 293 *\param closest_out Result position 294 */ 295 void closest_location_on_polygon( const CartVect& location, 296 const CartVect* vertices, 297 int num_vertices, 298 CartVect& closest_out ); 299 300 /**\brief find closest topological location on triangle 301 * 302 * Find closest location on linear triangle. 303 *\param location Input position to evaluate from 304 *\param vertices Array of three corner vertex coordinates. 305 *\param tolerance Tolerance to use when comparing to corners and edges 306 *\param closest_out Result position 307 *\param closest_topo Closest topological entity 308 * 0-2 : vertex index 309 * 3-5 : edge beginning at closest_topo - 3 310 * 6 : triangle interior 311 */ 312 void closest_location_on_tri( const CartVect& location, 313 const CartVect* vertices, 314 double tolerance, 315 CartVect& closest_out, 316 int& closest_topo ); 317 318 // Finds whether or not a box defined by the center and the half 319 // width intersects a trilinear hex defined by its eight vertices. 320 bool box_hex_overlap( const CartVect hexv[8], 321 const CartVect& box_center, 322 const CartVect& box_dims); 323 324 // Finds whether or not a box defined by the center and the half 325 // width intersects a linear tetrahedron defined by its four vertices. 326 bool box_tet_overlap( const CartVect tet_corners[4], 327 const CartVect& box_center, 328 const CartVect& box_dims); 329 330 // tests if 2 boxes overlap within a tolerance 331 // assume that data is valid, box_min1 << box_max1, and box_min2<< box_max2 332 // they are overlapping if they are overlapping in all directions 333 // for example in x direction: 334 // box_max2[0] + tolerance < box_min1[0] -- false 335 bool boxes_overlap( const CartVect & box_min1, const CartVect & box_max1, 336 const CartVect & box_min2, const CartVect & box_max2, double tolerance); 337 338 // see if boxes formed by 2 lists of "CartVect"s overlap 339 bool bounding_boxes_overlap (const CartVect * list1, int num1, const CartVect * list2, int num2, 340 double tolerance); 341 342 // see if boxes from vertices in 2d overlap (used in gnomonic planes right now) 343 bool bounding_boxes_overlap_2d (const double * list1, int num1, const double * list2, int num2, 344 double tolerance); 345 // point_in_trilinear_hex 346 // Tests if a point in xyz space is within a hex element defined with 347 // its eight vertex points forming a trilinear basis function. Computes 348 // the natural coordinates with respect to the hex of the xyz point 349 // and checks if each are between +/-1. If anyone is outside the range 350 // the function returns false, otherwise it returns true. 351 // 352 bool point_in_trilinear_hex(const CartVect *hex, 353 const CartVect& xyz, 354 double etol); 355 356 bool nat_coords_trilinear_hex( const CartVect* corner_coords, 357 const CartVect& x, 358 CartVect& xi, 359 double tol ); 360 } // namespace GeomUtil 361 362 } // namespace moab 363 364 #endif 365