Browse Source

MASSIVE performance improvement(barnes hut algorithm), better colours when converting to img

master
Stephen Downward 1 year ago
parent
commit
6b05c35d46
4 changed files with 352 additions and 77 deletions
  1. +2
    -2
      Makefile
  2. +235
    -0
      engine.c
  3. +36
    -0
      engine.h
  4. +79
    -75
      main.c

+ 2
- 2
Makefile View File

@@ -1,5 +1,5 @@
galaxysim: main.c qdbmp.c qdbmp.h
gcc -o galaxysim main.c qdbmp.c -lm -fopenmp
galaxysim: main.c qdbmp.c qdbmp.h engine.c engine.h
gcc -o galaxysim main.c qdbmp.c engine.c -lm -fopenmp

collide: collide.c
gcc -o collide collide.c

+ 235
- 0
engine.c View File

@@ -0,0 +1,235 @@
/* Implementation of the Barnes-Hut algorithm */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
#include <unistd.h>
#include <stdbool.h>
#include "engine.h"


#define THETA 0.5

void initTree(struct node_t *tree, double xmin, double xmax, double ymin, double ymax);
void calcForces(struct node_t *tree, RigidMass_t *mass);

//need a way to deal with masses that are on top of each other
void addToTree(struct node_t *tree, RigidMass_t *masses, int n, int rlev)
{
//printf("n=%d\r\n", n);
if(n == 2)
{
rlev ++; //debuggign only
}
RigidMass_t *masses_topleft;
int n_topleft = 0;

RigidMass_t *masses_topright;
int n_topright = 0;
RigidMass_t *masses_botleft;
int n_botleft = 0;
RigidMass_t *masses_botright;
int n_botright = 0;
masses_topleft = malloc(n * sizeof(RigidMass_t));
masses_topright = malloc(n * sizeof(RigidMass_t));
masses_botleft = malloc(n * sizeof(RigidMass_t));
masses_botright = malloc(n * sizeof(RigidMass_t));

tree->n += n;

double centerX = (tree->xmax + tree->xmin) / 2;
double centerY = (tree->ymax + tree->ymin) / 2;

// printf("Center of this node: (%f, %f)", centerX, centerY);
double totalMass = 0;
double weightedCenterX = 0;
double weightedCenterY = 0;
for(int i = 0; i < n; i++)
{
//add mass to tree stats
totalMass += masses[i].mass;
weightedCenterX += masses[i].mass * masses[i].x;
weightedCenterY += masses[i].mass * masses[i].y;

//categorize mass for child node
double curX = masses[i].x;
double curY = masses[i].y;
if(curX < centerX)
{
//Left side

if(curY < centerY)
{
//Top left
masses_topleft[n_topleft++] = masses[i];
}
else
{
//Bot left
masses_botleft[n_botleft++] = masses[i];
}
}
else
{
if(curY < centerY)
{
//Top right
masses_topright[n_topright++] = masses[i];
}
else
{
//Bot right
masses_botright[n_botright++] = masses[i];
}
}
}

tree->CoM_x = weightedCenterX / totalMass;
tree->CoM_y = weightedCenterY / totalMass;
tree->mass = totalMass;
//v printf("Distribution: %d, %d, %d, %d \r\n", n_topleft, n_topright, n_botleft, n_botright);
//printf("%d\r\n", rlev);
//recursion
if(rlev < 100) //can't go more than 100 levels deep.
{
if(n_topleft > 1)
{
// printf("Down an inception level\r\n");
tree->nodes[0] = malloc(sizeof(struct node_t));
initTree(tree->nodes[0], tree->xmin, centerX, tree->ymin, centerY);
addToTree(tree->nodes[0], masses_topleft, n_topleft, rlev + 1);
}

if(n_topright > 1)
{
tree->nodes[1] = malloc(sizeof(struct node_t));
initTree(tree->nodes[1], centerX, tree->xmax, tree->ymin, centerY);
addToTree(tree->nodes[1], masses_topright, n_topright, rlev + 1);
}

if(n_botleft > 1) //seems broken
{
tree->nodes[2] = malloc(sizeof(struct node_t));
initTree(tree->nodes[2], tree->xmin, centerX, centerY, tree->ymax);
addToTree(tree->nodes[2], masses_botleft, n_botleft, rlev + 1);
}

if(n_botright > 1)
{
tree->nodes[3] = malloc(sizeof(struct node_t));
initTree(tree->nodes[3], centerX, tree->xmax, centerY, tree->ymax);
addToTree(tree->nodes[3], masses_botright, n_botright, rlev + 1);
}
}

free(masses_topleft);
free(masses_topright);
free(masses_botleft);
free(masses_botright);
}

