From 0c18b7515c76693de1e3d45538bb04481d402f1b Mon Sep 17 00:00:00 2001 From: Thomas Steinreiter Date: Tue, 22 Nov 2016 17:13:28 +0100 Subject: [PATCH] replaced std::array with Vec3 --- bhtree_mpi/src/datastructures/Body.cpp | 33 ++-- bhtree_mpi/src/datastructures/Body.hpp | 17 +- bhtree_mpi/src/datastructures/Box.cpp | 181 ++++++-------------- bhtree_mpi/src/datastructures/Box.hpp | 15 +- bhtree_mpi/src/datastructures/Vec3.hpp | 166 ++++++++++++++++++ bhtree_mpi/src/simulation/MpiSimulation.cpp | 19 +- bhtree_mpi/src/simulation/MpiSimulation.hpp | 1 + 7 files changed, 262 insertions(+), 170 deletions(-) create mode 100644 bhtree_mpi/src/datastructures/Vec3.hpp diff --git a/bhtree_mpi/src/datastructures/Body.cpp b/bhtree_mpi/src/datastructures/Body.cpp index 9aa7971..73a9e35 100644 --- a/bhtree_mpi/src/datastructures/Body.cpp +++ b/bhtree_mpi/src/datastructures/Body.cpp @@ -7,9 +7,9 @@ #include namespace nbody { - Body::Body(const std::array& position_, - const std::array& velocity_, - const std::array& acceleration_, + Body::Body(const Vec3& position_, + const Vec3& velocity_, + const Vec3& acceleration_, double mass_, bool refinement_, std::size_t id_) : @@ -21,20 +21,15 @@ namespace nbody { refinement(refinement_) {} void Body::resetAcceleration() { - std::fill(std::begin(acceleration), std::end(acceleration), 0.0); + acceleration = Vec3{ 0.0 }; } //helper method for intergration Derivative Body::evaluate(double dt, const Derivative& d) const { - std::array tmp{}; - Derivative output; - - for (std::size_t i = 0; i < 3; i++) { - tmp[i] = velocity[i] + d.dv[i] * dt; - output.dx[i] = tmp[i]; - output.dv[i] = acceleration[i]; - } - return output; + return Derivative{ + velocity + d.dv * dt, + acceleration + }; } //integrating acceleration -> velocity -> position @@ -81,9 +76,9 @@ namespace nbody { } bool Body::isValid() const { - if (!std::all_of(std::begin(position), std::end(position), validDouble)) return false; - if (!std::all_of(std::begin(velocity), std::end(velocity), validDouble)) return false; - if (!std::all_of(std::begin(acceleration), std::end(acceleration), validDouble)) return false; + if (!position.IsValid()) return false; + if (!velocity.IsValid()) return false; + if (!acceleration.IsValid()) return false; if (!validDouble(mass)) return false; return mass >= 0.0; } @@ -111,9 +106,9 @@ namespace nbody { vx = vy = vz = 0.0; iss >> mass >> px >> py >> pz >> vx >> vy >> vz; Body b{ - {{px, py, pz}}, - {{vx, vy, vz}}, - {{0.0, 0.0, 0.0}}, + {px, py, pz}, + {vx, vy, vz}, + {0.0, 0.0, 0.0}, mass, false, id++ diff --git a/bhtree_mpi/src/datastructures/Body.hpp b/bhtree_mpi/src/datastructures/Body.hpp index d1e26b1..b28c4ae 100644 --- a/bhtree_mpi/src/datastructures/Body.hpp +++ b/bhtree_mpi/src/datastructures/Body.hpp @@ -4,28 +4,29 @@ #include #include #include +#include namespace nbody { static const double timestep = 1.0; struct Derivative { - std::array dx{}; - std::array dv{}; + Vec3 dx{}; + Vec3 dv{}; }; struct Body { Body() = default; - Body(const std::array& position_, - const std::array& velocity_, - const std::array& acceleration_, + Body(const Vec3& position_, + const Vec3& velocity_, + const Vec3& acceleration_, double mass_, bool refinement_, std::size_t id_); std::size_t id{0}; - std::array position{}; - std::array velocity{}; - std::array acceleration{}; + Vec3 position{}; + Vec3 velocity{}; + Vec3 acceleration{}; double mass{0.0}; bool refinement{false}; diff --git a/bhtree_mpi/src/datastructures/Box.cpp b/bhtree_mpi/src/datastructures/Box.cpp index 2181457..24b09a3 100644 --- a/bhtree_mpi/src/datastructures/Box.cpp +++ b/bhtree_mpi/src/datastructures/Box.cpp @@ -1,58 +1,37 @@ #include #include #include +#include #include "Box.hpp" namespace nbody { //extend box to form cube void Box::extendToCube() { - std::size_t longestSide = std::numeric_limits::max(); - double sidelength = 0.0; + const auto maxL = maxSidelength(); + const auto diag = max - min; + const auto middle = min + (diag / 2); - for (std::size_t i = 0; i < 3; i++) { - if (max[i] - min[i] >= sidelength) { - longestSide = i; - sidelength = max[i] - min[i]; - } - } - if (longestSide == std::numeric_limits::max()) { - return; - } - for (std::size_t i = 0; i < 3; i++) { - if (i != longestSide) { - double extend = (sidelength - (max[i] - min[i])) / 2.0; + const auto normalizedDiag = diag / Norm(diag); - min[i] -= extend; - max[i] += extend; - } - } + const auto halfCubeDiag = normalizedDiag * (maxL / 2); + + min = middle - halfCubeDiag; + max = middle + halfCubeDiag; } //extend for bodies void Box::extendForBodies(const std::vector& bodies) { for (const auto& body : bodies) { - for (std::size_t i = 0; i < 3; i++) { - min[i] = std::min(body.position[i], min[i]); - max[i] = std::max(body.position[i], max[i]); - } + min.MinComponents(body.position); + max.MaxComponents(body.position); } } //extract bodies within box std::vector Box::extractBodies(std::vector& bodies) const { - std::vector result; - auto it = std::begin(bodies); - - while (it != std::end(bodies)) { - if (it->position[0] >= min[0] && it->position[0] <= max[0] && - it->position[1] >= min[1] && it->position[1] <= max[1] && - it->position[2] >= min[2] && it->position[2] <= max[2]) { - result.push_back(*it); - it = bodies.erase(it); - } else { - it++; - } - } + auto i = std::stable_partition(std::begin(bodies), std::end(bodies), [&](const Body& b) {return !(min <= b.position && b.position <= max); }); + std::vector result{ i, std::end(bodies) }; + bodies.erase(i, std::end(bodies)); return result; } @@ -60,31 +39,19 @@ namespace nbody { std::vector Box::copyBodies(const std::vector& bodies) const { std::vector result; std::copy_if(std::begin(bodies), std::end(bodies), std::back_inserter(result), [&](const Body& body) { - return body.position[0] >= min[0] && body.position[0] <= max[0] && - body.position[1] >=min[1] && body.position[1] <= max[1] && - body.position[2] >=min[2] && body.position[2] <= max[2]; + return min <= body.position && body.position <= max; }); return result; } //check for body inside box bool isContained(const Body& body, const Box& box) { - for (std::size_t i = 0; i < 3; i++) { - if (body.position[i] < box.min[i] || body.position[i] > box.max[i]) { - return false; - } - } - return true; + return !(box.min > body.position || body.position > box.max); } //check for box inside box bool isContained(const Box& inner, const Box& outer) { - for (std::size_t i = 0; i < 3; i++) { - if (inner.min[i] < outer.min[i] || inner.max[i] > outer.max[i]) { - return false; - } - } - return true; + return !(outer.min > inner.min || inner.max > outer.max); } //box volume @@ -92,24 +59,16 @@ namespace nbody { if (!isValid()) { return -1.0; } - double result = 1.0; - - for (std::size_t i = 0; i < 3; i++) { - result *= max[i] - min[i]; - } - return result; + auto diagonal = max - min; + return diagonal.X * diagonal.Y * diagonal.Z; } double Box::maxSidelength() const { - double maxVal = 0.0; - if (!isValid()) { return -1.0; } - for (std::size_t i = 0; i < 3; i++) { - maxVal = std::max(max[i] - min[i], maxVal); - } - return maxVal; + auto diagonal = max - min; + return std::max({ diagonal.X, diagonal.Y, diagonal.Z }); } bool Box::isCorrectBox() const { @@ -121,53 +80,30 @@ namespace nbody { } bool Box::isValid() const { - for (std::size_t i = 0; i < 3; i++) { - if (max[i] < min[i]) { - return false; - } - } - return true; + return max > min; } void Box::printBB(std::size_t parallelId) const { - std::cout << parallelId << ": min "; - for (std::size_t i = 0; i < 3; i++) { - std::cout << ": " << min[i] << " "; - } - std::cout << parallelId << ": max "; - for (std::size_t i = 0; i < 3; i++) { - std::cout << max[i] << " "; - } + std::cout << parallelId << ": min :" << min; + std::cout << parallelId << ": max :" << max; std::cout << '\n'; } //distance from nearest box order to position - double Box::distanceToPosition(const std::array& position) const { - std::size_t inside = 0; - std::array nextPosition = position; - + double Box::distanceToPosition(const Vec3& position) const { if (!isValid()) { return std::numeric_limits::max(); } - for (std::size_t i = 0; i < 3; i++) { - if (nextPosition[i] < min[i]) { - nextPosition[i] = min[i]; - } else if (nextPosition[i] > max[i]) { - nextPosition[i] = max[i]; - } else { - inside++; - } - } - if (inside == 3) { + + bool isInside{ min <= position && position <= max }; + + if (isInside) { return 0.0; } else { - double dist = 0.0; - - for (std::size_t i = 0; i < 3; i++) { - dist += (nextPosition[i] - position[i]) * (nextPosition[i] - position[i]); - } - return sqrt(dist); + auto nextpos = position; + nextpos.ClampComponents(min, max); + return Distance(nextpos, position); } } @@ -178,6 +114,7 @@ namespace nbody { if (!isValid() || !box2.isValid()) { return std::numeric_limits::max(); } + for (std::size_t i = 0; i < 3; i++) { double elem; @@ -190,32 +127,27 @@ namespace nbody { } length += elem * elem; } - return sqrt(length); + return std::sqrt(length); } //determine octree subboxes std::vector Box::octreeSplit() const { - std::vector result; - if (!isValid()) { - return result; + return std::vector{}; } - for (std::size_t i = 0; i < 8; i++) { - Box current = *this; - for (std::size_t j = 0; j < 3; j++) { - double middle = current.min[j] + (current.max[j] - current.min[j]) / 2.0; + const auto diag = max - min; + const auto mid = min + (diag / 2); - //TODO(steinret): unroll the loop, construct explicitly for better readability - if ((i & (1 << j)) != 0) { - current.min[j] = middle; - } else { - current.max[j] = middle; - } - } - result.push_back(current); - } - return result; + return std::vector { + { { min.X, min.Y, min.Z }, { mid.X, mid.Y, mid.Z } }, + { { mid.X, min.Y, min.Z }, { max.X, mid.Y, mid.Z } }, + { { min.X, mid.Y, min.Z }, { mid.X, max.Y, mid.Z } }, + { { mid.X, mid.Y, min.Z }, { max.X, max.Y, mid.Z } }, + { { min.X, min.Y, mid.Z }, { mid.X, mid.Y, max.Z } }, + { { mid.X, min.Y, mid.Z }, { max.X, mid.Y, max.Z } }, + { { min.X, mid.Y, mid.Z }, { mid.X, max.Y, max.Z } }, + { { mid.X, mid.Y, mid.Z }, { max.X, max.Y, max.Z } }}; } //split box into two across longest side @@ -242,16 +174,11 @@ namespace nbody { } //check for position in box - bool Box::contained(const std::array& position) const { + bool Box::contained(const Vec3& position) const { if (!isValid()) { return false; } - for (std::size_t i = 0; i < 3; i++) { - if (position[i] < min[i] || position[i] > max[i]) { - return false; - } - } - return true; + return min <= position && position <= max; } //extend box by box @@ -259,18 +186,14 @@ namespace nbody { if (!extender.isValid()) { return; } - for (std::size_t i = 0; i < 3; i++) { - min[i] = std::min(min[i], extender.min[i]); - max[i] = std::max(max[i], extender.max[i]); - } + min.MinComponents(extender.min); + max.MaxComponents(extender.max); } //extend box by body void Box::extend(const Body& extender) { - for (std::size_t i = 0; i < 3; i++) { - min[i] = std::min(min[i], extender.position[i]); - max[i] = std::max(max[i], extender.position[i]); - } + min.MinComponents(extender.position); + max.MaxComponents(extender.position); } } // namespace nbody diff --git a/bhtree_mpi/src/datastructures/Box.hpp b/bhtree_mpi/src/datastructures/Box.hpp index e0c3dd3..d8e8964 100644 --- a/bhtree_mpi/src/datastructures/Box.hpp +++ b/bhtree_mpi/src/datastructures/Box.hpp @@ -5,18 +5,13 @@ #include #include #include +#include namespace nbody { struct Box { using flt_lim = std::numeric_limits; - std::array min{ { - flt_lim::max(), - flt_lim::max(), - flt_lim::max() }}; - std::array max{{ - flt_lim::min(), - flt_lim::min(), - flt_lim::min() }}; + Vec3 min{ flt_lim::max() }; + Vec3 max{ flt_lim::min() }; void extendToCube(); void extendForBodies(const std::vector& bodies); @@ -27,11 +22,11 @@ namespace nbody { bool isCorrectBox() const; bool isValid() const; void printBB(std::size_t parallelId) const; - double distanceToPosition(const std::array& position) const; + double distanceToPosition(const Vec3& position) const; double distanceToBox(const Box& box2) const; std::vector octreeSplit() const; std::vector splitLongestSide() const; - bool contained(const std::array& position) const; + bool contained(const Vec3& position) const; void extend(const Box& extender); void extend(const Body& extender); }; diff --git a/bhtree_mpi/src/datastructures/Vec3.hpp b/bhtree_mpi/src/datastructures/Vec3.hpp new file mode 100644 index 0000000..8b582e4 --- /dev/null +++ b/bhtree_mpi/src/datastructures/Vec3.hpp @@ -0,0 +1,166 @@ +#pragma once +#include +#include + +namespace { + template + T clamp(const T& v, const T& lo, const T& hi) { // use std::clamp in c++17 + return std::max(lo, std::min(v, hi)); + } +} + +struct Vec3 { + Vec3() = default; + explicit Vec3(double d) : X(d), Y(d), Z(d) {} + Vec3(double x, double y, double z) : X(x), Y(y), Z(z) {}; + + double X{0}; + double Y{0}; + double Z{0}; + + double operator[](std::size_t idx) const { + switch (idx) + { + case 0: + return X; + case 1: + return Y; + case 2: + return Z; + default: + std::abort(); + } + } + + double& operator[](std::size_t idx) { + switch (idx) + { + case 0: + return X; + case 1: + return Y; + case 2: + return Z; + default: + std::abort(); + } + } + + Vec3& operator+=(const Vec3& v) { + X += v.X; + Y += v.Y; + Z += v.Z; + return *this; + } + + Vec3& operator-=(const Vec3& v) { + X -= v.X; + Y -= v.Y; + Z -= v.Z; + return *this; + } + + Vec3& operator*=(double d) { + X *= d; + Y *= d; + Z *= d; + return *this; + } + + Vec3& operator*=(const Vec3& v) { + X *= v.X; + Y *= v.Y; + Z *= v.Z; + return *this; + } + + Vec3& operator/=(double d) { + X /= d; + Y /= d; + Z /= d; + return *this; + } + + bool IsValid() const { + return !std::isnan(X) && std::isfinite(X) && + !std::isnan(Y) && std::isfinite(Y) && + !std::isnan(Z) && std::isfinite(Z); + } + + void MaxComponents(const Vec3& other) { + X = std::max(X, other.X); + Y = std::max(Y, other.Y); + Z = std::max(Z, other.Z); + } + void MinComponents(const Vec3& other) { + X = std::min(X, other.X); + Y = std::min(Y, other.Y); + Z = std::min(Z, other.Z); + } + void ClampComponents(const Vec3& lo, const Vec3& hi) { + X = clamp(X, lo.X, hi.X); + Y = clamp(Y, lo.Y, hi.Y); + Z = clamp(Z, lo.Z, hi.Z); + } + +}; + +inline Vec3 operator-(Vec3 v1, const Vec3& v2) { + v1 -= v2; + return v1; +} + +inline Vec3 operator+(Vec3 v1, const Vec3& v2) { + v1 += v2; + return v1; +} + +inline Vec3 operator*(Vec3 v1, const Vec3& v2) { + v1 *= v2; + return v1; +} + +inline Vec3 operator*(Vec3 v, double d) { + v *= d; + return v; +} + +inline Vec3 operator/(Vec3 v, double d) { + v /= d; + return v; +} + +inline bool operator<(const Vec3& lhs, const Vec3& rhs) { + return lhs.X < rhs.X && + lhs.Y < rhs.Y && + lhs.Z < rhs.Z; +} +inline bool operator>(const Vec3& lhs, const Vec3& rhs) { + return lhs.X > rhs.X && + lhs.Y > rhs.Y && + lhs.Z > rhs.Z; +} +inline bool operator<=(const Vec3& lhs, const Vec3& rhs) { + return lhs.X <= rhs.X && + lhs.Y <= rhs.Y && + lhs.Z <= rhs.Z; +} +inline bool operator>=(const Vec3& lhs, const Vec3& rhs) { + return lhs.X >= rhs.X && + lhs.Y >= rhs.Y && + lhs.Z >= rhs.Z; +} +inline bool operator==(const Vec3& lhs, const Vec3& rhs) { + return lhs.X == rhs.X && + lhs.Y == rhs.Y && + lhs.Z == rhs.Z; +} + +inline double Norm(const Vec3& v) { return std::sqrt(v.X * v.X + v.Y * v.Y + v.Z * v.Z); } + +inline double Distance(const Vec3& v1, const Vec3& v2) { return Norm(v1 - v2); } + +inline std::ostream& operator<<(std::ostream& os, const Vec3 v) { + os << v.X << ":" << v.Y << ":" << v.Z; + return os; +} diff --git a/bhtree_mpi/src/simulation/MpiSimulation.cpp b/bhtree_mpi/src/simulation/MpiSimulation.cpp index cb3acb3..16dfb04 100644 --- a/bhtree_mpi/src/simulation/MpiSimulation.cpp +++ b/bhtree_mpi/src/simulation/MpiSimulation.cpp @@ -13,8 +13,18 @@ namespace nbody { MpiSimulation::MpiSimulation(const std::string& inputFile) { //create MPI datatypes for bodies and domain boxes - int bodyBlocklengths[6] = { 1, 3, 3, 3, 1, 1 }; - MPI_Datatype bodyDatatypes[6] = { MPI_UINT64_T, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_CXX_BOOL }; + int vec3Blocklengths[3] = { 1,1,1 }; + MPI_Datatype vec3Datatypes[3] = { MPI_DOUBLE, MPI_DOUBLE , MPI_DOUBLE }; + MPI_Aint vec3Offsets[3]{ + offsetof(Vec3, X), + offsetof(Vec3, Y), + offsetof(Vec3, Z) + }; + MPI_Type_create_struct(3, vec3Blocklengths, vec3Offsets, vec3Datatypes, &vec3Type); + MPI_Type_commit(&vec3Type); + + int bodyBlocklengths[6] = { 1, 1, 1, 1, 1, 1 }; + MPI_Datatype bodyDatatypes[6] = { MPI_UINT64_T, vec3Type, vec3Type, vec3Type, MPI_DOUBLE, MPI_CXX_BOOL }; MPI_Aint bodyOffsets[6]; bodyOffsets[0] = offsetof(Body, id); bodyOffsets[1] = offsetof(Body, position); @@ -25,8 +35,8 @@ namespace nbody { MPI_Type_create_struct(6, bodyBlocklengths, bodyOffsets, bodyDatatypes, &bodyType); MPI_Type_commit(&bodyType); - int boxBlocklengths[2] = { 3, 3 }; - MPI_Datatype boxDatatypes[2] = { MPI_DOUBLE, MPI_DOUBLE }; + int boxBlocklengths[2] = { 1, 1 }; + MPI_Datatype boxDatatypes[2] = { vec3Type, vec3Type }; MPI_Aint boxOffsets[2]; boxOffsets[0] = offsetof(Box, min); boxOffsets[1] = offsetof(Box, max); @@ -67,6 +77,7 @@ namespace nbody { //cleanup MPI types MPI_Type_free(&bodyType); MPI_Type_free(&boxType); + MPI_Type_free(&vec3Type); delete this->tree; tree = nullptr; } diff --git a/bhtree_mpi/src/simulation/MpiSimulation.hpp b/bhtree_mpi/src/simulation/MpiSimulation.hpp index b266b9e..191fb61 100644 --- a/bhtree_mpi/src/simulation/MpiSimulation.hpp +++ b/bhtree_mpi/src/simulation/MpiSimulation.hpp @@ -19,6 +19,7 @@ namespace nbody { //MPI simulation class MpiSimulation : public Simulation { protected: + MPI_Datatype vec3Type{ MPI_DATATYPE_NULL }; MPI_Datatype bodyType{MPI_DATATYPE_NULL}; MPI_Datatype boxType{MPI_DATATYPE_NULL}; std::vector domains; -- GitLab