ExternalVendorCode/trick
1
0
Fork 0
mirror of https://github.com/nasa/trick.git synced 2024-12-24 07:16:41 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Trick Smoothed Particle Hydrodynamics (TSPH) (#1170)

Browse Source 
* Create base simulation file/directory structure

* Create particle class in a Trick header file

* Create Fluid class with simulation initialization and update jobs in Trick header

* Create S_define for SIM_fluid

* Create S_overrides.mk for SIM_fluid

* Change Particle to a struct and initialize pos, velocity, and force arrays in constructor

* Store particles as an std::vector instead of an array

* Implement default_data method and dummy update_SPH for Fluid.cpp

* Add core CUDA SPH simulation code to Trick simulation

* Copy CPU-based SPH fluid sim code to Trick src directory

* Remove all code involving rigid body simulation

* Modify updateSPH routine to use CPU-based routines rather than GPU-based routines

* Temporarily comment out neighbor list logic

* Replace use of glm::vec3 with float[3] in sph.cpp and sph.h

* Move simulation parameters from sph.h to Trick Fluid.hh

* Remove particle struct and function prototypes from sph.h

* Move core simulation routines from sph.cpp to Fluid.cpp

* Remove neighbor list comments from Fluid.cpp

* Remove extraneous initSPH function

* Add calls to core simulation routines in Fluid::update_SPH()

* Refactor core simulation functions into private member functions of the Fluid class

* Remove static qualifier from simulation parameters in Fluid.hh

* Include <cmath> and use std scope resolution operater on pow

* Comment out isnan check in computeForces for now

* Initialize p_start and p_end in update_SPH

* Remove const qualifier from simulation paramters in Fluid class

* Remove n_particles member variable from Fluid class

* Create runtime.py to enable real-time synchronization and sim control panel

* Create input.py to exec realtime.py and limit simulation time

* Create simple python variable server client to unit test connection by accessing VISC member of fluid

* Create unit test to see if particle position is being updated correctly

* Create pointer to particles vector data for use in variable client

* Implement simple 2D GUI in order to render a single fluid particle

* Temporarily comment out unit tests

* Support 2D rendering of multiple particles in python variable client

* Scale size of simulation GUI based on BOUND variable from Trick variable server

* Comment out trick.stop(5) temporarily

* Set default value of BOUND to 400

* Fix sideways rendering bug by indexing field 1 position earlier

* Start variable server client from input.py rather than from command line

* Determine number of particles to render based on value of NUM_PARTICLES from Trick Variable Server. Also added TODO for correctly indexing into field.

* Reduce number of particles for 2D simulation

* Add CUDA runtime library flag and include directories

* Rename .cu/.cuh files to .cpp/.h

* Add sph_gpu.o as a library dependency of Fluid.cpp and call updateSPH_GPU in update_SPH

* Successfully compile and link Trick with CUDA in S_overrides.mk

* Call CUDA wrapper callVectorAdd() in update_SPH() to test Trick compilation/linking

* Create vectorAdd.cu to test CUDA kernel with Trick

* Change sph_gpu to a .cu file

* Pass pointer to Fluid object into updateSPH_GPU in order to access fluid parameters on GPU

* Refactor comptueDensityAndPressureGPU to use float[3] and fluid pointer to access sim parameters

* Refactor computeForcesGPU to use float[3] and fluid pointer to access sim parameters

* Refactor timeIntegrationGPU to use float[3] and fluid pointer to access sim parameters

* Refactor verletUpdatePosition to use float[3] and fluid pointer to access sim parameters

* Add base code for OpenGL particle renderer and TCP socket communication

* Create .gitignore(s) and delete built files

* Replace char* with char[] to remove warnings

* Receive and parse number of particles and particle positions

* Move connection initialization to stupComm

* Move TCP client to graphics folder and use it within OpenGL particle renderer

* Render particle position data from Trick

* Distinguish between trick and custom client

* Attempt to fix variable server client freezing

* Create input file that doesn't launch python VSC

* Remove extraneous call to initSPH

* Modify number of particles to 2048

* Create .gitignore for sim

* Create FluidServer class header

* Add Trick Header to FluidServer.hh and add clientfd member to FluidServer

* Send simple message from FluidServer in Trick to custom_client.cpp via TCP

* Add FluidServer to S_define

* Delete client directory

* Send particle data on initialization and decrease frequency of sendParticleData scheduled job

* Implement sendParticleData() to stream particle positions to client

* Implement client to receive particle position data from custom TCP fluid server

* Use custom_client interface in particle renderer of Trick variable server client

* Pass fluid pointer as parameter to FluidServer::sendParticleData(Fluid* fluid)

* Remove sim binary

* Remove TCP client usage from OpenGL particle renderer

* Rename particle renderer main.cc to opengl.cc

* Create openGLCaller, which spawns a thread that handles the OpenGL particle rendering

* Change openGLCaller signature

* Remove FluidServer usage from S_define

* Add make targets to OpenGL particle renderer source files and include OpenGL libraries in S_overrides.mk

* Temporarily set argc and argv in openGLCaller and change signature of openGLCaller

* Rename Particle to ParticleGL in sph.cpp and sph.h to avoid naming conflicts

* Hard code .obj file name in openGLMain rather than passing it as a command line arg

* Call openGLCaller inside of Fluid::default_data()

* Fix compilation bug with g++ and add .o files to nvcc device linking and TRICK_USER_LINK_LIBS

* Add 100_sphere.obj to root sim directory

* Pass pointer to Fluid object to openGLMain

* Make getParticlePositions() member function public

* Update particle positions with data from Trick sim

* Remove usage of Trickless SPH sim

* Include .o files in .gitignore

* Remove glm folder from repo

* Ignore local glm folder in .gitignore

* Copy Marching Cubes lookup tables from Cory Bloyd's implementation

* Create simple shaders to render the fluid mesh

* Create GridCell struct to store the isoValues and vertices of each voxel for Marching Cubes

* Add Trickless SPH sim back temporarily to test Marching Cubes functionality

* Implement generateCellMesh similar to Polygonise from Paul Borke
's 'Polygonising a Scalar Field'

* Implement vertex isovalue interpolation

* Create marching_cubes.h

* Modify OpenGL program to render fluid mesh rather than particles

* Create prototypes for updateIsoValues and initializeGridCells

* Implement function to initialize vector of GridCell structs

* Implement function to update isoValues of each vertex in each GridCell. (currently setting each isoValue to distance from origin)

* Support triangle indexing for more than one GridCell and do quick fix for interpolation bug

* Initialize and update both SPH particles and isoValues in grid for Marching Cubes

* Control OpenGL graphics refresh rate via usleep

* Set isoValues for vertices equal to distance from origin in order to demonstrate expanding sphere

* Fix bug with mesh vertex and face data failing to update

* Reduce number of particles to speed up sim testing

* Fix bug where only a quarter of vertex positions update

* Make particle depth greater than one to test marching cubes behavior for 3D fluid

* Redefine isoValues of gridCells to have a value equal to the number of particles within a given radius of the grid cell vertex

* Successfully update fluid mesh in response to particle positions

* Add spatial grid member and prototypes for neighbor list methods in Fluid header

* Implement buildSpatialGrid which assigns each particle to a grid cell

* Implement getCandidateNeighbors, which returns a vector of all particles within a given grid cell and up to its 26 neighboring cells

* Implement Marching Cubes isoValue update on the GPU

* Use GPU implementation of isoValue update

* Include files with build commands

* Write code for unoptimized GPU neighborlist in comments

* Add constants for spatial grid in Fluid.hh

* Use neighbor list in computeForces and computeDensityAndPressure procedures

* Keep count of the number of timesteps simulated for use in the marching cubes update

* Remove trickless SPH from unit testing

* Use Trick SPH particle data to update mesh from marching cubes

* Modify S_overrides.mk to compile and link files from Marching Cubes implementation

* Remove extraneous include

* Remove debugging code from opengl.cc

* Add back OpenGL program to render particles

* Implement ability to toggle between fluid mesh and particle modes

* Reverse surface normal direction for fluid mesh

* Add directory variables to S_overrides.mk

* Remove comment from S_overrides.mk

* Remove unused files

* Cleanup commented code

* Cleanup includes and remove unused variable

* Move marching cubes parameters to Fluid class

* Create variables to toggle GPU mode and neighbor list

* Fix malloc to include null terminated character of string and call free in loadObj

* Remove graphics window title

* Add GPU sim support for a number of particles that is not a multiple of NUM_THREADS

* Create initial conditions to demo

* Rename .cc files to .cpp

* Perform major structural refactoring on opengl.cpp

* Create refreshRate and mcUpdateFreq members in Fluid class

* Credit tables to Cory Bloyd

* Add Paul Bourke attribution

* Replace extern usage with header files

* Rename input files

* Remove unused Particle.cpp

* Move sphere models and move model file path to Fluid class

* Replace marching cubes triangle construction while loop with for loop

Co-authored-by: Scott Fennell <spfennell@gmail.com>

* Condense MC logic using loops

* Remove deprecated FluidServer from master

* Fix memory leak in iso_values.cu

* Change environment specific path to relative path in launch_vsc.py

* Write comment and TODO for CPU isoValue update code

* Add comments for demos and reduce usage of magic numbers in input.py

* Add enum for demos and fix concentric circle demo bug

* Fix style in input.py

* Replace scheduled job with flag for SPH initialization with initialization job

* Add shutdown job to free memory used by CUDA

* Add newlines to the end of files

* Delete deprecated FluidServer.cpp

Co-authored-by: Scott Fennell <spfennell@gmail.com>
This commit is contained in:
nateh7 2021-08-13 11:15:04 -05:00 committed by GitHub
parent f6ce855683
commit e537f2bae4
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
37 changed files with 35235 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

300
trick_sims/SIM_fluid/100_sphere.obj Normal file
View File

@ -0,0 +1,300 @@
# Blender v2.82 (sub 7) OBJ File: ''
# www.blender.org
o Sphere_Sphere.002
v -0.001779 1.010472 -0.485388
v -0.001779 0.086593 -1.102704
v -0.001779 -0.620514 -0.809811
v -0.001779 -0.837287 -0.485388
v 0.164261 1.010472 -0.447490
v 0.305023 0.793700 -0.739785
v 0.399077 0.469276 -0.935091
v 0.432105 0.086593 -1.003673
v 0.399077 -0.296091 -0.935091
v 0.305023 -0.620514 -0.739785
v 0.164261 -0.837287 -0.447490
v 0.297415 1.010472 -0.341303
v 0.551059 0.793700 -0.543578
v 0.720539 0.469276 -0.678733
v 0.780052 0.086593 -0.726194
v 0.720539 -0.296091 -0.678733
v 0.551059 -0.620514 -0.543578
v 0.297415 -0.837287 -0.341303
v 0.371310 1.010472 -0.187859
v 0.687599 0.793700 -0.260050
v 0.898937 0.469276 -0.308287
v 0.973149 0.086593 -0.325225
v 0.898937 -0.296091 -0.308287
v 0.687599 -0.620514 -0.260050
v 0.371310 -0.837287 -0.187859
v 0.371310 1.010472 -0.017549
v 0.687599 0.793700 0.054642
v 0.898937 0.469276 0.102878
v 0.973149 0.086593 0.119817
v 0.898937 -0.296091 0.102878
v 0.687599 -0.620514 0.054642
v 0.371310 -0.837287 -0.017549
v 0.297415 1.010472 0.135895
v 0.551059 0.793700 0.338170
v 0.720539 0.469276 0.473325
v 0.780052 0.086593 0.520786
v 0.720539 -0.296091 0.473325
v 0.551059 -0.620514 0.338170
v 0.297415 -0.837287 0.135895
v 0.164261 1.010472 0.242082
v 0.305023 0.793700 0.534377
v 0.399077 0.469276 0.729682
v 0.432104 0.086593 0.798265
v 0.399077 -0.296091 0.729682
v 0.305023 -0.620514 0.534377
v 0.164261 -0.837287 0.242082
v -0.001779 1.010472 0.279979
v -0.001779 0.793700 0.604402
v -0.001779 0.469276 0.821175
v -0.001779 0.086593 0.897296
v -0.001779 -0.296091 0.821175
v -0.001779 -0.620514 0.604402
v -0.001779 -0.837287 0.279979
v -0.167819 1.010472 0.242081
v -0.308581 0.793700 0.534377
v -0.402635 0.469276 0.729682
v -0.435663 0.086593 0.798264
v -0.402635 -0.296091 0.729682
v -0.308581 -0.620514 0.534377
v -0.167819 -0.837287 0.242082
v -0.300973 1.010472 0.135895
v -0.554617 0.793700 0.338169
v -0.724097 0.469276 0.473325
v -0.783611 0.086593 0.520785
v -0.724097 -0.296091 0.473325
v -0.554617 -0.620514 0.338169
v -0.300973 -0.837287 0.135895
v -0.001779 -0.913407 -0.102704
v -0.374868 1.010472 -0.017549
v -0.691157 0.793700 0.054642
v -0.902495 0.469276 0.102878
v -0.976707 0.086593 0.119816
v -0.902495 -0.296091 0.102878
v -0.691157 -0.620514 0.054642
v -0.374868 -0.837287 -0.017549
v -0.374868 1.010472 -0.187859
v -0.691157 0.793700 -0.260051
v -0.902495 0.469276 -0.308287
v -0.976707 0.086593 -0.325225
v -0.902495 -0.296091 -0.308287
v -0.691157 -0.620514 -0.260051
v -0.374868 -0.837287 -0.187859
v -0.300973 1.010472 -0.341303
v -0.554617 0.793700 -0.543578
v -0.724097 0.469276 -0.678734
v -0.783610 0.086593 -0.726194
v -0.724097 -0.296091 -0.678734
v -0.554617 -0.620514 -0.543578
v -0.300973 -0.837287 -0.341303
v -0.167819 1.010472 -0.447490
v -0.308581 0.793700 -0.739785
v -0.402635 0.469276 -0.935091
v -0.435662 0.086593 -1.003673
v -0.402635 -0.296091 -0.935091
v -0.308581 -0.620514 -0.739785
v -0.167819 -0.837287 -0.447490
v -0.001779 1.086593 -0.102704
v -0.001779 0.793700 -0.809811
v -0.001779 0.469276 -1.026583
v -0.001779 -0.296091 -1.026583
s off
f 98 5 6
f 3 11 4
f 100 8 9
f 98 7 99
f 1 97 5
f 68 4 11
f 100 10 3
f 99 8 2
f 9 15 16
f 7 13 14
f 5 97 12
f 68 11 18
f 9 17 10
f 7 15 8
f 6 12 13
f 10 18 11
f 68 18 25
f 17 23 24
f 14 22 15
f 12 20 13
f 17 25 18
f 16 22 23
f 13 21 14
f 12 97 19
f 21 29 22
f 19 27 20
f 25 31 32
f 23 29 30
f 21 27 28
f 19 97 26
f 68 25 32
f 23 31 24
f 31 39 32
f 30 36 37
f 28 34 35
f 26 97 33
f 68 32 39
f 30 38 31
f 28 36 29
f 26 34 27
f 34 42 35
f 33 97 40
f 68 39 46
f 37 45 38
f 36 42 43
f 33 41 34
f 38 46 39
f 36 44 37
f 68 46 53
f 44 52 45
f 43 49 50
f 41 47 48
f 45 53 46
f 43 51 44
f 41 49 42
f 40 97 47
f 47 55 48
f 52 60 53
f 50 58 51
f 48 56 49
f 47 97 54
f 68 53 60
f 51 59 52
f 49 57 50
f 59 67 60
f 57 65 58
f 56 62 63
f 54 97 61
f 68 60 67
f 58 66 59
f 56 64 57
f 55 61 62
f 61 97 69
f 68 67 75
f 66 73 74
f 63 72 64
f 62 69 70
f 66 75 67
f 64 73 65
f 62 71 63
f 73 81 74
f 71 79 72
f 69 77 70
f 74 82 75
f 72 80 73
f 70 78 71
f 69 97 76
f 68 75 82
f 81 89 82
f 79 87 80
f 77 85 78
f 76 97 83
f 68 82 89
f 81 87 88
f 78 86 79
f 76 84 77
f 86 94 87
f 84 92 85
f 83 97 90
f 68 89 96
f 87 95 88
f 85 93 86
f 84 90 91
f 88 96 89
f 68 96 4
f 94 3 95
f 92 2 93
f 91 1 98
f 95 4 96
f 94 2 100
f 91 99 92
f 90 97 1
f 98 1 5
f 3 10 11
f 100 2 8
f 98 6 7
f 100 9 10
f 99 7 8
f 9 8 15
f 7 6 13
f 9 16 17
f 7 14 15
f 6 5 12
f 10 17 18
f 17 16 23
f 14 21 22
f 12 19 20
f 17 24 25
f 16 15 22
f 13 20 21
f 21 28 29
f 19 26 27
f 25 24 31
f 23 22 29
f 21 20 27
f 23 30 31
f 31 38 39
f 30 29 36
f 28 27 34
f 30 37 38
f 28 35 36
f 26 33 34
f 34 41 42
f 37 44 45
f 36 35 42
f 33 40 41
f 38 45 46
f 36 43 44
f 44 51 52
f 43 42 49
f 41 40 47
f 45 52 53
f 43 50 51
f 41 48 49
f 47 54 55
f 52 59 60
f 50 57 58
f 48 55 56
f 51 58 59
f 49 56 57
f 59 66 67
f 57 64 65
f 56 55 62
f 58 65 66
f 56 63 64
f 55 54 61
f 66 65 73
f 63 71 72
f 62 61 69
f 66 74 75
f 64 72 73
f 62 70 71
f 73 80 81
f 71 78 79
f 69 76 77
f 74 81 82
f 72 79 80
f 70 77 78
f 81 88 89
f 79 86 87
f 77 84 85
f 81 80 87
f 78 85 86
f 76 83 84
f 86 93 94
f 84 91 92
f 87 94 95
f 85 92 93
f 84 83 90
f 88 95 96
f 94 100 3
f 92 99 2
f 91 90 1
f 95 3 4
f 94 93 2
f 91 98 99

