6
\$\begingroup\$

Link to previous question: Barnes-Hut N-body simulator


It still has issues I'm sure.

Node.h

#pragma once
#include <vector>

struct Body
{
    double posX, posY;          //position x and y
    double velX, velY;          //velocity x and y
    double AccelX, AccelY;      //force acting on object since last frame
    float mass;                 //mass of object
};

class Node
{
public:
    std::vector<Body*> Bodies;
    std::vector<Node*> Child;
    bool HasChildren;

    float posX, posY;
    float width, height;
    float TotalMass;
    float CenterOfMassx;
    float CenterOfMassy;
    unsigned int Depth;

    Node();
    Node(std::vector<Body*> pBodies, float pwidth, float pheight, float px = 0, float py = 0, unsigned int pDepth = 0);
    ~Node();
    void GenerateChildren();
    void SetParam(std::vector<Body*> pBodies, float pwidth, float pheight, float px = 0, float py = 0, unsigned int pDepth = 0);
    void Reset();
};

Node.cpp

#include "Node.h"
#include <memory>

#define _MAX_DEPTH 100

Node::Node(std::vector<Body*> pBodies, float pwidth, float pheight, float px, float py, unsigned int pDepth)
{
    SetParam(pBodies, pwidth, pheight, px, py, pDepth);
}

void Node::SetParam(std::vector<Body*> pBodies, float pwidth, float pheight, float px, float py, unsigned int pDepth)
{
    HasChildren = false;
    Depth = pDepth;

    unsigned int size = pBodies.size();

    Bodies = pBodies;
    posX = px;
    posY = py;
    width = pwidth;
    height = pheight;

    float mass = 0;
    double Centerx = 0;
    double Centery = 0;

    for (unsigned int i = 0; i < size; i++)
    {
        mass += pBodies[i]->mass;
        Centerx += pBodies[i]->posX;
        Centery += pBodies[i]->posY;
    }

    TotalMass = mass;

    if (size > 0)
    {
        CenterOfMassx = static_cast<float>(Centerx / size);
        CenterOfMassy = static_cast<float>(Centery / size);
    }
    else
    {
        CenterOfMassx = 0;
        CenterOfMassy = 0;
    }

    if (size > 1 && Depth < _MAX_DEPTH)
    {
        GenerateChildren();
    }
}

Node::Node()
{
    HasChildren = false;
    Depth = 0;
}

void Node::GenerateChildren()
{
    std::vector<Body*> q1Bodies, q2Bodies, q3Bodies, q4Bodies;

    for (unsigned int i = 0; i < Bodies.size(); i++)
    {
        if (Bodies[i]->posX < (posX + (width / 2)))     //if true, 1st or 3rd
        {
            if (Bodies[i]->posY < (posY + (height / 2)))    //1
            {
                q1Bodies.push_back(Bodies[i]);
            }
            else //3
            {
                q3Bodies.push_back(Bodies[i]);
            }
        }
        else                                            //2 or 4
        {
            if (Bodies[i]->posY < (posY + (height / 2)))    //2
            {
                q2Bodies.push_back(Bodies[i]);
            }
            else //4
            {
                q4Bodies.push_back(Bodies[i]);
            }
        }
    }

    Node* q1 = new Node(q1Bodies, width / 2, height / 2, posX, posY, Depth + 1);
    Node* q2 = new Node(q2Bodies, width / 2, height / 2, posX + (width / 2), posY, Depth + 1);
    Node* q3 = new Node(q3Bodies, width / 2, height / 2, posX, posY + (height / 2), Depth + 1);
    Node* q4 = new Node(q4Bodies, width / 2, height / 2, posX + (width / 2), posY + (height / 2), Depth + 1);

    Child.push_back(q1);
    Child.push_back(q2);
    Child.push_back(q3);
    Child.push_back(q4);

    HasChildren = true;
}

void Node::Reset()
{
    Bodies.clear();

    for (unsigned int i = 0; i < Child.size(); i++)
    {
        Child[i]->Reset();
        delete Child[i];
    }       

    Child.clear();

    HasChildren = false;
}

Node::~Node()
{

}

Main.cpp

#include "SFML/Graphics.hpp"
#include "Node.h"
#include <ctime>
#include <vector>

