Skip to content
Snippets Groups Projects
Commit 6eab07a8 authored by Paul Heinzlreiter's avatar Paul Heinzlreiter
Browse files

* force computation and body movement

parent d5d88aea
Branches
No related merge requests found
......@@ -107,6 +107,9 @@ namespace nbody {
double dxdt[3];
double dvdt[3];
if (this->refinement) {
return;
}
a = evaluate(0.0, start);
b = evaluate(timestep * 0.5, a);
c = evaluate(timestep * 0.5, b);
......@@ -119,21 +122,21 @@ namespace nbody {
}
}
void Body::accumulateForce(Body from) {
void Body::accumulateForceOnto(Body& to) {
double distance2 = 0.0;
distance2 += (this->position[0] - from.position[0]) * (this->position[0] - from.position[0]);
distance2 += (this->position[1] - from.position[1]) * (this->position[1] - from.position[1]);
distance2 += (this->position[2] - from.position[2]) * (this->position[2] - from.position[2]);
distance2 += (this->position[0] - to.position[0]) * (this->position[0] - to.position[0]);
distance2 += (this->position[1] - to.position[1]) * (this->position[1] - to.position[1]);
distance2 += (this->position[2] - to.position[2]) * (this->position[2] - to.position[2]);
if (fabs(distance2) < FLT_EPSILON) {
distance2 = FLT_EPSILON;
}
double distance = sqrt(distance2);
double mdist = -1.0 * ((this->mass * from.mass) / distance2);
double mdist = -1.0 * ((this->mass * to.mass) / distance2);
for (int i = 0; i < 3; i++) {
double vec = (this->position[i] - from.position[i]) / distance;
double vec = (this->position[i] - to.position[i]) / distance;
this->acceleration[i] += vec * mdist;
to.acceleration[i] += vec * mdist;
}
}
......
......@@ -26,7 +26,7 @@ namespace nbody {
virtual ~Body();
virtual void resetAcceleration();
virtual void integrate();
virtual void accumulateForce(Body from);
virtual void accumulateForceOnto(Body& from);
virtual void print();
};
}
......
#include <iostream>
#include <cmath>
#include "Box.hpp"
namespace nbody {
......@@ -136,6 +137,32 @@ namespace nbody {
return max;
}
double Box::distanceToPosition(double* position) {
int inside = 0;
double nextPosition[3] = {position[0], position[1], position[2]};
for (int i = 0; i < 3; i++) {
if (nextPosition[i] < this->min[i]) {
nextPosition[i] = this->min[i];
} else if (nextPosition[i] > this->max[i]) {
nextPosition[i] = this->max[i];
} else {
inside++;
}
}
if (inside == 3) {
return 0.0;
} else {
double dist = 0.0;
for (int i = 0; i < 3; i++) {
dist += (nextPosition[i] - position[i]) * (nextPosition[i] - position[i]);
}
return sqrt(dist);
}
}
bool Box::isCorrect() {
for (int i = 0; i < 3; i++) {
if (this->max[i] < this->min[i]) {
......
......@@ -29,6 +29,7 @@ namespace nbody {
virtual bool isCorrect();
virtual void print();
virtual bool overlapsSphere(double* sphereCenter, double sphereRadius);
virtual double distanceToPosition(double* position);
};
}
......
#include <cfloat>
#include <cstdlib>
#include <cmath>
#include <iostream>
#include "Node.hpp"
......@@ -149,7 +150,7 @@ namespace nbody {
}
} else {
for (Node* node = this->next; node != this->tree->nodes && node != NULL; node = node->afterSubtree) {
mass += node->representative.mass;
}
}
for (int i = 0; i < 3; i++) {
......@@ -169,4 +170,13 @@ namespace nbody {
it->print();
}
}
bool Node::sufficientForBody(Body body) {
double distance = 0.0;
for (int i = 0; i < 3; i++) {
distance += (this->representative.position[i] - body.position[i]) * (this->representative.position[i] - body.position[i]);
}
return sqrt(distance) > this->getL();
}
}
......@@ -17,6 +17,7 @@ namespace nbody {
protected:
Box bb;
vector<Body> bodies;
vector<Body> refinements;
Node* prev;
Node* next;
Node* nextSibling;
......@@ -41,6 +42,7 @@ namespace nbody {
virtual double getL();
virtual bool isCorrect();
virtual void print();
virtual bool sufficientForBody(Body body);
};
}
......
......@@ -51,12 +51,27 @@ namespace nbody {
return nodes;
}
void Tree::accumulateForceFor(Body& body) {
void Tree::accumulateForceOnto(Body& body) {
Node* n = this->nodes->next;
body.resetAcceleration();
while (n != this->nodes) {
if (n->sufficientForBody(body)) {
n->representative.accumulateForceOnto(body);
} else if (n->leaf) {
for (vector<Body>::iterator it = n->bodies.begin(); it != n->bodies.end(); it++) {
it->accumulateForceOnto(body);
}
}
n = n->afterSubtree;
}
}
void Tree::computeForces() {
for (Node* n = this->nodes->next; n != this->nodes; n = n->next) {
if (n->leaf) {
for (vector<Body>::iterator it = n->bodies.begin(); it != n->bodies.end(); it++) {
body.accumulateForce(*it);
this->accumulateForceOnto(*it);
}
}
}
......@@ -76,4 +91,33 @@ namespace nbody {
infile.close();
return result;
}
void Tree::moveBodies() {
for (Node* n = this->nodes->next; n != this->nodes; n = n->next) {
if (n->leaf) {
for (vector<Body>::iterator it = n->bodies.begin(); it != n->bodies.end(); it++) {
it->integrate();
}
}
}
}
vector<Body> Tree::extractBodies() {
vector<Body> result;
while (this->nodes->next != this->nodes) {
if (this->nodes->next->leaf) {
result.insert(result.end(), nodes->next->bodies.begin(), nodes->next->bodies.end());
}
Node* h = this->nodes->next;
this->nodes->next->unlink();
delete(h);
}
this->clean();
return result;
}
void Tree::rebuildTree() {
this->build(this->extractBodies());
}
}
......@@ -24,7 +24,11 @@ namespace nbody {
virtual void build(vector<Body> bodies) = 0;
virtual int numberOfChildren() = 0;
virtual unsigned long numberOfNodes();
virtual void accumulateForceFor(Body& body);
virtual void accumulateForceOnto(Body& body);
virtual void computeForces();
virtual void moveBodies();
virtual vector<Body> extractBodies();
virtual void rebuildTree();
virtual bool isCorrect();
};
......
......@@ -12,11 +12,11 @@ int main(int argc, char* argv[]) {
cout << bodies.size() << endl;
tree.build(bodies);
if (tree.isCorrect()) {
cout << "Tree correct" << endl;
} else {
if (!tree.isCorrect()) {
cout << "Tree not correct" << endl;
return -1;
}
} else {
std::cout << "You need to specify the body input file." << endl;
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment