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