void BodyAttraction(std::vector<Body*> &pBodies);                                                                       //Attracts each body to each other body in the given vector of pointers to body objects
void CalculateForceNode(Body* bi, Node* bj);                                                                            //Calculate force exerted on body from node
void CalculateForce(Body* bi, Body* bj);                                                                                //Calculate force exerted on eachother between two bodies
Body* CreateBody(float px, float py, float pmass, float pvx = 0, float pvy = 0);                                        //return a pointer to new body object defined on the heap with given paramiters
void DeleteBodies(std::vector<Body*> &pBodies);                                                                         //Deletes objects pointed to by given vector
void PollEvent(sf::RenderWindow* pTarget, sf::View* pSimView, float &pZoom, bool* pIsPaused);                           //Call all polled events for the sf::window
void PopulateBodyVectorDisk(std::vector<Body*> *pBodies, unsigned int pParticlesCount, float pWidth, float pHeight, float pMaxDist, float pMinDist, float pMinMass, float pMaxMass, float pGalaticCenterMass = 0);  //populate given vector with bodies with given paramiters in a disk formation
void PopulateBodyVectorUniform(std::vector<Body*> *pBodies, unsigned int pParticlesCount, float pWidth, float pHeight, float pMinMass, float pMaxMass);
void PopulateBodyVectorExplosion(std::vector<Body*> *pBodies, unsigned int pParticlesCount, float pWidth, float pHeight, float pMinMass, float pMaxMass, float ExplosionForce);
void Render(sf::RenderWindow* pTarget, std::vector<Body*> &pBodies, float pZoom);                   //Render given body objects to given screen
void SetView(sf::View* pView, sf::RenderWindow* pTarget, float pViewWidth, float pViewHeight);                          //set the window to the simulation view
void UpdateBodies(std::vector<Body*> &pBodies);                                                                         //Calculate velocity chance from the bodies force exerted since last update, update position based on velocity, reset force to 0
void DrawNode(Node* pNode, sf::RenderWindow* pTarget, float pZoom);                                                     //Draw a node to the screen, and all of its children (recursive)
void CheckNode(Node* pNode, Body* pBody);                                                                               //Checks if a node is sufficently far away for force calculation, if not recureses on nodes children
void OctreeBodyAttraction(std::vector<Body*> &pBodies, Node *pGlobalNode);                                              //Using a calculated oct-tree, calculate the force exerted on each object
void AttractToCenter(std::vector<Body*> &pBodies, float width, float height, float centerMass);                         //Attract each particle to the center of the simulation
void ResetAccel(std::vector<Body*> &pBodies);
void RepelFromCenter(std::vector<Body*> &pBodies, float width, float height, float centerMass);             

double const _PI = 3.14159265;      //Pi, used for calculations and rounded to 8 decimal places. 
double const NodeThresholdSqr = 0.5 * 0.5;  //Threshold for node calculations   
double const _GRAV_CONST = 0.1;     //the gravitational constant. This is the timestep between each frame. Lower for slower but more accurate simulations
float const Softener = 10;          //A softener used for the force calculations, 10 is a good amount
float const SimWidth = 327680;      //Width and height of simulation, needs to be large, particles outside of this range will not be included in the octree
float const SimHeight = 327680;                                                                                                                 

int main()
{
    unsigned int const NumParticles = 10000;        //Number of particles in simtulation, currently 2^15    

    float const DiskRadiusMax = 20000;              //Max and min distances objects will be from the galatic center
    float const DiskRadiusMin = 50;
    float const GalaticCenterMass = 1000000;        //The mass of the very large object simulating a black hole at the center of a galixy;
    float const ObjectMassMax = 2;                  //The max and min mass of the objects in the galixy
    float const ObjectMassMin = 1;  
    float const ViewWidth = 1920;                   //Width and height of view of the simulation for the screen. 
    float const ViewHeight = 1080;

    Node GlobalNode;
    std::vector<Body*> Bodies;                      //Container of all Bodies in simulation
    sf::Color ObjColor(255, 255, 255, 128);         //the defult colour of the objects
    sf::View SimulationView;
    sf::RenderWindow window(sf::VideoMode(1920, 1080), "N-Body simulation");

    bool IsPaused = false;                          //Contains the state of weather the simulation is paused or not
    float zoom = 1;                                 //The current amount of zoom in or out the user has inputed in total

    srand(static_cast<unsigned int> (time(NULL)));

    PopulateBodyVectorDisk(&Bodies, NumParticles, SimWidth, SimHeight, DiskRadiusMax, DiskRadiusMin, ObjectMassMin, ObjectMassMax, GalaticCenterMass);
    SetView(&SimulationView, &window, ViewWidth, ViewHeight);

    while (window.isOpen())
    {
        PollEvent(&window, &SimulationView, zoom, &IsPaused);   //These will always be done

        if (!IsPaused)  //These will not if the simulation is paused
        {       
            AttractToCenter(Bodies, SimWidth, SimHeight, GalaticCenterMass);
            UpdateBodies(Bodies);
            ResetAccel(Bodies);
            GlobalNode.Reset();
            GlobalNode.SetParam(Bodies, SimWidth, SimHeight);
            OctreeBodyAttraction(Bodies, &GlobalNode);
        }   

        Render(&window, Bodies, zoom);
    }

    DeleteBodies(Bodies);
}

void AttractToCenter(std::vector<Body*> &pBodies, float width, float height, float centerMass)
{
    Body* Temp = CreateBody(width / 2, height / 2, centerMass); //Create a body at the center of the simulation

    for (unsigned int i = 0; i < pBodies.size(); i++)
    {
        CalculateForce(pBodies[i], Temp);
    }

    delete Temp;
}

void RepelFromCenter(std::vector<Body*> &pBodies, float width, float height, float centerMass)
{
    Body* Temp = CreateBody(width / 2, height / 2, centerMass); //Create a body at the center of the simulation

    for (unsigned int i = 0; i < pBodies.size(); i++)
    {
        double vectorx = Temp->posX - pBodies[i]->posX;
        double vectory = Temp->posY - pBodies[i]->posY;

        double distSqr = vectorx * vectorx + vectory * vectory;

        double Dist = (sqrt(distSqr));

        double accel = (pBodies[i]->mass * Dist * _GRAV_CONST * 0.0001);

        pBodies[i]->AccelX -= vectorx * accel;
        pBodies[i]->AccelY -= vectory * accel;
    }

    delete Temp;
}

