1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214 | #include "moab/Core.hpp"
#include "moab/SpatialLocator.hpp"
#include "moab/Tree.hpp"
#include "moab/HomXform.hpp"
#include "moab/ScdInterface.hpp"
#include "moab/CartVect.hpp"
#include "moab/AdaptiveKDTree.hpp"
#include "moab/BVHTree.hpp"
#include "moab/ProgOptions.hpp"
#include "moab/CpuTimer.hpp"
#include "moab/ElemEvaluator.hpp"
#ifdef MOAB_HAVE_MPI
#include "moab_mpi.h"
#endif
#include <cstdlib>
#include <sstream>
using namespace moab;
ErrorCode test_locator( SpatialLocator& sl, int npoints, double rtol, double& cpu_time, double& percent_outside );
ErrorCode create_hex_mesh( Interface& mb, Range& elems, int n, int dim );
int main( int argc, char** argv )
{
#ifdef MOAB_HAVE_MPI
int fail = MPI_Init( &argc, &argv );
if( fail ) return fail;
#else
// silence the warning of parameters not used, in serial; there should be a smarter way :(
argv[0] = argv[argc - argc];<--- Same expression on both sides of '-'. [+]Finding the same expression on both sides of an operator is suspicious and might indicate a cut and paste or logic error. Please examine this code carefully to determine if it is correct.
#endif
int npoints = 100, dim = 3;
int dints = 1, dleafs = 1, ddeps = 1;
bool eval = false;
double rtol = 1.0e-10;
ProgOptions po( "tree_searching_perf options" );
po.addOpt< void >( ",e", "Use ElemEvaluator in tree search", &eval );
po.addOpt< int >( "ints,i", "Number of doublings of intervals on each side of scd mesh", &dints );
po.addOpt< int >( "leaf,l", "Number of doublings of maximum number of elements per leaf", &dleafs );
po.addOpt< int >( "max_depth,m", "Number of 5-intervals on maximum depth of tree", &ddeps );
po.addOpt< int >( "npoints,n", "Number of query points", &npoints );
po.addOpt< int >( "dim,d", "Dimension of the mesh", &dim );
po.addOpt< double >( "tol,t", "Relative tolerance of point search", &rtol );
// po.addOpt<void>( "print,p", "Print tree details", &print_tree);
po.parseCommandLine( argc, argv );
std::vector< int > ints, deps, leafs;
ints.push_back( 10 );
for( int i = 1; i < dints; i++ )
ints.push_back( 2 * ints[i - 1] );
deps.push_back( 30 );
for( int i = 1; i < ddeps; i++ )
deps.push_back( deps[i - 1] - 5 );
leafs.push_back( 6 );
for( int i = 1; i < dleafs; i++ )
leafs.push_back( 2 * leafs[i - 1] );
ErrorCode rval = MB_SUCCESS;<--- Variable 'rval' is assigned a value that is never used.
std::cout << "Tree_type"
<< " "
<< "Elems_per_leaf"
<< " "
<< "Tree_depth"
<< " "
<< "Ints_per_side"
<< " "
<< "N_elements"
<< " "
<< "search_time"
<< " "
<< "perc_outside"
<< " "
<< "initTime"
<< " "
<< "nodesVisited"
<< " "
<< "leavesVisited"
<< " "
<< "numTraversals"
<< " "
<< "leafObjectTests" << std::endl;
// outermost iteration: # elements
for( std::vector< int >::iterator int_it = ints.begin(); int_it != ints.end(); ++int_it )
{
Core mb;
Range elems;
rval = create_hex_mesh( mb, elems, *int_it, dim );
if( MB_SUCCESS != rval ) return rval;
// iteration: tree depth
for( std::vector< int >::iterator dep_it = deps.begin(); dep_it != deps.end(); ++dep_it )
{
// iteration: tree max elems/leaf
for( std::vector< int >::iterator leafs_it = leafs.begin(); leafs_it != leafs.end(); ++leafs_it )
{
// iteration: tree type
for( int tree_tp = 0; tree_tp < 2; tree_tp++ )
{
// create tree
Tree* tree;
if( 0 == tree_tp )
tree = new BVHTree( &mb );
else
tree = new AdaptiveKDTree( &mb );
ElemEvaluator* eeval = NULL;
if( eval )
{
// create an element evaluator
eeval = new ElemEvaluator( &mb );
rval = eeval->set_eval_set( *elems.begin() );
if( MB_SUCCESS != rval ) return rval;
}
std::ostringstream opts;
opts << "MAX_DEPTH=" << *dep_it << ";MAX_PER_LEAF=" << *leafs_it;
FileOptions fo( opts.str().c_str() );
rval = tree->parse_options( fo );
if( MB_SUCCESS != rval ) return rval;
SpatialLocator sl( &mb, elems, tree, eeval );
// call evaluation
double cpu_time, perc_outside;
rval = test_locator( sl, npoints, rtol, cpu_time, perc_outside );
if( MB_SUCCESS != rval ) return rval;
std::cout << ( tree_tp == 0 ? "BVH" : "KD" ) << " " << *leafs_it << " " << *dep_it << " " << *int_it
<< " " << ( *int_it ) * ( *int_it ) * ( *int_it ) << " " << cpu_time << " "
<< perc_outside << " ";
tree->tree_stats().output_all_stats();
if( eeval ) delete eeval;
} // tree_tp
} // max elems/leaf
} // max depth
} // # elements
#ifdef MOAB_HAVE_MPI
fail = MPI_Finalize();
if( fail ) return fail;
#endif
return 0;
}
ErrorCode test_locator( SpatialLocator& sl, int npoints, double rtol, double& cpu_time, double& percent_outside )
{
BoundBox box = sl.local_box();
CartVect box_del = box.bMax - box.bMin;
std::vector< CartVect > test_pts( npoints ), test_res( npoints );
std::vector< EntityHandle > ents( npoints );
int* is_in = new int[npoints];
double denom = 1.0 / (double)RAND_MAX;
for( int i = 0; i < npoints; i++ )
{
// generate a small number of random point to test
double rx = (double)rand() * denom, ry = (double)rand() * denom, rz = (double)rand() * denom;
test_pts[i] = box.bMin + CartVect( rx * box_del[0], ry * box_del[1], rz * box_del[2] );
}
CpuTimer ct;
// call spatial locator to locate points
ErrorCode rval =
sl.locate_points( test_pts[0].array(), npoints, &ents[0], test_res[0].array(), &is_in[0], rtol, 0.0 );
if( MB_SUCCESS != rval )
{
delete[] is_in;
return rval;
}
cpu_time = ct.time_elapsed();
int num_out = std::count( is_in, is_in + npoints, false );
percent_outside = ( (double)num_out ) / npoints;
delete[] is_in;
return rval;
}
ErrorCode create_hex_mesh( Interface& mb, Range& elems, int n, int dim )
{
ScdInterface* scdi;
ErrorCode rval = mb.query_interface( scdi );
if( MB_SUCCESS != rval ) return rval;
HomCoord high( n - 1, -1, -1 );
if( dim > 1 ) high[1] = n - 1;
if( dim > 2 ) high[2] = n - 1;
ScdBox* new_box;
rval = scdi->construct_box( HomCoord( 0, 0, 0 ), high, NULL, 0, new_box );
if( MB_SUCCESS != rval ) return rval;
rval = mb.release_interface( scdi );
if( MB_SUCCESS != rval ) return rval;
rval = mb.get_entities_by_dimension( 0, dim, elems );
if( MB_SUCCESS != rval ) return rval;
return rval;
}
|