void initTree(struct node_t *tree, double xmin, double xmax, double ymin, double ymax)
{
tree->mass = 0;
tree->n = 0;
tree->xmin = xmin;
tree->xmax = xmax;
tree->ymin = ymin;
tree->ymax = ymax;
tree->nodes[0] = 0;
tree->nodes[1] = 0;
tree->nodes[2] = 0;
tree->nodes[3] = 0;
}

void delTree(struct node_t *tree, bool isParent)
{
if(tree == 0) return;
//printf("delTree recursion!\r\n");
delTree(tree->nodes[0], false);
delTree(tree->nodes[1], false);
delTree(tree->nodes[2], false);
delTree(tree->nodes[3], false);
if(!isParent)
{
free(tree);
//printf("tree freed!\r\n");
}
}

void barnesIterate(RigidMass_t *masses, int n, int w, int h)
{
printf("Attempting to make a tree...\r\n");
struct node_t tree;
//quad tree encapsulates an area of 5w * 5h, centered on the screen view
initTree(&tree, -2*w, 3*w, -2*h, 3*h);
addToTree(&tree, masses, n, 0);
printf("Tree done. Calculating forces...\r\n");

#pragma omp parallel for
for(int i = 0; i < n; i++)
{
calcForces(&tree, &masses[i]);
}
printf("Deleting tree...\r\n");
delTree(&tree, true);
printf("Processing tick..\r\n");
for(int i = 0; i < n; i++)
{
masses[i].vx += masses[i].dvx;
masses[i].vy += masses[i].dvy;
masses[i].dvx = 0;
masses[i].dvy = 0;
masses[i].x += masses[i].vx * masses[i].tick_time;
masses[i].y += masses[i].vy * masses[i].tick_time;
}
}

const double G = 0.00667408;

void calcForces(struct node_t *tree, RigidMass_t *mass)
{
double dx = tree->CoM_x - mass->x;
double dy = tree->CoM_y - mass->y;
//done in many parts b/c we need dist_squared and dist when approxing anyway
double dist_squared = dx*dx + dy*dy;
double dist = sqrt(dist_squared);
double sd = (tree->xmax - tree->xmin) / dist;

bool approx = true;
//if all of the children are null ptrs, we can not go any lower
for(int i = 0; i < 4; i++)
{
if(tree->nodes[i] != 0) approx = false;
}

if(sd < THETA || approx) //we can approximate
{
//printf("s/d = %f\r\n", sd);
if(dist > 1)
{
double v = G * tree->mass / dist_squared * mass->tick_time;
mass->dvx += v * (dx) / dist;
mass->dvy += v * (dy) / dist;
}
}
else
{
//go down a level
for(int i = 0; i < 4; i++)
{
if(tree->nodes[i] != 0)
{
calcForces(tree->nodes[i], mass);
}
}
}
}

+ 36
- 0
engine.h View File

@@ -0,0 +1,36 @@
#ifndef ENGINE_H
#define ENGINE_H

typedef struct RigidMass
{
long mass;
double vx, vy;
double x, y;
double dvx, dvy;
double tick_time;
} RigidMass_t;

struct node_t
{
long mass;
int n; //num masses
double xmin;
double xmax;
double ymin;
double ymax;
double CoM_x;
double CoM_y;
//0 - top left
//1 - top right
//2 - bottom left
//3 - bottom right
struct node_t *nodes[4];
};

extern const double G;

void barnesIterate(RigidMass_t *masses, int n, int w, int h);

#endif

+ 79
- 75
main.c View File

@@ -3,52 +3,11 @@
#include <math.h>
#include <time.h>
#include <string.h>
#include <stdbool.h>
#include <unistd.h>
#include "qdbmp.h"
#include "engine.h"

typedef struct RigidMass
{
long mass;
double vx, vy;
double x, y;
double dvx, dvy;
double tick_time;
} RigidMass_t;

double G = 0.00667408;

//Calculates delta V on body, due to other.
void calcGravity(RigidMass_t *body, RigidMass_t other)
{
double dx = body->x - other.x;
double dy = body->y - other.y;
double dist_squared = dx * dx + dy * dy;

double v = -G * other.mass / dist_squared * body->tick_time;

double hyp = sqrt(dist_squared);
body->dvx += v * (dx) / hyp;
body->dvy += v * (dy) / hyp;
}

void processTick(RigidMass_t *body)
{
//technically, we could just remove dvx/dvy
//and work with vx/vy directly, since they don't change anything other than when
//processing the tick here
body->vx += body->dvx;
body->vy += body->dvy;
body->dvx = 0;
body->dvy = 0;
body->x += body->vx * body->tick_time;
body->y += body->vy * body->tick_time;
}