void ResetAccel(std::vector<Body*> &pBodies)
{
    for (unsigned int i = 0; i < pBodies.size(); i++)
    { 
        pBodies[i]->AccelX = 0;
        pBodies[i]->AccelY = 0;
    }
}

void BodyAttraction(std::vector<Body*> &pBodies)
{
    for (unsigned int i = 0; i < pBodies.size(); i++)
    {
        for (unsigned int j = 0; j < pBodies.size(); j++)
        {
            CalculateForce(pBodies.at(i), pBodies.at(j)); //for each body in pBodies: each other body in pBodies: Calculate attractive force exerted on the first body from the second one
        }
    }
}

void OctreeBodyAttraction(std::vector<Body*> &pBodies, Node *pGlobalNode)
{
    for (unsigned int i = 0; i < pBodies.size(); i++)
    {   
        CheckNode(pGlobalNode, pBodies[i]);
    }
}

void CheckNode(Node* pNode, Body* pBody)
{
    double diffX = (pNode->CenterOfMassx - pBody->posX);
    double diffY = (pNode->CenterOfMassy - pBody->posY);

    double distanceSqr = ((diffX) * (diffX) + (diffY) * (diffY));                       //Distance from the node to the object                          
    double widthSqr = pNode->width * pNode->width;

    if (widthSqr / distanceSqr < NodeThresholdSqr || pNode->HasChildren == false)       //if sufficently far away or has no children (external node) group node and calculate force
    {   
        CalculateForceNode(pBody, pNode);
    }
    else                                                                            //if not, repeat function with children
    {
        if (pNode->Child[0]->Bodies.size() > 0)
            CheckNode(pNode->Child[0], pBody);
        if (pNode->Child[1]->Bodies.size() > 0)
            CheckNode(pNode->Child[1], pBody);
        if (pNode->Child[2]->Bodies.size() > 0)
            CheckNode(pNode->Child[2], pBody);
        if (pNode->Child[3]->Bodies.size() > 0)
            CheckNode(pNode->Child[3], pBody);
    }
}

void CalculateForceNode(Body* bi, Node* bj)  //bi is being attracted to bj. 15 flops of calculation
{
    //vector from the body to the center of mass
    double vectorx = bj->CenterOfMassx - bi->posX;
    double vectory = bj->CenterOfMassy - bi->posY;

    //c^2 = a^2 + b^2 + softener^2
    double distSqr = vectorx * vectorx + vectory * vectory + Softener * Softener;

    // ivnDistCube = 1/distSqr^(3/2)
    double distSixth = distSqr * distSqr * distSqr;
    double invDistCube = 1.0f / (sqrt(distSixth));


    double Accel = (bj->TotalMass * invDistCube * _GRAV_CONST);

    bi->AccelX += vectorx * Accel;
    bi->AccelY += vectory * Accel;
}

void CalculateForce(Body* bi, Body* bj)  //bi is being attracted to bj. 15 flops of calculation
{
    //std::vector from i to j
    double vectorx = bj->posX - bi->posX;
    double vectory = bj->posY - bi->posY;

    //c^2 = a^2 + b^2 + softener^2
    double distSqr = vectorx * vectorx + vectory * vectory + Softener * Softener;

    // ivnDistCube = 1/distSqr^(3/2)
    double distSixth = distSqr * distSqr * distSqr;
    double invDistCube = 1.0f / (sqrt(distSixth));


    double Accel = (bj->mass * invDistCube * _GRAV_CONST);

    bi->AccelX += vectorx * Accel;
    bi->AccelY += vectory * Accel;
}

Body* CreateBody(float px, float py, float pmass, float pvx, float pvy)
{
    Body* Temp = new Body;

    Temp->posX = px;
    Temp->posY = py;
    Temp->mass = pmass;
    Temp->velX = pvx;
    Temp->velY = pvy;
    Temp->AccelX = Temp->AccelY = 0;

    return Temp;
}

void DeleteBodies(std::vector<Body*> &pBodies)
{
    for (unsigned int i = 0; i < pBodies.size(); i++)
    {
        delete pBodies[i];
    }

    pBodies.clear();
}

void PollEvent(sf::RenderWindow* pTarget, sf::View* pSimView, float &pZoom, bool* pIsPaused)
{
    sf::Event event;

    while (pTarget->pollEvent(event))
    {
        if (event.type == sf::Event::Closed)
            pTarget->close();
        if (event.type == sf::Event::KeyPressed)
        {
            if (event.key.code == sf::Keyboard::Space)
                *pIsPaused = !*pIsPaused;                   //toggle what is pointed to by IsPaused
        }
        if (event.type == sf::Event::MouseWheelScrolled)
        {
            pZoom *= 1 + (static_cast <float> (-event.mouseWheelScroll.delta) / 10);    //for each notch down, -10%, for each notch up, +10%
            pSimView->zoom(1 + (static_cast <float> (-event.mouseWheelScroll.delta) / 10));
        }
    }

    if (sf::Mouse::getPosition().x > (1920 - 20))
        pSimView->move(2 * pZoom, 0);
    if (sf::Mouse::getPosition().x < (0 + 20))
        pSimView->move(-2 * pZoom, 0);
    if (sf::Mouse::getPosition().y > (1080 - 20))
        pSimView->move(0, 2 * pZoom);
    if (sf::Mouse::getPosition().y < (0 + 20))
        pSimView->move(0, -2 * pZoom);

    pTarget->setView(*pSimView);
}

