No Description

engine.c 5.5KB

    /* 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); } } } }