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 OrientedBox.hpp
17 *\author Jason Kraftcheck (kraftche@cae.wisc.edu)
18 *\date 2006-07-18
19 */
20
21 #ifndef MB_ORIENTED_BOX_HPP
22 #define MB_ORIENTED_BOX_HPP
23
24 #include "moab/Forward.hpp"
25 #include "moab/CartVect.hpp"
26 #include "moab/Matrix3.hpp"
27
28 #include <iosfwd>
29
30 namespace moab {
31
32 #define MB_ORIENTED_BOX_UNIT_VECTORS 1
33 #define MB_ORIENTED_BOX_OUTER_RADIUS 1
34
35 class Range;
36
37
38 /**\brief Oriented bounding box
39 */
40 class OrientedBox
41 {
42 private:
43 void order_axes_by_length( double ax1_len, double ax2_len, double ax3_len ); //!< orders the box axes by the given lengths for each axis
44
45 public:
46 CartVect center; //!< Box center
47 Matrix3 axes; //!< Box axes, unit vectors sorted by extent of box along axis
48 #if MB_ORIENTED_BOX_UNIT_VECTORS
49 CartVect length; //!< distance from center to plane along each axis
50 #endif
51 #if MB_ORIENTED_BOX_OUTER_RADIUS
52 double radius; //!< outer radius (1/2 diagonal length) of box
53 #endif
54
OrientedBox()55 inline OrientedBox() : radius(0.0) {}
56
57 OrientedBox( const Matrix3& axes_mat, const CartVect& center );
58 OrientedBox( const CartVect axes_in[3], const CartVect& center );
59
60 inline double inner_radius() const; //!< radius of inscribed sphere
61 inline double outer_radius() const; //!< radius of circumscribed sphere
62 inline double outer_radius_squared(const double reps) const; //!< square of (radius+at least epsilon) of circumsphere
63 inline double inner_radius_squared(const double reps) const; //!< square of (radius-epsilon) of inscribed sphere
64 inline double volume() const; //!< volume of box
65 inline CartVect dimensions() const; //!< number of dimensions for which box is not flat
66 inline double area() const; //!< largest side area
67 inline CartVect axis( int index ) const; //!< get unit vector in direction of axis
68 inline CartVect scaled_axis( int index ) const; //!< get vector in direction of axis, scaled to its true length
69
70 /** Test if point is contained in box */
71 bool contained( const CartVect& point, double tolerance ) const;
72
73 //bool contained( const OrientedBox& other, double tolerance ) const;
74
75 /**\brief get tag handle for storing oriented box
76 *
77 * Get the handle for the tag with the specified name and
78 * check that the tag is appropriate for storing instances
79 * of OrientedBox. The resulting tag may be used to store
80 * instances of OrientedBox directly.
81 *
82 *\param handle_out The TagHandle, passed back to caller
83 *\param name The tag name
84 *\param create If true, tag will be created if it does not exist
85 */
86 static ErrorCode tag_handle( Tag& handle_out,
87 Interface* instance,
88 const char* name);
89
90 /**\brief Calculate an oriented box from a set of vertices */
91 static ErrorCode compute_from_vertices( OrientedBox& result,
92 Interface* instance,
93 const Range& vertices );
94
95 /**\brief Calculate an oriented box from a set of 2D elements */
96 static ErrorCode compute_from_2d_cells( OrientedBox& result,
97 Interface* instance,
98 const Range& elements );
99
100 /** Structure to hold temporary accumulated triangle data for
101 * calculating box orientation. See box_from_covariance_data
102 * to see how this is used to calculate the final covariance matrix
103 * and resulting box orientation.
104 */
105 struct CovarienceData {
CovarienceDatamoab::OrientedBox::CovarienceData106 CovarienceData() : area(0.0) {}
CovarienceDatamoab::OrientedBox::CovarienceData107 CovarienceData( const Matrix3& m, const CartVect& c, double a)
108 : matrix(m), center(c), area(a) {}
109 Matrix3 matrix; //!< Running sum for covariance matrix
110 CartVect center; //!< Sum of triangle centroids weighted by 2*triangle area
111 double area; //!< 2x the sum of the triangle areas
112 };
113
114 /** Calculate a CovarienceData struct from a list of triangles */
115 static ErrorCode covariance_data_from_tris( CovarienceData& result,
116 Interface* moab_instance,
117 const Range& elements );
118
119 /** Calculate an OrientedBox given an array of CovarienceData and
120 * the list of vertices the box is to bound.
121 */
122 static ErrorCode compute_from_covariance_data( OrientedBox& result,
123 Interface* moab_instance,
124 const CovarienceData* orient_array,
125 unsigned orient_array_length,
126 const Range& vertices );
127
128 /** Test for intersection of a ray (or line segment) with this box.
129 * Ray length limits are used to optimize Monte Carlo particle tracking.
130 *\param ray_start_point The base point of the ray
131 *\param ray_unit_direction The direction of the ray (must be unit length)
132 *\param distance_tolerance Tolerance to use in intersection checks
133 *\param nonnegative_ray_len Optional length of ray in forward direction
134 *\param negative_ray_len Optional length of ray in reverse direction
135 */
136 bool intersect_ray( const CartVect& ray_start_point,
137 const CartVect& ray_unit_direction,
138 const double distance_tolerance,
139 const double* nonnegatve_ray_len = 0,
140 const double* negative_ray_len = 0 ) const;
141
142 /**\brief Find closest position on/within box to input position.
143 *
144 * Find the closest position in the solid box to the input position.
145 * If the input position is on or within the box, then the output
146 * position will be the same as the input position. If the input
147 * position is outside the box, the outside position will be the
148 * closest point on the box boundary to the input position.
149 */
150 void closest_location_in_box( const CartVect& input_position,
151 CartVect& output_position ) const;
152
153 //! Construct a hexahedral element with the same shape as this box.
154 ErrorCode make_hex( EntityHandle& hex, Interface* instance );
155
156
157 /** Calculate an OrientedBox given a CovarienceData struct and
158 * the list of points the box is to bound.
159 */
160 static ErrorCode compute_from_covariance_data( OrientedBox& result,
161 Interface* moab_instance,
162 CovarienceData& orientation_data,
163 const Range& vertices );
164 };
165
166 std::ostream& operator<<( std::ostream&, const OrientedBox& );
167
inner_radius() const168 double OrientedBox::inner_radius() const
169 {
170 #if MB_ORIENTED_BOX_UNIT_VECTORS
171 return length[0];
172 #else
173 return axes.col(0).length();
174 #endif
175 }
176
outer_radius() const177 double OrientedBox::outer_radius() const
178 {
179 #if MB_ORIENTED_BOX_OUTER_RADIUS
180 return radius;
181 #elif MB_ORIENTED_BOX_UNIT_VECTORS
182 return length.length();
183 #else
184 return (axes.col(0) + axes.col(1) + axes.col(2)).length();
185 #endif
186 }
187
188 // Add at least epsilon to the radius, before squaring it.
outer_radius_squared(const double reps) const189 double OrientedBox::outer_radius_squared(const double reps) const
190 {
191 #if MB_ORIENTED_BOX_OUTER_RADIUS
192 return (radius+reps)*(radius+reps);
193 #elif MB_ORIENTED_BOX_UNIT_VECTORS
194 CartVect tmp(length[0]+reps,length[1]+reps,length[2]+reps);
195 return tmp % tmp;
196 #else
197 CartVect half_diag = axes.col(0) + axes.col(1) + axes.col(2);
198 half_diag += CartVect(reps,reps,reps);
199 return half_diag % half_diag;
200 #endif
201 }
202
203 // Subtract epsilon from the length of the shortest axis, before squaring it.
inner_radius_squared(const double reps) const204 double OrientedBox::inner_radius_squared(const double reps) const
205 {
206 #if MB_ORIENTED_BOX_UNIT_VECTORS
207 return (length[0]-reps) * (length[0]-reps);
208 #else
209 CartVect tmp = axes.col(0);
210 tmp -= CartVect(reps,reps,reps);
211 return (tmp % tmp);
212 #endif
213 }
214
volume() const215 double OrientedBox::volume() const
216 {
217 #if MB_ORIENTED_BOX_UNIT_VECTORS
218 return 8 * length[0] * length[1] * length[2];
219 #else
220 return fabs(8 * axes.col(0) % (axes.col(1) * axes.col(2)));
221 #endif
222 }
223
dimensions() const224 CartVect OrientedBox::dimensions() const
225 {
226 #if MB_ORIENTED_BOX_UNIT_VECTORS
227 return 2.0 * length;
228 #else
229 return 2.0 * CartVect( axes.col(0).length(), axes.col(1).length(), axes.col(2).length() );
230 #endif
231 }
232
area() const233 double OrientedBox::area() const
234 {
235 #if MB_ORIENTED_BOX_UNIT_VECTORS
236 return 4 * length[1] * length[2];
237 #else
238 return 4 * (axes.col(1) * axes.col(2)).length();
239 #endif
240 }
241
axis(int index) const242 CartVect OrientedBox::axis( int index ) const
243 {
244 return axes.col(index);
245 }
246
scaled_axis(int index) const247 CartVect OrientedBox::scaled_axis( int index ) const
248 {
249 #if MB_ORIENTED_BOX_UNIT_VECTORS
250 return length[index] * axes.col(index);
251 #else
252 return axes.col(index);
253 #endif
254 }
255
256 } // namespace moab
257
258 #endif
259