void PopulateBodyVectorDisk(std::vector<Body*> *pBodies, unsigned int pParticlesCount, float pWidth, float pHeight, float pMaxDist, float pMinDist, float pMinMass, float pMaxMass, float pGalaticCenterMass)
{
    for (unsigned int i = 0; i < pParticlesCount; i++)
    {
        float angle = static_cast <float> (rand()) / (static_cast <float> (RAND_MAX / (2 * static_cast <float> (_PI))));    //sets angle to random float range (0, 2 pi) 

        float distanceCoefficent = static_cast <float> (rand()) / static_cast <float> (RAND_MAX);
        float distance = pMinDist + ((pMaxDist - pMinDist) * (distanceCoefficent * distanceCoefficent));                    //Distance point will be from the galatic center, between MinDiskRadius and MaxDiskRadius

        float positionx = cos(angle) * distance + (pWidth / 2);                                                             //set positionx and positiony to be the point you get when you go in the direction of 'angle' till you have traveled 'distance' 
        float positiony = sin(angle) * distance + (pHeight / 2);

        float orbitalVelocity = sqrt((pGalaticCenterMass * static_cast <float> (_GRAV_CONST)) / distance);                  //Calculate the orbital velocity required to orbit the galatic centre   

        float velocityx = (sin(angle) * orbitalVelocity);
        float velocityy = (-cos(angle) * orbitalVelocity);

        float mass = pMinMass + static_cast <float> (rand() % static_cast <int> (pMaxMass - pMinMass));                     //random mass (int) in range (MinObjectMass, MaxObjectMass)

        pBodies->push_back(CreateBody(positionx, positiony, mass, velocityx, velocityy));                                   
    }
}

void PopulateBodyVectorExplosion(std::vector<Body*> *pBodies, unsigned int pParticlesCount, float pWidth, float pHeight, float pMinMass, float pMaxMass, float ExplosionForce)
{
    for (unsigned int i = 0; i < pParticlesCount; i++)
    {
        float angle = static_cast <float> (rand()) / (static_cast <float> (RAND_MAX / (2 * static_cast <float> (_PI))));    //sets angle to random float range (0, 2 pi) 
        float velocity = static_cast <float> (rand()) / static_cast <float> (RAND_MAX);

        float positionx = (pWidth / 2);                                                             //set positionx and positiony to be the point you get when you go in the direction of 'angle' till you have traveled 'distance' 
        float positiony = (pHeight / 2);

        float velocityx = (cos(angle) * ExplosionForce * (1 - (velocity * velocity)));
        float velocityy = (sin(angle) * ExplosionForce * (1 - (velocity * velocity)));

        float mass = pMinMass + static_cast <float> (rand() % static_cast <int> (pMaxMass - pMinMass));                     //random mass (int) in range (MinObjectMass, MaxObjectMass)

        pBodies->push_back(CreateBody(positionx, positiony, mass, velocityx, velocityy));
    }
}

void PopulateBodyVectorUniform(std::vector<Body*> *pBodies, unsigned int pParticlesCount, float pWidth, float pHeight, float pMinMass, float pMaxMass)
{
    for (unsigned int i = 0; i < pParticlesCount; i++)
    {
        float positionx = static_cast <float> (rand()) / static_cast <float> (RAND_MAX) * pWidth + (SimWidth / 2 - pWidth / 2);                         
        float positiony = static_cast <float> (rand()) / static_cast <float> (RAND_MAX) * pHeight + (SimHeight / 2 - pHeight / 2);

        float mass = pMinMass + static_cast <float> (rand() % static_cast <int> (pMaxMass - pMinMass));                     //random mass (int) in range (MinObjectMass, MaxObjectMass)

        pBodies->push_back(CreateBody(positionx, positiony, mass));
    }
}

void Render(sf::RenderWindow* pTarget, std::vector<Body*> &pBodies, float pZoom)
{
    pTarget->clear();

    sf::RectangleShape Temp;

    for (unsigned int i = 0; i < pBodies.size(); i++)
    {       
        if (pZoom > 1)
            Temp.setSize(sf::Vector2f(pBodies.at(i)->mass * pZoom, pBodies.at(i)->mass * pZoom));
        else
            Temp.setSize(sf::Vector2f(pBodies.at(i)->mass, pBodies.at(i)->mass));   

        float AccelCoefficient = static_cast<float> (sqrt(pBodies.at(i)->AccelX * pBodies.at(i)->AccelX + pBodies.at(i)->AccelY * pBodies.at(i)->AccelY) * (20000 * _GRAV_CONST));

        if (AccelCoefficient > 1)
            AccelCoefficient = 1;

        float Red, Green, Blue;

        Blue = 1 - (AccelCoefficient);

        if (AccelCoefficient < 0.2)
            Red = (AccelCoefficient)* 5;
        else
            Red = 1;

        if (AccelCoefficient < 0.5)
            Green = (AccelCoefficient)* 2;
        else
            Green = 1;

        Temp.setFillColor(sf::Color(static_cast<sf::Uint8> (Red * 255), static_cast<sf::Uint8> (Green * 255), static_cast<sf::Uint8> (Blue * 255), 128));

        Temp.setPosition(static_cast<float>(pBodies.at(i)->posX), (static_cast<float>(pBodies.at(i)->posY)));
        pTarget->draw(Temp);
    }

    //DrawNode(&GlobalNode, pTarget);

    pTarget->display();
}