29702
trick_sims/SIM_fluid/10k_sphere.obj Normal file
View File

File diff suppressed because it is too large Load Diff

2778
trick_sims/SIM_fluid/1k_sphere.obj Normal file
View File

File diff suppressed because it is too large Load Diff

88
trick_sims/SIM_fluid/Fluid_VSC.py Executable file
View File

@ -0,0 +1,88 @@
#!/usr/bin/python
import sys
import socket
from tkinter import *
# Process command line args
if ( len(sys.argv) == 2) :
trick_varserver_port = int(sys.argv[1])
else :
print("Usage: vsclient <port_number>")
sys.exit()
# Connect to the variable server
client_socket = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
client_socket.connect( ("localhost", trick_varserver_port) )
insock = client_socket.makefile("r")
# Request size of BOUND in order to determine window size
# Also request the number of particles for rendering
client_socket.send("trick.var_pause()\n".encode())
client_socket.send("trick.var_ascii()\n".encode())
client_socket.send("trick.var_add(\"dyn.fluid.BOUND\")\n".encode())
client_socket.send("trick.var_add(\"dyn.fluid.NUM_PARTICLES\")\n".encode())
client_socket.send("trick.var_unpause()\n".encode())
line = insock.readline()
field = line.split("\t")
#print(field)
BOUND = int(field[1])
NUM_PARTICLES = int(field[2])
HEIGHT = 2 * BOUND
WIDTH = HEIGHT
"""
# Unit test 0: request viscosity of the fluid
client_socket.send("trick.var_pause()\n".encode())
client_socket.send("trick.var_ascii()\n".encode())
client_socket.send("trick.var_add(\"dyn.fluid.VISC\")\n".encode())
client_socket.send("trick.var_unpause()\n".encode())
# Unit test 1: request the x position of a particle
client_socket.send("trick.var_pause()\n".encode())
client_socket.send("trick.var_ascii()\n".encode())
client_socket.send("trick.var_add(\"dyn.fluid.particlesArr[0].pos[0]\")\n".encode())
client_socket.send("trick.var_unpause()\n".encode())
"""
# Create the GUI
tk = Tk()
canvas = Canvas(tk, width=WIDTH, height=HEIGHT)
tk.title("SPH Display")
canvas.pack()
# Create oval objects to represent fluid particles
particle_radius = 2;
fluidParticles = []
for i in range(NUM_PARTICLES):
fluidParticles.append(canvas.create_oval(0, 0, particle_radius, particle_radius, fill="blue"))
client_socket.send("trick.var_pause()\n".encode())
client_socket.send("trick.var_ascii()\n".encode())
for i in range(NUM_PARTICLES):
client_socket.send("trick.var_add(\"dyn.fluid.particlesArr[{}].pos[0]\")\n".format(i).encode())
client_socket.send("trick.var_add(\"dyn.fluid.particlesArr[{}].pos[1]\")\n".format(i).encode())
client_socket.send("trick.var_unpause()\n".encode())
# Repeatedly read and process the responses from the variable server.
while (True):
line = insock.readline()
if line == '':
break
field = line.split("\t")
# Get position of particle and update it on the canvas
for i in range(NUM_PARTICLES):
# TODO: calculate correct offset to index into field (one particle is not rendering correctly)
y = float(field[2 * i])
x = float(field[2 * i + 1])
cx = x + BOUND
cy = HEIGHT - (y + BOUND)
canvas.coords(fluidParticles[i], cx - particle_radius, cy - particle_radius, cx + particle_radius, cy+ particle_radius)
#print(line)
# Update the Tk graphics
tk.update()
# Keep the window open, when the data stops.
tk.mainloop()

7
trick_sims/SIM_fluid/Modified_data/realtime.py Normal file
View File

@ -0,0 +1,7 @@
#trick.frame_log_on()
trick.real_time_enable()
trick.exec_set_software_frame(0.01)
#trick.itimer_enable()
trick.exec_set_enable_freeze(True)
trick.exec_set_freeze_command(True)
trick.sim_control_panel_set_enabled(True)

51
trick_sims/SIM_fluid/RUN_test/input.py Normal file
View File

@ -0,0 +1,51 @@
from enum import Enum
import math
exec(open("Modified_data/realtime.py").read())
dyn.fluid.BOUND = 300
dyn.fluid.ISO_RADIUS = 32
class Mode(Enum):
DEFAULT = -1
PARABOLOID = 0
COSINE = 1
CIRCLES = 2
mode = Mode.CIRCLES
if mode == Mode.PARABOLOID:
# paraboloid initial condition
for i in range(dyn.fluid.NUM_PARTICLES):
ROOT_NUM_PARTICLES = math.sqrt(dyn.fluid.NUM_PARTICLES)
tx = i / ROOT_NUM_PARTICLES - (ROOT_NUM_PARTICLES / 2)
tz = i % ROOT_NUM_PARTICLES - (ROOT_NUM_PARTICLES / 2)
dyn.fluid.particlesArr[i].pos[0] = tx * 4
dyn.fluid.particlesArr[i].pos[2] = tz * 4
dyn.fluid.particlesArr[i].pos[1] = tx * tx + tz * tz
elif mode == Mode.COSINE:
# cosine wave initial condition
for i in range(dyn.fluid.NUM_PARTICLES):
t = i / 2
dyn.fluid.particlesArr[i].pos[1] = 15 * math.cos(math.radians(t))
elif mode == Mode.CIRCLES:
# concentric circles initial conditions
dyn.fluid.NUM_PARTICLES = 800
n = int(dyn.fluid.NUM_PARTICLES)
for i in range(n):
theta = (360 / (n / 2)) * i
print(theta)
dyn.fluid.particlesArr[i].pos[0] = 100 * math.cos(math.radians(theta))
dyn.fluid.particlesArr[i].pos[1] = 100 * math.sin(math.radians(theta))
if i > n / 2 and i <= (7 / 8) * n:
theta = (360 / (3 / 8 * n)) * i
dyn.fluid.particlesArr[i].pos[0] = 70 * math.cos(math.radians(theta))
dyn.fluid.particlesArr[i].pos[1] = 70 * math.sin(math.radians(theta))
if i > (7 / 8) * n:
theta = (360 / (1 / 8 * n)) * i
dyn.fluid.particlesArr[i].pos[0] = 40 * math.cos(math.radians(theta))
dyn.fluid.particlesArr[i].pos[1] = 40 * math.sin(math.radians(theta))

17
trick_sims/SIM_fluid/RUN_test/launch_vsc.py Normal file
View File

@ -0,0 +1,17 @@
exec(open("Modified_data/realtime.py").read())
#==================================
# Start the variable server client.
#==================================
var_server_port = trick.var_server_get_port()
fluid_graphics_path = "Fluid_VSC.py"
if (os.path.isfile(fluid_graphics_path)) :
fluid_graphics_cmd = "python3 " + fluid_graphics_path + " " + str(var_server_port) + " &"
print(fluid_graphics_cmd)
os.system(fluid_graphics_cmd)
else :
print('Oops! Can\'t find ' + fluid_graphics_path)
trick.var_set_copy_mode(0)
trick.var_cycle(0.02)

25
trick_sims/SIM_fluid/S_define Normal file
View File

@ -0,0 +1,25 @@
/************************TRICK HEADER*************************
PURPOSE:
( Simulate a Lagrangian fluid using smoothed particle hydrodynamics. )
LIBRARY DEPENDENCIES:
((fluid/src/Fluid.cpp)
(fluid/src/Particle.cpp))
*************************************************************/
#include "sim_objects/default_trick_sys.sm"
##include "fluid/include/Fluid.hh"
class FluidSimObject : public Trick::SimObject {
public:
Fluid fluid;
FluidSimObject() {
("default_data") fluid.default_data() ;
("initialization") fluid.init_cuda() ;
(0.02, "scheduled") fluid.update_SPH() ;
("shutdown") fluid.shutdown_cuda() ;
}
};
FluidSimObject dyn;

42
trick_sims/SIM_fluid/S_overrides.mk Normal file
View File

@ -0,0 +1,42 @@
SIM_DIR = models/fluid/src
GRAPHICS_DIR = models/fluid/graphics/src
TRICK_CFLAGS += -Imodels -I/usr/local/cuda/include
TRICK_CXXFLAGS += -Imodels -I/usr/local/cuda/include
TRICK_USER_LINK_LIBS += $(SIM_DIR)/gpuCode.o $(SIM_DIR)/sph_gpu.o $(GRAPHICS_DIR)/iso_values.o $(GRAPHICS_DIR)/opengl_caller.o $(GRAPHICS_DIR)/opengl.o $(GRAPHICS_DIR)/gui.o $(GRAPHICS_DIR)/grid_cell.o $(GRAPHICS_DIR)/marching_cubes.o -L/usr/local/cuda/lib64 -lcudart -lcudadevrt -lglfw -lGL -lGLEW -pthread
$(S_MAIN): $(GRAPHICS_DIR)/opengl_caller.o $(GRAPHICS_DIR)/opengl.o $(GRAPHICS_DIR)/gui.o $(GRAPHICS_DIR)/grid_cell.o $(GRAPHICS_DIR)/marching_cubes.o $(GRAPHICS_DIR)/iso_values.o
$(GRAPHICS_DIR)/grid_cell.o: $(GRAPHICS_DIR)/grid_cell.cpp
g++ -c $(GRAPHICS_DIR)/grid_cell.cpp -o $(GRAPHICS_DIR)/grid_cell.o
$(GRAPHICS_DIR)/marching_cubes.o: $(GRAPHICS_DIR)/marching_cubes.cpp
g++ -c $(GRAPHICS_DIR)/marching_cubes.cpp -o $(GRAPHICS_DIR)/marching_cubes.o
$(GRAPHICS_DIR)/opengl_caller.o: $(GRAPHICS_DIR)/opengl_caller.cpp
g++ -c $(GRAPHICS_DIR)/opengl_caller.cpp -o $(GRAPHICS_DIR)/opengl_caller.o
$(GRAPHICS_DIR)/opengl.o: $(GRAPHICS_DIR)/opengl.cpp
g++ -c $(GRAPHICS_DIR)/opengl.cpp -o $(GRAPHICS_DIR)/opengl.o
$(GRAPHICS_DIR)/gui.o: $(GRAPHICS_DIR)/gui.cpp
g++ -c $(GRAPHICS_DIR)/gui.cpp -o $(GRAPHICS_DIR)/gui.o
$(S_MAIN): $(SIM_DIR)/gpuCode.o
$(SIM_DIR)/sph_gpu.o: $(SIM_DIR)/sph_gpu.cu
nvcc -arch=sm_35 -I. -dc $(SIM_DIR)/sph_gpu.cu -o $(SIM_DIR)/sph_gpu.o
$(GRAPHICS_DIR)/iso_values.o: $(GRAPHICS_DIR)/iso_values.cu
nvcc -arch=sm_35 -I. -dc $(GRAPHICS_DIR)/iso_values.cu -o $(GRAPHICS_DIR)/iso_values.o
$(SIM_DIR)/gpuCode.o: $(SIM_DIR)/sph_gpu.o
nvcc -arch=sm_35 -dlink $(SIM_DIR)/sph_gpu.o $(GRAPHICS_DIR)/iso_values.o $(GRAPHICS_DIR)/opengl_caller.o $(GRAPHICS_DIR)/opengl.o $(GRAPHICS_DIR)/gui.o $(GRAPHICS_DIR)/grid_cell.o $(GRAPHICS_DIR)/marching_cubes.o $(S_OBJECTS) -o $(SIM_DIR)/gpuCode.o
clean: my_clean
my_clean:
-rm -rf $(SIM_DIR)/sph_gpu.o $(SIM_DIR)/gpuCode.o $(GRAPHICS_DIR)/grid_cell.o $(GRAPHICS_DIR)/gui.o $(GRAPHICS_DIR)/iso_values.o $(GRAPHICS_DIR)/marching_cubes.o $(GRAPHICS_DIR)/opengl_caller.o $(GRAPHICS_DIR)/opengl.o

3
trick_sims/SIM_fluid/models/fluid/graphics/src/.gitignore vendored Normal file
View File

@ -0,0 +1,3 @@
sim
glm
*.o

1
trick_sims/SIM_fluid/models/fluid/graphics/src/cpu_build Executable file
View File

@ -0,0 +1 @@
g++ gui.cc sph.cpp opengl.cc marching_cubes.cpp grid_cell.cpp -lglfw -lGL -lGLEW -o sim -pthread

1
trick_sims/SIM_fluid/models/fluid/graphics/src/gpu_build Executable file
View File

@ -0,0 +1 @@
g++ -c gui.cc sph.cpp opengl.cc marching_cubes.cpp grid_cell.cpp; nvcc -c iso_values.cu; nvcc gui.o sph.o opengl.o marching_cubes.o grid_cell.o iso_values.o -lglfw -lGL -lGLEW -o mc_gpu

67
trick_sims/SIM_fluid/models/fluid/graphics/src/grid_cell.cpp Normal file
View File

