372Project/physics.c at main · spoduminister/372Project · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
#include "physics.h"
#include <math.h>
//should move by one physics step

float vec_dot(vec_t a, vec_t b)
{
    float product_sum = 0.0f;

    for(int d = 0; d < dims; d++)
    {
        product_sum += a.v[d] * b.v[d];
    }
    return product_sum;
}
vec_t vec_sub(vec_t a, vec_t b)
{
    vec_t diff;

    for(int d = 0; d < dims; d++)
    {
        diff.v[d] = a.v[d] - b.v[d];
    }
    return diff;
}
vec_t vec_add(vec_t a, vec_t b)
{
    vec_t sum;

    for(int d = 0; d < dims; d++)
    {
        sum.v[d] = a.v[d] + b.v[d];
    }
    return sum;
}

void physics_step(Particle * particle_list, int numParticles, Box * Boxes, float timestep)
{
    while(timestep > 0)
    {
        //TODO parallel find collision time step
        float collision_step = timestep;

        float tempcoll = 0;
        float dist;
	int index[2];

	// UNCOMMENT BELOW FOR PARALLEL MULTICORE ONLY
	#pragma acc parallel loop
        for (int i = 0; i < numParticles; i++) {
            int pnum = 0;
            vec_t min;
            vec_t max;

            for(int j = 0; j < dims; j++)
		    {
			    max.v[j] = fmax(particle_list[i].pos.v[j], particle_list[i].pos.v[j] + (particle_list[i].vel.v[j]*timestep)) + particle_list[i].radius;
			    min.v[j] = fmin(particle_list[i].pos.v[j], particle_list[i].pos.v[j] + (particle_list[i].vel.v[j]*timestep)) - particle_list[i].radius;
		    }
            int* colls = get_within_bounds(Boxes, 0, min, max, &pnum);

	    if (colls != NULL) {

		// UNCOMMENT BELOW FOR PARALLEL GPU ONLY
		// #pragma acc data copy(colls[0:pnum], particle_list[0:numParticles], index[0:2]) create(colls[0:pnum], particle_list[0:numParticles], index[0:2])
		// #pragma acc parallel loop private(collision_step)
                for (int j = 0; j < pnum; j++) {
                    if (colls[j] >= i) continue;

                    float squareR = pow((particle_list[i].radius + particle_list[colls[j]].radius),2); // Sum of particle radius' squared
                    vec_t posdiff = vec_sub(particle_list[i].pos, particle_list[colls[j]].pos); // vector difference in particle positions
                    vec_t veldiff = vec_sub(particle_list[i].vel, particle_list[colls[j]].vel); // vector difference in particle velocities

                    float pos_dot = vec_dot(posdiff,posdiff);
                    float pos_vel_dot = vec_dot(posdiff,veldiff);
                    float vel_dot = vec_dot(veldiff,veldiff);

                    float t = pow((pos_vel_dot/2.0),2.0) - (vel_dot * (pos_dot - squareR));

                    if (t < 0)continue;

                    t -=  pos_vel_dot;
                    t /= vel_dot;
                    if (t < collision_step) {
                        index[0] = i;
                        index[1] = colls[j];
                        collision_step = t;
                    }

                }
            }

        }

        //step collision timestep
        //TODO parallel step collision time step
	// UNCOMMENT DATA CLAUSE FOR PARALLEL GPU
	// #pragma acc data copy(particle_list[0:numParticles])
	#pragma acc parallel loop
        for(int i = 0; i < numParticles; i++)
        {
            for(int j = 0; j < dims; j++)
            {
                particle_list[i].pos.v[j] += particle_list[i].vel.v[j]* collision_step;
            }
        }

        //TODO actually resolve collisions
        //resolve collisions
        if (collision_step < timestep) {
            float v1 = (float)sqrt(pow((double)particle_list[index[0]].vel.v[0], 2) + pow((double)particle_list[index[0]].vel.v[1], 2));
            float v2 = (float)sqrt(pow((double)particle_list[index[1]].vel.v[0], 2) + pow((double)particle_list[index[1]].vel.v[1], 2));
            //final x direction of particle 1; taking a random direction and making it a unit vector
            float scalar = (float)rand() / (float)(RAND_MAX);
            float fvx1_dir = (scalar * -1) * particle_list[index[0]].vel.v[0] / (particle_list[index[0]].vel.v[0] + particle_list[index[0]].vel.v[1]);
            float fvy1_dir = ((1-scalar) * -1) * particle_list[index[0]].vel.v[1] / (particle_list[index[0]].vel.v[0] + particle_list[index[0]].vel.v[1]);
            float fvx2_dir = (scalar * -1) * particle_list[index[1]].vel.v[0] / (particle_list[index[0]].vel.v[0] + particle_list[index[0]].vel.v[1]);
            float fvy2_dir = ((1 - scalar) * -1) * particle_list[index[1]].vel.v[1] / (particle_list[index[0]].vel.v[0] + particle_list[index[0]].vel.v[1]);
            particle_list[index[0]].vel.v[0] = fvx1_dir * v1;
            particle_list[index[0]].vel.v[1] = fvy1_dir * v1;
            particle_list[index[1]].vel.v[0] = fvx2_dir * v2;
            particle_list[index[1]].vel.v[1] = fvy2_dir * v2;
        }

        timestep -= collision_step;
    }

    return;
}