void DrawNode(Node* pNode, sf::RenderWindow* pTarget, float pZoom)
{
    sf::RectangleShape Temp;
    Temp.setFillColor(sf::Color(0, 0, 0, 0));
    Temp.setOutlineThickness(pZoom);
    Temp.setOutlineColor(sf::Color(0, 255, 0, 16));
    Temp.setPosition(pNode->posX, pNode->posY);
    Temp.setSize(sf::Vector2f(pNode->width, pNode->height));

    pTarget->draw(Temp);

    if (pNode->HasChildren) //recercivly draw all children
    {
        DrawNode(pNode->Child[0], pTarget, pZoom);
        DrawNode(pNode->Child[1], pTarget, pZoom);
        DrawNode(pNode->Child[2], pTarget, pZoom);
        DrawNode(pNode->Child[3], pTarget, pZoom);
    }
}

void SetView(sf::View* pView, sf::RenderWindow* pTarget, float pViewWidth, float pViewHeight)
{
    pView->reset(sf::FloatRect(SimWidth / 2 - pViewWidth / 2, SimHeight / 2 - pViewHeight / 2, pViewWidth, pViewHeight));
    pView->setViewport(sf::FloatRect(0.f, 0.f, 1.f, 1.f));
    pTarget->setView(*pView);
}

void UpdateBodies(std::vector<Body*> &pBodies)
{
    for (unsigned int i = 0; i < pBodies.size(); i++)
    {
        if ((pBodies[i]->posX > SimWidth && pBodies[i]->velX > 0) || (pBodies[i]->posX < 0 && pBodies[i]->velX < 0))
        {
            pBodies[i]->velX = -pBodies[i]->velX;
        }

        if ((pBodies[i]->posY > SimHeight && pBodies[i]->velY > 0) || (pBodies[i]->posY < 0 && pBodies[i]->velY < 0))
        {
            pBodies[i]->velY = -pBodies[i]->velY;
        }

        pBodies.at(i)->velX += pBodies.at(i)->AccelX;       
        pBodies.at(i)->velY += pBodies.at(i)->AccelY;

        pBodies.at(i)->posX += pBodies.at(i)->velX;
        pBodies.at(i)->posY += pBodies.at(i)->velY;
    }
}
\$\endgroup\$
1
  • 2
    \$\begingroup\$ While the sentiments are very much appreciated, they don't quite belong in posts. Since this is a follow-up, it would be helpful to provide the changes that were made. This could also prevent this post from possibly looking like a duplicate of the previous one. \$\endgroup\$
    – Jamal
    Commented Jul 7, 2015 at 6:27

1 Answer 1

5
\$\begingroup\$

Here are some observations that may help you improve your code.

Refine your objects

The Node and Body objects both have exclusively public members, which is not very good design. Every piece of code that interacts with them then reaches inside and manipulates the object members directly which is prone to error and makes maintenance difficult. Also, the entire simulation should be encapsulated into its own object.

Only expose the minimally sufficient interface to an object

Rather than exposing all data members, it's better to encapsulate data members and functionality within a class and to only make public those functions and data members that are required for the interface.

Use const where practical

Many of the function calls don't alter one or more of their parameters. Unaltered parameters and functions that don't alter the underlying object should be declared const to make clear what the interface may and may not do. This helps the compiler as well as the human reading the code.

Name things consistently

Using a consistent naming scheme makes it easier to understand your code. I generally like to reserve uppercase words like Node and Body as class names and member functions and data are lowercase. Also, one common scheme is to put a trailing underscore after member data names. That makes is easy to spot which ones are member variables and which are not.

Use constructors

The CreateBody function in main should instead be a constructor for Body. That way you can either create one on the stack or on the heap using new. It also means that any Body will be well-formed.

Use references to pass complex objects

Most of the functions in main which take a std::vector<Body*> as the first argument should instead be declared to take a reference instead. Doing so allows the compiler to avoid creating a copy of the variable. It is also better than using a raw pointer because unlike a raw pointer, a reference cannot dereference nullptr.

Use "range-for" to simplify your code

If your compiler fully supports the C++11 standard, you can use a "range-for" to simplify your code. This is often better, in terms of both readability and performance. If your compiler does not support C++11, I would strongly encourage you to upgrade to one that does.

Order conditionals for performance

The code currently contains this line

if (widthSqr / distanceSqr < NodeThresholdSqr || pNode->HasChildren == false)

Due to short-circuit evaluation, if the first clause is true the second does not need to be evaluated. This can be used for performance improvement by ordering the simpler case first:

if (!pNode->HasChildren || widthSqr / distanceSqr < NodeThresholdSqr)

Note also that rather than comparing to false, the boolean not operator ! is used. There is likely no performance difference, but it makes the line shorter and easier to read.

Use all required #includes

The program uses sqrt but does not #include <cmath>. Carefully check to assure that your program uses all required #includes and that it does not include ones that are not needed.

Consider using a better random number generator

If you are using a compiler that supports at least C++11, consider using a better random number generator. In particular, instead of rand, you might want to look at std::uniform_real_distribution and friends in the <random> header.

Use existing C++ libraries where appropriate

The PopulateVectorDisk routine creates a random distance and random angle and then converts those into discrete orthogonal (x and y) components for both the position and velocity of each point. Rather than doing all of that using primitive function calls with sin and cos, it would be much easier to read and more efficient to use the std::complex routines.