//END PHYSICS ENGINE

//stackoverflow stuff
/* reverse: reverse string s in place */
void reverse(char s[])
{
int i, j;
@@ -103,14 +62,9 @@ char s[255]; //save filename
char l[255]; //load filename
/* END DEFAULT CONFIGURATION OPTIONS */


//booleans are cool
typedef int bool;
#define true 1
#define false 0

//Allocate this array later!
//Allocate these arrays later!
RigidMass_t *stars;
int *starsAtPixel;

void generateStars()
{
@@ -151,40 +105,85 @@ void generateStars()
}
}

void iterate()
void hsl_to_rgb(int h, double s, double l, char *rOut, char *gOut, char *bOut)
{
#pragma omp parallel for
for(int a = 0; a < nStars; a++)
if(h == 360) h = 0;
if(s > 1) s = 1;
if(l > 1) l = 1;
double c = (1 - abs(2 * l - 1)) * s;
double x = c * (1 - fabs(fmod(((double)h)/60, 2) - 1));
double m = l - c / 2;
double r, g, b;
int q = h / 60;
switch(q)
{
for(int b = 0; b < nStars; b++)
{
//If stars are closer than 0.1 units apart, we don't calculate the force on each other
//this is important since we're stepping through time. If we end up with two stars ridiculously
//close together, the math will act as if they were in that position for the entire tick,
//whereas they would likely only be there for a small fraction of that.
//we don't want any stars to be flung away at a fraction of the speed of light!
if(abs(stars[a].x - stars[b].x) > 0.1 && abs(stars[a].y - stars[b].y) > 0.1)
{
//printf("%d, %d \r\n", stars[a].dvx, stars[a].dvy);
calcGravity(stars + a, stars[b]);
}
}
case(0):
r = c;
g = x;
b = 0;
break;
case(1):
r = x;
g = c;
b = 0;
break;
case(2):
r = 0;
g = c;
b = x;
break;
case(3):
r = 0;
g = x;
b = c;
break;
case(4):
r = x;
g = 0;
b = c;
break;
case(5):
r = c;
g = 0;
b = x;
}

//'commit' the deltaV
for(int a = 0; a < nStars; a++)
{
processTick(stars + a);
//printf("%f, %f \r\n", stars[a].x, stars[a].y);
}
*rOut = 255 * (r + m);
*gOut = 255 * (g + m);
*bOut = 255 * (b + m);
}

BMP *out;
void addToImage()
{
memset(starsAtPixel, 0, sizeof(int) * w * h);
#define IMAGE_INC_AMOUNT 50
for(int a = 0; a < nStars; a++)
{
BMP_SetPixelRGB(out, stars[a].x, stars[a].y, 0, 150, 0);
int starX = stars[a].x;
int starY = stars[a].y;
if(starX > 0 && starX < w && starY > 0 && starY < h)
{
starsAtPixel[starX * h + starY]++;
}
}
for(int a = 0; a < nStars; a++)
{
//pixel colour depends on how many stars are at that location
int starX = stars[a].x;
int starY = stars[a].y;
if(starX > 0 && starX < w && starY > 0 && starY < h)
{
int nStarsAtLoc = starsAtPixel[starX * h + starY] * 10;
//if(nStarsAtLoc > 360) nStarsAtLoc = 360;
char curR = 0, curG = 0, curB = 0;
if(nStarsAtLoc > 50) nStarsAtLoc = 50;
curR = nStarsAtLoc * 5;
curG = curR;
curB = curR;
//hsl_to_rgb(nStarsAtLoc, 1, nStarsAtLoc / (double)100, &curR, &curG, &curB);
BMP_SetPixelRGB(out, starX, starY, curR, curG, curB);
}
}
}

@@ -324,8 +323,9 @@ int main(int argc, char **argv)
printf("\t m = %ldE8 kg\r\n", m);
printf("\t Autosave: %s-ID.sgs every %d iterations\r\n", s, a);

//allocate array
//allocate arrays
stars = malloc(sizeof(RigidMass_t) * nStars);
starsAtPixel = malloc(sizeof(int) * w * h);

srand(time(NULL));
generateStars();
@@ -365,7 +365,11 @@ int main(int argc, char **argv)
bool justStarted = true;
while(1)
{
iterate();
//iterate();
printf("Started gentree\n");
barnesIterate(stars, nStars, w, h);
printf("Gentree finished\n");
//exit(1);
addToImage();
//gives the next frame some iterations to build up
//without justStarted, the loaded iteration is rerendered, but only from


Loading…
Cancel
Save