Get smaller neighbors

dccrg.hpp

@@ -2923,6 +2923,11 @@ template <
 		return this->neighborhood_length;
 	}
 
+	// round away from zero, stackoverflow.com/a/2745086
+	static int divide_away_from_zero(const int x, const int y)
+	{
+		return x > 0 ? (x + y - 1) / y : (x - y + 1) / y;
+	}
 
 	/*!
 	Returns all neighbors of given cell that are at given offset from it.
@@ -2961,24 +2966,15 @@ template <
 		for (const auto& neighbor_i: this->neighbors_of.at(cell)) {
 			const auto& offsets = neighbor_i.second;
 			if (offsets[3] <= 1) {
-				if (offsets[0] == x and offsets[1] == y and offsets[2] == z) {
+				if (offsets[0] == x && offsets[1] == y && offsets[2] == z) {
 					return_neighbors.push_back(neighbor_i);
 				}
 			} else {
 				auto scaled_offsets = offsets;
 				for (size_t i = 0; i < 3; i++) {
-					if (std::abs(scaled_offsets[i]) <= 1) {
-						continue;
-					}
-					// round away from zero, stackoverflow.com/a/2745086
-					if (scaled_offsets[i] > 1) {
-						scaled_offsets[i] += scaled_offsets[3] - 1;
-					} else {
-						scaled_offsets[i] -= scaled_offsets[3] - 1;
-					}
-					scaled_offsets[i] /= scaled_offsets[3];
+					scaled_offsets[i] = divide_away_from_zero(scaled_offsets[i], scaled_offsets[3]);
 				}
-				if (scaled_offsets[0] == x and scaled_offsets[1] == y and scaled_offsets[2] == z) {
+				if (scaled_offsets[0] == x && scaled_offsets[1] == y && scaled_offsets[2] == z) {
 					return_neighbors.push_back(neighbor_i);
 				}
 			}
@@ -2987,6 +2983,69 @@ template <
 		return return_neighbors;
 	}
 
+	// Doesn't support denominator < 0
+	std::vector<
+		std::pair<
+			uint64_t,
+			std::array<int, 4>
+		>
+	> get_neighbors_of_at_offset(
+		const uint64_t cell,
+		const int x,
+		const int y,
+		const int z,
+		const int denom
+	) const {
+		if (denom == 1) {
+			return get_neighbors_of_at_offset(cell, x, y, z);
+		}
+
+		std::vector<std::pair<uint64_t, std::array<int, 4>>> return_neighbors;
+		if (
+			this->cell_process.count(cell) == 0 || 
+			this->cell_process.at(cell) != this->rank || 
+			(x == 0 and y == 0 and z == 0) || 
+			denom < 0
+		) {
+			return return_neighbors;
+		}
+
+		std::vector<int> xyz = {x, y, z};
+		for (size_t i = 0; i < 3; i++) {
+			xyz[i] = divide_away_from_zero(xyz[i], denom);
+		}
+
+		for (const auto& neighbor_i : this->get_neighbors_of_at_offset(cell, xyz[0], xyz[1], xyz[2])) {
+			const auto& offsets = neighbor_i.second;
+			if (offsets[3] == denom) {
+				if (offsets[0] == x && offsets[1] == y && offsets[2] == z) {
+					return_neighbors.push_back(neighbor_i);
+				}
+			} else if (offsets[3] < denom) {
+				// If we're looking at cells larger than denom, calculate our desired xyz offset in that size
+				xyz = {x, y, z};
+				int factor = offsets[3] > 0 ? denom/offsets[3] : -offsets[3] * denom;
+				for (size_t i = 0; i < 3; i++) {
+					xyz[i] = divide_away_from_zero(xyz[i], factor);
+				}
+
+				if (offsets[0] == xyz[0] && offsets[1] == xyz[1] && offsets[2] == xyz[2]) {
+					return_neighbors.push_back(neighbor_i);
+				}
+			} else { // offsets[3] > denom
+				// If we're looking at cells smaller than denom, calculate their offsets in denom size
+				auto scaled_offsets = offsets;
+				for (size_t i = 0; i < 3; i++) {
+					scaled_offsets[i] = divide_away_from_zero(scaled_offsets[i], scaled_offsets[3] / denom);
+				}
+				if (scaled_offsets[0] == x && scaled_offsets[1] == y && scaled_offsets[2] == z) {
+					return_neighbors.push_back(neighbor_i);
+				}
+			}
+		}
+
+		return return_neighbors;
+	}
 
 	/*!
 	Returns neighbors of given local cell that are on another process.