1/*
2 * ORB.hpp
3 *
4 * Created on: Mar 13, 2015
5 * Author: Pietro Incardona
6 */
7
8#ifndef ORB_HPP_
9#define ORB_HPP_
10
11#include "util/mathutil.hpp"
12
13/*! \brief this class is a functor for "for_each" algorithm
14 *
15 * This class is a functor for "for_each" algorithm. For each
16 * element of a range_c the operator() is called.
17 * Is mainly used to unroll the bisect cycle inside one ORB cycle
18 *
19 */
20
21template<unsigned int dim, typename ORB>
22struct bisect_unroll
23{
24 ORB & orb;
25
26 /*! \brief constructor
27 *
28 */
29 bisect_unroll(ORB & orb)
30 :orb(orb)
31 {};
32
33 //! It call the copy function for each property
34 template<typename T>
35 void operator()(T& t)
36 {
37 orb.template bisect<T::value>();
38 }
39};
40
41/*! \brief This structure use SFINAE to avoid instantiation of invalid code
42 *
43 * birect_unroll cannot be performed when i > dim (basically you will never cut
44 * on a dimension that does not exist)
45 *
46 * What is happening here is that at compile-time enable_if< (i < dim) >::type
47 * is evaluated, if the condition is true enable_if<true>::type is a valid expression
48 * the the second specialization will be used, if is false enable_if<false>::type
49 * is not a valid expression, and the first specialization is chosen
50 *
51 */
52
53template<unsigned int dim, unsigned int i, typename ORB, class Enable=void>
54struct do_when_dim_gr_i
55{
56 static void bisect_loop(bisect_unroll<dim,ORB> & bu)
57 {
58 }
59};
60
61template<unsigned int dim, unsigned int i, typename ORB>
62struct do_when_dim_gr_i<dim,i,ORB,typename boost::enable_if< boost::mpl::bool_<(i < dim)> >::type>
63{
64 static void bisect_loop(bisect_unroll<dim,ORB> & bu, typename boost::enable_if< boost::mpl::bool_<(i < dim)> >::type* dummy = 0)
65 {
66 boost::mpl::for_each< boost::mpl::range_c<int,0,i> >(bu);
67 }
68};
69
70
71/*! \brief ORB node
72 * \tparam type of space float, double, complex ...
73 *
74 * it store the Center of mass of the local calculated particles
75 * on one direction (only one direction, mean that one scalar is enough)
76 *
77 */
78
79template<typename T> class ORB_node : public aggregate<T>
80{
81public:
82
83 static const unsigned int CM = 0;
84};
85
86/*! \brief This class implement orthogonal recursive bisection
87 *
88 * \tparam dim Dimensionality of the ORB
89 * \tparam T type of the space
90 * \tparam loc_wg type of structure that store the local weight of particles
91 * \tparam loc_pos type of structure that store the local position of particles
92 * \tparam Box type of structure that contain the domain extension
93 * \tparam Tree type of structure that store the tree structure
94 *
95 */
96
97template<unsigned int dim, typename T, typename loc_wg=openfpm::vector<float>, typename loc_pos=openfpm::vector<Point<dim,T>> , typename Box=Box<dim,T>, template<typename,typename> class Tree=Graph_CSR_s>
98class ORB
99{
100 // Virtual cluster
101 Vcluster<> & v_cl;
102
103 // particle coordinate accumulator
104 openfpm::vector<T> cm;
105 // number of elements summed into cm
106 openfpm::vector<size_t> cm_cnt;
107
108 // Local particle vector
109 loc_pos & lp;
110
111 // Label id of the particle, the label identify where the "local" particles
112 // are in the local decomposition
113 openfpm::vector<size_t> lp_lbl;
114
115 size_t dim_div;
116
117 // Structure that store the orthogonal bisection tree
118 Tree<ORB_node<T>,no_edge> grp;
119
120 /*! \brief Calculate the local center of mass on direction dir
121 *
122 * WARNING: with CM we mean the sum of the particle coordinate over one direction
123 *
124 * \tparam dir Direction witch to calculate the center of mass
125 *
126 * \param start from where the last leafs start
127 */
128 template<unsigned int dir> void local_cm(size_t start)
129 {
130 typedef Point<dim,T> s;
131
132 // reset the counters and accumulators
133 cm.fill(0);
134 cm_cnt.fill(0);
135
136 // Calculate the local CM
137 auto it = lp.getIterator();
138
139 while (it.isNext())
140 {
141 // vector key
142 auto key = it.get();
143
144 size_t lbl = lp_lbl.get(key) - start;
145
146 // add the particle coordinate to the CM accumulator
147 cm.get(lbl) += lp.template get<s::x>(key)[dir];
148
149 // increase the counter
150 cm_cnt.get(lbl)++;
151
152 ++it;
153 }
154 }
155
156 /*! \brief It label the particles
157 *
158 * It identify where the particles are
159 *
160 * \param n level
161 *
162 */
163
164 template<unsigned int dir> inline void local_label()
165 {
166 typedef Point<dim,T> s;
167
168 // direction of the previous split
169 const size_t dir_cm = (dir == 0)?(dim-1):(dir-1);
170
171 // if it is the first iteration we just have to create the labels array with zero
172 if (grp.getNVertex() == 1)
173 {lp_lbl.resize(lp.size());lp_lbl.fill(0); return;}
174
175 // we check where (the node) the particles live
176
177 auto p_it = lp.getIterator();
178
179 while(p_it.isNext())
180 {
181 auto key = p_it.get();
182
183 // Get the label
184 size_t lbl = lp_lbl.get(key);
185
186 // each node has two child's (FOR SURE)
187 // No leaf are processed here
188 size_t n1 = grp.getChild(lbl,0);
189 size_t n2 = grp.getChild(lbl,1);
190
191 // get the splitting center of mass
192 T cm = grp.template vertex_p<ORB_node<T>::CM>(lbl);
193
194 // if the particle n-coordinate is smaller than the CM is inside the child n1
195 // otherwise is inside the child n2
196
197 if (lp.template get<s::x>(key)[dir_cm] < cm)
198 {lp_lbl.get(key) = n1;}
199 else
200 {lp_lbl.get(key) = n2;}
201
202 ++p_it;
203 }
204 }
205
206 /*! \brief Bisect the domains along one direction
207 *
208 * \tparam dir direction where to split
209 *
210 */
211
212 template<unsigned int dir> size_t bisect()
213 {
214 //
215 size_t start = grp.getNVertex();
216
217 // first label the local particles
218 local_label<dir>();
219
220 // Index from where start the first leaf
221 size_t n_node = (start + 1) / 2;
222
223 // Calculate the local cm
224 local_cm<dir>(start - n_node);
225
226 // reduce the local cm and cm_cnt
227 v_cl.sum(cm);
228 v_cl.sum(cm_cnt);
229 v_cl.execute();
230
231 // set the CM for the previous leaf (they are not anymore leaf)
232
233 for (size_t i = 0 ; i < n_node ; i++)
234 {
235 grp.template vertex_p<ORB_node<T>::CM>(start-n_node+i) = cm.get(i) / cm_cnt.get(i);
236 }
237
238 // append on the 2^(n-1) previous end node to the 2^n leaf
239
240 for (size_t i = start - n_node, s = 0 ; i < start ; i++, s++)
241 {
242 // Add 2 leaf nodes and connect them with the node
243 grp.addVertex();
244 grp.addVertex();
245 grp.template addEdge(i,start+2*s);
246 grp.template addEdge(i,start+2*s+1);
247 }
248
249 return 2*n_node;
250 }
251
252 // define the friend class bisect_unroll
253
254 friend class bisect_unroll<dim,ORB<dim,T,loc_wg,loc_pos,Box,Tree>>;
255
256public:
257
258 /*! \brief constructor
259 *
260 * \param dom Box domain
261 * \param n_sub number of sub-domain to create (it is approximated to the biggest 2^n number)
262 * \param lp Local position of the particles
263 *
264 */
265 ORB(Box dom, size_t n_sub, loc_pos & lp)
266 :v_cl(create_vcluster()),lp(lp)
267 {
268 typedef ORB<dim,T,loc_wg,loc_pos,Box,Tree> ORB_class;
269
270 dim_div = 0;
271
272 n_sub = openfpm::math::round_big_2(n_sub);
273 size_t nsub = log2(n_sub);
274
275 // number of center or mass needed
276 cm.resize(2 << nsub);
277 cm_cnt.resize(2 << nsub);
278
279 // Every time we divide from 0 to dim-1, calculate how many cycle from 0 to dim-1 we have
280 size_t dim_cycle = nsub / dim;
281
282 // Create the root node in the graph
283 grp.addVertex();
284
285 // unroll bisection cycle
286 bisect_unroll<dim,ORB_class> bu(*this);
287 for (size_t i = 0 ; i < dim_cycle ; i++)
288 {
289 boost::mpl::for_each< boost::mpl::range_c<int,0,dim> >(bu);
290 // bu is recreated several time internaly
291 }
292
293 // calculate and execute the remaining cycles
294 switch(nsub - dim_cycle * dim)
295 {
296 case 1:
297 do_when_dim_gr_i<dim,1,ORB_class>::bisect_loop(bu);
298 break;
299 case 2:
300 do_when_dim_gr_i<dim,2,ORB_class>::bisect_loop(bu);
301 break;
302 case 3:
303 do_when_dim_gr_i<dim,3,ORB_class>::bisect_loop(bu);
304 break;
305 case 4:
306 do_when_dim_gr_i<dim,4,ORB_class>::bisect_loop(bu);
307 break;
308 case 5:
309 do_when_dim_gr_i<dim,5,ORB_class>::bisect_loop(bu);
310 break;
311 case 6:
312 do_when_dim_gr_i<dim,6,ORB_class>::bisect_loop(bu);
313 break;
314 case 7:
315 do_when_dim_gr_i<dim,7,ORB_class>::bisect_loop(bu);
316 break;
317 default:
318 // To extend on higher dimension create other cases or create runtime
319 // version of local_cm and bisect
320 std::cerr << "Error: " << __FILE__ << ":" << __LINE__ << " ORB is not working for dimension bigger than 8";
321
322 }
323 }
324};
325
326
327#endif /* ORB_HPP_ */
328