Eliminate memory leaks

There are a number of places where new is called without any corresponding delete calls. This is the recipe for a memory leak. Rational use of destructors can help with this if you've already implemented constructors.

Use calls that are in the std namespace

This code calls sqrt, cos and sin but should instead be using std::sqrt, std::cos and std::sin.

Don't overuse static_cast

In a number of places in the code, many static_cast instances are used where they are not needed. It's useful to know the numeric promotion rules so that you don't use too many casts. Doing so can clutter up the code and make it hard to read and maintain.

Use simple data types where possible

The Child member of the Node class is a std::vector<> but it really doesn't need to be. It is always ether empty or has four values, so for efficiency, it should be declared as an array of 4 Nodes.

Avoid calculations where practical

The SetParam routine is often called with the same set of Bodys as already exist within the Node object. Rather than redundantly copying and recalculating the center of mass and total mass, the code can simply check to see if they're equal and skip all of that if they are. Similarly, in the Render function, if a coefficient is < 0.2 it is certain that it is also < 0.5 so there's not much point in testing again.

Putting it all together

With all of these suggestions put together, the resulting code is around 10% faster on my machine, easier to read and understand, and shorter. Here are the components:

main.cpp

#include "SFML/Graphics.hpp"
#include "Node.h"
#include "Sim.h"

//Width and height of simulation, needs to be large, 
//particles outside of this range will not be included in the octree
float const SimWidth = 327680;
float const SimHeight = 327680;

void PollEvent(sf::RenderWindow& pTarget, bool &pIsPaused, sf::View &pSimView, float zoom)
{
    sf::Event event;

    while (pTarget.pollEvent(event))
    {
        if (event.type == sf::Event::Closed)
            pTarget.close();
        if (event.type == sf::Event::KeyPressed)
        {
            if (event.key.code == sf::Keyboard::Space)
                pIsPaused = !pIsPaused;
        }
        if (event.type == sf::Event::MouseWheelScrolled)
        {
            zoom *= 1 + (static_cast <float> (-event.mouseWheelScroll.delta) / 10); 
            pSimView.zoom(1 + (static_cast <float> (-event.mouseWheelScroll.delta) / 10));
        }
    }

    if (sf::Mouse::getPosition().x > (1920 - 20))
        pSimView.move(2 * zoom, 0);
    if (sf::Mouse::getPosition().x < (0 + 20))
        pSimView.move(-2 * zoom, 0);
    if (sf::Mouse::getPosition().y > (1080 - 20))
        pSimView.move(0, 2 * zoom);
    if (sf::Mouse::getPosition().y < (0 + 20))
        pSimView.move(0, -2 * zoom);

    pTarget.setView(pSimView);
}

void SetView(sf::View &pView, sf::RenderWindow& pTarget, float pViewWidth, float pViewHeight)
{
    pView.reset(sf::FloatRect(SimWidth / 2 - pViewWidth / 2, SimHeight / 2 - pViewHeight / 2, pViewWidth, pViewHeight));
    pView.setViewport(sf::FloatRect(0.f, 0.f, 1.f, 1.f));
    pTarget.setView(pView);
}

int main()
{
    float const ViewWidth = 1920;
    float const ViewHeight = 1080;

    float const DiskRadiusMax = 20000;
    float const DiskRadiusMin = 50;
    float const GalacticCenterMass = 1000000;
    float const ObjectMassMax = 2;
    float const ObjectMassMin = 1;
    unsigned int const NumParticles = 10000;

    float zoom = 1;
    bool IsPaused = false;

    sf::RenderWindow window(sf::VideoMode(1920, 1080), "N-Body simulation");
    sf::View SimulationView;                        

    Sim sim(NumParticles, SimWidth, SimHeight, DiskRadiusMin, 
            DiskRadiusMax, ObjectMassMin, ObjectMassMax, GalacticCenterMass);
    SetView(SimulationView, window, ViewWidth, ViewHeight);

    while (window.isOpen())
    {
        PollEvent(window, IsPaused, SimulationView, zoom); 
        if (!IsPaused)  
            ++sim;
        sim.Render(window, zoom);
    }
}

Sim.h

#ifndef SIM_H
#define SIM_H
#include "SFML/Graphics.hpp"
#include "Node.h"
#include <vector>

class Sim {
public:
    Sim(unsigned int pParticlesCount, float pWidth, float pHeight, float pMaxDist, float pMinDist, float pMinMass, float pMaxMass, float pGalaticCenterMass);
    Sim &operator++() {
        attractToCenter();
        updateBodies();
        globalNode_.setParam(bodies_, simWidth, simHeight);
        octreeBodyAttraction();
        return *this;
    }   
    void Render(sf::RenderWindow& pTarget, float zoom) const;
    ~Sim();
private:    
    void attractToCenter() {
        for (auto &pB : bodies_)
            pB->calculateAccel(center_);
    }
    void updateBodies() {
        for (auto &pB : bodies_) {
            pB->update(simWidth, simHeight);
            pB->resetAccel();
        }
    }
    void octreeBodyAttraction() {
        for (auto &pB : bodies_)
            globalNode_.checkNode(*pB);
    }
    float simWidth, simHeight, galacticCenterMass;
    Body center_;
    std::vector<Body*> bodies_;
    Node globalNode_;
};
#endif // SIM_H

Node.h

#pragma once
#include <vector>
#include <cmath>

extern const float _GRAV_CONST;
class Node;

