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