From c9c7856fa9c0c8181fb9146e1ee81fe921a44bc7 Mon Sep 17 00:00:00 2001
From: "armindamon.riess" <armindamon.riess@uzh.ch>
Date: Wed, 7 Dec 2022 22:48:25 +0100
Subject: [PATCH] finish force calculation. in theory sim is done
---
lib/node.cpp | 63 +++++++++++++++++++++++++++++++++++++++++++++++-----
lib/tree.cpp | 1 +
2 files changed, 58 insertions(+), 6 deletions(-)
diff --git a/lib/node.cpp b/lib/node.cpp
index 108d063..432569a 100644
--- a/lib/node.cpp
+++ b/lib/node.cpp
@@ -93,10 +93,14 @@ void Node::update(List& localParticles, unsigned nLocalParticles) {
calculateCenterOfMass();
localParticles_ = localParticles;
nLocalParticles_ = nLocalParticles;
+
// for each local particle, determine which octant it is in
+ // get the mass in this node while we're at it
+ mass_ = 0;
List* particleDivision = new List[nLocalParticles_];
for (unsigned i=0; i<nLocalParticles_; ++i) {
particleDivision[getOctant(localParticles_[i])].push_back(i);
+ mass_ += masses_[localParticles_[i]];
}
isLeaf_ = true;
@@ -130,6 +134,8 @@ void Node::update(List& localParticles, unsigned nLocalParticles) {
}
double* Node::calculateForce(unsigned particle) {
+ double* force = new double[3];
+ for (unsigned i=0; i<3; ++i) force[i] = 0;
// calculate theta
double dx = center_[0] - particles_[3*particle + 0];
double dy = center_[1] - particles_[3*particle + 1];
@@ -137,22 +143,67 @@ double* Node::calculateForce(unsigned particle) {
double lambda = std::sqrt(dx*dx + dy*dy + dz*dz);
double theta = size_/lambda;
- // TODO
- // If angle too large, call this function for every child which is not a nullptr.
+ // If angle too large:
// If the node is a leaf, calculate force on particle from particles in this node.
// If the node is not a leaf, call this function for every child which is not a nullptr.
// If angle small enough, calculate force on particle from this node.
if (theta > theta0_) {
if (isLeaf_) {
- // TODO
+
+ for (unsigned i : localParticles_) {
+ if (i != particle) {
+
+ double dx = particles_[3*i + 0] - particles_[3*particle + 0];
+ double dy = particles_[3*i + 1] - particles_[3*particle + 1];
+ double dz = particles_[3*i + 2] - particles_[3*particle + 2];
+
+ double r = std::sqrt(dx*dx + dy*dy + dz*dz);
+ double r3 = r*r*r;
+ double m = masses_[i];
+
+ for (unsigned j=0; j<3; ++j) {
+ force[j] += m * (particles_[3*i + j] - particles_[3*particle + j]) / r3;
+ }
+ }
+ }
+
} else {
- // TODO
+
+ // store the forces from each child in an array so that we can delete them later
+ double* childForces[8];
+
+ for (unsigned i=0; i<8; ++i) {
+ // get force for each child that exists and add them together
+ if (children_[i] != nullptr) {
+
+ childForces[i] = children_[i]->calculateForce(particle);
+ for (unsigned j=0; j<3; ++j) {
+ force[j] += childForces[i][j];
+ }
+
+ } else {
+ childForces[i] = nullptr;
+ }
+ }
+
+ for (unsigned i=0; i<8; ++i) {
+ delete[] childForces[i];
+ }
+
}
} else {
- // TODO
+
+ double r = std::sqrt(dx*dx + dy*dy + dz*dz);
+ double r3 = r*r*r;
+ double m = mass_;
+
+ for (unsigned j=0; j<3; ++j) {
+ force[j] += m * (centerOfMass_[j] - particles_[3*particle + j]) / r3;
+ }
+
}
- return nullptr;
+ return force;
}
void Node::calculateCenterOfMass() {
diff --git a/lib/tree.cpp b/lib/tree.cpp
index eb27fa8..3904c6d 100644
--- a/lib/tree.cpp
+++ b/lib/tree.cpp
@@ -29,6 +29,7 @@ Tree::~Tree() {
void Tree::calculateForce(unsigned particle) {
double* force = root_->calculateForce(particle);
for (unsigned i=0; i < 3; ++i) forces_[3*particle + i] = force[i];
+ delete[] force;
}
void Tree::update() {
--
GitLab