commit 6e4dc27862e4114bd8141524b1d2258af3b77458
parent 73d0407a413799f3075c25516674e761756ba5d0
Author: Mikolaj Lenczewski <33129490+EnderRifter@users.noreply.github.com>
Date: Thu, 4 Jul 2019 21:23:49 +0100
Finished initial Barnes-Hut algorithm implementation.
Diffstat:
4 files changed, 91 insertions(+), 17 deletions(-)
diff --git a/StarSim/StarSimLib/Constants.cs b/StarSim/StarSimLib/Constants.cs
@@ -69,6 +69,13 @@ namespace StarSimLib
public const double TimeStep = SecondsPerTick * (SimulationRate / (float)FrameRate) * SimulationRate;
/// <summary>
+ /// The tolerance of the mass grouping approximation in the simulation. A
+ /// body is only accelerated when the ratio of the tree's width to the
+ /// distance (from the tree's center of mass to the body) is less than this.
+ /// </summary>
+ public const double TreeTheta = 0.5;
+
+ /// <summary>
/// The maximum radius within which <see cref="Body"/>s will be placed (1e18f).
/// </summary>
public const double UniverseSize = 1e18f;
diff --git a/StarSim/StarSimLib/Data Structures/Body.cs b/StarSim/StarSimLib/Data Structures/Body.cs
@@ -159,8 +159,9 @@ namespace StarSimLib.Data_Structures
/// <param name="otherBody">The other instance with which to collide.</param>
public void Collide(Body otherBody)
{
- mass += otherBody.Mass;
- velocity += otherBody.Velocity;
+ mass += otherBody.mass;
+ velocity += otherBody.velocity;
+ force += otherBody.force;
}
/// <summary>
@@ -186,7 +187,7 @@ namespace StarSimLib.Data_Structures
/// </summary>
/// <param name="octant">The octant to check.</param>
/// <returns>Whether this <see cref="Body"/> instance is inside the given <see cref="Octant"/> instance.</returns>
- public bool InOctant(Octant octant)
+ public bool IsInOctant(Octant octant)
{
return octant.ContainsPoint(position);
}
diff --git a/StarSim/StarSimLib/Data Structures/OctantTree.cs b/StarSim/StarSimLib/Data Structures/OctantTree.cs
@@ -77,16 +77,36 @@ namespace StarSimLib.Data_Structures
// this is an empty instance that has not yet had any bodies added to it.
body = newBody;
}
+ else if (IsExternal())
+ {
+ // this instance is 'external' and contains another body. figure out where the new body should go and
+ // create a new octant tree instance to hold the new body
+ foreach (OctantTree tree in childTrees)
+ {
+ if (body.IsInOctant(tree.octant))
+ {
+ tree.AddBody(body);
+ }
+ }
+
+ AddBody(newBody);
+ }
else if (!IsExternal())
{
// this instance already has a body to represent it, and it is not an 'external' tree instance, that is
// it has child trees of its own. figure out in which child tree the new body should be stored and update
- // any further child nodes with recursion
- }
- else if (IsExternal())
- {
- // this instance is 'external' and contains another body. figure out where the new body should go and
- // create a new octant tree instance to hold the new body. no recursion is necessary
+ // any further child nodes
+
+ // make the held body an aggregate body
+ body.Collide(newBody);
+
+ foreach (OctantTree tree in childTrees)
+ {
+ if (newBody.IsInOctant(tree.octant))
+ {
+ tree.AddBody(newBody);
+ }
+ }
}
}
@@ -134,5 +154,35 @@ namespace StarSimLib.Data_Structures
"The given specifier was outside of the valid range.");
}
}
+
+ /// <summary>
+ /// Recursively updates the forces on each <see cref="Body"/> instance held in this tree, with respect to the
+ /// given reference <see cref="Body"/> instance.
+ /// </summary>
+ /// <param name="referenceBody">The body instance against which force updates are made.</param>
+ public void UpdateForces(Body referenceBody)
+ {
+ if (IsExternal())
+ {
+ // since this tree instance is 'external' it has no children. we can treat it as a single body
+ if (body != referenceBody)
+ {
+ referenceBody?.AddForce(body);
+ }
+ }
+ else if (octant.Length / body.DistanceTo(referenceBody, out _) < Constants.TreeTheta)
+ {
+ // otherwise if the octant length divided by the distance to the body (the width to distance ratio) is
+ // within a defined tolerance, we consider the tree to be effectively a single massive body
+ referenceBody?.AddForce(body);
+ }
+ else
+ {
+ foreach (OctantTree tree in childTrees)
+ {
+ tree?.UpdateForces(referenceBody);
+ }
+ }
+ }
}
}
\ No newline at end of file
diff --git a/StarSim/StarSimLib/Physics/BodyUpdater.cs b/StarSim/StarSimLib/Physics/BodyUpdater.cs
@@ -17,6 +17,12 @@ namespace StarSimLib.Physics
public static class BodyUpdater
{
/// <summary>
+ /// Represents the main universe. Anything outside this octant of space is not updated when using the
+ /// <see cref="UpdateBodiesBarnesHut"/> update method.
+ /// </summary>
+ public static readonly Octant UniverseOctant = new Octant(new Vector4(), Constants.UniverseSize);
+
+ /// <summary>
/// Updates the positions of all the given <see cref="Body"/>s with O(n^2) time complexity, with the given time step.
/// </summary>
/// <param name="bodies">The collection of <see cref="Body"/>s whose positions to update.</param>
@@ -48,19 +54,29 @@ namespace StarSimLib.Physics
public static void UpdateBodiesBruteForce(IEnumerable<Body> bodies, double deltaTime)
{
IEnumerable<Body> bodyEnumerable = bodies as Body[] ?? bodies.ToArray();
- Vector4 forceVector = new Vector4();
+ OctantTree barnesHutTree = new OctantTree(UniverseOctant);
+ // we construct our barnes-hut octant tree
foreach (Body body in bodyEnumerable)
{
- // resets the force vector to avoid another instantiation and allocation
- forceVector.X = 0;
- forceVector.Y = 0;
- forceVector.Z = 0;
+ if (body.IsInOctant(UniverseOctant))
+ {
+ // only update a body's position if it is within the bounds of the universe
+ barnesHutTree.AddBody(body);
+ }
+ }
- // use LINQ expression as it is more concise; sum attraction vectors for all other bodies
- forceVector = bodyEnumerable.Where(b => b != body).Aggregate(forceVector, (current, b) => current + Body.GetForceBetween(body, b));
+ // we update the positions of the bodies in the populated tree
+ foreach (Body body in bodyEnumerable)
+ {
+ body.ResetForce();
- body.Update(forceVector, deltaTime);
+ if (body.IsInOctant(UniverseOctant))
+ {
+ barnesHutTree.UpdateForces(body);
+
+ body.Update(deltaTime);
+ }
}
}
}