@ -0,0 +1,67 @@
#include "grid_cell.h"
#include "iso_values.h"
void updateIsoValues(std::vector<GridCell>& gridCells, std::vector<float> particlePositions, float radius) {
updateIsoValuesGPUCaller(gridCells, particlePositions, radius);
// TODO: give user the ability to toggle between GPU and CPU marching cubes
/* updateIsoValues implementation on CPU
int numParticles = particlePositions.size() / 3;
for (int i = 0; i < gridCells.size(); i++) {
GridCell& current = gridCells[i];
for (int j = 0; j < 8; j++) {
current.isoValues[j] = 0;
glm::vec3 vertexPos(current.vertices[j]);
for (int k = 0; k < numParticles; k++) {
float x = particlePositions[3*k];
float y = particlePositions[3*k+1];
float z = particlePositions[3*k+2];
glm::vec3 particlePos(x, y, z);
if (glm::length(vertexPos - particlePos) < radius) {
current.isoValues[j]+=1;
}
}
}
}*/
}
void initializeGridCells(std::vector<GridCell>& gridCells, double bound, int amountPerDim) {
double min = -bound;
double max = bound;
double size = 2 * bound / amountPerDim;
for (double x = min; x < max; x += size) {
for (double y = min; y < max; y += size) {
for (double z = min; z < max; z += size) {
GridCell current;
// vertex ordering convention accoriding to Paul Bourke
current.vertices[0] = glm::vec4(x, y, z, 1);
current.vertices[1] = glm::vec4(x + size, y, z, 1);
current.vertices[2] = glm::vec4(x + size, y, z + size, 1);
current.vertices[3] = glm::vec4(x, y, z + size, 1);
current.vertices[4] = glm::vec4(x, y + size, z, 1);
current.vertices[5] = glm::vec4(x + size, y + size, z, 1);
current.vertices[6] = glm::vec4(x + size, y + size, z + size, 1);
current.vertices[7] = glm::vec4(x, y + size, z + size, 1);
current.isoValues[0] = 0;
current.isoValues[1] = 0;
current.isoValues[2] = 0;
current.isoValues[3] = 0;
current.isoValues[4] = 0;
current.isoValues[5] = 0;
current.isoValues[6] = 0;
current.isoValues[7] = 0;
gridCells.push_back(current);
}
}
}
}

12
trick_sims/SIM_fluid/models/fluid/graphics/src/grid_cell.h Normal file
View File

@ -0,0 +1,12 @@
#pragma once
#include "glm/glm.hpp"
#include <vector>
struct GridCell {
glm::vec4 vertices[8];
float isoValues[8] = {0, 0, 0, 0, 0, 0, 0, 0};
};
void updateIsoValues(std::vector<GridCell>& gridCells, std::vector<float> particlePostions, float radius);
void initializeGridCells(std::vector<GridCell>& gridCells, double bound, int amountPerDim);

191
trick_sims/SIM_fluid/models/fluid/graphics/src/gui.cpp Normal file
View File

@ -0,0 +1,191 @@
#include "gui.h"
#include "glm/gtc/matrix_access.hpp"
#include "glm/gtx/transform.hpp"
const float kNear = 0.1f;
const float kFar = 1000.0f;
const float kFov = 45.0f;
extern bool paused;
extern bool meshMode;
GUI::GUI(GLFWwindow* window, int view_width, int view_height, int preview_height)
:window_(window), preview_height_(preview_height)
{
glfwSetWindowUserPointer(window_, this);
glfwSetKeyCallback(window_, KeyCallback);
glfwSetCursorPosCallback(window_, MousePosCallback);
glfwSetMouseButtonCallback(window_, MouseButtonCallback);
glfwGetWindowSize(window_, &window_width_, &window_height_);
if (view_width < 0 || view_height < 0) {
view_width_ = window_width_;
view_height_ = window_height_;
} else {
view_width_ = view_width;
view_height_ = view_height;
}
float aspect_ = static_cast<float>(view_width_) / view_height_;
projection_matrix_ = glm::perspective((float)(kFov * (M_PI / 180.0f)), aspect_, kNear, kFar);
}
GUI::~GUI()
{
}
void GUI::keyCallback(int key, int scancode, int action, int mods)
{
#if 0
if (action != 2)
std::cerr << "Key: " << key << " action: " << action << std::endl;
#endif
if (key == GLFW_KEY_ESCAPE && action == GLFW_PRESS) {
glfwSetWindowShouldClose(window_, GL_TRUE);
return ;
}
if (mods == 0 && captureWASDUPDOWN(key, action))
return ;
if (key == GLFW_KEY_C && action != GLFW_RELEASE) {
fps_mode_ = !fps_mode_;
} else if (key == GLFW_KEY_P && action != GLFW_RELEASE) {
paused = !paused;
} else if (key == GLFW_KEY_M && action != GLFW_RELEASE) {
meshMode = !meshMode;
}
}
void GUI::mousePosCallback(double mouse_x, double mouse_y)
{
last_x_ = current_x_;
last_y_ = current_y_;
current_x_ = mouse_x;
current_y_ = window_height_ - mouse_y;
float delta_x = current_x_ - last_x_;
float delta_y = current_y_ - last_y_;
if (sqrt(delta_x * delta_x + delta_y * delta_y) < 1e-15)
return;
if (mouse_x > view_width_)
return;
glm::vec3 mouse_direction = glm::normalize(glm::vec3(delta_x, delta_y, 0.0f));
bool drag_camera = drag_state_ && current_button_ == GLFW_MOUSE_BUTTON_RIGHT;
if (drag_camera) {
glm::vec3 axis = glm::normalize(
orientation_ *
glm::vec3(mouse_direction.y, -mouse_direction.x, 0.0f)
);
orientation_ =
glm::mat3(glm::rotate(rotation_speed_, axis) * glm::mat4(orientation_));
tangent_ = glm::column(orientation_, 0);
up_ = glm::column(orientation_, 1);
look_ = glm::column(orientation_, 2);
}
}
void GUI::mouseButtonCallback(int button, int action, int mods)
{
if (current_x_ <= view_width_) {
drag_state_ = (action == GLFW_PRESS);
current_button_ = button;
return ;
}
}
void GUI::updateMatrices()
{
if (fps_mode_)
center_ = eye_ + camera_distance_ * look_;
else
eye_ = center_ - camera_distance_ * look_;
view_matrix_ = glm::lookAt(eye_, center_, up_);
light_position_ = glm::vec4(eye_, 1.0f);
aspect_ = static_cast<float>(view_width_) / view_height_;
projection_matrix_ =
glm::perspective((float)(kFov * (M_PI / 180.0f)), aspect_, kNear, kFar);
model_matrix_ = glm::mat4(1.0f);
}
MatrixPointers GUI::getMatrixPointers() const
{
MatrixPointers ret;
ret.projection = &projection_matrix_;
ret.model= &model_matrix_;
ret.view = &view_matrix_;
return ret;
}
bool GUI::captureWASDUPDOWN(int key, int action)
{
if (key == GLFW_KEY_W) {
if (fps_mode_)
eye_ += zoom_speed_ * look_;
else
camera_distance_ -= zoom_speed_;
return true;
} else if (key == GLFW_KEY_S) {
if (fps_mode_)
eye_ -= zoom_speed_ * look_;
else
camera_distance_ += zoom_speed_;
return true;
} else if (key == GLFW_KEY_A) {
if (fps_mode_)
eye_ -= pan_speed_ * tangent_;
else
center_ -= pan_speed_ * tangent_;
return true;
} else if (key == GLFW_KEY_D) {
if (fps_mode_)
eye_ += pan_speed_ * tangent_;
else
center_ += pan_speed_ * tangent_;
return true;
} else if (key == GLFW_KEY_DOWN) {
if (fps_mode_)
eye_ -= pan_speed_ * up_;
else
center_ -= pan_speed_ * up_;
return true;
} else if (key == GLFW_KEY_UP) {
if (fps_mode_)
eye_ += pan_speed_ * up_;
else
center_ += pan_speed_ * up_;
return true;
}
return false;
}
// Delegate to the actual GUI object.
void GUI::KeyCallback(GLFWwindow* window, int key, int scancode, int action, int mods)
{
GUI* gui = (GUI*)glfwGetWindowUserPointer(window);
gui->keyCallback(key, scancode, action, mods);
}
void GUI::MousePosCallback(GLFWwindow* window, double mouse_x, double mouse_y)
{
GUI* gui = (GUI*)glfwGetWindowUserPointer(window);
gui->mousePosCallback(mouse_x, mouse_y);
}
void GUI::MouseButtonCallback(GLFWwindow* window, int button, int action, int mods)
{
GUI* gui = (GUI*)glfwGetWindowUserPointer(window);
gui->mouseButtonCallback(button, action, mods);
}

60
trick_sims/SIM_fluid/models/fluid/graphics/src/gui.h Normal file
View File

@ -0,0 +1,60 @@
#include "glm/glm.hpp"
#include "glm/gtc/matrix_transform.hpp"
#include <GLFW/glfw3.h>
struct MatrixPointers {
const glm::mat4 *projection, *model, *view;
};
class GUI {
public:
GUI(GLFWwindow*, int view_width = -1, int view_height = -1, int preview_height = -1);
~GUI();
void keyCallback(int key, int scancode, int action, int mods);
void mousePosCallback(double mouse_x, double mouse_y);
void mouseButtonCallback(int button, int action, int mods);
void updateMatrices();
MatrixPointers getMatrixPointers() const;
static void KeyCallback(GLFWwindow* window, int key, int scancode, int action, int mods);
static void MousePosCallback(GLFWwindow* window, double mouse_x, double mouse_y);
static void MouseButtonCallback(GLFWwindow* window, int button, int action, int mods);
private:
GLFWwindow* window_;
int window_width_, window_height_;
int view_width_, view_height_;
int preview_height_;
bool drag_state_ = false;
bool fps_mode_ = true;
int current_button_ = -1;
float last_x_ = 0.0f, last_y_ = 0.0f, current_x_ = 0.0f, current_y_ = 0.0f;
float camera_distance_ = 300.0;
float camera_depth_ = 0.f;
float pan_speed_ = 5.f;
float rotation_speed_ = 0.02f;
float zoom_speed_ = 5.f;
float aspect_;
glm::vec3 eye_ = glm::vec3(0.0f, camera_depth_, camera_distance_);
glm::vec3 up_ = glm::vec3(0.0f, 1.0f, 0.0f);
glm::vec3 look_ = glm::vec3(0.0f, 0.0f, -1.0f);
glm::vec3 tangent_ = glm::cross(look_, up_);
glm::vec3 center_ = eye_ - camera_distance_ * look_;
glm::mat3 orientation_ = glm::mat3(tangent_, up_, look_);
glm::vec4 light_position_;
glm::mat4 view_matrix_ = glm::lookAt(eye_, center_, up_);
glm::mat4 projection_matrix_;
glm::mat4 model_matrix_ = glm::mat4(1.0f);
bool captureWASDUPDOWN(int key, int action);
};

65
trick_sims/SIM_fluid/models/fluid/graphics/src/iso_values.cu Normal file
View File

@ -0,0 +1,65 @@
#include "cuda_runtime.h"
#include "grid_cell.h"
#include <stdio.h>
#define NUM_THREADS 1024
__global__ void updateIsoValuesGPU(GridCell* gridCells, float* particlePositions, int numParticles, int numCells, float radius) {
int tid = threadIdx.x;
int blockSize = numCells / NUM_THREADS;
int g_start = tid * blockSize;
int g_end = (tid + 1) * blockSize;
for (int i = g_start; i < g_end; i++) {
for (int j = 0; j < 8; j++) {
gridCells[i].isoValues[j] = 0;
glm::vec3 vertexPos(gridCells[i].vertices[j]);
//printf("(%f, %f, %f): %d\n", vertexPos.x, vertexPos.y, vertexPos.z, tid);
for (int k = 0; k < numParticles; k++) {
float x = particlePositions[3*k];
float y = particlePositions[3*k+1];
float z = particlePositions[3*k+2];
glm::vec3 particlePos(x, y, z);
if (glm::length(vertexPos - particlePos) < radius) {
gridCells[i].isoValues[j]+=1;
}
}
}
}
}
void updateIsoValuesGPUCaller(std::vector<GridCell>& gridCells, std::vector<float> particlePositions, float radius) {
GridCell* d_gridCells;
float* d_particlePositions;
cudaMalloc(&d_gridCells, gridCells.size() * sizeof(GridCell));
cudaMalloc(&d_particlePositions, particlePositions.size() * sizeof(float));
cudaMemcpy(d_gridCells, gridCells.data(), gridCells.size() * sizeof(GridCell), cudaMemcpyHostToDevice);
cudaMemcpy(d_particlePositions, particlePositions.data(), particlePositions.size() * sizeof(float), cudaMemcpyHostToDevice);
updateIsoValuesGPU<<<1, NUM_THREADS>>>(d_gridCells, d_particlePositions, particlePositions.size() / 3, gridCells.size(), radius);
cudaDeviceSynchronize();
cudaMemcpy(gridCells.data(), d_gridCells, gridCells.size() * sizeof(GridCell), cudaMemcpyDeviceToHost);
cudaFree(d_gridCells);
cudaFree(d_particlePositions);
}

6
trick_sims/SIM_fluid/models/fluid/graphics/src/iso_values.h Normal file
View File

@ -0,0 +1,6 @@
#pragma once
#include <vector>
#include <stdio.h>
void updateIsoValuesGPUCaller(std::vector<GridCell>& gridCells, std::vector<float> particlePositions, float radius);

52
trick_sims/SIM_fluid/models/fluid/graphics/src/marching_cubes.cpp Normal file
View File

@ -0,0 +1,52 @@
/* Adapted from Paul Bourke's "Polygonising a scalar field" */
#include "marching_cubes.h"
#include "mc_tables.h"
glm::vec4 vlerp(float isoLevel, glm::vec4 vectorA, glm::vec4 vectorB, float valueA, float valueB) {
// TODO: debug vlerp
float t = (isoLevel - std::min(valueA, valueB)) / std::abs(valueA - valueB);
t = 0.5;
return vectorA * t + vectorB * (1 - t);
}
void generateCellMesh(GridCell &grid, double isoLevel, std::vector<glm::uvec3> &triangles, std::vector<glm::vec4> &vertices)
{
int geometryCaseIndex = 0;
int startTriangleIdx = vertices.size();
for (int i = 0; i < 8; i++) {
if (grid.isoValues[i] > isoLevel) {
geometryCaseIndex |= (1 << i);
}
}
// no edges in grid cell contain a triangle vertex
if(edgeTable[geometryCaseIndex] == 0)
return;
static const int EDGE_CONNECTIONS[12][2] = {{0,1}, {1,2}, {2,3}, {3, 0}, {4,5}, {5,6}, {6,7}, {7, 4}, {0, 4}, {1, 5}, {2, 6}, {3, 7}};
glm::vec4 origin = glm::vec4(0, 0, 0, 1);
// On edges that contain a triangle vertex, interpolate that triangle vertex based on the isovalues of the endpoint vertices of that edge
for (int i = 0; i < 12; i++) {
if (edgeTable[geometryCaseIndex] & (1 << i)) {
int vertexOne = EDGE_CONNECTIONS[i][0];
int vertexTwo = EDGE_CONNECTIONS[i][1];
vertices.push_back(vlerp(isoLevel, grid.vertices[vertexOne], grid.vertices[vertexTwo], grid.isoValues[vertexOne], grid.isoValues[vertexTwo]));
} else {
vertices.push_back(origin);
}
}
// Add triangle indices to triangles
for(int i = 0; triTable[geometryCaseIndex][i] != -1; i += 3) {
int idxOne = triTable[geometryCaseIndex][i];
int idxTwo = triTable[geometryCaseIndex][i+1];
int idxThree = triTable[geometryCaseIndex][i+2];
triangles.push_back(glm::uvec3(idxOne + startTriangleIdx, idxTwo + startTriangleIdx, idxThree + startTriangleIdx));
}
}

6
trick_sims/SIM_fluid/models/fluid/graphics/src/marching_cubes.h Normal file
View File

