1 #ifndef EXAMPLES_MPI_DOMAIN_PARTITION_HPP
2 #define EXAMPLES_MPI_DOMAIN_PARTITION_HPP
3 
4 #include <vector>
5 #include <utility>
6 
7 #include <boost/geometry.hpp>
8 #include <boost/geometry/geometries/adapted/boost_array.hpp>
9 #include <boost/geometry/geometries/box.hpp>
10 #include <boost/geometry/index/rtree.hpp>
11 
12 BOOST_GEOMETRY_REGISTER_BOOST_ARRAY_CS(cs::cartesian)
13 
14 template <int NDIM>
15 class domain_partition {
16     public:
17         typedef boost::array<ptrdiff_t, NDIM>      point;
18         typedef boost::geometry::model::box<point> box;
19         typedef std::pair<box, int>                process;
20 
domain_partition(point lo,point hi,int num_processes)21         domain_partition(point lo, point hi, int num_processes) {
22             split(box(lo, hi), num_processes);
23 
24             for(int i = 0; i < num_processes; ++i)
25                 rtree.insert( std::make_pair(subdomains[i], i) );
26         }
27 
index(point p) const28         std::pair<int, ptrdiff_t> index(point p) const {
29             namespace bgi = boost::geometry::index;
30 
31             for(const process &v : rtree | bgi::adaptors::queried(bgi::intersects(p)) )
32             {
33                 return std::make_pair(v.second, local_index(v.first, p));
34             }
35 
36             // Unreachable:
37             return std::make_pair(0, 0l);
38         }
39 
size(size_t process) const40         size_t size(size_t process) const {
41             if (process >= subdomains.size()) return 0;
42 
43             point lo = subdomains[process].min_corner();
44             point hi = subdomains[process].max_corner();
45 
46             size_t v = 1;
47 
48             for(int i = 0; i < NDIM; ++i)
49                 v *= hi[i] - lo[i] + 1;
50 
51             return v;
52         }
53 
domain(size_t process) const54         box domain(size_t process) const {
55             if (process < subdomains.size())
56                 return subdomains[process];
57             else {
58                 point lo;
59                 point hi;
60                 for(int i = 0; i < NDIM; ++i) {
61                     lo[i] = 0;
62                     hi[i] = -1;
63                 }
64                 return box(lo, hi);
65             }
66         }
67     private:
68         std::vector<box> subdomains;
69 
70         boost::geometry::index::rtree<
71             process,
72             boost::geometry::index::quadratic<16>
73             > rtree;
74 
local_index(box domain,point p)75         static ptrdiff_t local_index(box domain, point p) {
76             point lo = domain.min_corner();
77             point hi = domain.max_corner();
78 
79             ptrdiff_t stride = 1, idx = 0;
80             for(int i = 0; i < NDIM; ++i) {
81                 idx += (p[i] - lo[i]) * stride;
82                 stride *= hi[i] - lo[i] + 1;
83             }
84 
85             return idx;
86         }
87 
split(box domain,int np)88         void split(box domain, int np) {
89             if (np == 1) {
90                 subdomains.push_back(domain);
91                 return;
92             }
93 
94             point lo = domain.min_corner();
95             point hi = domain.max_corner();
96 
97             // Get longest dimension of the domain
98             int wd = 0;
99             for(int i = 1; i < NDIM; ++i)
100                 if (hi[i] - lo[i] > hi[wd] - lo[wd]) wd = i;
101 
102             ptrdiff_t mid = lo[wd] + (hi[wd] - lo[wd]) * (np / 2) / np;
103 
104             box sd1 = domain;
105             box sd2 = domain;
106 
107             sd1.max_corner()[wd] = mid;
108             sd2.min_corner()[wd] = mid + 1;
109 
110             split(sd1, np / 2);
111             split(sd2, np - np / 2);
112         }
113 };
114 
115 #endif
116