Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
GeomQueryTool.hpp
Go to the documentation of this file.
1
#ifndef MOAB_GEOM_QUERY_TOOL_HPP
2
#define MOAB_GEOM_QUERY_TOOL_HPP
3
4
#ifdef _MSC_VER
/* windows */
5
#define _USE_MATH_DEFINES
// For M_PI
6
#endif
7
8
#include "
MBTagConventions.hpp
"
9
#include "
moab/CartVect.hpp
"
10
#include "
moab/Range.hpp
"
11
#include "
moab/Core.hpp
"
12
#include "
moab/GeomUtil.hpp
"
13
#include "
moab/FileOptions.hpp
"
14
#include "
moab/EntityHandle.hpp
"
15
#include "
moab/GeomTopoTool.hpp
"
16
#include "
moab/OrientedBoxTreeTool.hpp
"
17
18
#include <vector>
19
#include <map>
20
#include <string>
21
#include <cassert>
22
23
namespace
moab
24
{
25
26
/** \class GeomQueryTool
27
*
28
* \brief Tool for querying different aspects of geometric topology sets in MOAB
29
*
30
* Given the conventions established in GeomTopoTool for representing
31
* geometric topology through a hierarchy of meshsets, this tool provides a
32
* set of methods to query different geometric aspects of those geometric
33
* topology sets including:
34
*
35
* * measures of surface area and volume
36
* * distance to a bounding surface from a point within a volume
37
* * test the inclusion of a point within a volume
38
* * find the angle of a surface at a point
39
*
40
* A feature of these queries is that there is some support for overlapping
41
* geometries.
42
*
43
*/
44
45
class
GeomQueryTool
46
{
47
public
:
48
/* \class RayHistory
49
*
50
* In many use cases, it is useful to track some of the history of a ray as
51
* it passes through a geometry, particularly a geometry represented by
52
* facets. For example, given round-off erorr in ray-triangle tests
53
* (GeomUtil::ray_tri_intersect) used as part of a test for ray-surface
54
* intersection, it is possible for subsequent queries to oscillate between
55
* adjacent surfaces. This class stores information about history of a ray
56
* that can be used to test for such circumstances so that they can be
57
* accommodated.
58
*
59
*/
60
61
class
RayHistory
62
{
63
64
public
:
65
/**
66
* Clear this entire history-- logically equivalent to creating a new history,
67
* but probably more efficient.
68
*/
69
void
reset
();
70
71
/**
72
* Clear the history up to the most recent intersection. This should be
73
* called when a ray changes direction at the site of a surface crossing,
74
* a situation that most commonly occurs at a reflecting boundary.
75
*/
76
void
reset_to_last_intersection
();
77
78
/**
79
* Remove the most recent intersection. This allows a subsequent call
80
* along the same ray to return the same intersection.
81
*/
82
void
rollback_last_intersection
();
83
84
/**
85
* Get the last intersection in the RayHistory. This will return a null
86
* EntityHandle (0) if the history is empty.
87
*/
88
ErrorCode
get_last_intersection
(
EntityHandle
& last_facet_hit )
const
;
89
90
/**
91
* @return the number of surface crossings currently represented by this ray history
92
*/
93
int
size
()
const
94
{
95
return
prev_facets
.size();
96
}
97
98
/**
99
* @return Boolean indicating if this entity is in the RayHistory
100
*/
101
bool
in_history
(
EntityHandle
ent )
const
;
102
103
/**
104
* Add entity to the RayHistory
105
*/
106
void
add_entity
(
EntityHandle
ent );
107
108
private
:
109
std::vector< EntityHandle >
prev_facets
;
110
111
friend
class
GeomQueryTool
;
112
};
113
114
// Constructors
115
116
GeomQueryTool
(
Interface
* impl,
117
bool
find_geoments =
false
,
118
EntityHandle
modelRootSet = 0,
119
bool
p_rootSets_vector =
true
,
120
bool
restore_rootSets =
true
,
121
bool
trace_counting =
false
,
122
double
overlap_thickness = 0.,
123
double
numerical_precision = 0.001 );
124
125
GeomQueryTool
(
GeomTopoTool
* geomtopotool,
126
bool
trace_counting =
false
,
127
double
overlap_thickness = 0.,
128
double
numerical_precision = 0.001 );
129
130
// Destructor
131
~GeomQueryTool
();
132
133
ErrorCode
initialize
();
134
135
/**\brief find the next surface crossing from a given point in a given direction
136
*
137
* This is the primary method to enable ray tracing through a geometry.
138
* Given a volume and a ray, it determines the distance to the nearest intersection
139
* with a bounding surface of that volume and returns that distance and the
140
* EntityHandle indicating on which surface that intersection occurs.
141
* The caller can compute the location of the intersection by adding the
142
* distance to the ray.
143
*
144
* When a series of calls to this function are made along the same ray (e.g. for
145
* the purpose of tracking a ray through several volumes), the optional history
146
* argument should be given. The history prevents previously intersected facets
147
* from being intersected again. A single history should be used as long as a
148
* ray is proceeding forward without changing direction. This situation is
149
* sometimes referred to as "streaming."
150
*
151
* If a ray changes direction at an intersection site, the caller should call
152
* reset_to_last_intersection() on the history object before the next ray fire.
153
*
154
* @param volume The volume at which to fire the ray
155
* @param ray_start An array of x,y,z coordinates from which to start the ray.
156
* @param ray_dir An array of x,y,z coordinates indicating the direction of the ray.
157
* Must be of unit length.
158
* @param next_surf Output parameter indicating the next surface intersected by the ray.
159
* If no intersection is found, will be set to 0.
160
* @param next_surf_dist Output parameter indicating distance to next_surf. If next_surf is
161
* 0, this value is undefined and should not be used.
162
* @param history Optional RayHistory object. If provided, the facets in the history are
163
* assumed to not intersect with the given ray. The facet intersected
164
* by this query will also be added to the history.
165
* @param dist_limit Optional distance limit. If provided and > 0, no intersections at a
166
* distance further than this value will be returned.
167
* @param ray_orientation Optional ray orientation. If provided determines intersections
168
* along the normal provided, e.g. if -1 allows intersections back along the
169
* the ray direction, Default is 1, i.e. exit intersections
170
* @param stats Optional TrvStats object used to measure performance of underlying OBB
171
* ray-firing query. See OrientedBoxTreeTool.hpp for details.
172
*
173
*/
174
ErrorCode
ray_fire
(
const
EntityHandle
volume,
175
const
double
ray_start[3],
176
const
double
ray_dir[3],
177
EntityHandle
& next_surf,
178
double
& next_surf_dist,
179
RayHistory
* history = NULL,
180
double
dist_limit = 0,
181
int
ray_orientation = 1,
182
OrientedBoxTreeTool::TrvStats
* stats = NULL );
183
184
/**\brief Test if a point is inside or outside a volume
185
*
186
* This method finds the point on the boundary of the volume that is nearest
187
* the test point (x,y,z). If that point is "close" to a surface, a boundary test
188
* is performed based on the normal of the surface at that point and the
189
* optional ray direction (u,v,w).
190
* @param volume The volume to test
191
* @param xyz The location to test for volume containment
192
* @param result Set to 0 if xyz it outside volume, 1 if inside, and -1 if on boundary.
193
* @param Optional direction to use for underlying ray fire query. Used to ensure
194
* consistent results when a ray direction is known. If NULL or {0,0,0} is
195
* given, a random direction will be used.
196
* @param history Optional RayHistory object to pass to underlying ray fire query.
197
* The history is not modified by this call.
198
*/
199
ErrorCode
point_in_volume
(
const
EntityHandle
volume,
200
const
double
xyz[3],
201
int
& result,
202
const
double
* uvw = NULL,
203
const
RayHistory
* history = NULL );
204
205
/**\brief Robust test if a point is inside or outside a volume using unit sphere area method
206
*
207
* This test may be more robust that the standard point_in_volume, but is much slower.
208
* It does not detect 'on boundary' situations as point_in_volume does.
209
* @param volume The volume to test
210
* @param xyz The location to test for volume containment
211
* @param result Set to 0 if xyz it outside volume, 1 if inside.
212
*/
213
ErrorCode
point_in_volume_slow
(
const
EntityHandle
volume,
const
double
xyz[3],
int
& result );
214
215
/**\brief Find volume for a given location.
216
*
217
* Determines which volume contains the point if possible. If no volume is found,
218
* a null EntityHandle is returned along with a MB_ENTITY_NOT_FOUND ErrorCode.
219
* @param xyz The location to test
220
* @param volume Set to volume handle containing the location if found
221
* @param dir Optional direction to use for underlying ray fire query. Used to ensure
222
* consistent results when a ray direction is known. If NULL or {0,0,0} is
223
* given, a random direction will be used.
224
*/
225
ErrorCode
find_volume
(
const
double
xyz[3],
EntityHandle
& volume,
const
double
* dir = NULL );
226
227
/**\brief Find volume for a given location using loop. (slow)
228
*
229
* Loops over all volumes in the model, checking for point containment
230
* @param xyz The location to test
231
* @param volume Set to volume handle containing the location if found
232
* @param dir Optional direction to use for underlying ray fire query. Used to ensure
233
* consistent results when a ray direction is known. If NULL or {0,0,0} is
234
* given, a random direction will be used.
235
*/
236
ErrorCode
find_volume_slow
(
const
double
xyz[3],
EntityHandle
& volume,
const
double
* dir = NULL );
237
238
/**\brief Test if a point is inside or outsize a volume's axis-aligned bounding box
239
*
240
* This is used as a fast way to determine whether a point is inside or outside of a volume
241
* before performing a point_in_volume test which involves firing a ray.
242
* @param volume The volume to test
243
* @param point The location to test for bounding box containment
244
* @param inside set to 0 if point is outside the box, 1 if inside
245
*/
246
ErrorCode
point_in_box
(
const
EntityHandle
volume,
const
double
point[3],
int
& inside );
247
248
/** \brief Given a ray starting at a surface of a volume, check whether the ray enters or exits
249
* the volume
250
*
251
* This function is most useful for rays that change directions at a surface crossing.
252
* It can be used to check whether a direction change redirects the ray back into the
253
* originating volume.
254
*
255
* @param volume The volume to test
256
* @param surface A surface on volume
257
* @param xyz A point location on surface
258
* @param uvw A (unit) direction vector
259
* @param result Set to 1 if ray is entering volume, or 0 if it is leaving
260
* @param history Optional ray history object from a previous call to ray_fire. If present and
261
* non-empty, the history is used to look up the surface facet at which the ray begins. Absent
262
* a history, the facet nearest to xyz will be looked up. The history should always be provided
263
* if available, as it avoids the computational expense of a nearest-facet query.
264
*/
265
ErrorCode
test_volume_boundary
(
const
EntityHandle
volume,
266
const
EntityHandle
surface,
267
const
double
xyz[3],
268
const
double
uvw[3],
269
int
& result,
270
const
RayHistory
* history = NULL );
271
272
/**\brief Find the distance to the point on the boundary of the volume closest to the test point
273
*
274
* @param volume Volume to query
275
* @param point Coordinates of test point
276
* @param result Set to the minimum distance from point to a surface in volume
277
*/
278
ErrorCode
closest_to_location
(
EntityHandle
volume,
279
const
double
point[3],
280
double
& result,
281
EntityHandle
* closest_surface = 0 );
282
283
/** Calculate the volume contained in a 'volume' */
284
ErrorCode
measure_volume
(
EntityHandle
volume,
double
& result );
285
286
/** Calculate sum of area of triangles */
287
ErrorCode
measure_area
(
EntityHandle
surface,
double
& result );
288
289
/** Get the normal to a given surface at the point on the surface closest to a given point
290
*
291
* This method first identifies which facet contains this point and then
292
* calculates the unit outward normal of that facet. The facet of the
293
* provided volume that is nearest the provided point is used for this
294
* calculation. The search for that facet can be circumvented by providing
295
* a RayHistory, in which case the last facet of the history will be used.
296
*
297
* @param surf Surface on which to get normal
298
* @param xyz Point on surf
299
* @param angle Set to coordinates of surface normal nearest xyz
300
* @param history Optional ray history from a previous call to ray_fire().
301
* If present and non-empty, return the normal
302
* of the most recently intersected facet, ignoring xyz.
303
*/
304
ErrorCode
get_normal
(
EntityHandle
surf,
const
double
xyz[3],
double
angle
[3],
const
RayHistory
* history = NULL );
305
306
private
:
307
/**\brief determine the point membership when the point is effectively on the boundary
308
*
309
* Called by point_in_volume when the point is with tolerance of the boundary. Compares the
310
* ray direction with the surface normal to determine a volume membership.
311
*/
312
ErrorCode
boundary_case
(
EntityHandle
volume,
313
int
& result,
314
double
u,
315
double
v,
316
double
w,
317
EntityHandle
facet,
318
EntityHandle
surface );
319
320
/** get the solid angle projected by a facet on a unit sphere around a point
321
* - used by point_in_volume_slow
322
*/
323
ErrorCode
poly_solid_angle
(
EntityHandle
face
,
const
CartVect
& point,
double
& area );
324
325
/**\brief State object used in calls to ray_fire()
326
*
327
* Storage for the "history" of a ray. This represents the surface facets
328
* that the ray is known to have crossed, which cannot be crossed again
329
* as long as the ray does not change direction. It is intended to be used
330
* with a series of consecutive calls to ray_fire(), in which a ray passes
331
* over potentially many surfaces.
332
*/
333
334
public
:
335
/*
336
Overlap Thickness:
337
This tolerance is the maximum distance across an overlap. It should be zero
338
unless the geometry has overlaps. The overlap thickness is set using the dagmc
339
card. Overlaps must be small enough not to significantly affect physics.
340
Performance: increasing tolerance decreases performance
341
Robustness: increasing tolerance increases robustness
342
Knowledge: user must have intuition of overlap thickness
343
*/
344
345
/** Attempt to set a new overlap thickness tolerance, first checking for sanity */
346
347
void
set_overlap_thickness
(
double
new_overlap_thickness );
348
349
/*
350
Numerical Precision:
351
This tolerance is used for obb.intersect_ray, finding neighborhood of
352
adjacent triangle for edge/node intersections, and error in advancing
353
geometric position of particle (x' ~= x + d*u). When determining the
354
neighborhood of adjacent triangles for edge/node intersections, the facet
355
based model is expected to be watertight.
356
Performance: increasing tolerance decreases performance (but not very much)
357
Robustness: increasing tolerance increases robustness
358
Knowledge: user should not change this tolerance
359
*/
360
361
/** Attempt to set a new numerical precision , first checking for sanity
362
* Use of this function is discouraged.
363
*/
364
void
set_numerical_precision
(
double
new_precision );
365
366
void
set_verbosity
(
bool
value )
367
{
368
verbose
= value;
369
}
370
371
bool
get_verbosity
()
const
372
{
373
return
verbose
;
374
}
375
376
double
get_numerical_precision
()
377
{
378
return
numericalPrecision
;
379
}
380
381
double
get_overlap_thickness
()
382
{
383
return
overlapThickness
;
384
}
385
386
GeomTopoTool
*
gttool
()
387
{
388
return
geomTopoTool
;
389
}
390
391
Interface
*
moab_instance
()
392
{
393
return
MBI
;
394
}
395
396
private
:
397
GeomTopoTool
*
geomTopoTool
;
398
bool
verbose
;
399
bool
owns_gtt
;
400
Interface
*
MBI
;
401
OrientedBoxTreeTool
*
obbTreeTool
;
402
bool
counting
;
403
long
long
int
n_pt_in_vol_calls
;
404
long
long
int
n_ray_fire_calls
;
405
double
overlapThickness
,
numericalPrecision
;
406
Tag
senseTag
;
407
};
408
409
}
// namespace moab
410
411
#endif
src
moab
GeomQueryTool.hpp
Generated on Tue Jul 30 2024 02:18:01 for Mesh Oriented datABase by
1.9.1.