@ -0,0 +1,6 @@
#pragma once
#include "grid_cell.h"
#include "glm/glm.hpp"
#include <vector>
void generateCellMesh(GridCell &grid, double isoLevel, std::vector<glm::uvec3> &triangles, std::vector<glm::vec4> &vertices);

293
trick_sims/SIM_fluid/models/fluid/graphics/src/mc_tables.h Normal file
View File

@ -0,0 +1,293 @@
/* Tables credit of Cory Bloyd (corysama@yahoo.com) */
/* http://paulbourke.net/geometry/polygonise/marchingsource.cpp */
int edgeTable[256]={
0x0 , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c,
0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac,
0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c,
0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc,
0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c,
0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc ,
0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc,
0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c,
0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc,
0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c,
0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460,
0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0,
0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c,
0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230,
0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,
0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190,
0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0 };
int triTable[256][16] =
{{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 8, 3, 9, 8, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 2, 10, 0, 2, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 8, 3, 2, 10, 8, 10, 9, 8, -1, -1, -1, -1, -1, -1, -1},
{3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 11, 2, 8, 11, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 9, 0, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 11, 2, 1, 9, 11, 9, 8, 11, -1, -1, -1, -1, -1, -1, -1},
{3, 10, 1, 11, 10, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 10, 1, 0, 8, 10, 8, 11, 10, -1, -1, -1, -1, -1, -1, -1},
{3, 9, 0, 3, 11, 9, 11, 10, 9, -1, -1, -1, -1, -1, -1, -1},
{9, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 3, 0, 7, 3, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 1, 9, 4, 7, 1, 7, 3, 1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 4, 7, 3, 0, 4, 1, 2, 10, -1, -1, -1, -1, -1, -1, -1},
{9, 2, 10, 9, 0, 2, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
{2, 10, 9, 2, 9, 7, 2, 7, 3, 7, 9, 4, -1, -1, -1, -1},
{8, 4, 7, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{11, 4, 7, 11, 2, 4, 2, 0, 4, -1, -1, -1, -1, -1, -1, -1},
{9, 0, 1, 8, 4, 7, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
{4, 7, 11, 9, 4, 11, 9, 11, 2, 9, 2, 1, -1, -1, -1, -1},
{3, 10, 1, 3, 11, 10, 7, 8, 4, -1, -1, -1, -1, -1, -1, -1},
{1, 11, 10, 1, 4, 11, 1, 0, 4, 7, 11, 4, -1, -1, -1, -1},
{4, 7, 8, 9, 0, 11, 9, 11, 10, 11, 0, 3, -1, -1, -1, -1},
{4, 7, 11, 4, 11, 9, 9, 11, 10, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 4, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 5, 4, 1, 5, 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{8, 5, 4, 8, 3, 5, 3, 1, 5, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 8, 1, 2, 10, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
{5, 2, 10, 5, 4, 2, 4, 0, 2, -1, -1, -1, -1, -1, -1, -1},
{2, 10, 5, 3, 2, 5, 3, 5, 4, 3, 4, 8, -1, -1, -1, -1},
{9, 5, 4, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 11, 2, 0, 8, 11, 4, 9, 5, -1, -1, -1, -1, -1, -1, -1},
{0, 5, 4, 0, 1, 5, 2, 3, 11, -1, -1, -1, -1, -1, -1, -1},
{2, 1, 5, 2, 5, 8, 2, 8, 11, 4, 8, 5, -1, -1, -1, -1},
{10, 3, 11, 10, 1, 3, 9, 5, 4, -1, -1, -1, -1, -1, -1, -1},
{4, 9, 5, 0, 8, 1, 8, 10, 1, 8, 11, 10, -1, -1, -1, -1},
{5, 4, 0, 5, 0, 11, 5, 11, 10, 11, 0, 3, -1, -1, -1, -1},
{5, 4, 8, 5, 8, 10, 10, 8, 11, -1, -1, -1, -1, -1, -1, -1},
{9, 7, 8, 5, 7, 9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 3, 0, 9, 5, 3, 5, 7, 3, -1, -1, -1, -1, -1, -1, -1},
{0, 7, 8, 0, 1, 7, 1, 5, 7, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 7, 8, 9, 5, 7, 10, 1, 2, -1, -1, -1, -1, -1, -1, -1},
{10, 1, 2, 9, 5, 0, 5, 3, 0, 5, 7, 3, -1, -1, -1, -1},
{8, 0, 2, 8, 2, 5, 8, 5, 7, 10, 5, 2, -1, -1, -1, -1},
{2, 10, 5, 2, 5, 3, 3, 5, 7, -1, -1, -1, -1, -1, -1, -1},
{7, 9, 5, 7, 8, 9, 3, 11, 2, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 7, 9, 7, 2, 9, 2, 0, 2, 7, 11, -1, -1, -1, -1},
{2, 3, 11, 0, 1, 8, 1, 7, 8, 1, 5, 7, -1, -1, -1, -1},
{11, 2, 1, 11, 1, 7, 7, 1, 5, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 8, 8, 5, 7, 10, 1, 3, 10, 3, 11, -1, -1, -1, -1},
{5, 7, 0, 5, 0, 9, 7, 11, 0, 1, 0, 10, 11, 10, 0, -1},
{11, 10, 0, 11, 0, 3, 10, 5, 0, 8, 0, 7, 5, 7, 0, -1},
{11, 10, 5, 7, 11, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 0, 1, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 8, 3, 1, 9, 8, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
{1, 6, 5, 2, 6, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 6, 5, 1, 2, 6, 3, 0, 8, -1, -1, -1, -1, -1, -1, -1},
{9, 6, 5, 9, 0, 6, 0, 2, 6, -1, -1, -1, -1, -1, -1, -1},
{5, 9, 8, 5, 8, 2, 5, 2, 6, 3, 2, 8, -1, -1, -1, -1},
{2, 3, 11, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{11, 0, 8, 11, 2, 0, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1, -1, -1, -1},
{5, 10, 6, 1, 9, 2, 9, 11, 2, 9, 8, 11, -1, -1, -1, -1},
{6, 3, 11, 6, 5, 3, 5, 1, 3, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 11, 0, 11, 5, 0, 5, 1, 5, 11, 6, -1, -1, -1, -1},
{3, 11, 6, 0, 3, 6, 0, 6, 5, 0, 5, 9, -1, -1, -1, -1},
{6, 5, 9, 6, 9, 11, 11, 9, 8, -1, -1, -1, -1, -1, -1, -1},
{5, 10, 6, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 3, 0, 4, 7, 3, 6, 5, 10, -1, -1, -1, -1, -1, -1, -1},
{1, 9, 0, 5, 10, 6, 8, 4, 7, -1, -1, -1, -1, -1, -1, -1},
{10, 6, 5, 1, 9, 7, 1, 7, 3, 7, 9, 4, -1, -1, -1, -1},
{6, 1, 2, 6, 5, 1, 4, 7, 8, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 5, 5, 2, 6, 3, 0, 4, 3, 4, 7, -1, -1, -1, -1},
{8, 4, 7, 9, 0, 5, 0, 6, 5, 0, 2, 6, -1, -1, -1, -1},
{7, 3, 9, 7, 9, 4, 3, 2, 9, 5, 9, 6, 2, 6, 9, -1},
{3, 11, 2, 7, 8, 4, 10, 6, 5, -1, -1, -1, -1, -1, -1, -1},
{5, 10, 6, 4, 7, 2, 4, 2, 0, 2, 7, 11, -1, -1, -1, -1},
{0, 1, 9, 4, 7, 8, 2, 3, 11, 5, 10, 6, -1, -1, -1, -1},
{9, 2, 1, 9, 11, 2, 9, 4, 11, 7, 11, 4, 5, 10, 6, -1},
{8, 4, 7, 3, 11, 5, 3, 5, 1, 5, 11, 6, -1, -1, -1, -1},
{5, 1, 11, 5, 11, 6, 1, 0, 11, 7, 11, 4, 0, 4, 11, -1},
{0, 5, 9, 0, 6, 5, 0, 3, 6, 11, 6, 3, 8, 4, 7, -1},
{6, 5, 9, 6, 9, 11, 4, 7, 9, 7, 11, 9, -1, -1, -1, -1},
{10, 4, 9, 6, 4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 10, 6, 4, 9, 10, 0, 8, 3, -1, -1, -1, -1, -1, -1, -1},
{10, 0, 1, 10, 6, 0, 6, 4, 0, -1, -1, -1, -1, -1, -1, -1},
{8, 3, 1, 8, 1, 6, 8, 6, 4, 6, 1, 10, -1, -1, -1, -1},
{1, 4, 9, 1, 2, 4, 2, 6, 4, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 8, 1, 2, 9, 2, 4, 9, 2, 6, 4, -1, -1, -1, -1},
{0, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{8, 3, 2, 8, 2, 4, 4, 2, 6, -1, -1, -1, -1, -1, -1, -1},
{10, 4, 9, 10, 6, 4, 11, 2, 3, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 2, 2, 8, 11, 4, 9, 10, 4, 10, 6, -1, -1, -1, -1},
{3, 11, 2, 0, 1, 6, 0, 6, 4, 6, 1, 10, -1, -1, -1, -1},
{6, 4, 1, 6, 1, 10, 4, 8, 1, 2, 1, 11, 8, 11, 1, -1},
{9, 6, 4, 9, 3, 6, 9, 1, 3, 11, 6, 3, -1, -1, -1, -1},
{8, 11, 1, 8, 1, 0, 11, 6, 1, 9, 1, 4, 6, 4, 1, -1},
{3, 11, 6, 3, 6, 0, 0, 6, 4, -1, -1, -1, -1, -1, -1, -1},
{6, 4, 8, 11, 6, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{7, 10, 6, 7, 8, 10, 8, 9, 10, -1, -1, -1, -1, -1, -1, -1},
{0, 7, 3, 0, 10, 7, 0, 9, 10, 6, 7, 10, -1, -1, -1, -1},
{10, 6, 7, 1, 10, 7, 1, 7, 8, 1, 8, 0, -1, -1, -1, -1},
{10, 6, 7, 10, 7, 1, 1, 7, 3, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 6, 1, 6, 8, 1, 8, 9, 8, 6, 7, -1, -1, -1, -1},
{2, 6, 9, 2, 9, 1, 6, 7, 9, 0, 9, 3, 7, 3, 9, -1},
{7, 8, 0, 7, 0, 6, 6, 0, 2, -1, -1, -1, -1, -1, -1, -1},
{7, 3, 2, 6, 7, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 11, 10, 6, 8, 10, 8, 9, 8, 6, 7, -1, -1, -1, -1},
{2, 0, 7, 2, 7, 11, 0, 9, 7, 6, 7, 10, 9, 10, 7, -1},
{1, 8, 0, 1, 7, 8, 1, 10, 7, 6, 7, 10, 2, 3, 11, -1},
{11, 2, 1, 11, 1, 7, 10, 6, 1, 6, 7, 1, -1, -1, -1, -1},
{8, 9, 6, 8, 6, 7, 9, 1, 6, 11, 6, 3, 1, 3, 6, -1},
{0, 9, 1, 11, 6, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{7, 8, 0, 7, 0, 6, 3, 11, 0, 11, 6, 0, -1, -1, -1, -1},
{7, 11, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 8, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 9, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{8, 1, 9, 8, 3, 1, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
{10, 1, 2, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 3, 0, 8, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
{2, 9, 0, 2, 10, 9, 6, 11, 7, -1, -1, -1, -1, -1, -1, -1},
{6, 11, 7, 2, 10, 3, 10, 8, 3, 10, 9, 8, -1, -1, -1, -1},
{7, 2, 3, 6, 2, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{7, 0, 8, 7, 6, 0, 6, 2, 0, -1, -1, -1, -1, -1, -1, -1},
{2, 7, 6, 2, 3, 7, 0, 1, 9, -1, -1, -1, -1, -1, -1, -1},
{1, 6, 2, 1, 8, 6, 1, 9, 8, 8, 7, 6, -1, -1, -1, -1},
{10, 7, 6, 10, 1, 7, 1, 3, 7, -1, -1, -1, -1, -1, -1, -1},
{10, 7, 6, 1, 7, 10, 1, 8, 7, 1, 0, 8, -1, -1, -1, -1},
{0, 3, 7, 0, 7, 10, 0, 10, 9, 6, 10, 7, -1, -1, -1, -1},
{7, 6, 10, 7, 10, 8, 8, 10, 9, -1, -1, -1, -1, -1, -1, -1},
{6, 8, 4, 11, 8, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 6, 11, 3, 0, 6, 0, 4, 6, -1, -1, -1, -1, -1, -1, -1},
{8, 6, 11, 8, 4, 6, 9, 0, 1, -1, -1, -1, -1, -1, -1, -1},
{9, 4, 6, 9, 6, 3, 9, 3, 1, 11, 3, 6, -1, -1, -1, -1},
{6, 8, 4, 6, 11, 8, 2, 10, 1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 3, 0, 11, 0, 6, 11, 0, 4, 6, -1, -1, -1, -1},
{4, 11, 8, 4, 6, 11, 0, 2, 9, 2, 10, 9, -1, -1, -1, -1},
{10, 9, 3, 10, 3, 2, 9, 4, 3, 11, 3, 6, 4, 6, 3, -1},
{8, 2, 3, 8, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1},
{0, 4, 2, 4, 6, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 9, 0, 2, 3, 4, 2, 4, 6, 4, 3, 8, -1, -1, -1, -1},
{1, 9, 4, 1, 4, 2, 2, 4, 6, -1, -1, -1, -1, -1, -1, -1},
{8, 1, 3, 8, 6, 1, 8, 4, 6, 6, 10, 1, -1, -1, -1, -1},
{10, 1, 0, 10, 0, 6, 6, 0, 4, -1, -1, -1, -1, -1, -1, -1},
{4, 6, 3, 4, 3, 8, 6, 10, 3, 0, 3, 9, 10, 9, 3, -1},
{10, 9, 4, 6, 10, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 9, 5, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 4, 9, 5, 11, 7, 6, -1, -1, -1, -1, -1, -1, -1},
{5, 0, 1, 5, 4, 0, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
{11, 7, 6, 8, 3, 4, 3, 5, 4, 3, 1, 5, -1, -1, -1, -1},
{9, 5, 4, 10, 1, 2, 7, 6, 11, -1, -1, -1, -1, -1, -1, -1},
{6, 11, 7, 1, 2, 10, 0, 8, 3, 4, 9, 5, -1, -1, -1, -1},
{7, 6, 11, 5, 4, 10, 4, 2, 10, 4, 0, 2, -1, -1, -1, -1},
{3, 4, 8, 3, 5, 4, 3, 2, 5, 10, 5, 2, 11, 7, 6, -1},
{7, 2, 3, 7, 6, 2, 5, 4, 9, -1, -1, -1, -1, -1, -1, -1},
{9, 5, 4, 0, 8, 6, 0, 6, 2, 6, 8, 7, -1, -1, -1, -1},
{3, 6, 2, 3, 7, 6, 1, 5, 0, 5, 4, 0, -1, -1, -1, -1},
{6, 2, 8, 6, 8, 7, 2, 1, 8, 4, 8, 5, 1, 5, 8, -1},
{9, 5, 4, 10, 1, 6, 1, 7, 6, 1, 3, 7, -1, -1, -1, -1},
{1, 6, 10, 1, 7, 6, 1, 0, 7, 8, 7, 0, 9, 5, 4, -1},
{4, 0, 10, 4, 10, 5, 0, 3, 10, 6, 10, 7, 3, 7, 10, -1},
{7, 6, 10, 7, 10, 8, 5, 4, 10, 4, 8, 10, -1, -1, -1, -1},
{6, 9, 5, 6, 11, 9, 11, 8, 9, -1, -1, -1, -1, -1, -1, -1},
{3, 6, 11, 0, 6, 3, 0, 5, 6, 0, 9, 5, -1, -1, -1, -1},
{0, 11, 8, 0, 5, 11, 0, 1, 5, 5, 6, 11, -1, -1, -1, -1},
{6, 11, 3, 6, 3, 5, 5, 3, 1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 10, 9, 5, 11, 9, 11, 8, 11, 5, 6, -1, -1, -1, -1},
{0, 11, 3, 0, 6, 11, 0, 9, 6, 5, 6, 9, 1, 2, 10, -1},
{11, 8, 5, 11, 5, 6, 8, 0, 5, 10, 5, 2, 0, 2, 5, -1},
{6, 11, 3, 6, 3, 5, 2, 10, 3, 10, 5, 3, -1, -1, -1, -1},
{5, 8, 9, 5, 2, 8, 5, 6, 2, 3, 8, 2, -1, -1, -1, -1},
{9, 5, 6, 9, 6, 0, 0, 6, 2, -1, -1, -1, -1, -1, -1, -1},
{1, 5, 8, 1, 8, 0, 5, 6, 8, 3, 8, 2, 6, 2, 8, -1},
{1, 5, 6, 2, 1, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 6, 1, 6, 10, 3, 8, 6, 5, 6, 9, 8, 9, 6, -1},
{10, 1, 0, 10, 0, 6, 9, 5, 0, 5, 6, 0, -1, -1, -1, -1},
{0, 3, 8, 5, 6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{10, 5, 6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{11, 5, 10, 7, 5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{11, 5, 10, 11, 7, 5, 8, 3, 0, -1, -1, -1, -1, -1, -1, -1},
{5, 11, 7, 5, 10, 11, 1, 9, 0, -1, -1, -1, -1, -1, -1, -1},
{10, 7, 5, 10, 11, 7, 9, 8, 1, 8, 3, 1, -1, -1, -1, -1},
{11, 1, 2, 11, 7, 1, 7, 5, 1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 1, 2, 7, 1, 7, 5, 7, 2, 11, -1, -1, -1, -1},
{9, 7, 5, 9, 2, 7, 9, 0, 2, 2, 11, 7, -1, -1, -1, -1},
{7, 5, 2, 7, 2, 11, 5, 9, 2, 3, 2, 8, 9, 8, 2, -1},
{2, 5, 10, 2, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1},
{8, 2, 0, 8, 5, 2, 8, 7, 5, 10, 2, 5, -1, -1, -1, -1},
{9, 0, 1, 5, 10, 3, 5, 3, 7, 3, 10, 2, -1, -1, -1, -1},
{9, 8, 2, 9, 2, 1, 8, 7, 2, 10, 2, 5, 7, 5, 2, -1},
{1, 3, 5, 3, 7, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 7, 0, 7, 1, 1, 7, 5, -1, -1, -1, -1, -1, -1, -1},
{9, 0, 3, 9, 3, 5, 5, 3, 7, -1, -1, -1, -1, -1, -1, -1},
{9, 8, 7, 5, 9, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{5, 8, 4, 5, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1},
{5, 0, 4, 5, 11, 0, 5, 10, 11, 11, 3, 0, -1, -1, -1, -1},
{0, 1, 9, 8, 4, 10, 8, 10, 11, 10, 4, 5, -1, -1, -1, -1},
{10, 11, 4, 10, 4, 5, 11, 3, 4, 9, 4, 1, 3, 1, 4, -1},
{2, 5, 1, 2, 8, 5, 2, 11, 8, 4, 5, 8, -1, -1, -1, -1},
{0, 4, 11, 0, 11, 3, 4, 5, 11, 2, 11, 1, 5, 1, 11, -1},
{0, 2, 5, 0, 5, 9, 2, 11, 5, 4, 5, 8, 11, 8, 5, -1},
{9, 4, 5, 2, 11, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 5, 10, 3, 5, 2, 3, 4, 5, 3, 8, 4, -1, -1, -1, -1},
{5, 10, 2, 5, 2, 4, 4, 2, 0, -1, -1, -1, -1, -1, -1, -1},
{3, 10, 2, 3, 5, 10, 3, 8, 5, 4, 5, 8, 0, 1, 9, -1},
{5, 10, 2, 5, 2, 4, 1, 9, 2, 9, 4, 2, -1, -1, -1, -1},
{8, 4, 5, 8, 5, 3, 3, 5, 1, -1, -1, -1, -1, -1, -1, -1},
{0, 4, 5, 1, 0, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{8, 4, 5, 8, 5, 3, 9, 0, 5, 0, 3, 5, -1, -1, -1, -1},
{9, 4, 5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 11, 7, 4, 9, 11, 9, 10, 11, -1, -1, -1, -1, -1, -1, -1},
{0, 8, 3, 4, 9, 7, 9, 11, 7, 9, 10, 11, -1, -1, -1, -1},
{1, 10, 11, 1, 11, 4, 1, 4, 0, 7, 4, 11, -1, -1, -1, -1},
{3, 1, 4, 3, 4, 8, 1, 10, 4, 7, 4, 11, 10, 11, 4, -1},
{4, 11, 7, 9, 11, 4, 9, 2, 11, 9, 1, 2, -1, -1, -1, -1},
{9, 7, 4, 9, 11, 7, 9, 1, 11, 2, 11, 1, 0, 8, 3, -1},
{11, 7, 4, 11, 4, 2, 2, 4, 0, -1, -1, -1, -1, -1, -1, -1},
{11, 7, 4, 11, 4, 2, 8, 3, 4, 3, 2, 4, -1, -1, -1, -1},
{2, 9, 10, 2, 7, 9, 2, 3, 7, 7, 4, 9, -1, -1, -1, -1},
{9, 10, 7, 9, 7, 4, 10, 2, 7, 8, 7, 0, 2, 0, 7, -1},
{3, 7, 10, 3, 10, 2, 7, 4, 10, 1, 10, 0, 4, 0, 10, -1},
{1, 10, 2, 8, 7, 4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 9, 1, 4, 1, 7, 7, 1, 3, -1, -1, -1, -1, -1, -1, -1},
{4, 9, 1, 4, 1, 7, 0, 8, 1, 8, 7, 1, -1, -1, -1, -1},
{4, 0, 3, 7, 4, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{4, 8, 7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{9, 10, 8, 10, 11, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 9, 3, 9, 11, 11, 9, 10, -1, -1, -1, -1, -1, -1, -1},
{0, 1, 10, 0, 10, 8, 8, 10, 11, -1, -1, -1, -1, -1, -1, -1},
{3, 1, 10, 11, 3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 2, 11, 1, 11, 9, 9, 11, 8, -1, -1, -1, -1, -1, -1, -1},
{3, 0, 9, 3, 9, 11, 1, 2, 9, 2, 11, 9, -1, -1, -1, -1},
{0, 2, 11, 8, 0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{3, 2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 8, 2, 8, 10, 10, 8, 9, -1, -1, -1, -1, -1, -1, -1},
{9, 10, 2, 0, 9, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{2, 3, 8, 2, 8, 10, 0, 1, 8, 1, 10, 8, -1, -1, -1, -1},
{1, 10, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{1, 3, 8, 9, 1, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 9, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}};

523
trick_sims/SIM_fluid/models/fluid/graphics/src/opengl.cpp Normal file
View File

@ -0,0 +1,523 @@
#include <GL/glew.h>
#include <fstream>
#include <iostream>
#include <cstring>
#include <vector>
#include <unistd.h>
#include "../../include/Fluid.hh"
#include "gui.h"
#include "grid_cell.h"
#include "marching_cubes.h"
#include "glm/gtx/io.hpp"
const char* particle_vertex_shader =
#include "shaders/particle.vert"
;
const char* particle_geometry_shader =
#include "shaders/particle.geom"
;
const char* particle_fragment_shader =
#include "shaders/particle.frag"
;
const char* mesh_vertex_shader =
#include "shaders/mesh.vert"
;
const char* mesh_geometry_shader =
#include "shaders/mesh.geom"
;
const char* mesh_fragment_shader =
#include "shaders/mesh.frag"
;
int window_width = 960;
int window_height = 720;
int main_view_width = 960;
int main_view_height = 720;
int preview_width = window_width - main_view_width;
int preview_height = preview_width / 4 * 3;
int preview_bar_width = preview_width;
int preview_bar_height = main_view_height;
const std::string window_title = "Animation";
bool paused = true;
bool meshMode = false;
enum {kVertexBufferMesh, kIndexBufferMesh, kNumVbosMesh};
enum {kVertexBufferParticle, kParticleIndex, kIndexBufferParticle, kNumVbosParticle};
GLuint mesh_buffer_objects[kNumVbosMesh];
GLuint particle_buffer_objects[kNumVbosParticle];
GLuint meshVAO;
GLint mesh_projection_matrix_location = 0;
GLint mesh_view_matrix_location = 0;
GLint mesh_model_matrix_location = 0;
GLint mesh_light_position_location = 0;
GLuint mesh_program_id = 0;
std::vector<GridCell> gridCells;
std::vector<glm::vec4> mesh_vertices;
std::vector<glm::uvec3> mesh_faces;
GLuint particleVAO;
GLint particle_projection_matrix_location = 0;
GLint particle_view_matrix_location = 0;
GLint particle_model_matrix_location = 0;
GLint particle_light_position_location = 0;
GLuint particle_program_id = 0;
std::vector<glm::vec4> particle_vertices;
std::vector<glm::uvec3> particle_faces;
std::vector<int> particle_idx;
GLuint sphTBO;
GLuint sph_tbo_texture;
GLint sph_tbo_pos;
std::vector<float> particlePositions;
MatrixPointers mats;
glm::vec4 light_position;
const int ONE_SECOND = 1000000;
GLFWwindow *window;
void ErrorCallback(int error, const char* description) {
std::cerr << "GLFW Error: " << description << "\n";
}
GLFWwindow* init_glefw()
{
if (!glfwInit())
exit(EXIT_FAILURE);
glfwSetErrorCallback(ErrorCallback);
glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 4);
glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 1);
glfwWindowHint(GLFW_OPENGL_FORWARD_COMPAT, GL_TRUE);
glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE);
glfwWindowHint(GLFW_RESIZABLE, GL_FALSE); // Disable resizing, for simplicity
glfwWindowHint(GLFW_SAMPLES, 4);
auto ret = glfwCreateWindow(window_width, window_height, window_title.data(), nullptr, nullptr);
glfwMakeContextCurrent(ret);
glewExperimental = GL_TRUE;
glewInit();
glGetError(); // clear GLEW's error for it
glfwSwapInterval(1);
const GLubyte* renderer = glGetString(GL_RENDERER);
const GLubyte* version = glGetString(GL_VERSION);
std::cout << "Renderer: " << renderer << "\n";
std::cout << "OpenGL version supported:" << version << "\n";
return ret;
}
void
loadObj(const std::string& file, std::vector<glm::vec4>& vertices,
std::vector<glm::uvec3>& indices, int cur_vertex)
{
std::ifstream obj_file(file);
std::string line;
while (getline(obj_file, line)) {
if (line.at(0) == 'v') {
const char* buf = line.c_str();
char* ptr = (char*)malloc(strlen(buf) + 1);
strcpy(ptr, buf);
char* token = strtok(ptr, " ");
std::vector<float> line_vertices;
int i = 0;
while (token != NULL) {
if (i > 0) {
line_vertices.push_back(atof(token));
}
token = strtok(NULL, " ");
i++;
}
vertices.push_back(glm::vec4(line_vertices[0], line_vertices[1], line_vertices[2], 1.0f));
free(ptr);
}
else if (line.at(0) == 'f') {
const char* buf = line.c_str();
char* ptr = (char*)malloc(strlen(buf) + 1);
strcpy(ptr, buf);
char* token = strtok(ptr, " ");
std::vector<int> line_indicies;
int i = 0;
while (token != NULL) {
if (i > 0) {
line_indicies.push_back(atoi(token));
}
token = strtok(NULL, " ");
i++;
}
indices.push_back(glm::uvec3(cur_vertex + line_indicies[0] - 1, cur_vertex + line_indicies[1] - 1, cur_vertex + line_indicies[2] - 1));
free(ptr);
}
}
}
void meshGLInit() {
/* Initialize mesh OpenGL Program */
// Setup VAO for mesh program
glGenVertexArrays(1, &meshVAO);
// Bind VAO for mesh program
glBindVertexArray(meshVAO);
// Generate VBOs for mesh program
glGenBuffers(kNumVbosMesh, &mesh_buffer_objects[0]);
// Setup particle vertex data in VBO
glBindBuffer(GL_ARRAY_BUFFER, mesh_buffer_objects[kVertexBufferMesh]);
// Describe vertex data to OpenGL
glBufferData(GL_ARRAY_BUFFER,
sizeof(float) * mesh_vertices.size() * 4, mesh_vertices.data(),
GL_DYNAMIC_DRAW);
glVertexAttribPointer(0, 4, GL_FLOAT, GL_FALSE, 0, 0);
glEnableVertexAttribArray(0);
// Setup element array buffer.
glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, mesh_buffer_objects[kIndexBufferMesh]);
// Describe elemnt array buffer to OpenGL
glBufferData(GL_ELEMENT_ARRAY_BUFFER,
sizeof(uint32_t) * mesh_faces.size() * 3,
mesh_faces.data(), GL_DYNAMIC_DRAW);
// Setup mesh vertex shader.
GLuint mesh_vertex_shader_id = 0;
mesh_vertex_shader_id = glCreateShader(GL_VERTEX_SHADER);
glShaderSource(mesh_vertex_shader_id, 1, &mesh_vertex_shader, nullptr);
glCompileShader(mesh_vertex_shader_id);
// Setup mesh geometry shader.
GLuint mesh_geometry_shader_id = 0;
mesh_geometry_shader_id = glCreateShader(GL_GEOMETRY_SHADER);
glShaderSource(mesh_geometry_shader_id, 1, &mesh_geometry_shader, nullptr);
glCompileShader(mesh_geometry_shader_id);
// Setup mesh fragment shader.
GLuint mesh_fragment_shader_id = 0;
mesh_fragment_shader_id = glCreateShader(GL_FRAGMENT_SHADER);
glShaderSource(mesh_fragment_shader_id, 1, &mesh_fragment_shader, nullptr);
glCompileShader(mesh_fragment_shader_id);
// Setup program for the particles, and get its locations.
mesh_program_id = glCreateProgram();
glAttachShader(mesh_program_id, mesh_vertex_shader_id);
glAttachShader(mesh_program_id, mesh_geometry_shader_id);
glAttachShader(mesh_program_id, mesh_fragment_shader_id);
// Bind attributes.
glBindAttribLocation(mesh_program_id, 0, "vertex_position");
glBindFragDataLocation(mesh_program_id, 0, "fragment_color");
glLinkProgram(mesh_program_id);
// Get the mesh uniform locations.
mesh_projection_matrix_location =
glGetUniformLocation(mesh_program_id, "projection");
mesh_view_matrix_location =
glGetUniformLocation(mesh_program_id, "view");
mesh_model_matrix_location =
glGetUniformLocation(mesh_program_id, "model");
mesh_light_position_location =
glGetUniformLocation(mesh_program_id, "light_position");
}
void particleGLInit() {
/* Initialize particle OpenGL Program */
// Setup VAO for particle program
glGenVertexArrays(1, &particleVAO);
// Bind VAO for particle program
glBindVertexArray(particleVAO);
// Generate VBOs for particle program
glGenBuffers(kNumVbosParticle, &particle_buffer_objects[0]);
// Setup particle vertex data in VBO
glBindBuffer(GL_ARRAY_BUFFER, particle_buffer_objects[kVertexBufferParticle]);
// Describe vertex data to OpenGL
glBufferData(GL_ARRAY_BUFFER,
sizeof(float) * particle_vertices.size() * 4, particle_vertices.data(),
GL_DYNAMIC_DRAW);
glVertexAttribPointer(0, 4, GL_FLOAT, GL_FALSE, 0, 0);
glEnableVertexAttribArray(0);
// Setup particle index datas in VBO
glBindBuffer(GL_ARRAY_BUFFER, particle_buffer_objects[kParticleIndex]);
// Describe particle index data to OpenGL
glBufferData(GL_ARRAY_BUFFER, sizeof(int) * particle_idx.size(), particle_idx.data(), GL_STATIC_DRAW);
glVertexAttribIPointer(1, 1, GL_INT, 0, 0);
glEnableVertexAttribArray(1);
// Setup element array buffer.
glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, particle_buffer_objects[kIndexBufferParticle]);
// Describe element array buffer to OpenGL
glBufferData(GL_ELEMENT_ARRAY_BUFFER,
sizeof(uint32_t) * particle_faces.size() * 3,
particle_faces.data(), GL_DYNAMIC_DRAW);
// Setup particle vertex shader.
GLuint particle_vertex_shader_id = 0;
particle_vertex_shader_id = glCreateShader(GL_VERTEX_SHADER);
glShaderSource(particle_vertex_shader_id, 1, &particle_vertex_shader, nullptr);
glCompileShader(particle_vertex_shader_id);
// Setup particle geometry shader.
GLuint particle_geometry_shader_id = 0;
particle_geometry_shader_id = glCreateShader(GL_GEOMETRY_SHADER);
glShaderSource(particle_geometry_shader_id, 1, &particle_geometry_shader, nullptr);
glCompileShader(particle_geometry_shader_id);
// Setup particle fragment shader.
GLuint particle_fragment_shader_id = 0;
particle_fragment_shader_id = glCreateShader(GL_FRAGMENT_SHADER);
glShaderSource(particle_fragment_shader_id, 1, &particle_fragment_shader, nullptr);
glCompileShader(particle_fragment_shader_id);
// Setup program for the particles, and get its locations.
particle_program_id = glCreateProgram();
glAttachShader(particle_program_id, particle_vertex_shader_id);
glAttachShader(particle_program_id, particle_geometry_shader_id);
glAttachShader(particle_program_id, particle_fragment_shader_id);
// Bind attributes.
glBindAttribLocation(particle_program_id, 0, "vertex_position");
glBindAttribLocation(particle_program_id, 1, "idx");
glBindFragDataLocation(particle_program_id, 0, "fragment_color");
glLinkProgram(particle_program_id);
// Get the particle uniform locations.
particle_projection_matrix_location =
glGetUniformLocation(particle_program_id, "projection");
particle_view_matrix_location =
glGetUniformLocation(particle_program_id, "view");
particle_model_matrix_location =
glGetUniformLocation(particle_program_id, "model");
particle_light_position_location =
glGetUniformLocation(particle_program_id, "light_position");
}
void particleTBOInit() {
/* Setup TBO for particle position data*/
glGenBuffers(1, &sphTBO);
glBindBuffer(GL_TEXTURE_BUFFER, sphTBO);
glBufferData(GL_TEXTURE_BUFFER, sizeof(float) * particlePositions.size(), particlePositions.data(), GL_DYNAMIC_DRAW);
glGenTextures(1, &sph_tbo_texture);
glBindBuffer(GL_TEXTURE_BUFFER, 0);
sph_tbo_pos = glGetUniformLocation(particle_program_id, "sph_tbo_pos");
}
void setupParticleProgram() {
// Switch VAO
glBindVertexArray(particleVAO);
// Switch Program
glUseProgram(particle_program_id);
// Setup TBO
glBindBuffer(GL_TEXTURE_BUFFER, sphTBO);
glBufferData(GL_TEXTURE_BUFFER, sizeof(float) * particlePositions.size(), particlePositions.data(), GL_DYNAMIC_DRAW);
// Pass uniforms
glUniformMatrix4fv(particle_projection_matrix_location, 1, GL_FALSE,
&((*mats.projection)[0][0]));
glUniformMatrix4fv(particle_view_matrix_location, 1, GL_FALSE,
&((*mats.view)[0][0]));
glUniformMatrix4fv(particle_model_matrix_location, 1, GL_FALSE,
&((*mats.model)[0][0]));
glUniform4fv(particle_light_position_location, 1, &light_position[0]);
glActiveTexture(GL_TEXTURE0);
glBindTexture(GL_TEXTURE_BUFFER, sph_tbo_texture);
glTexBuffer(GL_TEXTURE_BUFFER, GL_R32F, sphTBO);
glUniform1i(sph_tbo_pos, 0);
}
void setupMeshProgram() {
// Switch VAO
glBindVertexArray(meshVAO);
// Switch Program
glUseProgram(mesh_program_id);
// Pass uniforms
glUniformMatrix4fv(mesh_projection_matrix_location, 1, GL_FALSE,
&((*mats.projection)[0][0]));
glUniformMatrix4fv(mesh_view_matrix_location, 1, GL_FALSE,
&((*mats.view)[0][0]));
glUniformMatrix4fv(mesh_model_matrix_location, 1, GL_FALSE,
&((*mats.model)[0][0]));
glUniform4fv(mesh_light_position_location, 1, &light_position[0]);
}
void updateMesh(Fluid* fluid) {
updateIsoValues(gridCells, particlePositions, fluid->ISO_RADIUS);
mesh_vertices.clear();
mesh_faces.clear();
gridCells.clear();
initializeGridCells(gridCells, fluid->BOUND, fluid->MC_GRID_DIM);
updateIsoValues(gridCells, particlePositions, fluid->ISO_RADIUS);
for (int i = 0; i < gridCells.size(); i++) {
generateCellMesh(gridCells[i], fluid->PARTICLES_WITHIN_VERTEX, mesh_faces, mesh_vertices);
}
glBindVertexArray(meshVAO);
glBindBuffer(GL_ARRAY_BUFFER, mesh_buffer_objects[kVertexBufferMesh]);
glBufferData(GL_ARRAY_BUFFER, 4 * sizeof(float) * mesh_vertices.size(), mesh_vertices.data(), GL_DYNAMIC_DRAW);
glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, mesh_buffer_objects[kIndexBufferMesh]);
glBufferData(GL_ELEMENT_ARRAY_BUFFER,
sizeof(uint32_t) * mesh_faces.size() * 3,
mesh_faces.data(), GL_DYNAMIC_DRAW);
}
void initGeometry(Fluid* fluid) {
loadObj(fluid->obj_file, particle_vertices, particle_faces, 0);
int num_vert = particle_vertices.size();
for (int i = 1; i < fluid->NUM_PARTICLES; i++) {
loadObj(fluid->obj_file, particle_vertices, particle_faces, i * num_vert);
}
/* Assign index to each particle model */
for(int i = 0; i < fluid->NUM_PARTICLES; i++) {
for (int p = 0; p < num_vert; p++) {
particle_idx.push_back(i);
}
}
particlePositions = fluid->getParticlePositions();
initializeGridCells(gridCells, fluid->BOUND, fluid->MC_GRID_DIM);
updateIsoValues(gridCells, particlePositions, fluid->ISO_RADIUS);
for (int i = 0; i < gridCells.size(); i++) {
generateCellMesh(gridCells[i], fluid->PARTICLES_WITHIN_VERTEX, mesh_faces, mesh_vertices);
}
glm::vec4 light_position = glm::vec4(0.0f, 0.0f, fluid->BOUND + 10.0f, 1.0f);
}
void configureViewport() {
glfwGetFramebufferSize(window, &window_width, &window_height);
glViewport(0, 0, window_width, window_height);
glClearColor(0.0f, 0.0f, 0.0f, 0.0f);
glEnable(GL_DEPTH_TEST);
glEnable(GL_MULTISAMPLE);
glEnable(GL_BLEND);
glEnable(GL_CULL_FACE);
glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
glDepthFunc(GL_LESS);
glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
glCullFace(GL_BACK);
}
int openGLMain(Fluid* fluid)
{
window = init_glefw();
GUI gui(window, window_width, window_height, preview_height);
initGeometry(fluid);
meshGLInit();
particleGLInit();
particleTBOInit();
while (!glfwWindowShouldClose(window)) {
configureViewport();
gui.updateMatrices();
mats = gui.getMatrixPointers();
particlePositions = fluid->getParticlePositions();
if (meshMode) {
if (fluid->timeSteps % fluid->mcUpdateFreq == 0)
updateMesh(fluid);
setupMeshProgram();
glDrawElements(GL_TRIANGLES, 3 * mesh_faces.size(), GL_UNSIGNED_INT, 0);
} else {
setupParticleProgram();
glDrawElements(GL_TRIANGLES, 3 * particle_faces.size(), GL_UNSIGNED_INT, 0);
}
glfwPollEvents();
glfwSwapBuffers(window);
usleep(ONE_SECOND / fluid->refreshRate);
}
glfwDestroyWindow(window);
glfwTerminate();
exit(EXIT_SUCCESS);
}

1
trick_sims/SIM_fluid/models/fluid/graphics/src/opengl.h Normal file
View File

@ -0,0 +1 @@
int openGLMain(Fluid* fluid);

10
trick_sims/SIM_fluid/models/fluid/graphics/src/opengl_caller.cpp Normal file
View File

@ -0,0 +1,10 @@
#include <thread>
#include "../../include/Fluid.hh"
#include "opengl.h"
void openGLCaller(Fluid* fluid) {
{
std::thread openGLThread(openGLMain, fluid);
openGLThread.detach();
}
}

1
trick_sims/SIM_fluid/models/fluid/graphics/src/opengl_caller.h Normal file
View File

@ -0,0 +1 @@
void openGLCaller(Fluid* fluid);

15
trick_sims/SIM_fluid/models/fluid/graphics/src/shaders/mesh.frag Normal file
View File

@ -0,0 +1,15 @@
R"zzz(
#version 330 core
out vec4 fragment_color;
in vec4 face_normal;
in vec4 light_direction;
in vec4 world_position;
void main() {
float dot_nl = clamp(dot(vec3(light_direction), vec3(face_normal)), 0, 1);
vec3 ambient = vec3(0.2, 0.6, 1) * .1;
vec3 diffuse = dot_nl * vec3(0.2, 0.6, 1);
fragment_color = vec4(diffuse + ambient, 1.0);
}
)zzz"

40
trick_sims/SIM_fluid/models/fluid/graphics/src/shaders/mesh.geom Normal file
View File

@ -0,0 +1,40 @@
R"zzz(#version 330 core
layout (triangles) in;
layout (triangle_strip, max_vertices = 3) out;
uniform mat4 projection;
uniform mat4 model;
uniform mat4 view;
uniform vec4 light_position;
in vec4 vs_light_direction[];
in vec4 vs_camera_direction[];
in vec4 vs_normal[];
in vec2 vs_uv[];
out vec4 face_normal;
out vec4 light_direction;
out vec4 camera_direction;
out vec4 world_position;
out vec4 vertex_normal;
out vec2 uv_coords;
in int vs_idx[];
flat out int idx;
void main() {
int n = 0;
vec3 a = gl_in[0].gl_Position.xyz;
vec3 b = gl_in[1].gl_Position.xyz;
vec3 c = gl_in[2].gl_Position.xyz;
vec3 u = normalize(b - a);
vec3 v = normalize(c - a);
face_normal = normalize(vec4(normalize(cross(u, v)), 0.0));
for (n = 0; n < gl_in.length(); n++) {
light_direction = normalize(vs_light_direction[n]);
camera_direction = normalize(vs_camera_direction[n]);
world_position = gl_in[n].gl_Position;
vertex_normal = vs_normal[n];
uv_coords = vs_uv[n];
idx = vs_idx[n];
gl_Position = projection * view * model * gl_in[n].gl_Position;
EmitVertex();
}
EndPrimitive();
}
)zzz"

20
trick_sims/SIM_fluid/models/fluid/graphics/src/shaders/mesh.vert Normal file
View File

@ -0,0 +1,20 @@
R"zzz(#version 330 core
in vec4 vertex_position;
uniform vec4 light_position;
uniform mat4 projection;
uniform mat4 model;
uniform mat4 view;
out vec4 vs_light_direction;
out int vs_idx;
void main() {
vec3 vertex_pos = vertex_position.xyz;
gl_Position = vec4(vertex_pos, 1);
vs_light_direction = light_position - gl_Position;
}
)zzz"

15
trick_sims/SIM_fluid/models/fluid/graphics/src/shaders/particle.frag Normal file
View File

@ -0,0 +1,15 @@
R"zzz(
#version 330 core
out vec4 fragment_color;
in vec4 face_normal;
in vec4 light_direction;
in vec4 world_position;
void main() {
float dot_nl = clamp(dot(vec3(light_direction), vec3(face_normal)), 0, 1);
vec3 ambient = vec3(0.2, 0.6, 1) * .1;
vec3 diffuse = dot_nl * vec3(0.2, 0.6, 1);
fragment_color = vec4(diffuse + ambient, 1.0);
}
)zzz"

40
trick_sims/SIM_fluid/models/fluid/graphics/src/shaders/particle.geom Normal file
View File

@ -0,0 +1,40 @@
R"zzz(#version 330 core
layout (triangles) in;
layout (triangle_strip, max_vertices = 3) out;
uniform mat4 projection;
uniform mat4 model;
uniform mat4 view;
uniform vec4 light_position;
in vec4 vs_light_direction[];
in vec4 vs_camera_direction[];
in vec4 vs_normal[];
in vec2 vs_uv[];
out vec4 face_normal;
out vec4 light_direction;
out vec4 camera_direction;
out vec4 world_position;
out vec4 vertex_normal;
out vec2 uv_coords;
in int vs_idx[];
flat out int idx;
void main() {
int n = 0;
vec3 a = gl_in[0].gl_Position.xyz;
vec3 b = gl_in[1].gl_Position.xyz;
vec3 c = gl_in[2].gl_Position.xyz;
vec3 u = normalize(b - a);
vec3 v = normalize(c - a);
face_normal = normalize(vec4(normalize(cross(u, v)), 0.0));
for (n = 0; n < gl_in.length(); n++) {
light_direction = normalize(vs_light_direction[n]);
camera_direction = normalize(vs_camera_direction[n]);
world_position = gl_in[n].gl_Position;
vertex_normal = vs_normal[n];
uv_coords = vs_uv[n];
idx = vs_idx[n];
gl_Position = projection * view * model * gl_in[n].gl_Position;
EmitVertex();
}
EndPrimitive();
}
)zzz"

27
trick_sims/SIM_fluid/models/fluid/graphics/src/shaders/particle.vert Normal file
View File

@ -0,0 +1,27 @@
R"zzz(#version 330 core
in vec4 vertex_position;
in int idx;
uniform samplerBuffer sph_tbo_pos;
uniform vec4 light_position;
uniform mat4 projection;
uniform mat4 model;
uniform mat4 view;
out vec4 vs_light_direction;
out int vs_idx;
void main() {
vs_idx = idx;
vec3 vertex_pos = vertex_position.xyz;
float x = texelFetch(sph_tbo_pos, 3 * idx).r;
float y = texelFetch(sph_tbo_pos, 3 * idx + 1).r;
float z = texelFetch(sph_tbo_pos, 3 * idx + 2).r;
vec3 particle_pos = vec3(x, y, z);
gl_Position = vec4(1 * vertex_pos + particle_pos, 1);
vs_light_direction = light_position - gl_Position;
}
)zzz"

90
trick_sims/SIM_fluid/models/fluid/include/Fluid.hh Normal file
View File

@ -0,0 +1,90 @@
/********************************* TRICK HEADER *******************************
PURPOSE: (Simulate a fluid using smoothed particle hydrodynamics (SPH).)
**************************************************************************/
#ifndef _fluid_hh_
#define _fluid_hh_
#include "Particle.hh"
#include <vector>
#include <unordered_map>
#include <cmath>
class Fluid {
private:
void computeDensityAndPressure(int p_start, int p_end);
void computeForces(int p_start, int p_end);
void verletUpdatePosition(int p_start, int p_end);
void timeIntegration(int p_start, int p_end);
void buildSpatialGrid();
std::vector<Particle> getNeighbors(int grid_x, int grid_y, int grid_z);
bool checkBounds(int gridX, int gridY, int gridZ);
public:
Fluid(){};
std::unordered_map<int, std::vector<Particle>> spatialGrid;
std::vector<float> getParticlePositions();
std::vector<Particle> particles;
int default_data();
int init_cuda();
int update_SPH();
int shutdown_cuda();
float G_STRENGTH = 600;
// external graviational force
float G[3] = {0.f, G_STRENGTH * -9.8f, 0.f};
float REST_DENS = 1000.f;
float GAS_CONST = 100.f;
// kernel radius
float H = 16.f;
float HSQ = H * H;
// assume all particles have the same mass
float MASS = 65.f;
// viscosity constant
float VISC = 700.f;
// timestep size
float DT = 0.001f;
// smoothing kernels and their gradients
float POLY6 = 315.f / (64.f * M_PI * std::pow(H, 9.f));
float SPIKY_GRAD = -45.f / (M_PI * std::pow(H, 6.f));
float VISC_LAP = 45.f / (M_PI * std::pow(H, 6.f));
// number of particles on the edge of the square
int EDGE_NUM_PARTICLES = 20;
int PARTICLE_DEPTH = 4;
int NUM_PARTICLES = PARTICLE_DEPTH * EDGE_NUM_PARTICLES * EDGE_NUM_PARTICLES;
float PARTICLE_DIST = 5;
Particle* particlesArr;
// boundary parameters
float EPS = H; // boundary epsilon
float BOUND_DAMPING = -.5f;
float BOUND = 200;
// parameters for spatial grid
int CELL_SIZE = 2 * H;
int CELLS_PER_DIM = 2 * BOUND / CELL_SIZE;
// number of timesteps simulated
int timeSteps = 0;
// parameters for marching cubes
int MC_GRID_DIM = 16;
int ISO_RADIUS = H;
int PARTICLES_WITHIN_VERTEX = 1; // number of particles within MC GridCell vertex;
// simulation modes
bool gpuMode = true;
bool cpuNeighborList = true;
int refreshRate = 60; // set refresh rate for graphics
int mcUpdateFreq = 1; // update the fluid mesh every multiple of mcUpdateFreq frames
std::string obj_file = "100_sphere.obj";
};
#endif

22
trick_sims/SIM_fluid/models/fluid/include/Particle.hh Normal file
View File

@ -0,0 +1,22 @@
/********************************* TRICK HEADER *******************************
PURPOSE: (Simulate fluid particles in SPH)
*******************************************************************************/
#ifndef _particle_hh_
#define _particle_hh_
struct Particle {
Particle(float x_, float y_, float z_) :
pos {x_, y_, z_},
velocity {0.f, 0.f, 0.f},
force {0.f, 0.f, 0.f},
rho(0),
pressure(0.f) {}
float pos[3];
float velocity[3];
float force[3];
float rho;
float pressure;
};
#endif

1
trick_sims/SIM_fluid/models/fluid/src/.gitignore vendored Normal file
View File

@ -0,0 +1 @@
*.o

248
trick_sims/SIM_fluid/models/fluid/src/Fluid.cpp Normal file
View File

@ -0,0 +1,248 @@
/********************************* TRICK HEADER *******************************
PURPOSE: (Simulate a fluid using smoothed particle hydrodynamics (SPH).)
LIBRARY DEPENDENCY:
((Fluid.o))
*******************************************************************************/
#include "../include/Fluid.hh"
#include "../include/Particle.hh"
#include "../graphics/src/opengl_caller.h"
#include "sph_gpu.h"
#include <iostream>
#include <iterator>
int Fluid::default_data() {
for (int i = 0; i < EDGE_NUM_PARTICLES; i++) {
for (int j = 0; j < EDGE_NUM_PARTICLES; j++) {
for (int k = 0; k < PARTICLE_DEPTH; k++) {
particles.push_back(Particle(PARTICLE_DIST * i, PARTICLE_DIST * j, PARTICLE_DIST * k));
}
}
}
particlesArr = &particles[0];
openGLCaller(this);
return 0;
}
int Fluid::init_cuda() {
initSPH_GPU(particles, this);
return 0;
}
void Fluid::buildSpatialGrid() {
spatialGrid.clear();
for (auto& pi : particles) {
int grid_x = CELLS_PER_DIM * ((pi.pos[0] + BOUND) / (2 * BOUND));
int grid_y = CELLS_PER_DIM * ((pi.pos[1] + BOUND) / (2 * BOUND));
int grid_z = CELLS_PER_DIM * ((pi.pos[2] + BOUND) / (2 * BOUND));
int grid_key = grid_x + grid_y * CELLS_PER_DIM + grid_z * CELLS_PER_DIM * CELLS_PER_DIM;
if (spatialGrid.find(grid_key) == spatialGrid.end()) {
std::vector<Particle> single_particle;
single_particle.push_back(pi);
spatialGrid[grid_key] = single_particle;
}
else {
spatialGrid[grid_key].push_back(pi);
}
}
}
bool Fluid::checkBounds(int x, int y, int z) {
return x >= 0 && x < CELLS_PER_DIM&& y >= 0 && y < CELLS_PER_DIM&& z >= 0 && z < CELLS_PER_DIM;
}
std::vector<Particle> Fluid::getNeighbors(int grid_x, int grid_y, int grid_z) {
std::vector<Particle> neighbors;
for (int x = grid_x - 1; x <= grid_x + 1; x++) {
for (int y = grid_y - 1; y <= grid_y + 1; y++) {
for (int z = grid_z - 1; z <= grid_z + 1; z++) {
int grid_key = x + y * CELLS_PER_DIM + z * CELLS_PER_DIM * CELLS_PER_DIM;
if (checkBounds(x, y, z) && spatialGrid.find(grid_key) != spatialGrid.end()) {
std::vector<Particle> cell_particles = spatialGrid[grid_key];
neighbors.insert(neighbors.end(), cell_particles.begin(), cell_particles.end());
}
}
}
}
return neighbors;
}
int Fluid::update_SPH() {
if (gpuMode) {
// GPU simulation
updateSPH_GPU(particles, this);
} else {
// CPU simulation
if (cpuNeighborList) {
buildSpatialGrid();
}
int p_start = 0;
int p_end = particles.size();
verletUpdatePosition(p_start, p_end);
computeDensityAndPressure(p_start, p_end);
computeForces(p_start, p_end);
timeIntegration(p_start, p_end);
}
timeSteps++;
return 0;
}
void Fluid::computeDensityAndPressure(int p_start, int p_end) {
for (int i = p_start; i < p_end; i++) {
Particle& pi = particles[i];
pi.rho = 0;
int grid_x = CELLS_PER_DIM * ((pi.pos[0] + BOUND) / (2 * BOUND));
int grid_y = CELLS_PER_DIM * ((pi.pos[1] + BOUND) / (2 * BOUND));
int grid_z = CELLS_PER_DIM * ((pi.pos[2] + BOUND) / (2 * BOUND));
std::vector<Particle> candidate_neighbors;
if(cpuNeighborList) {
candidate_neighbors = getNeighbors(grid_x, grid_y, grid_z);
} else {
candidate_neighbors = particles;
}
for (auto& pj : candidate_neighbors) {
float rij[3] = {pj.pos[0] - pi.pos[0], pj.pos[1] - pi.pos[1], pj.pos[2] - pi.pos[2]};
float r = std::sqrt(rij[0] * rij[0] + rij[1] * rij[1] + rij[2] * rij[2]);
if (r >= 0 && r <= H) {
//printf("size of neighbors: %lu\n", candidate_neighbors.size());
pi.rho += MASS * POLY6 * pow(HSQ - r * r, 3);
}
}
pi.pressure = GAS_CONST * (pi.rho - REST_DENS);
}
}
void Fluid::computeForces(int p_start, int p_end) {
for (int i = p_start; i < p_end; i++) {
Particle& pi = particles[i];
float pressure_force[3] = {0, 0, 0};
float viscosity_force[3] = {0, 0, 0};
int grid_x = CELLS_PER_DIM * ((pi.pos[0] + BOUND) / (2 * BOUND));
int grid_y = CELLS_PER_DIM * ((pi.pos[1] + BOUND) / (2 * BOUND));
int grid_z = CELLS_PER_DIM * ((pi.pos[2] + BOUND) / (2 * BOUND));
std::vector<Particle> candidate_neighbors;
if (cpuNeighborList) {
candidate_neighbors = getNeighbors(grid_x, grid_y, grid_z);
} else {
candidate_neighbors = particles;
}
for (auto& pj : candidate_neighbors) {
if (&pi != &pj) {
float rij[3] = {pj.pos[0] - pi.pos[0], pj.pos[1] - pi.pos[1], pj.pos[2] - pi.pos[2]};
float r = std::sqrt(rij[0] * rij[0] + rij[1] * rij[1] + rij[2] * rij[2]);
float rij_hat[3] = {rij[0] / r, rij[1] / r, rij[2] / r};
if (r > 0 && r <= H) {
pressure_force[0] -= rij_hat[0] * MASS * (pi.pressure + pj.pressure) / (2 * pj.rho) * SPIKY_GRAD * pow(H - r, 2);
pressure_force[1] -= rij_hat[1] * MASS * (pi.pressure + pj.pressure) / (2 * pj.rho) * SPIKY_GRAD * pow(H - r, 2);
pressure_force[2] -= rij_hat[2] * MASS * (pi.pressure + pj.pressure) / (2 * pj.rho) * SPIKY_GRAD * pow(H - r, 2);
viscosity_force[0] += VISC * MASS * ((pj.velocity[0] - pi.velocity[0]) / pj.rho) * VISC_LAP * (H - r);
viscosity_force[1] += VISC * MASS * ((pj.velocity[1] - pi.velocity[1]) / pj.rho) * VISC_LAP * (H - r);
viscosity_force[2] += VISC * MASS * ((pj.velocity[2] - pi.velocity[2]) / pj.rho) * VISC_LAP * (H - r);
}
}
}
float gravity_force[3] = {pi.rho * G[0], pi.rho * G[1], pi.rho * G[2]};
pi.force[0] = viscosity_force[0] + pressure_force[0] + gravity_force[0];
pi.force[1] = viscosity_force[1] + pressure_force[1] + gravity_force[1];
pi.force[2] = viscosity_force[2] + pressure_force[2] + gravity_force[2];
if (std::isnan(pi.force[0]) || std::isnan(pi.force[1]) || std::isnan(pi.force[2])) {
pi.force[0] = gravity_force[0];
pi.force[1] = gravity_force[1];
pi.force[2] = gravity_force[2];
}
}
}
void Fluid::verletUpdatePosition(int p_start, int p_end) {
for (int i = p_start; i < p_end; i++) {
Particle& pi = particles[i];
//pi.pos += DT * pi.velocity;
}
}
// Velocity Verlet time integrator
void Fluid::timeIntegration(int p_start, int p_end) {
for (int i = p_start; i < p_end; i++) {
Particle& pi = particles[i];
pi.pos[0] += DT * pi.velocity[0];
pi.pos[1] += DT * pi.velocity[1];
pi.pos[2] += DT * pi.velocity[2];
pi.velocity[0] += DT * pi.force[0] / pi.rho;
pi.velocity[1] += DT * pi.force[1] / pi.rho;
pi.velocity[2] += DT * pi.force[2] / pi.rho;
if (pi.pos[2] - EPS < -BOUND) {
pi.velocity[2] *= BOUND_DAMPING;
pi.pos[2] = -BOUND + EPS;
}
if (pi.pos[2] + EPS > BOUND) {
pi.velocity[2] *= BOUND_DAMPING;
pi.pos[2] = BOUND - EPS;
}
if (pi.pos[1] - EPS < -BOUND) {
pi.velocity[1] *= BOUND_DAMPING;
pi.pos[1] = -BOUND + EPS;
}
if (pi.pos[1] + EPS > BOUND) {
pi.velocity[1] *= BOUND_DAMPING;
pi.pos[1] = BOUND - EPS;
}
if (pi.pos[0] - EPS < -BOUND) {
pi.velocity[0] *= BOUND_DAMPING;
pi.pos[0] = -BOUND + EPS;
}
if (pi.pos[0] + EPS > BOUND) {
pi.velocity[0] *= BOUND_DAMPING;
pi.pos[0] = BOUND - EPS;
}
}
}
std::vector<float> Fluid::getParticlePositions() {
std::vector<float> positions;
for (auto& pi : particles) {
positions.push_back(pi.pos[0]);
positions.push_back(pi.pos[1]);
positions.push_back(pi.pos[2]);
}
return positions;
}
int Fluid::shutdown_cuda() {
shutdownSPH_GPU(particles, this);
return 0;
}

402
trick_sims/SIM_fluid/models/fluid/src/sph_gpu.cu Normal file
View File

@ -0,0 +1,402 @@
#include "cuda_runtime.h"
#include "sph_gpu.h"
#define NUM_THREADS 1024
Particle* d_particles;
Particle** d_spatial_grid;
int* d_cell_counts;
Fluid* d_fluid;
int* d_n;
__global__ void computeDensityAndPressureGPU(Particle* particles, int* n, Fluid* fluid) {
int tid = threadIdx.x;
// assuming n is a multiple of NUM_THREADS
int block_size = *n / NUM_THREADS;
int p_start = tid * block_size;
int p_end = (tid + 1) * block_size;
for (int i = p_start; i < p_end; i++) {
Particle& pi = particles[i];
pi.rho = 0;
for (int j = 0; j < *n; j++) {
//Particle &pj = candidate_neighbors[j];
Particle& pj = particles[j];
float rij[3] = {pj.pos[0] - pi.pos[0], pj.pos[1] - pi.pos[1], pj.pos[2] - pi.pos[2]};
float r = sqrt(rij[0] * rij[0] + rij[1] * rij[1] + rij[2] * rij[2]);
if (r >= 0 && r <= fluid->H) {
pi.rho += fluid->MASS * fluid->POLY6 * pow(fluid->HSQ - r * r, 3.f);
}
}
/*
int grid_x = fluid->CELLS_PER_DIM * ((pi.pos[0] + fluid->BOUND) / (2 * fluid->BOUND));
int grid_y = fluid->CELLS_PER_DIM * ((pi.pos[1] + fluid->BOUND) / (2 * fluid->BOUND));
int grid_z = fluid->CELLS_PER_DIM * ((pi.pos[2] + fluid->BOUND) / (2 * fluid->BOUND));
for (int x = grid_x - 1; x <= grid_x + 1; x++) {
for (int y = grid_y - 1; y <= grid_y + 1; y++) {
for (int z = grid_z - 1; z <= grid_z + 1; z++) {
if (x >= 0 && x < fluid->CELLS_PER_DIM && y >= 0 && y < fluid->CELLS_PER_DIM && z >= 0 && z < fluid->CELLS_PER_DIM) {
int grid_idx = x + y * fluid->CELLS_PER_DIM + z * fluid->CELLS_PER_DIM * fluid->CELLS_PER_DIM;
int num_neighbors = cell_counts[grid_idx];
if (num_neighbors == 0) {
continue;
}
Particle* candidate_neighbors = spatial_grid[grid_idx];
for (int j = 0; j < num_neighbors; j++) {
Particle &pj = candidate_neighbors[j];
//Particle& pj = particles[j];
float rij[3] = {pj.pos[0] - pi.pos[0], pj.pos[1] - pi.pos[1], pj.pos[2] - pi.pos[2]};
float r = sqrt(rij[0] * rij[0] + rij[1] * rij[1] + rij[2] * rij[2]);
if (r >= 0 && r <= fluid->H) {
pi.rho += fluid->MASS * fluid->POLY6 * pow(fluid->HSQ - r * r, 3.f);
}
}
}
}
}
}
*/
pi.pressure = fluid->GAS_CONST * (pi.rho - fluid->REST_DENS);
}
if (*n % NUM_THREADS != 0 && tid < (*n - NUM_THREADS * block_size)) {
int leftover_start = (*n / NUM_THREADS) * NUM_THREADS;
p_start = tid + leftover_start;
p_end = (tid + 1) + leftover_start;
for (int i = p_start; i < p_end; i++) {
Particle& pi = particles[i];
pi.rho = 0;
for (int j = 0; j < *n; j++) {
//Particle &pj = candidate_neighbors[j];
Particle& pj = particles[j];
float rij[3] = {pj.pos[0] - pi.pos[0], pj.pos[1] - pi.pos[1], pj.pos[2] - pi.pos[2]};
float r = sqrt(rij[0] * rij[0] + rij[1] * rij[1] + rij[2] * rij[2]);
if (r >= 0 && r <= fluid->H) {
pi.rho += fluid->MASS * fluid->POLY6 * pow(fluid->HSQ - r * r, 3.f);
}
}
pi.pressure = fluid->GAS_CONST * (pi.rho - fluid->REST_DENS);
}
}
}
__global__ void computeForcesGPU(Particle* particles, int* n, Fluid* fluid) {
int tid = threadIdx.x;
// assuming n is a multiple of NUM_THREADS
int block_size = *n / NUM_THREADS;
int p_start = tid * block_size;
int p_end = (tid + 1) * block_size;
//for(auto &pi : particles) {
for (int i = p_start; i < p_end; i++) {
Particle& pi = particles[i];
float pressure_force[3] = {0, 0, 0};
float viscosity_force[3] = {0, 0, 0};
//Particle* candidate_neighbors = all_neighbors[i];
for (int j = 0; j < *n; j++) {
//Particle& pj = candidate_neighbors[j];
Particle& pj = particles[j];
if (&pi != &pj) {
float rij[3] = {pj.pos[0] - pi.pos[0], pj.pos[1] - pi.pos[1], pj.pos[2] - pi.pos[2]};
float r = std::sqrt(rij[0] * rij[0] + rij[1] * rij[1] + rij[2] * rij[2]);
float rij_hat[3] = {rij[0] / r, rij[1] / r, rij[2] / r};
if (r > 0 && r <= fluid->H) {
pressure_force[0] -= rij_hat[0] * fluid->MASS * (pi.pressure + pj.pressure) / (2 * pj.rho) * fluid->SPIKY_GRAD * pow(fluid->H - r, 2.f);
pressure_force[1] -= rij_hat[1] * fluid->MASS * (pi.pressure + pj.pressure) / (2 * pj.rho) * fluid->SPIKY_GRAD * pow(fluid->H - r, 2.f);
pressure_force[2] -= rij_hat[2] * fluid->MASS * (pi.pressure + pj.pressure) / (2 * pj.rho) * fluid->SPIKY_GRAD * pow(fluid->H - r, 2.f);
viscosity_force[0] += fluid->VISC * fluid->MASS * ((pj.velocity[0] - pi.velocity[0]) / pj.rho) * fluid->VISC_LAP * (fluid->H - r);
viscosity_force[1] += fluid->VISC * fluid->MASS * ((pj.velocity[1] - pi.velocity[1]) / pj.rho) * fluid->VISC_LAP * (fluid->H - r);
viscosity_force[2] += fluid->VISC * fluid->MASS * ((pj.velocity[2] - pi.velocity[2]) / pj.rho) * fluid->VISC_LAP * (fluid->H - r);
}
}
}
float G[3] = {0.f, fluid->G_STRENGTH * -9.8f, 0.f};
float gravity_force[3] = {pi.rho * G[0], pi.rho * G[1], pi.rho * G[2]};
pi.force[0] = viscosity_force[0] + pressure_force[0] + gravity_force[0];
pi.force[1] = viscosity_force[1] + pressure_force[1] + gravity_force[1];
pi.force[2] = viscosity_force[2] + pressure_force[2] + gravity_force[2];
}
if (*n % NUM_THREADS != 0 && tid < (*n - NUM_THREADS * block_size)) {
int leftover_start = (*n / NUM_THREADS) * NUM_THREADS;
p_start = tid + leftover_start;
p_end = (tid + 1) + leftover_start;
//for(auto &pi : particles) {
for (int i = p_start; i < p_end; i++) {
Particle& pi = particles[i];
float pressure_force[3] = {0, 0, 0};
float viscosity_force[3] = {0, 0, 0};
//Particle* candidate_neighbors = all_neighbors[i];
for (int j = 0; j < *n; j++) {
//Particle& pj = candidate_neighbors[j];
Particle& pj = particles[j];
if (&pi != &pj) {
float rij[3] = {pj.pos[0] - pi.pos[0], pj.pos[1] - pi.pos[1], pj.pos[2] - pi.pos[2]};
float r = std::sqrt(rij[0] * rij[0] + rij[1] * rij[1] + rij[2] * rij[2]);
float rij_hat[3] = {rij[0] / r, rij[1] / r, rij[2] / r};
if (r > 0 && r <= fluid->H) {
pressure_force[0] -= rij_hat[0] * fluid->MASS * (pi.pressure + pj.pressure) / (2 * pj.rho) * fluid->SPIKY_GRAD * pow(fluid->H - r, 2.f);
pressure_force[1] -= rij_hat[1] * fluid->MASS * (pi.pressure + pj.pressure) / (2 * pj.rho) * fluid->SPIKY_GRAD * pow(fluid->H - r, 2.f);
pressure_force[2] -= rij_hat[2] * fluid->MASS * (pi.pressure + pj.pressure) / (2 * pj.rho) * fluid->SPIKY_GRAD * pow(fluid->H - r, 2.f);
viscosity_force[0] += fluid->VISC * fluid->MASS * ((pj.velocity[0] - pi.velocity[0]) / pj.rho) * fluid->VISC_LAP * (fluid->H - r);
viscosity_force[1] += fluid->VISC * fluid->MASS * ((pj.velocity[1] - pi.velocity[1]) / pj.rho) * fluid->VISC_LAP * (fluid->H - r);
viscosity_force[2] += fluid->VISC * fluid->MASS * ((pj.velocity[2] - pi.velocity[2]) / pj.rho) * fluid->VISC_LAP * (fluid->H - r);
}
}
}
float G[3] = {0.f, fluid->G_STRENGTH * -9.8f, 0.f};
float gravity_force[3] = {pi.rho * G[0], pi.rho * G[1], pi.rho * G[2]};
pi.force[0] = viscosity_force[0] + pressure_force[0] + gravity_force[0];
pi.force[1] = viscosity_force[1] + pressure_force[1] + gravity_force[1];
pi.force[2] = viscosity_force[2] + pressure_force[2] + gravity_force[2];
}
}
}
__global__ void verletUpdatePosition(Particle* particles, int* n, Fluid* fluid) {
int tid = threadIdx.x;
// assuming n is a multiple of NUM_THREADS
int block_size = *n / NUM_THREADS;
int p_start = tid * block_size;
int p_end = (tid + 1) * block_size;
for (int i = p_start; i < p_end; i++) {
Particle& pi = particles[i];
pi.pos[0] += fluid->DT * pi.velocity[0];
pi.pos[1] += fluid->DT * pi.velocity[1];
pi.pos[2] += fluid->DT * pi.velocity[2];
}
if (*n % NUM_THREADS != 0 && tid < (*n - NUM_THREADS * block_size)) {
int leftover_start = (*n / NUM_THREADS) * NUM_THREADS;
p_start = tid + leftover_start;
p_end = (tid + 1) + leftover_start;
for (int i = p_start; i < p_end; i++) {
Particle& pi = particles[i];
pi.pos[0] += fluid->DT * pi.velocity[0];
pi.pos[1] += fluid->DT * pi.velocity[1];
pi.pos[2] += fluid->DT * pi.velocity[2];
}
}
}
__global__ void timeIntegrationGPU(Particle* particles, int* n, Fluid* fluid) {
int tid = threadIdx.x;
// assuming n is a multiple of NUM_THREADS
int block_size = *n / NUM_THREADS;
int p_start = tid * block_size;
int p_end = (tid + 1) * block_size;
for (int i = p_start; i < p_end; i++) {
Particle& pi = particles[i];
//pi.pos += DT * pi.velocity;
pi.velocity[0] += fluid->DT * pi.force[0] / pi.rho;
pi.velocity[1] += fluid->DT * pi.force[1] / pi.rho;
pi.velocity[2] += fluid->DT * pi.force[2] / pi.rho;
if (pi.pos[2] - fluid->EPS < -fluid->BOUND) {
pi.velocity[2] *= fluid->BOUND_DAMPING;
pi.pos[2] = -fluid->BOUND + fluid->EPS;
}
if (pi.pos[2] + fluid->EPS > fluid->BOUND) {
pi.velocity[2] *= fluid->BOUND_DAMPING;
pi.pos[2] = fluid->BOUND - fluid->EPS;
}
if (pi.pos[1] - fluid->EPS < -fluid->BOUND) {
pi.velocity[1] *= fluid->BOUND_DAMPING;
pi.pos[1] = -fluid->BOUND + fluid->EPS;
}
if (pi.pos[1] + fluid->EPS > fluid->BOUND) {
pi.velocity[1] *= fluid->BOUND_DAMPING;
pi.pos[1] = fluid->BOUND - fluid->EPS;
}
if (pi.pos[0] - fluid->EPS < -fluid->BOUND) {
pi.velocity[0] *= fluid->BOUND_DAMPING;
pi.pos[0] = -fluid->BOUND + fluid->EPS;
}
if (pi.pos[0] + fluid->EPS > fluid->BOUND) {
pi.velocity[0] *= fluid->BOUND_DAMPING;
pi.pos[0] = fluid->BOUND - fluid->EPS;
}
}
if (*n % NUM_THREADS != 0 && tid < (*n - NUM_THREADS * block_size)) {
int leftover_start = (*n / NUM_THREADS) * NUM_THREADS;
p_start = tid + leftover_start;
p_end = (tid + 1) + leftover_start;
for (int i = p_start; i < p_end; i++) {
Particle& pi = particles[i];
//pi.pos += DT * pi.velocity;
pi.velocity[0] += fluid->DT * pi.force[0] / pi.rho;
pi.velocity[1] += fluid->DT * pi.force[1] / pi.rho;
pi.velocity[2] += fluid->DT * pi.force[2] / pi.rho;
if (pi.pos[2] - fluid->EPS < -fluid->BOUND) {
pi.velocity[2] *= fluid->BOUND_DAMPING;
pi.pos[2] = -fluid->BOUND + fluid->EPS;
}
if (pi.pos[2] + fluid->EPS > fluid->BOUND) {
pi.velocity[2] *= fluid->BOUND_DAMPING;
pi.pos[2] = fluid->BOUND - fluid->EPS;
}
if (pi.pos[1] - fluid->EPS < -fluid->BOUND) {
pi.velocity[1] *= fluid->BOUND_DAMPING;
pi.pos[1] = -fluid->BOUND + fluid->EPS;
}
if (pi.pos[1] + fluid->EPS > fluid->BOUND) {
pi.velocity[1] *= fluid->BOUND_DAMPING;
pi.pos[1] = fluid->BOUND - fluid->EPS;
}
if (pi.pos[0] - fluid->EPS < -fluid->BOUND) {
pi.velocity[0] *= fluid->BOUND_DAMPING;
pi.pos[0] = -fluid->BOUND + fluid->EPS;
}
if (pi.pos[0] + fluid->EPS > fluid->BOUND) {
pi.velocity[0] *= fluid->BOUND_DAMPING;
pi.pos[0] = fluid->BOUND - fluid->EPS;
}
}
}
}
void initSPH_GPU(std::vector<Particle>& particles, Fluid* fluid) {
int n = fluid->NUM_PARTICLES;
cudaMalloc(&d_particles, n * sizeof(Particle));
cudaMalloc(&d_fluid, sizeof(Fluid));
cudaMalloc(&d_n, sizeof(int));
cudaMemcpy(d_fluid, fluid, sizeof(Fluid), cudaMemcpyHostToDevice);
cudaMemcpy(d_particles, particles.data(), n * sizeof(Particle), cudaMemcpyHostToDevice);
cudaMemcpy(d_n, &n, sizeof(int), cudaMemcpyHostToDevice);
/*
int num_cells = std::pow(fluid->CELLS_PER_DIM, 3);
cudaMallocManaged(&d_spatial_grid, num_cells * sizeof(Particle*));
cudaMallocManaged(&d_cell_counts, num_cells * sizeof(int));
for (int i = 0; i < num_cells; i++)
{
if (fluid->spatialGrid.find(i) != fluid->spatialGrid.end()) {
std::vector<Particle> cellParticles = fluid->spatialGrid[i];
Particle* array;
cudaMallocManaged(&array, cellParticles.size() * sizeof(Particle));
d_cell_counts[i] = cellParticles.size();
for(int j = 0; j < cellParticles.size(); j++) {
array[j] = cellParticles[j];
}
d_spatial_grid[i] = array;
} else {
cudaMalloc(&d_spatial_grid[i], sizeof(Particle));
d_cell_counts[i] = n;
}
}
} else {
for (int i = 0; i < num_cells; i++)
{
if (fluid->spatialGrid.find(i) != fluid->spatialGrid.end()) {
std::vector<Particle> cellParticles = fluid->spatialGrid[i];
Particle* array;
cudaMallocManaged(&array, cellParticles.size() * sizeof(Particle));
d_cell_counts[i] = cellParticles.size();
for(int j = 0; j < cellParticles.size(); j++) {
array[j] = cellParticles[j];
}
cudaFree(d_spatial_grid[i]);
d_spatial_grid[i] = array;
} else {
d_cell_counts[i] = n;
}
}*/
}
void updateSPH_GPU(std::vector<Particle>& particles, Fluid* fluid) {
int n = fluid->NUM_PARTICLES;
verletUpdatePosition<<<1, NUM_THREADS>>>(d_particles, d_n, d_fluid);
computeDensityAndPressureGPU << <1, NUM_THREADS >> > (d_particles, d_n, d_fluid);
computeForcesGPU << <1, NUM_THREADS >> > (d_particles, d_n, d_fluid);
timeIntegrationGPU << <1, NUM_THREADS >> > (d_particles, d_n, d_fluid);
cudaDeviceSynchronize();
cudaMemcpy(particles.data(), d_particles, n * sizeof(Particle), cudaMemcpyDeviceToHost);
/*
for (int i = 0; i < num_cells; i++)
{
cudaFree(&d_spatial_grid[i]);
}
*/
}
void shutdownSPH_GPU(std::vector<Particle>& particles, Fluid* fluid) {
cudaFree(d_particles);
cudaFree(d_fluid);
cudaFree(d_n);
}

13
trick_sims/SIM_fluid/models/fluid/src/sph_gpu.h Normal file
View File

@ -0,0 +1,13 @@
#include <vector>
#include <stdio.h>
#include "../include/Particle.hh"
#include "../include/Fluid.hh"
#ifndef SPH_GPU_H
#define SPH_GPU_H
void initSPH_GPU(std::vector<Particle>& particles, Fluid* fluid);
void updateSPH_GPU(std::vector<Particle>& particles, Fluid* fluid);
void shutdownSPH_GPU(std::vector<Particle>& particles, Fluid* fluid);
#endif