class Body
{
public:
    Body(float px, float py, float pmass, float pvx=0, float pvy=0) :
        posX_(px), 
        posY_(py),
        velX_(pvx),
        velY_(pvy),
        accelX_(0),
        accelY_(0),
        mass_(pmass)
    {}

    // updates accel based on other body
    void calculateAccel(const Body& bj);  
    // updates accel based on other node
    void calculateAccel(const Node& bj);
    void resetAccel() { accelX_ = accelY_ = 0; }
    void update(float SimWidth, float SimHeight);
    float dist2(const Node& bj) const;  
    float mass() const { return mass_; }
    float posX() const { return posX_; }
    float posY() const { return posY_; }
    float accel() const { return std::sqrt(accelX_*accelX_ + accelY_*accelY_); }

private:
    void calculateAccel(double x, double y, double mass);

    float posX_, posY_;           //position x and y
    float velX_, velY_;           //velocity x and y
    double accelX_, accelY_;      //accel acting on object since last frame
    float mass_;                 //mass of object
};

class Node
{
public:
    Node() : 
        child_{nullptr, nullptr, nullptr, nullptr},
        depth_{0},
        hasChildren_{false}
    {}
    Node(unsigned int pDepth) : 
        child_{nullptr, nullptr, nullptr, nullptr},
        depth_{pDepth},
        hasChildren_{false}
    {}
    void setParam(const std::vector<Body*> &pBodies, float pwidth, float pheight, float px = 0, float py = 0);
    bool hasChildren() const { return hasChildren_; }
    ~Node() { 
        hasChildren_ = false;
        delete child_[0];
        delete child_[1];
        delete child_[2];
        delete child_[3];

        child_[0] = nullptr;
        child_[1] = nullptr;
        child_[2] = nullptr;
        child_[3] = nullptr;
    }

    void checkNode(Body& pBody);
    friend class Body;
private:
    void generateChildren();

    Node* child_[4];
    float posX_, posY_;
    float width_, height_;
    std::vector<Body*> bodies_;
    float totalMass_;
    float centerOfMassx_;
    float centerOfMassy_;
    unsigned int depth_;
    bool isUsed_;            //For testing, delete this later
    bool hasChildren_;
};

Sim.cpp

#include "SFML/Graphics.hpp"
#include "Node.h"
#include "Sim.h"
#include <random>
#include <complex>
#include <cmath>

const float TWO_PI = 2*3.14159265;
extern const float _GRAV_CONST = 0.1;

Sim::~Sim()
{
    for (auto &pB : bodies_)
        delete pB;
    bodies_.clear();
}

Sim::Sim(unsigned int pParticlesCount, float pWidth, float pHeight, 
        float pMinDist, float pMaxDist, float pMinMass, float pMaxMass, 
        float pGalacticCenterMass) : 
    simWidth(pWidth),
    simHeight(pHeight),
    galacticCenterMass(pGalacticCenterMass),
    //Create a body at the center of the simulation
    center_(simWidth / 2, simHeight / 2, galacticCenterMass)
{
    bodies_.reserve(pParticlesCount);
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<float> dist(0,1); //pMinDist, pMaxDist);
    std::uniform_real_distribution<float> angle(0, TWO_PI);
    std::uniform_real_distribution<float> mass(pMinMass, pMaxMass);
    const std::complex<float> ctr(simWidth/2, simHeight/2);

    for (unsigned int i = 0; i < pParticlesCount; i++)
    {
        float theta = angle(gen);
        float distance = dist(gen);
        // beta distribution with alpha=0.5, beta=1.0
        distance = distance*distance * (pMaxDist - pMinDist) + pMinDist;
        float orbitalVelocity = std::sqrt(galacticCenterMass * _GRAV_CONST / distance);
        std::complex<float> pos = std::polar(distance, theta)+ctr;
        std::complex<float> vel = std::polar(orbitalVelocity, theta-TWO_PI/4);

        bodies_.push_back(new Body(pos.real(), pos.imag(), mass(gen), vel.real(), vel.imag()));                                   
    }
}

void Sim::Render(sf::RenderWindow& pTarget, float zoom) const
{
    pTarget.clear();

    sf::RectangleShape temp;

    for (auto &pB : bodies_)
    {       
        auto a = zoom > 1 ? zoom*pB->mass() : pB->mass();
        temp.setSize(sf::Vector2f(a, a));

        float accelCoefficient = pB->accel() * 20000 * _GRAV_CONST;

        float red, green, blue;

        if (accelCoefficient > 1) {
            red = 5*255; 
            green = 2*255; 
            blue = 0*255;
        } else {
            blue = 1 - accelCoefficient;
            if (accelCoefficient < 0.2) {
                red = accelCoefficient * 5*255;
                green = accelCoefficient * 2*255;
            }
            else {
                red = 1*255;
                if (accelCoefficient < 0.5)
                    green = accelCoefficient * 2*255;
                else
                    green = 1*255;
            }
        }

        temp.setFillColor(sf::Color(red, green, blue, 128));
        temp.setPosition(pB->posX(), pB->posY());
        pTarget.draw(temp);
    }

    pTarget.display();
}

Node.cpp

#include "Node.h"
#include <cmath>

unsigned const _MAX_DEPTH = 100;
double const _NODE_THRESHOLD_SQ = 0.5 * 0.5;
float const Softener = 10;
const float Softener2 = Softener * Softener;

