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