void Node::generateChildren()
{
    std::vector<Body*> q1Bodies, q2Bodies, q3Bodies, q4Bodies;

    float halfWidth = posX_ + width_ / 2;
    float halfHeight = posY_ + height_ / 2;

    for (const auto &pB : bodies_)
    {
        if (pB->posX() < halfWidth)
        {
            if (pB->posY() < halfHeight)
                q1Bodies.push_back(pB);
            else 
                q3Bodies.push_back(pB);
        }
        else
        {
            if (pB->posY() < halfHeight)
                q2Bodies.push_back(pB);
            else
                q4Bodies.push_back(pB);
        }
    }

    if (child_[0] == nullptr) {
        child_[0] = new Node();
        child_[1] = new Node();
        child_[2] = new Node();
        child_[3] = new Node();
    }
    child_[0]->depth_ = child_[1]->depth_ = child_[2]->depth_ = child_[3]->depth_ = depth_ + 1;

    child_[0]->setParam(q1Bodies, width_ / 2, height_ / 2, posX_, posY_);
    child_[1]->setParam(q2Bodies, width_ / 2, height_ / 2, halfWidth, posY_);
    child_[2]->setParam(q3Bodies, width_ / 2, height_ / 2, posX_, halfHeight);
    child_[3]->setParam(q4Bodies, width_ / 2, height_ / 2, halfWidth, halfHeight);

    hasChildren_ = true;
}

void Node::setParam(const std::vector<Body*> &pBodies, float pwidth, float pheight, float px, float py)
{
    if (bodies_ != pBodies) {
        bodies_ = pBodies;
        totalMass_ = 0;
        centerOfMassx_ = 0;
        centerOfMassy_ = 0;
        for (const auto &pB : bodies_)
        {
            totalMass_ += pB->mass();
            centerOfMassx_ += pB->posX();
            centerOfMassy_ += pB->posY();
        }
        unsigned int size = bodies_.size();
        centerOfMassx_ /= size;
        centerOfMassy_ /= size;
    }
    posX_ = px;
    posY_ = py;
    width_ = pwidth;
    height_ = pheight;
    hasChildren_ = false;
    if (bodies_.size() > 1 && depth_ < _MAX_DEPTH)
    {
        generateChildren();
    }
}

void Node::checkNode(Body& pBody)
{
    if (bodies_.size() == 0)
        return;
    if (!hasChildren_ ||
        (width_*width_ / pBody.dist2(*this) < _NODE_THRESHOLD_SQ))
    {   
        pBody.calculateAccel(*this);
        isUsed_ = true;
    }
    else
    {
        child_[0]->checkNode(pBody);
        child_[1]->checkNode(pBody);
        child_[2]->checkNode(pBody);
        child_[3]->checkNode(pBody);
    }
}

void Body::calculateAccel(const Body& bj)
{
    calculateAccel(bj.posX_, bj.posY_, bj.mass_);
}

void Body::calculateAccel(const Node& bj)
{
    calculateAccel(bj.centerOfMassx_, bj.centerOfMassx_, bj.totalMass_);
}

void Body::calculateAccel(double x, double y, double mass)
{
    double vectorx = x - posX_;
    double vectory = y - posY_;

    //c^2 = a^2 + b^2 + softener^2
    double distSqr = vectorx * vectorx + vectory * vectory + Softener2;
    double distSixth = distSqr * distSqr * distSqr;
    double accel = mass / std::sqrt(distSixth) * _GRAV_CONST;

    accelX_ += vectorx * accel;
    accelY_ += vectory * accel;
}

void Body::update(float SimWidth, float SimHeight) {
    if ((posX_ > SimWidth && velX_ > 0) || (posX_ < 0 && velX_ < 0))
    {
        velX_ = -velX_;
    }
    if ((posY_ > SimHeight && velY_ > 0) || (posY_ < 0 && velY_ < 0))
    {
        velY_ = -velY_;
    }

    velX_ += accelX_;
    velY_ += accelY_;

    posX_ += velX_;
    posY_ += velY_;
}

float Body::dist2(const Node& bj) const
{
    float vectorx = bj.centerOfMassx_ - posX_;
    float vectory = bj.centerOfMassy_ - posY_;
    return vectorx * vectorx + vectory * vectory;
}

Further work

Use more efficient data structures

The setParam and generateChildren functions are not as efficient as they could be. Specifically, rather than recreating the four quadrant vectors every time, it would yield better performance to have a single list and to have the quadrant be an annotation for each Body member.

Use parallel threads

The C++11 and C++14 standards include simpler and more portable threading mechanisms. Because this code is largely parallel, it would not be too difficult to make it multithreaded and significantly boost performance.

Use native code

At the expense of portability, another approach might be to use native instructions to speed the code. The use of assembly language (especially SIMD) instruuctions would likely be quite effective in speeding the code.

Profile your code

Measuring your code using a profiler is one of the most effective ways to discover which routines are slow and to measure the effects of various things one might try to improve the performance of the code.

\$\endgroup\$
2
  • \$\begingroup\$ How is the best way to impliment multithreading? Is there a libary I need to use or is it part of the standard? \$\endgroup\$ Commented Jul 8, 2015 at 0:35
  • \$\begingroup\$ std::thread is part of the C++11 standard. That's how I'd do it. \$\endgroup\$
    – Edward
    Commented Jul 8, 2015 at 0:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.