diff --git a/.gitignore b/.gitignore index f3219dd9..50136996 100644 --- a/.gitignore +++ b/.gitignore @@ -1,8 +1,9 @@ -Thanda.xcodeproj -Thanda.build +MFluidSolver.xcodeproj +MFluidSolver.build DerivedData CMakeScripts CMakeFiles CMakeCache.txt cmake_install.cmake CMakeLists.txt.user +build diff --git a/CMakeLists.txt b/CMakeLists.txt index bed4e57a..969c2eca 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,23 +1,43 @@ -################## -# Thanda # -################## +######################## +# MFluidSolver # +######################## -# credit - base CMake config : Yining Karl Li , edited CMake config: Akshay Shah & Debanshu Singh +# Credit - base CMake config: Yining Karl Li, edited CMake config: Akshay Shah & Debanshu Singh -#name your project -project(Thanda) +# Name your project +project(MFluidSolver) cmake_minimum_required(VERSION 2.8) -# set creates a variable -set (NUPARU "nuparu") -# add include and src directories to path +# Versioning - https://cmake.org/cmake-tutorial/ +set(MFluidSolver_VERSION_MAJOR 1) +set(MFluidSolver_VERSION_MINOR 0) +configure_file( + "${PROJECT_SOURCE_DIR}/src/MFluidSolverConfig.hpp.in" + "${PROJECT_BINARY_DIR}/MFluidSolverConfig.hpp" +) +include_directories("${PROJECT_BINARY_DIR}") + +# Copy GLSL, JSON, and Texture files +file(GLOB CONFIG_FILE "${PROJECT_SOURCE_DIR}/config/*") +file(COPY ${CONFIG_FILE} DESTINATION "${PROJECT_BINARY_DIR}/") +file(GLOB GLSL_SRC "${PROJECT_SOURCE_DIR}/glsl/*.glsl") +file(COPY ${GLSL_SRC} DESTINATION "${PROJECT_BINARY_DIR}/glsl/") +file(GLOB JSON_SRC "${PROJECT_SOURCE_DIR}/scene/*.json") +file(COPY ${JSON_SRC} DESTINATION "${PROJECT_BINARY_DIR}/scene/") +file(GLOB TEXTURE_SRC "${PROJECT_SOURCE_DIR}/texture/*") +file(COPY ${TEXTURE_SRC} DESTINATION "${PROJECT_BINARY_DIR}/texture/") + +# Set creates a variable +set(NUPARU "nuparu") +# Add include and src directories to path include_directories( + ${PROJECT_SOURCE_DIR}/src ${NUPARU}/include ${NUPARU}/src ) # Add path for pre-compiled libraries here (we will later link them with our compiled source) -# Add Nuparu library to path for OSX, linux and windows +# Add Nuparu library to path for OSX, Linux and Windows if(${CMAKE_SYSTEM_NAME} MATCHES "Darwin") set(CMAKE_LIBRARY_PATH ${CMAKE_LIBRARY_PATH} ${NUPARU}/lib/osx) elseif(${CMAKE_SYSTEM_NAME} MATCHES "Linux") @@ -26,17 +46,25 @@ elseif(WIN32) set(CMAKE_LIBRARY_PATH ${CMAKE_LIBRARY_PATH} ${NUPARU}/lib/win) endif() -#set include variables -set(GLFW_INCLUDE_DIR ${NUPARU}/include ) -set(GLEW_INCLUDE_DIR ${NUPARU}/include ) +# Set include variables set(GLFW_LIBRARY_DIR ${CMAKE_LIBRARY_PATH}) set(GLEW_LIBRARY_DIR ${CMAKE_LIBRARY_PATH}) +set(GLFW_INCLUDE_DIR ${NUPARU}/include) +set(GLEW_INCLUDE_DIR ${NUPARU}/include) # Use find_package & find_library to link with -find_package(OPENGL REQUIRED) +find_package(Boost COMPONENTS filesystem system REQUIRED) +find_package(OpenGL REQUIRED) find_package(GLEW) find_library(GLFW_LIBRARY "glfw3" HINTS ${GLFW_LIBRARY_DIR}) find_library(JSONCPP "jsoncpp") +find_library(LIBZ "z") +find_library(OPENVDB "openvdb") +find_library(OPENEXR_HALF NAMES Half PATHS ${CMAKE_LIBRARY_PATH}) +find_library(PARTIO "partio") +find_library(TBB "tbb") + +include_directories(${Boost_INCLUDE_DIRS}) add_definitions( -DTW_STATIC @@ -44,9 +72,12 @@ add_definitions( -DTW_NO_DIRECT3D -DGLEW_STATIC -D_CRT_SECURE_NO_WARNINGS + -DBOOST_TEST_DYN_LINK ) -set(CORE_LIBS ${GLFW_LIBRARY} ${GLUT_LIBRARY} ${GLEW_LIBRARY} ${JSONCPP} ${OPENGL_LIBRARY} ) +set(CORE_LIBS ${GLFW_LIBRARY} ${GLUT_LIBRARY} ${GLEW_LIBRARY} ${OPENGL_LIBRARY} + ${Boost_FILESYSTEM_LIBRARY} ${Boost_SYSTEM_LIBRARY} + ${JSONCPP} ${LIBZ} ${OPENVDB} ${OPENEXR_HALF} ${PARTIO} ${TBB}) # OSX-specific hacks/fixes if(${CMAKE_SYSTEM_NAME} MATCHES "Darwin") @@ -62,19 +93,25 @@ if(${CMAKE_SYSTEM_NAME} MATCHES "Linux") set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -lX11 -lXxf86vm -lXrandr -lpthread -lXi") endif() -# set compiler flags for c++11 +# Set compiler flags for c++11 if(${CMAKE_SYSTEM_NAME} MATCHES "Linux" OR ${CMAKE_SYSTEM_NAME} MATCHES "Darwin") set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11 -O3 -m64 -msse2 -w") elseif(WIN32) set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11") endif() +# Add -O0 to remove optimizations when using gcc +IF(CMAKE_COMPILER_IS_GNUCC) + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -O0") + set(CMAKE_C_FLAGS_DEBUG "${CMAKE_C_FLAGS_DEBUG} -O0") +ENDIF(CMAKE_COMPILER_IS_GNUCC) + if(MSVC) - set(COMPILER_FLAGS + set(COMPILER_FLAGS CMAKE_CXX_FLAGS CMAKE_CXX_FLAGS_DEBUG CMAKE_CXX_FLAGS_RELEASE - CMAKE_C_FLAGS + CMAKE_C_FLAGS CMAKE_C_FLAGS_DEBUG CMAKE_C_FLAGS_RELEASE ) @@ -82,13 +119,66 @@ endif() # Add source files you want to compile (.cpp) set(CORE_SRC - src/main.cpp - src/camera/camera.cpp - src/viewer/viewer.cpp + nuparu/src/stb_image/stb_image_write.c + src/fluidSolver/sphSolver/iiSphSolver.cpp + src/fluidSolver/sphSolver/kernelFunctions.cpp + src/fluidSolver/sphSolver/neighborSearch.cpp + src/fluidSolver/sphSolver/sphGrid.cpp + src/fluidSolver/sphSolver/sphIndexSortedUniformGrid.cpp + src/fluidSolver/sphSolver/sphParticle.cpp + src/fluidSolver/sphSolver/sphSolver.cpp + src/fluidSolver/sphSolver/sphUniformGrid.cpp + src/fluidSolver/sphSolver/sphZIndexSortedUniformGrid.cpp src/fluidSolver/fluidSolver.cpp - src/scene/scene.cpp + src/fluidSolver/particle.cpp + src/fluidSolver/zCurve.cpp + src/geom/cube.cpp + src/geom/drawable.cpp src/geom/geom.cpp + src/geom/transform.cpp + src/scene/camera.cpp + src/scene/scene.cpp + src/viewer/input.cpp + src/viewer/particleShaderProgram.cpp + src/viewer/shaderProgram.cpp + src/viewer/viewer.cpp + src/main.cpp + src/utils.cpp +) + +add_executable(MFluidSolver ${CORE_SRC}) +target_link_libraries(MFluidSolver ${CORE_LIBS}) + +##### Testing ##### +add_subdirectory(unittest) + +# KernelFunctions Library +set(MKERNELFUNCTIONS_LIB + src/fluidSolver/sphSolver/kernelFunctions.cpp +) +add_library(MKernelFunctions ${MKERNELFUNCTIONS_LIB}) + +# NeighborSearch Library +set(MNEIGHBORSEARCH_LIB + src/fluidSolver/sphSolver/neighborSearch.cpp + src/fluidSolver/sphSolver/sphGrid.cpp + src/fluidSolver/sphSolver/sphIndexSortedUniformGrid.cpp + src/fluidSolver/sphSolver/sphParticle.cpp + src/fluidSolver/sphSolver/sphUniformGrid.cpp + src/fluidSolver/sphSolver/sphZIndexSortedUniformGrid.cpp + src/fluidSolver/particle.cpp + src/fluidSolver/zCurve.cpp + src/utils.cpp +) +add_library(MNeighborSearch ${MNEIGHBORSEARCH_LIB}) +target_link_libraries(MNeighborSearch ${CORE_LIBS}) + +# ZCurve Library +set(MZCURVE_LIB + src/fluidSolver/zCurve.cpp ) +add_library(MZCurve ${MZCURVE_LIB}) -add_executable(Thanda ${CORE_SRC}) -target_link_libraries(Thanda ${CORE_LIBS}) +# Tests +enable_testing() +add_test(NAME MTest COMMAND MTest) \ No newline at end of file diff --git a/README.md b/README.md index 54283261..51309411 100644 --- a/README.md +++ b/README.md @@ -1,16 +1,151 @@ -# CIS563-FluidSolver -(Credit : CIS565 README) +# MFluidSolver # -Fluid Solver Submission guidelines: +IISPH Demo: [Vimeo](https://vimeo.com/groups/cis563/videos/164437486) +![IISPH Cube in Cube Snapshot](images/cube-iisph-sprint1_0.png?raw=true "IISPH Cube in Cube snapshot.") -- If you have modified any of the CMakeLists.txt files at all (aside from the list of CORE_SRC), you must test that your project can build. Beware of any build issues. +### How Does It Work ### +This fluid solver uses GLFW, GLM, and GLEW. It reads a JSON configuration file +and a JSON scene file, creating a fluid container and a scene container. The +fluid generates particles once, then all the particles are subject to fluid +forces (currently none). The scene is initially paused. The camera is controlled +by a standard arcball interface described below. The solver currently set to +display the neighbor particles of the first particle, which can be changed +randomly. The solver can also export the backend grid for the fluid particles to +the VDB format. The solver runs both SESPH and IISPH. -- Open a GitHub pull request so that we can see that you have finished. The title should be "Submission: YOUR NAME". +### Dependencies ### -- In the body of the pull request, include a link to your repository. +- Boost (Filesytem, System, Unit Test Framework) +- OpenVDB 2.3.0 -- Submit on canvas with a direct link to your pull request on GitHub +### Build Instructions (*nix) ### +By default, the program is set to default debugging mode as set in +`src/MFluidSolverConfig.hpp.in`. There are many other options there including +logging settings and defaults, but note that defaults will be overwritten by +the configuration file `config/config.json`. + mkdir build + cd build + cmake .. + make + ./MFluidSolver -And you're done! \ No newline at end of file +### Unit Test Instructions (*nix) ### +This will run all tests in the unittest folder. Test results will be printed +to stdout, and details are found in `build/Testing/Temporary/LastTest.log`. +This does include a brief performance test on neighbor search. + + make test + +### How To Use ### +Press P to toggle pausing. Press Press N to randomly pick a particle to show +its neighbors. Press E to export a VDB file of the graph to `export.vdb` (if +OpenVDB is enabled in `MFluidSolverConfig`). Press R to reseed the scene +using the fluid source. Press right to step frame by frame (0.001 seconds). +Press S to write out a screenshot to `screenshot_%05d.tga`. + +All SPH parameters, as well as several general program parameters, can be found +in `config/config.json`. The scene file, which contains dimensions, particle +parameters, camera settings, and particle spawn information, can be edited from +`scene/scene.json`. Parameter names are case sensitive, but options are not. If +you change parameters outside of the build folder, you will need to rerun +`cmake ..`. You can avoid this issue by editing the files in the build folder. +Note the following case-insensitive string options available for several +parameters: + +neighborSearchType +- naive +- uniform +- indexSortedUniform +- indexSortedUniformInsertionSort +- zIndexSortedUniform +- zIndexSortedUniformInsertionSort + +spawnMethod +- jittered +- poissonDisk +- uniform + +visualization +- type + - density (set maxDensityDifference) + - index + - neighbors + - none + - particle (set particleId) + - pressure (set maxPressure) + - velocity (set maxVelocity) + - velocityDir +- maxDensityDifference +- maxPressure +- maxVelocity +- particleId +- velocityColor (8-bit rgb array) + +The `numUpdates` parameter can be used to set the number of frames to simulate +in order to get accurate and comparable performance ratings. + +#### Arcball #### +The camera uses the arcball interface and can be controlled with the left and +right mouse buttons and the scroll wheel. + +The left mouse button is handled by an arcball interface. Clicking and dragging +across the center of the GL widget rotates the camera around its focus point in +the X and Y directions. Clicking and dragging in a circular motion around the +outer edges of the GL widget rotates the camera around its focus point in the Z +direction. The movements of the camera should follow the movements of the mouse +intuitively. + +The scroll wheel zooms the camera. Scrolling up zooms out and scrolling down +zooms in. + +The right mouse button pans the camera. This changes both the camera's and its +focus point's position. Right-click and drag to move the camera in the opposite +direction of the mouse. It will look like you're dragging the scene with you. + +### IISPH Simulation Performance ### +CPU: Intel i7-4800MQ limited to 2.7GHz
+GPU: NVIDIA GeForce GTX 765M
+Compiler: GCC 5.3.0
+Measured in seconds per frame. Note that this includes time spent updating the +viewer and logging input events. + +![Graph of the simulation performance for different neighbor search types](images/nsPerfGraph2.png?raw=true "Index sorting improves performance without forced insertion sort.") + +- The performance gain from index sorting is a signficant improvement over a +standard uniform grid, especially when insertion sort is not forced. +- The fact that the uniform grid uses the Y axis last (i.e.: change in the +Y-coordinate causes the greatest index change) may play a role in making uniform +grids comparable to Z-indexed grids, as fluids tend to stay within the XZ-plane. +It would be interesting to explore this further. +- Z-indexing may be more beneficial in cases where the XZ-plane is large (i.e.: +when there are larger index changes in the uniform grid), as it is in the Cube +in Cube scene. +- Index sorted grids sort every frame with std::sort, which is a hybrid of +introsort and insertion sort in [GCC](https://gcc.gnu.org/onlinedocs/libstdc++/libstdc++-html-USERS-4.4/a01027.html). +This likely explains the deficiencies of using insertion sort only. +- Insertion sort may also be less effective in IISPH where timesteps can be +larger than in traditional SESPH, so particles travel much farther and therefore +change indices much more frequently. + +### Missing Required Features ### + +- Test Partio point export +- OpenVDB level set export +- Optional: Update derivatives +- Optional: Ihmsen boundary conditions +- Optional: Cache neighbor distances +- Optional: Better seeding (poisson distribution) +- Optional: Test bicubic spline kernel +- Optional: Solid body coupling +- Optional: Headless mode +- Optional: Refactor IISPH because update() is way too hard to read +- Optional: Flyaway prevention + +### Extra Features ### +None + +### Credits ### +- [Nuparu Dependency and Foundation Libraries](https://github.com/betajippity/Nuparu) +- DDS File Loader adapted from [opengl-tutorial.org](http://www.opengl-tutorial.org/beginners-tutorials/tutorial-5-a-textured-cube/) [(WTF Public License)](http://www.wtfpl.net/). \ No newline at end of file diff --git a/config/config-sesph.json b/config/config-sesph.json new file mode 100644 index 00000000..0836f570 --- /dev/null +++ b/config/config-sesph.json @@ -0,0 +1,28 @@ +{ + "sceneJSON": "scene/scene.json", + "wireVShader": "glsl/wire.vert.glsl", + "wireFShader": "glsl/wire.frag.glsl", + "particleVShader": "glsl/particle.vert.glsl", + "particleFShader": "glsl/particle.frag.glsl", + "particleTexture": "texture/particle.dds", + "sph": { + "kStiffness": 1000, + "muViscosity": 0.028, + "mMass": 0.125, + "dRestDensity": 1000, + "dtTimestep": 0.0006, + "kernelRadius": 0.100001, + "neighborSearchType": "indexSortedUniformInsertionSort" + }, + "visualization": { + "type": "velocity", + "maxDensityDifference": 150, + "maxPressure": 40000, + "maxVelocity": 4, + "particleId": 444, + "velocityColor": [0, 160, 204] + }, + "numUpdates": -1, + "autoRender": false, + "renderSkip": 10 +} \ No newline at end of file diff --git a/config/config.json b/config/config.json new file mode 100644 index 00000000..9f7fa6e7 --- /dev/null +++ b/config/config.json @@ -0,0 +1,29 @@ +{ + "sceneJSON": "scene/scene.json", + "wireVShader": "glsl/wire.vert.glsl", + "wireFShader": "glsl/wire.frag.glsl", + "particleVShader": "glsl/particle.vert.glsl", + "particleFShader": "glsl/particle.frag.glsl", + "particleTexture": "texture/particle.dds", + "sph": { + "kStiffness": 1000, + "muViscosity": 1e-6, + "mMass": 0.125, + "dRestDensity": 1000, + "dtTimestep": 0.004, + "kernelRadius": 0.100001, + "neighborSearchType": "indexSortedUniformInsertionSort", + "maxPressureSolveIterations": 1e6 + }, + "visualization": { + "type": "pressure", + "maxDensityDifference": 150, + "maxPressure": 40000, + "maxVelocity": 4, + "particleId": 444, + "velocityColor": [0, 160, 204] + }, + "numUpdates": -1, + "autoRender": false, + "renderSkip": 10 +} \ No newline at end of file diff --git a/glsl/particle.frag.glsl b/glsl/particle.frag.glsl new file mode 100644 index 00000000..310dbe69 --- /dev/null +++ b/glsl/particle.frag.glsl @@ -0,0 +1,26 @@ +#version 330 core + +in vec3 afs_Color; + +out vec4 out_Color; + +void main() +{ + out_Color = vec4(afs_Color, 1.f); +} + +/* Use points until I figure out how to draw lines as a triangle strip +uniform sampler2D u_BillboardTextureSampler; + +in vec3 afs_Color; +in vec2 afs_UV; + +out vec4 out_Color; + +void main() +{ + vec4 tex_Color = texture(u_BillboardTextureSampler, afs_UV); + vec4 in_Color = vec4(afs_Color, 1.f); + out_Color = tex_Color * in_Color; +} +*/ \ No newline at end of file diff --git a/glsl/particle.vert.glsl b/glsl/particle.vert.glsl new file mode 100644 index 00000000..78e0a8ee --- /dev/null +++ b/glsl/particle.vert.glsl @@ -0,0 +1,42 @@ +#version 330 core + +uniform mat4 u_Model; +uniform mat4 u_ViewProjection; + +in vec3 avs_Color; +in vec3 avs_Position; + +out vec3 afs_Color; + +void main() +{ + afs_Color = avs_Color; + gl_Position = u_ViewProjection * vec4(avs_Position, 1.f); +} + +/* Use points until I figure out how to draw lines as a triangle strip +uniform mat4 u_Model; +uniform mat4 u_ViewProjection; +uniform vec3 u_CameraRight; +uniform vec3 u_CameraUp; +uniform float u_ParticleSize; + +in vec3 avs_Color; +in vec3 avs_Position; +in vec3 avs_Billboard; + +out vec3 afs_Color; +out vec2 afs_UV; + +void main() +{ + afs_Color = avs_Color; + + vec3 modelPosition = avs_Position + + u_CameraRight * avs_Billboard.x * u_ParticleSize + + u_CameraUp * avs_Billboard.y * u_ParticleSize; + gl_Position = u_ViewProjection * vec4(modelPosition, 1.f); + + afs_UV = avs_Billboard.xy + vec2(0.5, 0.5); +} +*/ diff --git a/glsl/wire.frag.glsl b/glsl/wire.frag.glsl new file mode 100644 index 00000000..3243517f --- /dev/null +++ b/glsl/wire.frag.glsl @@ -0,0 +1,10 @@ +#version 330 core + +in vec3 afs_Color; + +out vec4 out_Color; + +void main() +{ + out_Color = vec4(afs_Color, 1.f); +} diff --git a/glsl/wire.vert.glsl b/glsl/wire.vert.glsl new file mode 100644 index 00000000..9dacb2f2 --- /dev/null +++ b/glsl/wire.vert.glsl @@ -0,0 +1,16 @@ +#version 330 core + +uniform mat4 u_Model; +uniform mat4 u_ViewProjection; + +in vec3 avs_Color; +in vec3 avs_Position; + +out vec3 afs_Color; + +void main() +{ + afs_Color = avs_Color; + vec4 modelPosition = u_Model * vec4(avs_Position, 1.f); + gl_Position = u_ViewProjection * modelPosition; +} diff --git a/images/cube-iisph-sprint1_0.png b/images/cube-iisph-sprint1_0.png new file mode 100644 index 00000000..d10e65d0 Binary files /dev/null and b/images/cube-iisph-sprint1_0.png differ diff --git a/images/nsPerfGraph0.png b/images/nsPerfGraph0.png new file mode 100644 index 00000000..5d3ab7f2 Binary files /dev/null and b/images/nsPerfGraph0.png differ diff --git a/images/nsPerfGraph1.png b/images/nsPerfGraph1.png new file mode 100644 index 00000000..37177f36 Binary files /dev/null and b/images/nsPerfGraph1.png differ diff --git a/images/nsPerfGraph2.png b/images/nsPerfGraph2.png new file mode 100644 index 00000000..062932c0 Binary files /dev/null and b/images/nsPerfGraph2.png differ diff --git a/nuparu/bin/linux/makecircle b/nuparu/bin/linux/makecircle new file mode 100644 index 00000000..c2e31c8d Binary files /dev/null and b/nuparu/bin/linux/makecircle differ diff --git a/nuparu/bin/linux/makeline b/nuparu/bin/linux/makeline new file mode 100644 index 00000000..1bc176da Binary files /dev/null and b/nuparu/bin/linux/makeline differ diff --git a/nuparu/bin/linux/partattr b/nuparu/bin/linux/partattr new file mode 100644 index 00000000..2aef7e37 Binary files /dev/null and b/nuparu/bin/linux/partattr differ diff --git a/nuparu/bin/linux/partconv b/nuparu/bin/linux/partconv new file mode 100644 index 00000000..a03230aa Binary files /dev/null and b/nuparu/bin/linux/partconv differ diff --git a/nuparu/bin/linux/partinfo b/nuparu/bin/linux/partinfo new file mode 100644 index 00000000..86c27058 Binary files /dev/null and b/nuparu/bin/linux/partinfo differ diff --git a/nuparu/bin/linux/partview b/nuparu/bin/linux/partview new file mode 100644 index 00000000..91084d0b Binary files /dev/null and b/nuparu/bin/linux/partview differ diff --git a/nuparu/bin/linux/test b/nuparu/bin/linux/test new file mode 100644 index 00000000..e8a6439b Binary files /dev/null and b/nuparu/bin/linux/test differ diff --git a/nuparu/bin/linux/testcache b/nuparu/bin/linux/testcache new file mode 100644 index 00000000..9b628721 Binary files /dev/null and b/nuparu/bin/linux/testcache differ diff --git a/nuparu/bin/linux/testiterator b/nuparu/bin/linux/testiterator new file mode 100644 index 00000000..c8d2a792 Binary files /dev/null and b/nuparu/bin/linux/testiterator differ diff --git a/nuparu/bin/linux/testkdtree b/nuparu/bin/linux/testkdtree new file mode 100644 index 00000000..42424c7c Binary files /dev/null and b/nuparu/bin/linux/testkdtree differ diff --git a/nuparu/bin/linux/teststr b/nuparu/bin/linux/teststr new file mode 100644 index 00000000..a749a42f Binary files /dev/null and b/nuparu/bin/linux/teststr differ diff --git a/nuparu/lib/linux/libopenvdb.so b/nuparu/lib/linux/libopenvdb.so new file mode 120000 index 00000000..d72ac950 Binary files /dev/null and b/nuparu/lib/linux/libopenvdb.so differ diff --git a/nuparu/lib/linux/libopenvdb.so.2.3 b/nuparu/lib/linux/libopenvdb.so.2.3 new file mode 120000 index 00000000..d72ac950 Binary files /dev/null and b/nuparu/lib/linux/libopenvdb.so.2.3 differ diff --git a/nuparu/lib/linux/libopenvdb.so.2.3.0 b/nuparu/lib/linux/libopenvdb.so.2.3.0 new file mode 100644 index 00000000..7d079afa Binary files /dev/null and b/nuparu/lib/linux/libopenvdb.so.2.3.0 differ diff --git a/nuparu/lib/linux/libpartio.a b/nuparu/lib/linux/libpartio.a new file mode 100644 index 00000000..2f1c5c25 Binary files /dev/null and b/nuparu/lib/linux/libpartio.a differ diff --git a/scene/scene-new.json b/scene/scene-new.json new file mode 100644 index 00000000..ee176cc8 --- /dev/null +++ b/scene/scene-new.json @@ -0,0 +1,26 @@ +{ + "containerDim": { + "scaleX": 2.0, + "scaleY": 2.0, + "scaleZ": 0.5 + }, + "particleDim": { + "scaleX": 1.0, + "scaleY": 1.5, + "scaleZ": 0.5, + "posX": -0.5, + "posY": 0.25, + "posZ": 0.0 + }, + "camera": { + "eyeX": 0.0, + "eyeY": 0.0, + "eyeZ": 3.5, + "refX": 0.0, + "refY": 0.0, + "refZ": 0.0 + }, + "particleSeparation": 0.05, + "maxParticles": 10000, + "spawnMethod": "jittered" +} \ No newline at end of file diff --git a/scene/scene.json b/scene/scene.json new file mode 100644 index 00000000..bf666c07 --- /dev/null +++ b/scene/scene.json @@ -0,0 +1,23 @@ +{ + "containerDim": { + "scaleX": 2.0, + "scaleY": 2.0, + "scaleZ": 2.0 + }, + "particleDim": { + "scaleX": 1.0, + "scaleY": 1.0, + "scaleZ": 1.0 + }, + "camera": { + "eyeX": 0.0, + "eyeY": 0.0, + "eyeZ": 3.5, + "refX": 0.0, + "refY": 0.0, + "refZ": 0.0 + }, + "particleSeparation": 0.05, + "maxParticles": 10000, + "spawnMethod": "jittered" +} \ No newline at end of file diff --git a/src/MFluidSolverConfig.hpp.in b/src/MFluidSolverConfig.hpp.in new file mode 100644 index 00000000..21eb0a3e --- /dev/null +++ b/src/MFluidSolverConfig.hpp.in @@ -0,0 +1,106 @@ +// Copyright 2016 Robert Zhou +// +// MFluidSolverConfig.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_CONFIG_HPP_ +#define MFLUIDSOLVER_CONFIG_HPP_ + +// Versioning +#define MFluidSolver_VERSION_MAJOR @MFluidSolver_VERSION_MAJOR@ +#define MFluidSolver_VERSION_MINOR @MFluidSolver_VERSION_MINOR@ + +// Logging/Debug +#define MFluidSolver_LOG_LEVEL 2 +#define MFluidSolver_OPENGL_DEBUG 1 +#define MFluidSolver_PARTICLE_STATS 0 +#define MFluidSolver_PARTICLE_STATS_FILES 0 +#define MFluidSolver_USE_ASSERTS 0 + +#define MFluidSolver_RECORD_PERFORMANCE 1 + +// Libraries +#define MFluidSolver_USE_OPENVDB 1 +#define MFluidSolver_USE_PARTIO 1 +#define MFluidSolver_USE_TBB 1 + +// Misc +#define MFluidSolver_ZCURVE_CACHING 1 + +// Defaults +#define MFluidSolver_DEFAULT_CONFIG_FILE "config.json" + +#define MFluidSolver_DEFAULT_SCENE_FILE "scene/scene.json" +#define MFluidSolver_DEFAULT_WIRE_VERT_FILE "glsl/wire.vert.glsl" +#define MFluidSolver_DEFAULT_WIRE_FRAG_FILE "glsl/wire.frag.glsl" +#define MFluidSolver_DEFAULT_PARTICLE_VERT_FILE "glsl/particle.vert.glsl" +#define MFluidSolver_DEFAULT_PARTICLE_FRAG_FILE "glsl/particle.frag.glsl" +#define MFluidSolver_DEFAULT_PARTICLE_TEX_FILE "texture/particle.dds" + +#define MFluidSolver_DEFAULT_AUTORENDER false +#define MFluidSolver_DEFAULT_BG_COLOR_R 0.3f +#define MFluidSolver_DEFAULT_BG_COLOR_G 0.3f +#define MFluidSolver_DEFAULT_BG_COLOR_B 0.3f +#define MFluidSolver_DEFAULT_GEOMETRY_COLOR_R 1 +#define MFluidSolver_DEFAULT_GEOMETRY_COLOR_G 1 +#define MFluidSolver_DEFAULT_GEOMETRY_COLOR_B 1 +#define MFluidSolver_DEFAULT_GRAVITY -9.8 +#define MFluidSolver_DEFAULT_LIMIT_UPDATES false +#define MFluidSolver_DEFAULT_MAX_PARTICLES 10000 +#define MFluidSolver_DEFAULT_MAX_UPDATES 1000 +#define MFluidSolver_DEFAULT_PARTICLE_BOUNCE 0.1f +#define MFluidSolver_DEFAULT_PARTICLE_COLOR_R 0 +#define MFluidSolver_DEFAULT_PARTICLE_COLOR_G 0 +#define MFluidSolver_DEFAULT_PARTICLE_COLOR_B 1 +#define MFluidSolver_DEFAULT_PARTICLE_SEPARATION 0.1f +#define MFluidSolver_DEFAULT_POINT_SIZE 5 +#define MFluidSolver_DEFAULT_SCENE_CAMERA_EYEX 0.f +#define MFluidSolver_DEFAULT_SCENE_CAMERA_EYEY 0.f +#define MFluidSolver_DEFAULT_SCENE_CAMERA_EYEZ 1.f +#define MFluidSolver_DEFAULT_SCENE_CAMERA_REFX 0.f +#define MFluidSolver_DEFAULT_SCENE_CAMERA_REFY 0.f +#define MFluidSolver_DEFAULT_SCENE_CAMERA_REFZ 0.f +#define MFluidSolver_DEFAULT_SCENE_CONTAINER_SCALEX 1.f +#define MFluidSolver_DEFAULT_SCENE_CONTAINER_SCALEY 1.f +#define MFluidSolver_DEFAULT_SCENE_CONTAINER_SCALEZ 1.f +#define MFluidSolver_DEFAULT_SCENE_SOURCE_SCALEX 0.5f +#define MFluidSolver_DEFAULT_SCENE_SOURCE_SCALEY 0.5f +#define MFluidSolver_DEFAULT_SCENE_SOURCE_SCALEZ 0.5f +#define MFluidSolver_DEFAULT_SCENE_SOURCE_POSX 0.f +#define MFluidSolver_DEFAULT_SCENE_SOURCE_POSY 0.f +#define MFluidSolver_DEFAULT_SCENE_SOURCE_POSZ 0.f +#define MFluidSolver_DEFAULT_SEARCH_RADIUS 0.1f +#define MFluidSolver_DEFAULT_SPAWNMETHOD ParticleSpawnMethod::Jittered +#define MFluidSolver_DEFAULT_SPAWNMETHODSTRING "jittered" +#define MFluidSolver_DEFAULT_SPH_KERNELRADIUS 0.100001 +#define MFluidSolver_DEFAULT_SPH_MASS 0.125 +#define MFluidSolver_DEFAULT_SPH_NEIGHBORSEARCHTYPE_STRING "indexSortedUniformInsertionSort" +#define MFluidSolver_DEFAULT_SPH_NEIGHBORSEARCHTYPE NeighborSearchType::IndexSortedUniformGridWithInsertion +#define MFluidSolver_DEFAULT_SPH_NEIGHBORSEARCH_USE_ZCURVE true +#define MFluidSolver_DEFAULT_SPH_RESTDENSITY 1000 +#define MFluidSolver_DEFAULT_SPH_STIFFNESS 1000 +#define MFluidSolver_DEFAULT_SPH_TIMESTEP 0.001 +#define MFluidSolver_DEFAULT_SPH_VISCOSITY 1 +#define MFluidSolver_DEFAULT_RENDERSKIP 0 +#define MFluidSolver_DEFAULT_UPDATE_STEP 0.001 +#define MFluidSolver_DEFAULT_VISUALIZATION FluidVisualizationType::Velocity +#define MFluidSolver_DEFAULT_VISUALIZATION_STRING "velocity" +#define MFluidSolver_DEFAULT_VISUALIZATION_MAXDENSITYDIFFERENCE 200 +#define MFluidSolver_DEFAULT_VISUALIZATION_MAXPRESSURE 40000 +#define MFluidSolver_DEFAULT_VISUALIZATION_MAXVELOCITY 4 +#define MFluidSolver_DEFAULT_VISUALIZATION_VELOCITYCOLOR_R 0 +#define MFluidSolver_DEFAULT_VISUALIZATION_VELOCITYCOLOR_G 160 +#define MFluidSolver_DEFAULT_VISUALIZATION_VELOCITYCOLOR_B 204 +#define MFluidSolver_DEFAULT_WINDOW_WIDTH 1024 +#define MFluidSolver_DEFAULT_WINDOW_HEIGHT 768 + +// Do not touch +#define MFluidSolver_LOG_NONE 1000 +#define MFluidSolver_LOG_FATAL 5 +#define MFluidSolver_LOG_ERROR 4 +#define MFluidSolver_LOG_WARN 3 +#define MFluidSolver_LOG_INFO 2 +#define MFluidSolver_LOG_DEBUG 1 +#define MFluidSolver_LOG_TRACE 0 + +#endif // MFLUIDSOLVER_CONFIG_HPP_ diff --git a/src/camera/camera.cpp b/src/camera/camera.cpp deleted file mode 100644 index fb246b12..00000000 --- a/src/camera/camera.cpp +++ /dev/null @@ -1,6 +0,0 @@ -// -// camera.cpp -// Thanda -// - -#include "camera.hpp" diff --git a/src/camera/camera.hpp b/src/camera/camera.hpp deleted file mode 100644 index 59ad2b12..00000000 --- a/src/camera/camera.hpp +++ /dev/null @@ -1,10 +0,0 @@ -// -// camera.hpp -// Thanda -// - -#ifndef camera_hpp -#define camera_hpp - - -#endif /* camera_hpp */ diff --git a/src/fluidSolver/fluidSolver.cpp b/src/fluidSolver/fluidSolver.cpp index 9c9a1663..1eb60526 100644 --- a/src/fluidSolver/fluidSolver.cpp +++ b/src/fluidSolver/fluidSolver.cpp @@ -1,6 +1,93 @@ +// Copyright 2016 Robert Zhou // // fluidSolver.cpp -// Thanda - +// MFluidSolver #include "fluidSolver.hpp" + +#include + +#include "utils.hpp" + +FluidSolver::FluidSolver() + : _maxParticles(MFluidSolver_DEFAULT_MAX_PARTICLES), + _gravity(MFluidSolver_DEFAULT_GRAVITY), + _particleSeparation(MFluidSolver_DEFAULT_PARTICLE_SEPARATION), + visualizationType(MFluidSolver_DEFAULT_VISUALIZATION), + visualizationMaxPressure(MFluidSolver_DEFAULT_VISUALIZATION_MAXPRESSURE), + visualizationMaxVelocity(MFluidSolver_DEFAULT_VISUALIZATION_MAXVELOCITY), + visualizationVelocityColor( + static_cast( + MFluidSolver_DEFAULT_VISUALIZATION_VELOCITYCOLOR_R / 255.f), + static_cast( + MFluidSolver_DEFAULT_VISUALIZATION_VELOCITYCOLOR_G / 255.f), + static_cast( + MFluidSolver_DEFAULT_VISUALIZATION_VELOCITYCOLOR_B / 255.f)), + _fixedTimestep(MFluidSolver_DEFAULT_UPDATE_STEP), + computeTime(0), numUpdates(0), firstRun(true), endedSimulation(false), + maxUpdates(MFluidSolver_DEFAULT_MAX_UPDATES), + limitNumUpdates(MFluidSolver_DEFAULT_LIMIT_UPDATES), + fluidSource(nullptr), fluidContainer(nullptr) { +} + +FluidSolver::~FluidSolver() { +} + +// Solver +void FluidSolver::update(double deltaT) { +} +void FluidSolver::updateStep() { + update(_fixedTimestep); +} + +// Getters +float FluidSolver::gravity() const { + return _gravity; +} +int FluidSolver::maxParticles() const { + return _maxParticles; +} +float FluidSolver::particleSeparation() const { + return _particleSeparation; +} + +// Setters +void FluidSolver::setGravity(float g) { + _gravity = g; +} +void FluidSolver::setMaxParticles(int mp) { + _maxParticles = mp; +} +void FluidSolver::setParticleSeparation(float ps) { + _particleSeparation = ps; +} +void FluidSolver::setFixedTimestep(float ft) { + _fixedTimestep = ft; +} + +// Simulation End +void FluidSolver::endSimulation() { + if (!endedSimulation) { + endTime = tbb::tick_count::now(); + endedSimulation = true; + } +} + +// Misc +void FluidSolver::printPerformanceStats() { + if (numUpdates > 0) { + computeTime = (endTime - startTime).seconds(); + std::cout << "PERF: Overall simulation ran for " << + MUtils::toHMS(computeTime) << " (" << computeTime << " seconds) over " + << numUpdates << " frames" << std::endl; + std::cout << "PERF: Average is " << + (computeTime / numUpdates) << " seconds per frame" << std::endl; + std::cout << "PERF: Simulated " << + (numUpdates * _fixedTimestep) << " seconds of fluid simulation" << + std::endl; + } +} + +unsigned int FluidSolver::updateNumber() { + return numUpdates; +} diff --git a/src/fluidSolver/fluidSolver.hpp b/src/fluidSolver/fluidSolver.hpp index 6429c4e3..42cb3fbe 100644 --- a/src/fluidSolver/fluidSolver.hpp +++ b/src/fluidSolver/fluidSolver.hpp @@ -1,9 +1,87 @@ +// Copyright 2016 Robert Zhou // // fluidSolver.hpp -// Thanda +// MFluidSolver -#ifndef fluidSolver_hpp -#define fluidSolver_hpp +#ifndef MFLUIDSOLVER_FLUIDSOLVER_HPP_ +#define MFLUIDSOLVER_FLUIDSOLVER_HPP_ +#include -#endif /* fluidSolver_hpp */ +#include "MFluidSolverConfig.hpp" + +#include "geom/geom.hpp" +#include "particle.hpp" + +enum FluidVisualizationType { + Density, Index, Neighbors, None, Particle, Pressure, Velocity, VelocityDir +}; + +class FluidSolver { + public: + FluidSolver(); + ~FluidSolver(); + + // Solver + virtual void update(double deltaT); + void updateStep(); + + // Particles + virtual void addParticleAt(const glm::vec3 &position) = 0; + virtual unsigned int numParticles() const = 0; + + // Getters + float gravity() const; + int maxParticles() const; + float particleSeparation() const; + + // Setters + void setFixedTimestep(float ft); + void setGravity(float g); + virtual void setMaxParticles(int mp); + virtual void setParticleSeparation(float ps); + + // Simulation End + inline bool checkIfEnded(); + inline bool hasEndedSimulation(); + void endSimulation(); + + // Misc + virtual void printPerformanceStats(); + virtual void sceneLoaded() = 0; + unsigned int updateNumber(); + + // Containers + Geometry *fluidSource; + Geometry *fluidContainer; + + protected: + // Helpers + inline void logTimestep(); + + // Parameters + float _gravity; + float _particleSeparation; + int _maxParticles; + float _fixedTimestep; + + // Visualization + FluidVisualizationType visualizationType; + float visualizationMaxDensityDifference; + float visualizationMaxPressure; + float visualizationMaxVelocity; + glm::vec3 visualizationVelocityColor; + unsigned int targetParticle; + + // Performance measuring + double computeTime; + tbb::tick_count startTime, endTime; + unsigned int numUpdates; + unsigned int maxUpdates; + bool limitNumUpdates; + bool firstRun, endedSimulation; +}; + +#include "fluidSolver/fluidSolver.inline.hpp" + +#endif // MFLUIDSOLVER_FLUIDSOLVER_HPP_ diff --git a/src/fluidSolver/fluidSolver.inline.hpp b/src/fluidSolver/fluidSolver.inline.hpp new file mode 100644 index 00000000..1266cf65 --- /dev/null +++ b/src/fluidSolver/fluidSolver.inline.hpp @@ -0,0 +1,39 @@ +// Copyright 2016 Robert Zhou +// +// fluidSolver.inline.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_FLUIDSOLVER_INLINE_HPP_ +#define MFLUIDSOLVER_FLUIDSOLVER_INLINE_HPP_ + +// Simulation End +inline bool FluidSolver::checkIfEnded() { + // Record times and updates + if (firstRun) { + startTime = tbb::tick_count::now(); + firstRun = false; + } + if (endedSimulation) { + return true; + } else if (limitNumUpdates && numUpdates >= maxUpdates) { + endTime = tbb::tick_count::now(); + endedSimulation = true; + return true; + } else { + ++numUpdates; + } + return false; +} + +inline bool FluidSolver::hasEndedSimulation() { + return endedSimulation; +} + +// Helpers +inline void FluidSolver::logTimestep() { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_DEBUG + std::cout << "DEBUG: Updating by " << deltaT << " seconds" << std::endl; + #endif +} + +#endif // MFLUIDSOLVER_FLUIDSOLVER_INLINE_HPP_ diff --git a/src/fluidSolver/particle.cpp b/src/fluidSolver/particle.cpp new file mode 100644 index 00000000..aa3933c8 --- /dev/null +++ b/src/fluidSolver/particle.cpp @@ -0,0 +1,31 @@ +// Copyright 2016 Robert Zhou +// +// particle.cpp +// MFluidSolver + +#include "particle.hpp" + +#if MFluidSolver_PARTICLE_STATS +unsigned int Particle::lastParticleID = 1; +#endif + +// Constructors / Destructors +Particle::Particle(float mass, const glm::vec3 &position) + : color( + MFluidSolver_DEFAULT_PARTICLE_COLOR_R, + MFluidSolver_DEFAULT_PARTICLE_COLOR_G, + MFluidSolver_DEFAULT_PARTICLE_COLOR_B), + _mass(mass), _position(position), _velocity(glm::vec3(0)) { + #if MFluidSolver_PARTICLE_STATS + // Note assumes atomic access to lastParticleID + ID = lastParticleID; + ++lastParticleID; + #endif +} + +// Operators +bool Particle::operator==(const Particle &b) const { + return b.mass() == _mass && + b.position() == _position && + b.velocity() == _velocity; +} diff --git a/src/fluidSolver/particle.hpp b/src/fluidSolver/particle.hpp new file mode 100644 index 00000000..f1816c76 --- /dev/null +++ b/src/fluidSolver/particle.hpp @@ -0,0 +1,48 @@ +// Copyright 2016 Robert Zhou +// +// particle.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_FLUIDSOLVER_PARTICLE_HPP_ +#define MFLUIDSOLVER_FLUIDSOLVER_PARTICLE_HPP_ + +#include + +#include "MFluidSolverConfig.hpp" + +class Particle { + public: + // Constructors / Destructors + Particle() : Particle(1, glm::vec3(0)) {} + explicit Particle(const glm::vec3 &position) : Particle(1, position) {} + Particle(float mass, const glm::vec3 &position); + + // Operators + virtual bool operator==(const Particle &b) const; + + // Update + inline void setPosition(const glm::vec3 &position); + inline void reverseVelocity(const glm::ivec3 &directions, + float bounceCoefficient = MFluidSolver_DEFAULT_PARTICLE_BOUNCE); + inline void stopVelocity(const glm::ivec3 &directions); + + // Properties + glm::vec3 color; + inline float mass() const; + inline glm::vec3 position() const; + inline glm::vec3 velocity() const; + + #if MFluidSolver_PARTICLE_STATS + unsigned int ID; + static unsigned int lastParticleID; + #endif + + protected: + float _mass; + glm::vec3 _position; + glm::vec3 _velocity; +}; + +#include "fluidSolver/particle.inline.hpp" + +#endif // MFLUIDSOLVER_FLUIDSOLVER_PARTICLE_HPP_ diff --git a/src/fluidSolver/particle.inline.hpp b/src/fluidSolver/particle.inline.hpp new file mode 100644 index 00000000..1d24ab99 --- /dev/null +++ b/src/fluidSolver/particle.inline.hpp @@ -0,0 +1,40 @@ +// Copyright 2016 Robert Zhou +// +// particle.inline.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_FLUIDSOLVER_PARTICLE_INLINE_HPP_ +#define MFLUIDSOLVER_FLUIDSOLVER_PARTICLE_INLINE_HPP_ + +// Update +inline void Particle::setPosition(const glm::vec3 &position) { + _position = position; +} + +inline void Particle::reverseVelocity( + const glm::ivec3 &directions, float bounceCoefficient) { + if (directions.x != 0) _velocity.x *= -1 * bounceCoefficient; + if (directions.y != 0) _velocity.y *= -1 * bounceCoefficient; + if (directions.z != 0) _velocity.z *= -1 * bounceCoefficient; +} + +inline void Particle::stopVelocity(const glm::ivec3 &directions) { + if (directions.x != 0) _velocity.x = 0; + if (directions.y != 0) _velocity.y = 0; + if (directions.z != 0) _velocity.z = 0; +} + +// Properties +inline float Particle::mass() const { + return _mass; +} + +inline glm::vec3 Particle::position() const { + return _position; +} + +inline glm::vec3 Particle::velocity() const { + return _velocity; +} + +#endif // MFLUIDSOLVER_FLUIDSOLVER_PARTICLE_INLINE_HPP_ diff --git a/src/fluidSolver/sphSolver/iiSphSolver.cpp b/src/fluidSolver/sphSolver/iiSphSolver.cpp new file mode 100644 index 00000000..0fd88374 --- /dev/null +++ b/src/fluidSolver/sphSolver/iiSphSolver.cpp @@ -0,0 +1,489 @@ +// Copyright 2016 Robert Zhou +// +// iiSphSolver.cpp +// MFluidSolver + +#include "iiSphSolver.hpp" + +#include +#include +#include +#include + +#if MFluidSolver_PARTICLE_STATS_FILES +#include +#include +#endif + +#include + +#include "utils.hpp" + +IISPHSolver::IISPHSolver() + : maxIterations(1e6) { +} + +// TODO: Clean this up +void IISPHSolver::loadConfig(const std::string &file) { + if (checkInited()) return; + SPHSolver::loadConfig(file); + + // Parse JSON file + Json::Reader reader; + Json::Value root; + std::ifstream sceneStream(file, std::ifstream::binary); + bool success = reader.parse(sceneStream, root, false); + if (!success) { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_ERROR + std::cerr << "ERROR: Failed to parse config file " << file << std::endl; + #endif + return; + } + + maxIterations = root["sph"].get( + "maxPressureSolveIterations", 1e6).asInt(); + + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + std::cout << "INFO: - Max Pressure Solve Iterations: " << + maxIterations << std::endl; + #endif +} + +void IISPHSolver::update(double deltaT) { + if (checkIfEnded()) return; + + // NOTE: TIMESTEP IS OVERWRITTEN HERE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + deltaT = _fixedTimestep; + const float deltaTF = static_cast(deltaT); + const float deltaTF2 = static_cast(deltaT * deltaT); + + logTimestep(); + prepNeighborSearch(); + runNeighborSearch(); + + #if MFluidSolver_PARTICLE_STATS + averagePosition = glm::vec3(0); + averageVelocityMagnitude = 0; + averageDensity = 0; + averagePressure = 0; + averageNonPressureForceMagnitude = 0; + averagePressureForceMagnitude = 0; + averageAdvectionDiagonal = 0; + averageAdvectionDisplacementMagnitude = 0; + averageDensityIntermediate = 0; + averageSumPressureDisplacementFromNeighborsMagnitude = 0; + averageVelocityIntermediateMagnitude = 0; + averageNumNeighbors = 0; + + #if MFluidSolver_PARTICLE_STATS_FILES + // Output file names + std::string frameString = MUtils::zeroPad(numUpdates, 6); + std::ostringstream pressureForceMagnitudeSS; + pressureForceMagnitudeSS << "particlestats/pfm-" << frameString << ".txt"; + std::ostringstream pressureSS; + pressureSS << "particlestats/psr-" << frameString << ".txt"; + std::ostringstream velocityMagnitudeSS; + velocityMagnitudeSS << "particlestats/vel-" << frameString << ".txt"; + + // Output files + std::ofstream pressureForceMagnitudeOutFile; + pressureForceMagnitudeOutFile.open( + pressureForceMagnitudeSS.str(), + std::ofstream::out | std::ofstream::trunc); + std::ofstream pressureOutFile; + pressureOutFile.open( + pressureSS.str(), + std::ofstream::out | std::ofstream::trunc); + std::ofstream velocityMagnitudeOutFile; + velocityMagnitudeOutFile.open( + velocityMagnitudeSS.str(), + std::ofstream::out | std::ofstream::trunc); + + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_ERROR + if (!pressureForceMagnitudeOutFile.is_open()) { + std::cerr << "ERROR: Could not open file: " << + pressureForceMagnitudeSS.str() << std::endl; + } + if (!pressureOutFile.is_open()) { + std::cerr << "ERROR: Could not open file: " << + pressureSS.str() << std::endl; + } + if (!velocityMagnitudeOutFile.is_open()) { + std::cerr << "ERROR: Could not open file: " << + velocityMagnitudeSS.str() << std::endl; + } + #endif + #endif + #endif + + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_WARN + numFlyaways = 0; + #endif + + // Procedure: Predict Advection + iter_all_sphparticles_start + calculateDensity(&p); + calculateNonPressureForce(&p); + + #if MFluidSolver_PARTICLE_STATS + // Check if flyaway is affecting other particles + if (p.flyaway && !p.neighbors()->empty()) { + for (SPHParticle *q : *(p.neighbors())) { + std::cout << "STAT: Particle ID " << q->ID << + " is now neighbors with flyaway ID " << p.ID << std::endl; + } + } + #endif + + // Flyaway particle detection (no neighbors, too much free will) + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_WARN || MFluidSolver_PARTICLE_STATS + if (p.neighbors()->empty()) { + #if MFluidSolver_PARTICLE_STATS + std::lock_guard guard(statsCoutMutex); + p.flyaway = true; + const glm::vec3 position = p.position(); + const glm::vec3 velocity = p.velocity(); + const glm::vec3 nonPressureForce = p.nonPressureForce(); + const glm::vec3 pressureForce = p.pressureForce(); + const glm::vec3 advectionDisplacementEstimate = + p.advectionDisplacementEstimate(); + const glm::vec3 sumPressureDisplacementFromNeighbors = + p.sumPressureDisplacementFromNeighbors(); + const glm::vec3 velocityIntermediate = p.velocityIntermediate(); + + std::cout << "WARN: Flyaway particle ID " << + p.ID << " detected! Printing information:" << + std::endl + << "WARN: - Position: (" << + position.x << ", " << + position.y << ", " << + position.z << ") with magnitude " << + glm::length(position) << std::endl + << "WARN: - Velocity: (" << + velocity.x << ", " << + velocity.y << ", " << + velocity.z << ") with magnitude " << + glm::length(velocity) << std::endl + << "WARN: - Density: " << p.density() << std::endl + << "WARN: - Pressure: " << p.pressure() << std::endl + << "WARN: - Non-pressure Force: (" << + nonPressureForce.x << ", " << + nonPressureForce.y << ", " << + nonPressureForce.z << ") with magnitude " << + glm::length(nonPressureForce) << std::endl + << "WARN: - Pressure Force: (" << + pressureForce.x << ", " << + pressureForce.y << ", " << + pressureForce.z << ") with magnitude " << + glm::length(pressureForce) << std::endl + << "WARN: - Advection Diagonal: " << + p.advectionDiagonal() << std::endl + << "WARN: - Advection Displacement Estimate: (" << + advectionDisplacementEstimate.x << ", " << + advectionDisplacementEstimate.y << ", " << + advectionDisplacementEstimate.z << + ") with magnitude " << + glm::length(advectionDisplacementEstimate) << + std::endl + << "WARN: - Density Intermediate: " << + p.densityIntermediate() << std::endl + << "WARN: - Sum Pressure Displacement From Neighbors: (" << + sumPressureDisplacementFromNeighbors.x << ", " << + sumPressureDisplacementFromNeighbors.y << ", " << + sumPressureDisplacementFromNeighbors.z << + ") with magnitude " << + glm::length(sumPressureDisplacementFromNeighbors) << + std::endl + << "WARN: - Velocity Intermediate: (" << + velocityIntermediate.x << ", " << + velocityIntermediate.y << ", " << + velocityIntermediate.z << ") with magnitude " << + glm::length(velocityIntermediate) << std::endl; + #endif + ++numFlyaways; + } + #endif + + // Particle statistics (TODO: Save to file for complex statistics) + #if MFluidSolver_PARTICLE_STATS + averagePosition += p.position(); + averageVelocityMagnitude += glm::length(p.velocity()); + averageDensity += p.density(); + averagePressure += p.pressure(); + averageNonPressureForceMagnitude += glm::length(p.nonPressureForce()); + averagePressureForceMagnitude += glm::length(p.pressureForce()); + averageAdvectionDiagonal += p.advectionDiagonal(); + averageAdvectionDisplacementMagnitude += + glm::length(p.advectionDisplacementEstimate()); + averageDensityIntermediate += p.densityIntermediate(); + averageSumPressureDisplacementFromNeighborsMagnitude += + glm::length(p.sumPressureDisplacementFromNeighbors()); + averageVelocityIntermediateMagnitude += + glm::length(p.velocityIntermediate()); + averageNumNeighbors += p.neighbors()->size(); + + #if MFluidSolver_PARTICLE_STATS_FILES + { + std::lock_guard guard(statsCoutMutex); + if (pressureForceMagnitudeOutFile.is_open()) { + pressureForceMagnitudeOutFile << + glm::length(p.pressureForce()) << std::endl; + } + if (pressureOutFile.is_open()) { + pressureOutFile << p.pressure() << std::endl; + } + if (velocityMagnitudeOutFile.is_open()) { + velocityMagnitudeOutFile << glm::length(p.velocity()) << std::endl; + } + } + #endif + #endif + + // Calculate intermediate velocity v^{adv}_i + p.setVelocityIntermediate( + p.velocity() + p.nonPressureForce() / p.mass() * deltaTF); + + // Calculate advection displacement estimate d_{ii} + glm::vec3 advectionDisplacementEstimate(0); + const float pDensity2 = p.density() * p.density(); + for (SPHParticle *n : *(p.neighbors())) { + advectionDisplacementEstimate -= n->mass() / pDensity2 * + kernelFunctions.computeSpikyGradient(p.position() - n->position()); + } + p.setAdvectionDisplacementEstimate( + advectionDisplacementEstimate * deltaTF2); + iter_all_sphparticles_end + + iter_all_sphparticles_start + // Set intermediate density and advection diagonal a_{ii} + float advectionDiagonal = 0; + float advectionDensity = 0; + for (SPHParticle *n : *(p.neighbors())) { + // Calculate spikyGradientIJ + const glm::vec3 spikyGradientFromNeighbor = + kernelFunctions.computeSpikyGradient(p.position() - n->position()); + + // Calculate dJI + glm::vec3 displacementToNeighbor = + -deltaTF2 * p.mass() / (p.density() * p.density()) * + kernelFunctions.computeSpikyGradient(n->position() - p.position()); + + advectionDensity += n->mass() * + glm::dot( + p.velocityIntermediate() - n->velocityIntermediate(), + spikyGradientFromNeighbor); + + advectionDiagonal += n->mass() * + glm::dot( + p.advectionDisplacementEstimate() - displacementToNeighbor, + spikyGradientFromNeighbor); + } + p.setDensityIntermediate(p.density() + advectionDensity * deltaTF); + p.setAdvectionDiagonal(advectionDiagonal); + + // Set iteration=0 pressure p^0_i + p.setPressure(0.5f * p.pressure()); + iter_all_sphparticles_end + + // Print average statistics + #if MFluidSolver_PARTICLE_STATS + const unsigned int _numParticles = _particles.size(); + averagePosition /= _numParticles; + averageVelocityMagnitude /= _numParticles; + averageDensity /= _numParticles; + averagePressure /= _numParticles; + averageNonPressureForceMagnitude /= _numParticles; + averagePressureForceMagnitude /= _numParticles; + averageAdvectionDiagonal /= _numParticles; + averageAdvectionDisplacementMagnitude /= _numParticles; + averageDensityIntermediate /= _numParticles; + averageSumPressureDisplacementFromNeighborsMagnitude /= _numParticles; + averageVelocityIntermediateMagnitude /= _numParticles; + averageNumNeighbors /= _numParticles; + std::cout << "STAT: Printing particle statistics for frame " << + numUpdates << std::endl + << "STAT: Average position: (" << + averagePosition.x << ", " << + averagePosition.y << ", " << + averagePosition.z << ")" << std::endl + << "STAT: Average velocity magnitude: " << + averageVelocityMagnitude << std::endl + << "STAT: Average density: " << + averageDensity << std::endl + << "STAT: Average pressure: " << + averagePressure << std::endl + << "STAT: Average non-pressure force magnitude: " << + averageNonPressureForceMagnitude << std::endl + << "STAT: Average pressure force magnitude: " << + averagePressureForceMagnitude << std::endl + << "STAT: Average advection diagonal: " << + averageAdvectionDiagonal << std::endl + << "STAT: Average advection displacement magnitude: " << + averageAdvectionDisplacementMagnitude << std::endl + << "STAT: Average density intermediate: " << + averageDensityIntermediate << std::endl + << "STAT: Average sum pressure displacement from neighbors magnitude: " << + averageSumPressureDisplacementFromNeighborsMagnitude << std::endl + << "STAT: Average velocity intermediate magnitude: " << + averageVelocityIntermediateMagnitude << std::endl + << "STAT: Average number of neighbors: " << + averageNumNeighbors << std::endl; + #endif + // Print number of flyaways (use WARN channel if enabled) + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_WARN + // WARN, >0, any STAT + if (numFlyaways > 0) { + std::cout << "WARN: We have " << numFlyaways << + " flyaway particle(s) that could potentially crash the simulation" << + std::endl; + } + #if MFluidSolver_PARTICLE_STATS + else { + // WARN, =0, STAT + std::cout << "STAT: We have 0 flyaway particles!" << std::endl; + } + #endif + #elif MFluidSolver_PARTICLE_STATS + // !WARN, any #, STAT + std::cout << "STAT: We have " << numFlyaways << + " flyaway particle(s)!" << std::endl; + #endif + + // Procedure: Pressure Solve + unsigned int iteration = 0; + unsigned int validParticles; + float averageDensity = 0; + const float tolerance = 0.01f * kernelRadius; + std::vector nextPressures(_particles.size()); + while (((averageDensity - dRestDensity) > tolerance || iteration < 2) + && iteration < 1e6) { + averageDensity = 0; + validParticles = _particles.size(); + + iter_all_sphparticles_start + // Calculate pressure displacement due to neighbors + glm::vec3 pressureDisplacementFromNeighbors(0); + for (SPHParticle *n : *(p.neighbors())) { + pressureDisplacementFromNeighbors -= + n->mass() * n->pressure() / (n->density() * n->density()) * + kernelFunctions.computeSpikyGradient(p.position() - n->position()); + } + p.setSumPressureDisplacementFromNeighbors( + pressureDisplacementFromNeighbors * deltaTF2); + iter_all_sphparticles_end + + #if MFluidSolver_USE_TBB + averageDensity = + tbb::parallel_reduce( + tbb::blocked_range(0, _particles.size()), 0.f, + [&](const tbb::blocked_range &r, float partialDensitySum) { + for (unsigned int i = r.begin(); i != r.end(); ++i) { + #else + for (unsigned int i = 0; i < _particles.size(); ++i) { + #endif + SPHParticle &p = _particles[i]; + if (MUtils::fequal(p.advectionDiagonal(), 0.f, 0.001f)) { + // We want to disable the pressure force and reset its values + nextPressures[i] = 0; + --validParticles; + continue; + } + + // Calculate next pressure + constexpr float omega = 0.5f; + float densityDifferenceBySelf = 0; + float densityDifferenceByNeighborsPressure = 0; + for (SPHParticle *n : *(p.neighbors())) { + // Calculate spikyGradientIJ + glm::vec3 spikyGradientFromNeighbor = + kernelFunctions.computeSpikyGradient( + p.position() - n->position()); + + // Calculate dJI + glm::vec3 displacementToNeighbor = + -deltaTF2 * p.mass() / (p.density() * p.density()) * + kernelFunctions.computeSpikyGradient( + n->position() - p.position()); + + // Calculate density difference components + densityDifferenceBySelf += n->mass() * + glm::dot(p.advectionDisplacementEstimate() - + displacementToNeighbor, spikyGradientFromNeighbor); + + densityDifferenceByNeighborsPressure += n->mass() * glm::dot( + p.sumPressureDisplacementFromNeighbors() - + n->advectionDisplacementEstimate() * n->pressure() - + n->sumPressureDisplacementFromNeighbors() + + displacementToNeighbor * p.pressure(), spikyGradientFromNeighbor); + } + // Sum for average density + // Note that density depends on old pressure + const float tempDensityEstimate = p.densityIntermediate() + + p.pressure() * densityDifferenceBySelf + + densityDifferenceByNeighborsPressure; + // TODO: Check if we should do this + // p.setDensity(tempDensityEstimate); + + // Update average density + #if MFluidSolver_USE_TBB + if (tempDensityEstimate > dRestDensity) { + partialDensitySum += tempDensityEstimate; + } else { + partialDensitySum += dRestDensity; + } + #else + if (tempDensityEstimate > dRestDensity) { + averageDensity += tempDensityEstimate; + } else { + averageDensity += dRestDensity; + } + #endif + + // Calculate new pressure + float newPressure = dRestDensity - p.densityIntermediate() - + densityDifferenceByNeighborsPressure; + newPressure *= omega / p.advectionDiagonal(); + newPressure += (1.f - omega) * p.pressure(); + nextPressures[i] = newPressure; + #if MFluidSolver_USE_TBB + } // end particle for + return partialDensitySum; + }, // end lambda + std::plus() + ); // end parallel_reduce + #else + } // end particle for + #endif + + // Update particles with next iteration pressure + iter_all_sphparticles_start + p.setPressure(nextPressures[i]); + iter_all_sphparticles_end + + averageDensity /= static_cast(validParticles); + ++iteration; + } // end while + + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_TRACE + std::cout << "TRACE: Iterated pressure solve " << + iteration << " times" << std::endl; + std::cout << "TRACE: Density error: " << + (averageDensity - dRestDensity) << std::endl; + #endif + + // Procedure: Iteration + iter_all_sphparticles_start + calculatePressureForce(&p); + iter_all_sphparticles_end + + iter_all_sphparticles_start + // Update + glm::vec3 newVel = p.velocityIntermediate() + + p.pressureForce() / p.mass() * deltaTF; + glm::vec3 newPos = p.position() + newVel * deltaTF; + p.update(newVel, newPos); + + enforceBounds(&p); + visualizeParticle(&p); + iter_all_sphparticles_end +} diff --git a/src/fluidSolver/sphSolver/iiSphSolver.hpp b/src/fluidSolver/sphSolver/iiSphSolver.hpp new file mode 100644 index 00000000..dc623bcb --- /dev/null +++ b/src/fluidSolver/sphSolver/iiSphSolver.hpp @@ -0,0 +1,48 @@ +// Copyright 2016 Robert Zhou +// +// iiSphSolver.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_SPHSOLVER_IISPHSOLVER_HPP_ +#define MFLUIDSOLVER_SPHSOLVER_IISPHSOLVER_HPP_ + +#include + +#include "sphSolver.hpp" + +#if MFluidSolver_PARTICLE_STATS +#include +#endif + +class IISPHSolver : public SPHSolver { + public: + IISPHSolver(); + + virtual void update(double deltaT); + virtual void loadConfig(const std::string &file); + + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_WARN + unsigned int numFlyaways; + #endif + + protected: + unsigned int maxIterations; + + #if MFluidSolver_PARTICLE_STATS + std::mutex statsCoutMutex; + glm::vec3 averagePosition; + float averageVelocityMagnitude; + float averageDensity; + float averagePressure; + float averageNonPressureForceMagnitude; + float averagePressureForceMagnitude; + float averageAdvectionDiagonal; + float averageAdvectionDisplacementMagnitude; + float averageDensityIntermediate; + float averageSumPressureDisplacementFromNeighborsMagnitude; + float averageVelocityIntermediateMagnitude; + float averageNumNeighbors; + #endif +}; + +#endif // MFLUIDSOLVER_SPHSOLVER_IISPHSOLVER_HPP_ diff --git a/src/fluidSolver/sphSolver/kernelFunctions.cpp b/src/fluidSolver/sphSolver/kernelFunctions.cpp new file mode 100644 index 00000000..260eb759 --- /dev/null +++ b/src/fluidSolver/sphSolver/kernelFunctions.cpp @@ -0,0 +1,20 @@ +// Copyright 2016 Robert Zhou +// +// kernelFunctions.cpp +// MFluidSolver + +#include "kernelFunctions.hpp" + +KernelFunctions::KernelFunctions() { + setKernelRadius(0.1); +} + +void KernelFunctions::setKernelRadius(const double &h) { + _h = h; + _h2 = h * h; + _h3 = _h2 * h; + _h4 = _h2 * _h2; + _h5 = _h3 * _h2; + _h6 = _h3 * _h3; + _h9 = _h6 * _h3; +} diff --git a/src/fluidSolver/sphSolver/kernelFunctions.hpp b/src/fluidSolver/sphSolver/kernelFunctions.hpp new file mode 100644 index 00000000..5ba56e96 --- /dev/null +++ b/src/fluidSolver/sphSolver/kernelFunctions.hpp @@ -0,0 +1,43 @@ +// Copyright 2016 Robert Zhou +// +// kernelFunctions.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_SPHSOLVER_KERNELFUNCTIONS_HPP_ +#define MFLUIDSOLVER_SPHSOLVER_KERNELFUNCTIONS_HPP_ + +#include + +class KernelFunctions { + public: + KernelFunctions(); + + void setKernelRadius(const double &h); + + // Bicubic spline functions according to Ihmsen thesis + inline double computeBicubicSpline(const glm::vec3 &r); + inline double computeBicubicSplineFirstDerivativeDivHHalf(const double q); + inline double computeBicubicSplineFirstDerivativeDivHWhole(const double invQ); + inline glm::vec3 computeBicubicSplineGradient(const glm::vec3 &r); + inline double computeBicubicSplineLaplacian(const glm::vec3 &r); + + // Standard functions according to Muller 2003 + inline double computePoly6(const glm::vec3 &r); + inline double computeSpiky(const glm::vec3 &r); + inline glm::vec3 computeSpikyGradient(const glm::vec3 &r); + inline double computeViscous(const glm::vec3 &r); + inline double computeViscousLaplacian(const glm::vec3 &r); + + private: + double _h; + double _h2; + double _h3; + double _h4; + double _h5; + double _h6; + double _h9; +}; + +#include "kernelFunctions.inline.hpp" + +#endif // MFLUIDSOLVER_SPHSOLVER_KERNELFUNCTIONS_HPP_ diff --git a/src/fluidSolver/sphSolver/kernelFunctions.inline.hpp b/src/fluidSolver/sphSolver/kernelFunctions.inline.hpp new file mode 100644 index 00000000..70454e13 --- /dev/null +++ b/src/fluidSolver/sphSolver/kernelFunctions.inline.hpp @@ -0,0 +1,111 @@ +// Copyright 2016 Robert Zhou +// +// kernelFunctions.inline.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_SPHSOLVER_KERNELFUNCTIONS_INLINE_HPP_ +#define MFLUIDSOLVER_SPHSOLVER_KERNELFUNCTIONS_INLINE_HPP_ + +#define GLM_FORCE_RADIANS +#include + +// Super long constants +// constexpr double POLY6_CONST = 1.5666814710608447114749495456981882512; +// constexpr double SPIKY_CONST = 4.7746482927568600730665129011754308610; +// constexpr double EIGHT_DIV_PI = 2.5464790894703253723021402139602297925; +// constexpr double VISCOUS_CONST = 2.3873241463784300365332564505877154305; +// constexpr double FORTYFIVE_DIV_PI = 14.323944878270580219199538703526292583; +// constexpr double FORTYEIGHT_DIV_PI = 15.278874536821952233812841283761378755; + +// Shorter constants +constexpr double POLY6_CONST = 1.56668147106084471; +constexpr double SPIKY_CONST = 4.77464829275686007; +constexpr double EIGHT_DIV_PI = 2.54647908947032537; +constexpr double VISCOUS_CONST = 2.38732414637843004; +constexpr double FORTYFIVE_DIV_PI = 14.323944878270580; +constexpr double FORTYEIGHT_DIV_PI = 15.278874536821952; + +// Bicubic spline functions according to Ihmsen thesis +inline double KernelFunctions::computeBicubicSpline(const glm::vec3 &r) { + double q = glm::length(r) / _h; + if (q <= 0.5) { + return EIGHT_DIV_PI / _h3 * 6 * (q * q * q - q * q) + 1; + } else { + // Assumes q <= 1 + double invQ = 1 - q; + return EIGHT_DIV_PI / _h3 * 2 * invQ * invQ * invQ; + } +} + +inline double KernelFunctions::computeBicubicSplineFirstDerivativeDivHHalf( + const double q) { + return FORTYEIGHT_DIV_PI / _h4 * (3 * q * q - 2 * q); +} + +inline double KernelFunctions::computeBicubicSplineFirstDerivativeDivHWhole( + const double invQ) { + return FORTYEIGHT_DIV_PI / _h4 * (-1 * invQ * invQ); +} + +inline glm::vec3 KernelFunctions::computeBicubicSplineGradient( + const glm::vec3 &r) { + double rLen = glm::length(r); + double q = rLen / _h; + if (q <= 0.5) { + return static_cast( + computeBicubicSplineFirstDerivativeDivHHalf(q) / rLen) * r; + } else { + return static_cast( + computeBicubicSplineFirstDerivativeDivHWhole(1 - q) / rLen) * r; + } +} + +inline double KernelFunctions::computeBicubicSplineLaplacian( + const glm::vec3 &r) { + double rLen = glm::length(r); + double q = rLen / _h; + if (q <= 0.5) { + return FORTYEIGHT_DIV_PI * (6 * q - 2) / _h5 + + computeBicubicSplineFirstDerivativeDivHHalf(q) * 2 / rLen; + } else { + return FORTYEIGHT_DIV_PI * (2 - 2 * q) / _h5 + + computeBicubicSplineFirstDerivativeDivHWhole(1 - q) * 2 / rLen; + } +} + +// Standard functions according to Muller 2003 +inline double KernelFunctions::computePoly6(const glm::vec3 &r) { + double _r2 = glm::length2(r); + if (_r2 > _h2) return 0; + double h2MinusR2 = _h2 - _r2; + return POLY6_CONST * h2MinusR2 * h2MinusR2 * h2MinusR2 / _h9; +} + +inline double KernelFunctions::computeSpiky(const glm::vec3 &r) { + double _r = glm::length(r); + double hMinusR = _h - _r; + return SPIKY_CONST * hMinusR * hMinusR * hMinusR / _h6; +} + +inline glm::vec3 KernelFunctions::computeSpikyGradient(const glm::vec3 &r) { + double _r = glm::length(r); + if (_r == 0) return glm::vec3(0); + double hMinusR = _h - _r; + return static_cast( + -1 * FORTYFIVE_DIV_PI * hMinusR * hMinusR / _r /_h6) * r; +} + +inline double KernelFunctions::computeViscous(const glm::vec3 &r) { + double _r = glm::length(r); + return VISCOUS_CONST * + (-0.5 * _r * _r * _r / _h3 + + _r * _r / _h2 + + 0.5 * _h / _r - 1) / _h3; +} + +inline double KernelFunctions::computeViscousLaplacian(const glm::vec3 &r) { + double _r = glm::length(r); + return FORTYFIVE_DIV_PI * (_h - _r) / _h6; +} + +#endif // MFLUIDSOLVER_SPHSOLVER_KERNELFUNCTIONS_INLINE_HPP_ diff --git a/src/fluidSolver/sphSolver/neighborSearch.cpp b/src/fluidSolver/sphSolver/neighborSearch.cpp new file mode 100644 index 00000000..d1a6fbed --- /dev/null +++ b/src/fluidSolver/sphSolver/neighborSearch.cpp @@ -0,0 +1,122 @@ +// Copyright 2016 Robert Zhou +// +// neighborSearch.cpp +// MFluidSolver + +#include "neighborSearch.hpp" + +#include +#include + +#define GLM_FORCE_RADIANS +#include + +NeighborSearch::NeighborSearch(float r) { + setSearchRadius(r); +} + +void NeighborSearch::setSearchRadius(float r) { + searchRadius = r; + searchRadius2 = r * r; +} + +void NeighborSearch::printPerformanceStats() { +} + +void NaiveNeighborSearch::findNeighbors(SPHParticle *p) { + const glm::vec3 pPos = p->position(); + std::vector *neighbors = p->neighbors(); + neighbors->clear(); + + // Filter neighbors by kernel radius + for (SPHParticle *n : particleList) { + if (n != p) { + if (glm::distance2(n->position(), pPos) < searchRadius2) { + neighbors->push_back(n); + } + } + } +} + +void NaiveNeighborSearch::addParticle(SPHParticle *p) { + particleList.push_back(p); +} + +void NaiveNeighborSearch::clear() { + particleList.clear(); +} + +GridNeighborSearch::GridNeighborSearch(float r) +: NeighborSearch(r), grid(nullptr) { +} + +GridNeighborSearch::~GridNeighborSearch() { + delete grid; +} + +void GridNeighborSearch::findNeighbors(SPHParticle *p) { + grid->getNeighbors(p); + + const glm::vec3 pPos = p->position(); + std::vector *neighbors = p->neighbors(); + + // Filter neighbors by kernel radius. Also check for itself + unsigned int numCandidates = neighbors->size(); + unsigned int idx = 0; + unsigned int count = 0; + SPHParticle *n; + while (count < numCandidates) { + n = (*neighbors)[idx]; + if (glm::distance2(n->position(), pPos) >= searchRadius2 || n == p) { + neighbors->erase(neighbors->begin() + idx); + } else { + ++idx; + } + count++; + } +} + +void GridNeighborSearch::addParticle(SPHParticle *p) { + grid->addParticle(p); +} + +void GridNeighborSearch::updateParticle(SPHParticle *p) { + grid->updateParticle(p); +} + +void GridNeighborSearch::clear() { + grid->clear(); +} + +void GridNeighborSearch::printDiagnostics() { + grid->printDiagnostics(); +} + +#if MFluidSolver_USE_OPENVDB +void GridNeighborSearch::exportVDB() { + std::string filename = "export.vdb"; + std::string exportname = "MExport"; + grid->exportVDB(filename, exportname); +} +#endif + +UniformGridNeighborSearch::UniformGridNeighborSearch( + float r, const glm::vec3 &gridMin, const glm::vec3 &gridMax, float cellSize) + : GridNeighborSearch(r) { + grid = new SPHUniformGrid(gridMin, gridMax, cellSize); +} + +IndexSortedUniformGridNeighborSearch::IndexSortedUniformGridNeighborSearch( + float r, const glm::vec3 &gridMin, const glm::vec3 &gridMax, float cellSize, + std::vector *master, bool useZCurve) + : GridNeighborSearch(r) { + // Store a pointer to the grid so we don't need to cast + if (useZCurve) { + isuGrid = + new SPHZIndexSortedUniformGrid(gridMin, gridMax, cellSize, master); + } else { + isuGrid = + new SPHIndexSortedUniformGrid(gridMin, gridMax, cellSize, master); + } + grid = isuGrid; +} diff --git a/src/fluidSolver/sphSolver/neighborSearch.hpp b/src/fluidSolver/sphSolver/neighborSearch.hpp new file mode 100644 index 00000000..b6dc33bd --- /dev/null +++ b/src/fluidSolver/sphSolver/neighborSearch.hpp @@ -0,0 +1,88 @@ +// Copyright 2016 Robert Zhou +// +// neighborSearch.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_SPHSOLVER_NEIGHBORSEARCH_HPP_ +#define MFLUIDSOLVER_SPHSOLVER_NEIGHBORSEARCH_HPP_ + +#include + +#include "MFluidSolverConfig.hpp" + +#include "sphParticle.hpp" +#include "sphUniformGrid.hpp" +#include "sphZIndexSortedUniformGrid.hpp" + +enum NeighborSearchType { + Naive, UniformGrid, + IndexSortedUniformGrid, IndexSortedUniformGridWithInsertion, + ZIndexSortedUniformGrid, ZIndexSortedUniformGridWithInsertion}; + +class NeighborSearch { + public: + NeighborSearch() : NeighborSearch(MFluidSolver_DEFAULT_SEARCH_RADIUS) {} + explicit NeighborSearch(float r); + void setSearchRadius(float r); + void printPerformanceStats(); + virtual void findNeighbors(SPHParticle *p) = 0; + virtual void addParticle(SPHParticle *p) = 0; + virtual void updateParticle(SPHParticle *p) = 0; + virtual void clear() = 0; + + protected: + float searchRadius; + float searchRadius2; +}; + +class NaiveNeighborSearch : public NeighborSearch { + public: + NaiveNeighborSearch() + : NaiveNeighborSearch(MFluidSolver_DEFAULT_SEARCH_RADIUS) {} + explicit NaiveNeighborSearch(float r) : NeighborSearch(r) {} + virtual void findNeighbors(SPHParticle *p); + virtual void addParticle(SPHParticle *p); + virtual void updateParticle(SPHParticle *p) {} // Do nothing + virtual void clear(); + + private: + std::vector particleList; +}; + +class GridNeighborSearch : public NeighborSearch { + public: + explicit GridNeighborSearch(float r); + ~GridNeighborSearch(); + virtual void findNeighbors(SPHParticle *p); + virtual void addParticle(SPHParticle *p); + virtual void updateParticle(SPHParticle *p); + virtual void clear(); + virtual void printDiagnostics(); + + #if MFluidSolver_USE_OPENVDB + virtual void exportVDB(); + #endif + + SPHGrid *grid; +}; + +class UniformGridNeighborSearch : public GridNeighborSearch { + public: + UniformGridNeighborSearch( + float r, + const glm::vec3 &gridMin, const glm::vec3 &gridMax, + float cellSize); +}; + +class IndexSortedUniformGridNeighborSearch : public GridNeighborSearch { + public: + IndexSortedUniformGridNeighborSearch( + float r, + const glm::vec3 &gridMin, const glm::vec3 &gridMax, float cellSize, + std::vector *master, + bool useZCurve = MFluidSolver_DEFAULT_SPH_NEIGHBORSEARCH_USE_ZCURVE); + + SPHIndexSortedUniformGrid *isuGrid; +}; + +#endif // MFLUIDSOLVER_SPHSOLVER_NEIGHBORSEARCH_HPP_ diff --git a/src/fluidSolver/sphSolver/sphGrid.cpp b/src/fluidSolver/sphSolver/sphGrid.cpp new file mode 100644 index 00000000..8606a9da --- /dev/null +++ b/src/fluidSolver/sphSolver/sphGrid.cpp @@ -0,0 +1,62 @@ +// Copyright 2016 Robert Zhou +// +// sphGrid.cpp +// MFluidSolver + +#include "sphGrid.hpp" + +#include + +SPHGrid::SPHGrid( + const glm::vec3 &minBounds, const glm::vec3 &maxBounds, float cellSize) + : minBounds(minBounds), maxBounds(maxBounds), cellSize(cellSize) { + glm::vec3 gridSize = maxBounds - minBounds; + // Check grid bounds for sanity + if (gridSize.x < 0 || gridSize.y < 0 || gridSize.z < 0) { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_FATAL + std::cerr << "FATAL: Grid size (" << + gridSize.x << ", " << + gridSize.y << ", " << + gridSize.z << ") is invalid" << std::endl; + #endif + throw SPHGridNegativeSizeException(); + } + + // Calculate grid info + cellBounds = (gridSize + 0.500001f) / cellSize; + numCells = cellBounds.x * cellBounds.y * cellBounds.z; +} + +glm::ivec3 SPHGrid::getGridCoordinates(const glm::vec3 &pt) { + #if MFluidSolver_USE_ASSERTS + assert(pt.x >= minBounds.x && pt.y >= minBounds.y && pt.z >= minBounds.z); + assert(pt.x <= maxBounds.x && pt.y <= maxBounds.y && pt.z <= maxBounds.z); + #endif + return glm::ivec3(static_cast((pt.x - minBounds.x) / cellSize), + static_cast((pt.y - minBounds.y) / cellSize), + static_cast((pt.z - minBounds.z) / cellSize)); +} + +uint32_t SPHGrid::getIndex(const glm::vec3 &pt) { + return getIndex(getGridCoordinates(pt)); +} + +uint32_t SPHGrid::getIndex(const glm::ivec3 &c) { + #if MFluidSolver_USE_ASSERTS + assert(c.x >= 0 && c.y >= 0 && c.z >= 0); + assert(c.x < cellBounds.x && c.y < cellBounds.y && c.z < cellBounds.z); + #endif + + return getIndex(c.x, c.y, c.z); +} + +uint32_t SPHGrid::getIndex( + unsigned int x, unsigned int y, unsigned int z) { + #if MFluidSolver_USE_ASSERTS + uint32_t idx = y * cellBounds.z * cellBounds.x + z * cellBounds.x + x; + assert(idx >= 0 && idx < numCells); + return idx; + #else + return y * cellBounds.z * cellBounds.x + z * cellBounds.x + x; + #endif +} diff --git a/src/fluidSolver/sphSolver/sphGrid.hpp b/src/fluidSolver/sphSolver/sphGrid.hpp new file mode 100644 index 00000000..eada567e --- /dev/null +++ b/src/fluidSolver/sphSolver/sphGrid.hpp @@ -0,0 +1,56 @@ +// Copyright 2016 Robert Zhou +// +// sphGrid.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_FLUIDSOLVER_SPHGRID_HPP_ +#define MFLUIDSOLVER_FLUIDSOLVER_SPHGRID_HPP_ + +#include +#include +#include +#include + +#include "MFluidSolverConfig.hpp" + +#include + +#include "sphParticle.hpp" + +struct SPHGridNegativeSizeException : std::exception { + const char *what() const noexcept { + return "SPH grid size cannot be negative.\n"; + }; +}; + +class SPHGrid { + public: + SPHGrid(const glm::vec3 &minBounds, const glm::vec3 &maxBounds, + float cellSize); + + virtual void getNeighbors(SPHParticle *p) = 0; + + virtual glm::ivec3 getGridCoordinates(const glm::vec3 &pt); + virtual uint32_t getIndex(const glm::vec3 &pt); + virtual uint32_t getIndex(const glm::ivec3 &c); + virtual uint32_t getIndex(unsigned int x, unsigned int y, unsigned int z); + + virtual void addParticle(SPHParticle *p) = 0; + virtual void updateParticle(SPHParticle *p) = 0; + virtual void clear() = 0; + virtual void printDiagnostics() = 0; + + #if MFluidSolver_USE_OPENVDB + virtual void exportVDB(const std::string &file, + const std::string &gridName) = 0; + #endif + + protected: + glm::vec3 minBounds; + glm::vec3 maxBounds; + glm::ivec3 cellBounds; + uint32_t numCells; + float cellSize; +}; + +#endif // MFLUIDSOLVER_FLUIDSOLVER_SPHGRID_HPP_ diff --git a/src/fluidSolver/sphSolver/sphIndexSortedUniformGrid.cpp b/src/fluidSolver/sphSolver/sphIndexSortedUniformGrid.cpp new file mode 100644 index 00000000..353b017c --- /dev/null +++ b/src/fluidSolver/sphSolver/sphIndexSortedUniformGrid.cpp @@ -0,0 +1,168 @@ +// Copyright 2016 Robert Zhou +// +// sphIndexSortedUniformGrid.cpp +// MFluidSolver + +#include "sphIndexSortedUniformGrid.hpp" + +#include +#include +#include + +#if MFluidSolver_USE_ASSERTS +#include +#endif +#if MFluidSolver_USE_TBB +#include +#endif + +#include "utils.hpp" + +SPHIndexSortedUniformGrid::SPHIndexSortedUniformGrid( + const glm::vec3 &minBounds, const glm::vec3 &maxBounds, + float cellSize, std::vector *master) + : SPHGrid(minBounds, maxBounds, cellSize), master(master) { + cells.resize(numCells); +} + +SPHIndexSortedUniformGrid::~SPHIndexSortedUniformGrid() { +} + +void SPHIndexSortedUniformGrid::getNeighbors(SPHParticle *p) { + glm::ivec3 pC = getGridCoordinates(p->position()); + std::vector *neighbors = p->neighbors(); + neighbors->clear(); + + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_TRACE + std::cout << "TRACE: Getting neighbors for cell (" << + pC.x << ", " << pC.y << ", " << pC.z << ")" << std::endl; + #endif + + // Search neighboring grid cells + glm::ivec3 coords; + for (int i = -1; i < 2; ++i) { + for (int j = -1; j < 2; ++j) { + for (int k = -1; k < 2; ++k) { + coords = glm::ivec3(pC.x + i, pC.y + j, pC.z + k); + // Note: we check kernel radius and self in neighbor search + if (coords.x >= 0 && coords.y >= 0 && coords.z >= 0 && + coords.x < cellBounds.x && + coords.y < cellBounds.y && + coords.z < cellBounds.z) { + uint32_t index = getIndex(coords); + + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_TRACE + unsigned int numNeighbors = 0; + std::cout << "TRACE: Accessing cell (" << + coords.x << ", " << coords.y << ", " << coords.z << ")" << + std::endl; + #endif + + // Add all particles in index sorted list + SPHParticle *c = cells[index]; + if (c != nullptr) { + do { + neighbors->push_back(c); + + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_TRACE + ++numNeighbors; + #endif + + if (c == endParticle) { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_TRACE + std::cout << "TRACE: Reached end particle" << std::endl; + #endif + break; + } + ++c; + } while (index == c->index); + } + + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_TRACE + std::cout << "TRACE: Found " << + numNeighbors << " neighbor(s)" << std::endl; + #endif + } + } // end for k + } // end for j + } // end for i +} + +void SPHIndexSortedUniformGrid::addParticle(SPHParticle *p) { +} + +void SPHIndexSortedUniformGrid::updateParticle(SPHParticle *p) { +} + +void SPHIndexSortedUniformGrid::clear() { + for (unsigned int i = 0; i < cells.size(); ++i) { + cells[i] = nullptr; + } +} + +void SPHIndexSortedUniformGrid::printDiagnostics() { + std::cout << "No diagnostics available for index sorted uniform grids." << + std::endl; +} + +void SPHIndexSortedUniformGrid::resetAndFillCells(bool initialSort) { + clear(); + if (master->empty()) return; + updateParticleIndices(); + sortParticles(initialSort); + insertSortedParticleListToGrid(); +} + +void SPHIndexSortedUniformGrid::updateParticleIndices() { + for (SPHParticle &p : *master) { + p.index = getIndex(p.position()); + } +} + +void SPHIndexSortedUniformGrid::insertSortedParticleListToGrid() { + #if MFluidSolver_USE_ASSERTS + assert((*master)[0].index >= 0 && (*master)[0].index < cells.size()); + #endif + cells[(*master)[0].index] = &((*master)[0]); + + // Set pointer in each cell to first particle with associated cell index + #if MFluidSolver_USE_TBB + tbb::parallel_for((size_t)1, master->size(), + [&](size_t i) { + #else + for (size_t i = 1; i < master->size(); ++i) { + #endif + #if MFluidSolver_USE_ASSERTS + assert((*master)[i - 1].index >= 0 && + (*master)[i - 1].index < cells.size()); + assert((*master)[i].index >= 0 && + (*master)[i].index < cells.size()); + #endif + if ((*master)[i].index != (*master)[i - 1].index) { + cells[(*master)[i].index] = &((*master)[i]); + } + #if MFluidSolver_USE_TBB + } + ); + #else + } + #endif +} + +void SPHIndexSortedUniformGrid::sortParticles(bool initialSort) { + if (initialSort) { + std::sort(master->begin(), master->end(), SPHParticle::indexCompare); + } else { + // Insertion sort is faster for minimal changes in cell indices + MUtils::insertionSort(master->begin(), master->end(), + SPHParticle::indexCompare); + } + endParticle = &(master->back()); +} + +#if MFluidSolver_USE_OPENVDB +void SPHIndexSortedUniformGrid::exportVDB( + const std::string &file, const std::string &gridName) { + // TODO: export VDB +} +#endif diff --git a/src/fluidSolver/sphSolver/sphIndexSortedUniformGrid.hpp b/src/fluidSolver/sphSolver/sphIndexSortedUniformGrid.hpp new file mode 100644 index 00000000..40a5b669 --- /dev/null +++ b/src/fluidSolver/sphSolver/sphIndexSortedUniformGrid.hpp @@ -0,0 +1,46 @@ +// Copyright 2016 Robert Zhou +// +// sphIndexSortedUniformGrid.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_FLUIDSOLVER_SPHINDEXSORTEDUNIFORMGRID_HPP_ +#define MFLUIDSOLVER_FLUIDSOLVER_SPHINDEXSORTEDUNIFORMGRID_HPP_ + +#include +#include + +#include "sphGrid.hpp" + +// Grid type XZY (X elements together, then Z, then Y) +class SPHIndexSortedUniformGrid : public SPHGrid { + public: + SPHIndexSortedUniformGrid( + const glm::vec3 &minBounds, const glm::vec3 &maxBounds, + float cellSize, std::vector *master); + ~SPHIndexSortedUniformGrid(); + + // SPHGrid methods + virtual void getNeighbors(SPHParticle *p); + + virtual void addParticle(SPHParticle *p); + virtual void updateParticle(SPHParticle *p); + virtual void clear(); + virtual void printDiagnostics(); + + #if MFluidSolver_USE_OPENVDB + virtual void exportVDB(const std::string &file, const std::string &gridName); + #endif + + // Index sort methods + virtual void resetAndFillCells(bool initialSort = false); + virtual void updateParticleIndices(); + virtual void insertSortedParticleListToGrid(); + void sortParticles(bool initialSort = false); + + protected: + std::vector *master; + std::vector cells; + SPHParticle *endParticle; +}; + +#endif // MFLUIDSOLVER_FLUIDSOLVER_SPHINDEXSORTEDUNIFORMGRID_HPP_ diff --git a/src/fluidSolver/sphSolver/sphParticle.cpp b/src/fluidSolver/sphSolver/sphParticle.cpp new file mode 100644 index 00000000..16fe7f86 --- /dev/null +++ b/src/fluidSolver/sphSolver/sphParticle.cpp @@ -0,0 +1,22 @@ +// Copyright 2016 Robert Zhou +// +// sphParticle.cpp +// MFluidSolver + +#include "sphParticle.hpp" + +SPHParticle::SPHParticle(float mass, const glm::vec3 &position) + : Particle(mass, position), + _density(1000.f), _pressure(0.f), + _nonPressureForce(0), _pressureForce(0), _oldPosition(position), + _advectionDiagonal(0.f), _advectionDisplacementEstimate(0), + _sumPressureDisplacementFromNeighbors(0), + _densityIntermediate(1000.f) { + #if MFluidSolver_PARTICLE_STATS + flyaway = false; + #endif +} + +bool SPHParticle::indexCompare(const SPHParticle &a, const SPHParticle &b) { + return a.index < b.index; +} diff --git a/src/fluidSolver/sphSolver/sphParticle.hpp b/src/fluidSolver/sphSolver/sphParticle.hpp new file mode 100644 index 00000000..75e59b08 --- /dev/null +++ b/src/fluidSolver/sphSolver/sphParticle.hpp @@ -0,0 +1,77 @@ +// Copyright 2016 Robert Zhou +// +// sphParticle.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_FLUIDSOLVER_SPHSOLVER_SPHPARTICLE_HPP_ +#define MFLUIDSOLVER_FLUIDSOLVER_SPHSOLVER_SPHPARTICLE_HPP_ + +#include + +#include "fluidSolver/particle.hpp" + +class SPHParticle : public Particle { + public: + SPHParticle() : SPHParticle(1, glm::vec3(0)) {} + explicit SPHParticle(const glm::vec3 &position) : SPHParticle(1, position) {} + SPHParticle(float mass, const glm::vec3 &position); + + virtual inline void update(const glm::vec3 &newVel, const glm::vec3 &newPos); + + // Properties + inline float density() const; + inline float pressure() const; + inline glm::vec3 force() const; + inline glm::vec3 nonPressureForce() const; + inline glm::vec3 oldPosition() const; + inline glm::vec3 pressureForce() const; + + inline void setDensity(float density); + inline void setNonPressureForce(const glm::vec3 &force); + inline void setPressure(float pressure); + inline void setPressureForce(const glm::vec3 &force); + inline void setVelocity(const glm::vec3 &velocity); + + inline std::vector *neighbors(); + inline void clearNeighbors(); + + // IISPH Properties + inline float advectionDiagonal() const; + inline glm::vec3 advectionDisplacementEstimate() const; + inline float densityIntermediate() const; + inline glm::vec3 sumPressureDisplacementFromNeighbors() const; + inline glm::vec3 velocityIntermediate() const; + + inline void setAdvectionDiagonal(float a); + inline void setAdvectionDisplacementEstimate(const glm::vec3 &d); + inline void setDensityIntermediate(float d); + inline void setSumPressureDisplacementFromNeighbors(const glm::vec3 &d); + inline void setVelocityIntermediate(const glm::vec3 &v); + + // Index + uint32_t index; + static bool indexCompare(const SPHParticle &a, const SPHParticle &b); + + #if MFluidSolver_PARTICLE_STATS + bool flyaway; + #endif + + protected: + float _density; + float _pressure; + glm::vec3 _nonPressureForce; + glm::vec3 _pressureForce; + glm::vec3 _oldPosition; + std::vector _neighbors; + + // IISPH + float _advectionDiagonal; + glm::vec3 _advectionDisplacementEstimate; + float _densityIntermediate; + glm::vec3 _sumPressureDisplacementFromNeighbors; + glm::vec3 _velocityIntermediate; +}; + +#include "sphParticle.inline.hpp" + +#endif // MFLUIDSOLVER_FLUIDSOLVER_SPHSOLVER_SPHPARTICLE_HPP_ diff --git a/src/fluidSolver/sphSolver/sphParticle.inline.hpp b/src/fluidSolver/sphSolver/sphParticle.inline.hpp new file mode 100644 index 00000000..2014bb95 --- /dev/null +++ b/src/fluidSolver/sphSolver/sphParticle.inline.hpp @@ -0,0 +1,116 @@ +// Copyright 2016 Robert Zhou +// +// sphParticle.inline.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_FLUIDSOLVER_SPHSOLVER_SPHPARTICLE_INLINE_HPP_ +#define MFLUIDSOLVER_FLUIDSOLVER_SPHSOLVER_SPHPARTICLE_INLINE_HPP_ + +#include + +inline void SPHParticle::update( + const glm::vec3 &newVel, const glm::vec3 &newPos) { + _oldPosition = _position; + _velocity = newVel; + _position = newPos; +} + +// Getters +inline float SPHParticle::density() const { + return _density; +} + +inline glm::vec3 SPHParticle::force() const { + return _nonPressureForce + _pressureForce; +} + +inline glm::vec3 SPHParticle::nonPressureForce() const { + return _nonPressureForce; +} + +inline glm::vec3 SPHParticle::oldPosition() const { + return _oldPosition; +} + +inline float SPHParticle::pressure() const { + return _pressure; +} + +inline glm::vec3 SPHParticle::pressureForce() const { + return _pressureForce; +} + +// Setters +inline void SPHParticle::setDensity(float density) { + _density = density; +} + +inline void SPHParticle::setNonPressureForce(const glm::vec3 &force) { + _nonPressureForce = force; +} + +inline void SPHParticle::setPressure(float pressure) { + _pressure = pressure; +} + +inline void SPHParticle::setPressureForce(const glm::vec3 &force) { + _pressureForce = force; +} + +inline void SPHParticle::setVelocity(const glm::vec3 &velocity) { + _velocity = velocity; +} + +// Neighbors +inline std::vector *SPHParticle::neighbors() { + return &_neighbors; +} + +inline void SPHParticle::clearNeighbors() { + _neighbors.clear(); +} + +// IISPH Getters +inline float SPHParticle::advectionDiagonal() const { + return _advectionDiagonal; +} + +inline glm::vec3 SPHParticle::advectionDisplacementEstimate() const { + return _advectionDisplacementEstimate; +} + +inline float SPHParticle::densityIntermediate() const { + return _densityIntermediate; +} + +inline glm::vec3 SPHParticle::sumPressureDisplacementFromNeighbors() const { + return _sumPressureDisplacementFromNeighbors; +} + +inline glm::vec3 SPHParticle::velocityIntermediate() const { + return _velocityIntermediate; +} + +// IISPH Setters +inline void SPHParticle::setAdvectionDiagonal(float a) { + _advectionDiagonal = a; +} + +inline void SPHParticle::setAdvectionDisplacementEstimate(const glm::vec3 &d) { + _advectionDisplacementEstimate = d; +} + +inline void SPHParticle::setDensityIntermediate(float d) { + _densityIntermediate = d; +} + +inline void SPHParticle::setSumPressureDisplacementFromNeighbors( + const glm::vec3 &d) { + _sumPressureDisplacementFromNeighbors = d; +} + +inline void SPHParticle::setVelocityIntermediate(const glm::vec3 &v) { + _velocityIntermediate = v; +} + +#endif // MFLUIDSOLVER_FLUIDSOLVER_SPHSOLVER_SPHPARTICLE_INLINE_HPP_ diff --git a/src/fluidSolver/sphSolver/sphSolver.cpp b/src/fluidSolver/sphSolver/sphSolver.cpp new file mode 100644 index 00000000..97e560f9 --- /dev/null +++ b/src/fluidSolver/sphSolver/sphSolver.cpp @@ -0,0 +1,555 @@ +// Copyright 2016 Robert Zhou +// +// sphSolver.cpp +// MFluidSolver + +#include "sphSolver.hpp" + +#include +#include +#include +#include +#include + +#if MFluidSolver_PARTICLE_STATS +#include +#endif + +#include + +#include "utils.hpp" + +// Constructor / Destructor +SPHSolver::SPHSolver() + : nSearch(nullptr), inited(false) { + setDefaultConfig(); +} + +SPHSolver::~SPHSolver() { + delete nSearch; +} + +// Configuration +void SPHSolver::init(const glm::vec3 &gridMin, const glm::vec3 &gridMax) { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + std::cout << "INFO: Initializing SPH Solver" << std::endl; + #endif + + kernelFunctions.setKernelRadius(kernelRadius); + + // Initialize neighbor search + // Note: Neighbor search assumes grid cell size is equal to kernelRadius + switch (nSearchType) { + case NeighborSearchType::Naive: { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + std::cout << "INFO: SPH Solver using naive neighbor search" << std::endl; + #endif + nSearch = new NaiveNeighborSearch(); + break; + } + case NeighborSearchType::IndexSortedUniformGrid: { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + std::cout << + "INFO: SPH Solver using index sorted uniform grid neighbor search" << + std::endl; + #endif + nSearch = new IndexSortedUniformGridNeighborSearch(kernelRadius, + gridMin, gridMax, + kernelRadius, + &_particles, false); + break; + } + case NeighborSearchType::ZIndexSortedUniformGrid: { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + std::cout << "INFO: SPH Solver using " << + "Z-curve index sorted uniform grid neighbor search" << std::endl; + #endif + nSearch = new IndexSortedUniformGridNeighborSearch(kernelRadius, + gridMin, gridMax, + kernelRadius, + &_particles, true); + break; + } + case NeighborSearchType::ZIndexSortedUniformGridWithInsertion: { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + std::cout << "INFO: SPH Solver using " << + "Z-curve index sorted uniform grid neighbor search " << + "with insertion sort" << std::endl; + #endif + nSearch = new IndexSortedUniformGridNeighborSearch(kernelRadius, + gridMin, gridMax, + kernelRadius, + &_particles, true); + break; + } + case NeighborSearchType::UniformGrid: { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + std::cout << "INFO: SPH Solver using uniform grid neighbor search" << + std::endl; + #endif + nSearch = new UniformGridNeighborSearch(kernelRadius, + gridMin, gridMax, kernelRadius); + break; + } + case NeighborSearchType::IndexSortedUniformGridWithInsertion: { + default: + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + std::cout << + "INFO: SPH Solver using " << + "index sorted uniform grid neighbor search with insertion sort" << + std::endl; + #endif + nSearch = new IndexSortedUniformGridNeighborSearch(kernelRadius, + gridMin, gridMax, + kernelRadius, + &_particles, false); + break; + } + } + + // Print useful parameter info + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + switch (visualizationType) { + case FluidVisualizationType::Density: { + std::cout << + "INFO: Now visualizing density difference from rest density" << + std::endl; + break; + } + case FluidVisualizationType::Index: { + std::cout << + "INFO: Now visualizing index within particle array" << + std::endl; + break; + } + case FluidVisualizationType::Neighbors: { + std::cout << "INFO: Now visualizing neighbors" << std::endl; + break; + } + case FluidVisualizationType::Particle: { + std::cout << + "INFO: Now highlighting particle ID " << targetParticle << std::endl; + break; + } + case FluidVisualizationType::Pressure: { + std::cout << "INFO: Now visualizing pressure" << std::endl; + break; + } + case FluidVisualizationType::Velocity: { + std::cout << "INFO: Now visualizing scalar velocity" << std::endl; + break; + } + case FluidVisualizationType::VelocityDir: { + std::cout << + "INFO: Now visualizing normalized velocity direction" << std::endl; + break; + } + case FluidVisualizationType::None: + default: { + std::cout << "INFO: Visualization disabled" << std::endl; + break; + } + } + if (muViscosity <= 0) { + std::cout << "INFO: As per config, viscosity is disabled" << std::endl; + } + #endif + + inited = true; +} + +void SPHSolver::setDefaultConfig() { + // Set default SPH parameters + if (checkInited()) return; + kStiffness = MFluidSolver_DEFAULT_SPH_STIFFNESS; + muViscosity = MFluidSolver_DEFAULT_SPH_VISCOSITY; + mMass = MFluidSolver_DEFAULT_SPH_MASS; + dRestDensity = MFluidSolver_DEFAULT_SPH_RESTDENSITY; + dtTimestep = MFluidSolver_DEFAULT_SPH_TIMESTEP; + kernelRadius = MFluidSolver_DEFAULT_SPH_KERNELRADIUS; + nSearchType = MFluidSolver_DEFAULT_SPH_NEIGHBORSEARCHTYPE; +} + +void SPHSolver::loadConfig(const std::string &file) { + if (checkInited()) return; + + logTimestep(); + + // Parse JSON file + Json::Reader reader; + Json::Value root; + std::ifstream sceneStream(file, std::ifstream::binary); + bool success = reader.parse(sceneStream, root, false); + if (!success) { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_ERROR + std::cerr << "ERROR: Failed to parse config file " << file << std::endl; + #endif + return; + } + + // Read SPH Parameters + kStiffness = root["sph"].get( + "kStiffness", MFluidSolver_DEFAULT_SPH_STIFFNESS).asFloat(); + muViscosity = root["sph"].get( + "muViscosity", MFluidSolver_DEFAULT_SPH_VISCOSITY).asFloat(); + mMass = root["sph"].get( + "mMass", MFluidSolver_DEFAULT_SPH_MASS).asFloat(); + dRestDensity = root["sph"].get( + "dRestDensity", MFluidSolver_DEFAULT_SPH_RESTDENSITY).asFloat(); + dtTimestep = root["sph"].get( + "dtTimestep", MFluidSolver_DEFAULT_SPH_TIMESTEP).asFloat(); + kernelRadius = root["sph"].get( + "kernelRadius", MFluidSolver_DEFAULT_SPH_KERNELRADIUS).asFloat(); + setFixedTimestep(dtTimestep); + + // Read simulation time limits + int tempMaxUpdates = root.get("numUpdates", 0).asInt(); + if (tempMaxUpdates <= 0) { + limitNumUpdates = false; + } else { + maxUpdates = tempMaxUpdates; + limitNumUpdates = true; + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + std::cout << + "INFO: Limiting number of updates to " << maxUpdates << std::endl; + #endif + } + + // Read neighbor search config + std::string neighborSearchTypeString = root["sph"].get( + "neighborSearchType", + MFluidSolver_DEFAULT_SPH_NEIGHBORSEARCHTYPE_STRING).asString(); + MUtils::toLowerInplace(&neighborSearchTypeString); + if (neighborSearchTypeString == "naive") { + nSearchType = NeighborSearchType::Naive; + } else if (neighborSearchTypeString == "uniform") { + nSearchType = NeighborSearchType::UniformGrid; + } else if (neighborSearchTypeString == "indexsorteduniform") { + nSearchType = NeighborSearchType::IndexSortedUniformGrid; + } else if (neighborSearchTypeString == "indexsorteduniforminsertionsort") { + nSearchType = NeighborSearchType::IndexSortedUniformGridWithInsertion; + } else if (neighborSearchTypeString == "zindexsorteduniform") { + nSearchType = NeighborSearchType::ZIndexSortedUniformGrid; + } else if (neighborSearchTypeString == "zindexsorteduniforminsertionsort") { + nSearchType = NeighborSearchType::ZIndexSortedUniformGridWithInsertion; + } + + // Read visualization config + std::string visualizationTypeString = root["visualization"].get( + "type", MFluidSolver_DEFAULT_VISUALIZATION_STRING).asString(); + MUtils::toLowerInplace(&visualizationTypeString); + if (visualizationTypeString == "density") { + visualizationType = FluidVisualizationType::Density; + } else if (visualizationTypeString == "index") { + visualizationType = FluidVisualizationType::Index; + } else if (visualizationTypeString == "neighbors") { + visualizationType = FluidVisualizationType::Neighbors; + } else if (visualizationTypeString == "none") { + visualizationType = FluidVisualizationType::None; + } else if (visualizationTypeString == "particle") { + visualizationType = FluidVisualizationType::Particle; + } else if (visualizationTypeString == "pressure") { + visualizationType = FluidVisualizationType::Pressure; + } else if (visualizationTypeString == "velocity") { + visualizationType = FluidVisualizationType::Velocity; + } else if (visualizationTypeString == "velocitydir") { + visualizationType = FluidVisualizationType::VelocityDir; + } + visualizationMaxDensityDifference = + root["visualization"].get( + "maxDensityDifference", + MFluidSolver_DEFAULT_VISUALIZATION_MAXDENSITYDIFFERENCE).asFloat(); + visualizationMaxPressure = + root["visualization"].get( + "maxPressure", + MFluidSolver_DEFAULT_VISUALIZATION_MAXPRESSURE).asFloat(); + visualizationMaxVelocity = + root["visualization"].get( + "maxVelocity", + MFluidSolver_DEFAULT_VISUALIZATION_MAXVELOCITY).asFloat(); + targetParticle = + root["visualization"].get( + "particleId", + 0).asInt(); + const Json::Value velocityColorArray = + root["visualization"]["velocityColor"]; + if (!velocityColorArray.isNull()) { + for (unsigned int i = 0; i < velocityColorArray.size() && i < 3; ++i) { + visualizationVelocityColor[i] = + static_cast(velocityColorArray[i].asInt()) / 255.f; + if (visualizationVelocityColor[i] > 1.f) + visualizationVelocityColor[i] = 1.f; + } + } + + // Print parameters + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + std::cout << "INFO: - Stiffness: " << kStiffness << std::endl; + std::cout << "INFO: - Viscosity: " << muViscosity << std::endl; + std::cout << "INFO: - Mass: " << mMass << std::endl; + std::cout << "INFO: - Rest Density: " << dRestDensity << std::endl; + std::cout << "INFO: - Timestep: " << dtTimestep << std::endl; + std::cout << "INFO: - Kernel Radius: " << kernelRadius << std::endl; + #endif +} + +// Solver! +void SPHSolver::update(double deltaT) { + if (checkIfEnded()) return; + + // NOTE: TIMESTEP IS OVERWRITTEN HERE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + deltaT = _fixedTimestep; + const float deltaTF = static_cast(deltaT); + + logTimestep(); + prepNeighborSearch(); + runNeighborSearch(); + + // Compute density and pressure + iter_all_sphparticles_start + calculateDensity(&p); + + // Pressure + double densityRatio = p.density() / dRestDensity; + double pressureTemp = kStiffness * dRestDensity / 7.f * ( + densityRatio * densityRatio * densityRatio * densityRatio * + densityRatio * densityRatio * densityRatio - 1.f); + if (pressureTemp < 0) pressureTemp = 0; + p.setPressure(pressureTemp); + iter_all_sphparticles_end + + // Compute forces + iter_all_sphparticles_start + calculateNonPressureForce(&p); + calculatePressureForce(&p); + iter_all_sphparticles_end + + // Compute velocity and position, check bounds + iter_all_sphparticles_start + // Update + glm::vec3 newVel = p.velocity() + p.force() / p.mass() * deltaTF; + glm::vec3 newPos = p.position() + newVel * deltaTF; + p.update(newVel, newPos); + + enforceBounds(&p); + + largestIndex = _particles.back().index; + visualizeParticle(&p); + iter_all_sphparticles_end +} + +// Particles +void SPHSolver::addParticleAt(const glm::vec3 &position) { + if (_particles.size() < _maxParticles) { + SPHParticle pStack(mMass, position); + _particles.push_back(pStack); + SPHParticle *pList = &(_particles.back()); + pList->index = _particles.size() - 1; + nSearch->addParticle(pList); + } +} + +unsigned int SPHSolver::numParticles() const { + return _particles.size(); +} + +std::vector &SPHSolver::particles() { + return _particles; +} + +void SPHSolver::setMaxParticles(int mp) { + FluidSolver::setMaxParticles(mp); + if (mp <= 0) { + // Remove all particles as requested + _particles.clear(); + } else if (mp > _particles.size()) { + // Prune excess particles + while (_particles.size() > mp) { + _particles.erase(_particles.end() - 1); + } + } +} + +void SPHSolver::setParticleSeparation(float ps) { + FluidSolver::setParticleSeparation(ps); +} + +// Neighbor Search +void SPHSolver::prepNeighborSearchAfterSceneLoad() { + // Initially sort particles with standard sort + if (nSearchType == NeighborSearchType::IndexSortedUniformGridWithInsertion || + nSearchType == NeighborSearchType::ZIndexSortedUniformGridWithInsertion) { + IndexSortedUniformGridNeighborSearch *isugSearch = + static_cast(nSearch); + isugSearch->isuGrid->resetAndFillCells(true); + } +} + +// Neighbor Visualization +void SPHSolver::visualizeParticleNeighbors(SPHParticle *target) { + if (visualizationType != FluidVisualizationType::Neighbors) return; + + // Reset color + for (SPHParticle &p : _particles) { + p.color = glm::vec3(0, 0, 1); + } + + // Redo neighbor search + prepNeighborSearch(); + target->clearNeighbors(); + nSearch->findNeighbors(target); + + // Color particles + target->color = glm::vec3(1, 1, 0); + for (SPHParticle *n : *(target->neighbors())) { + n->color = glm::vec3(1, 0, 0); + } +} + +void SPHSolver::visualizeRandomParticlesNeighbors() { + if (!_particles.empty()) { + visualizeParticleNeighbors(&(_particles[rand() % _particles.size()])); + } + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_WARN + else { + std::cout << + "WARN: No particles available for neighbor visualization" << std::endl; + } + #endif +} + +// Misc +void SPHSolver::initVisualization() { + if (visualizationType == FluidVisualizationType::Neighbors) { + if (!_particles.empty()) { + visualizeParticleNeighbors(&(_particles[0])); + } + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_WARN + else { + std::cout << + "WARN: No particles available for neighbor visualization" << std::endl; + } + #endif + } else if (visualizationType == FluidVisualizationType::Particle) { + if (_particles.size() > targetParticle) { + _particles[targetParticle].color = glm::vec3(1, 1, 0); + } + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_WARN + else { + std::cout << "WARN: Target particle " << targetParticle << + " not available for visualization" << std::endl; + } + #endif + } + + largestIndex = _particles.back().index; + iter_all_sphparticles_start + visualizeParticle(&p); + iter_all_sphparticles_end +} + +void SPHSolver::printPerformanceStats() { + FluidSolver::printPerformanceStats(); + nSearch->printPerformanceStats(); +} + +void SPHSolver::sceneLoaded() { + initVisualization(); + prepNeighborSearchAfterSceneLoad(); + + // Gather stats concerning proximity (could lead to explosive pressures) + #if MFluidSolver_PARTICLE_STATS + prepNeighborSearch(); + runNeighborSearch(); + if (!_particles.empty()) { + float minDist = std::numeric_limits::infinity(); + float distSum = 0; + float minDistSum = 0; + unsigned int numPairs = 0; + + for (SPHParticle &p : _particles) { + float minDistParticle = std::numeric_limits::infinity(); + if (!p.neighbors()->empty()) { + for (SPHParticle *n : *(p.neighbors())) { + ++numPairs; + const float dist = glm::length(p.position() - n->position()); + if (dist < minDistParticle) { + minDistParticle = dist; + } + distSum += dist; + } // end for neighbors + + minDistSum += minDistParticle; + if (minDistParticle < minDist) { + minDist = minDistParticle; + } + } + } // end for particles + + distSum /= static_cast(numPairs); + minDistSum /= static_cast(_particles.size()); + + std::cout << "STAT: Average particle distance to neighbor: " << + distSum << std::endl; + std::cout << "STAT: Average particle distance to nearest neighbor: " << + minDistSum << std::endl; + std::cout << "STAT: Smallest particle distance to neighbor: " << + minDist << std::endl; + } + #endif +} + +#if MFluidSolver_USE_PARTIO +void SPHSolver::exportBgeo() { + // http://www.disneyanimation.com/technology/partio.html + // http://wdas.github.io/partio/doxygen/ + Partio::ParticlesDataMutable *particleData = Partio::create(); + Partio::ParticlesDataMutable::iterator iterator = + particleData->addParticles(_particles.size()); + Partio::ParticleAttribute posAttr = + particleData->addAttribute( + "position", Partio::ParticleAttributeType::VECTOR, 3); + Partio::ParticleAccessor posAcc(posAttr); + iterator.addAccessor(posAcc); + + unsigned int idx = 0; + for (Partio::ParticlesDataMutable::iterator it = + particleData->begin(); it != particleData->end(); ++it) { + float *vectorData = posAcc.raw(it); + const glm::vec3 position = _particles.at(idx).position(); + for (unsigned int i = 0; i < 3; ++i) { + vectorData[i] = position[i]; + } + ++idx; + } + + Partio::write("test.bgeo", *particleData); + particleData->release(); +} +#endif + +#if MFluidSolver_USE_OPENVDB +void SPHSolver::exportVDB() { + if (nSearchType == NeighborSearchType::UniformGrid) { + static_cast(nSearch)->exportVDB(); + } else { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_WARN + std::cout << + "WARN: Cannot export VDB unless using grid neighbor search" << std::endl; + #endif + } +} +#endif + +// Initialization +bool SPHSolver::checkInited() { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_ERROR + if (inited) { + std::cerr << + "ERROR: Already inited SPH Solver. Config will not be applied" << + std::endl; + } + #endif + + return inited; +} diff --git a/src/fluidSolver/sphSolver/sphSolver.hpp b/src/fluidSolver/sphSolver/sphSolver.hpp new file mode 100644 index 00000000..7bf3fb6c --- /dev/null +++ b/src/fluidSolver/sphSolver/sphSolver.hpp @@ -0,0 +1,100 @@ +// Copyright 2016 Robert Zhou +// +// sphSolver.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_SPHSOLVER_HPP_ +#define MFLUIDSOLVER_SPHSOLVER_HPP_ + +#include +#include + +#include "MFluidSolverConfig.hpp" + +#if MFluidSolver_USE_PARTIO +#include +#endif +#if MFluidSolver_USE_TBB +#include +#include +#endif + +#include "fluidSolver/fluidSolver.hpp" +#include "kernelFunctions.hpp" +#include "neighborSearch.hpp" +#include "sphParticle.hpp" + +class SPHSolver : public FluidSolver { + public: + // Constructor / Destructor + SPHSolver(); + ~SPHSolver(); + + // Configuration + void init(const glm::vec3 &gridMin, const glm::vec3 &gridMax); + void setDefaultConfig(); + void loadConfig(const std::string &file); + + // Solver! + virtual void update(double deltaT); + + // Particles + virtual void addParticleAt(const glm::vec3 &position); + virtual unsigned int numParticles() const; + virtual std::vector &particles(); + virtual void setMaxParticles(int mp); + virtual void setParticleSeparation(float ps); + + // Neighbor Visualization + virtual void visualizeParticleNeighbors(SPHParticle *target); + virtual void visualizeRandomParticlesNeighbors(); + + // Misc + virtual void initVisualization(); + virtual void prepNeighborSearchAfterSceneLoad(); + virtual void printPerformanceStats(); + virtual void sceneLoaded(); + + #if MFluidSolver_USE_PARTIO + virtual void exportBgeo(); + #endif + #if MFluidSolver_USE_OPENVDB + virtual void exportVDB(); + #endif + + protected: + // Initialization + virtual bool checkInited(); + bool inited; + + // Neighbor Search + virtual inline void prepNeighborSearch(); + virtual inline void runNeighborSearch(); + + // Helpers + inline void calculateDensity(SPHParticle *p); + virtual inline void calculateNonPressureForce(SPHParticle *p); + virtual inline void calculatePressureForce(SPHParticle *p); + inline void enforceBounds(SPHParticle *p); + virtual inline void visualizeParticle(SPHParticle *p); + + // SPH Required Objects + NeighborSearchType nSearchType; + KernelFunctions kernelFunctions; + NeighborSearch *nSearch; + std::vector _particles; + + // Configuration + float kStiffness; + float muViscosity; + float kernelRadius; + float mMass; + float dRestDensity; + float dtTimestep; + + unsigned int largestIndex; +}; + +#include "sphSolver.inline.hpp" + +#endif // MFLUIDSOLVER_SPHSOLVER_HPP_ diff --git a/src/fluidSolver/sphSolver/sphSolver.inline.hpp b/src/fluidSolver/sphSolver/sphSolver.inline.hpp new file mode 100644 index 00000000..ddb986ab --- /dev/null +++ b/src/fluidSolver/sphSolver/sphSolver.inline.hpp @@ -0,0 +1,174 @@ +// Copyright 2016 Robert Zhou +// +// sphSolver.inline.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_SPHSOLVER_INLINE_HPP_ +#define MFLUIDSOLVER_SPHSOLVER_INLINE_HPP_ + +#if MFluidSolver_USE_TBB + #define iter_all_sphparticles_start \ + tbb::parallel_for((size_t)0, _particles.size(), \ + [&](size_t i) { \ + SPHParticle &p = _particles[i]; + #define iter_all_sphparticles_end }); +#else + #define iter_all_sphparticles_start \ + for (size_t i = 0; i < _particles.size(); ++i) { \ + SPHParticle &p = _particles[i]; + #define iter_all_sphparticles_end } +#endif + +inline void SPHSolver::calculateDensity(SPHParticle *p) { + float densitySum = p->mass() * kernelFunctions.computePoly6(glm::vec3(0)); + if (p->neighbors()->size() > 0) { + for (SPHParticle *n : *(p->neighbors())) { + densitySum += n->mass() * + kernelFunctions.computePoly6(p->position() - n->position()); + } + } + p->setDensity(densitySum); +} + +inline void SPHSolver::calculateNonPressureForce(SPHParticle *p) { + glm::vec3 viscosityForce(0); + glm::vec3 gravityForce = glm::vec3(0, p->mass() * _gravity, 0); + + if (muViscosity > 0) { + for (SPHParticle *n : *(p->neighbors())) { + viscosityForce += glm::vec3( + n->mass() * (n->velocity() - p->velocity()) * + static_cast(kernelFunctions.computeViscousLaplacian( + p->position() - n->position() ))); + } + } + viscosityForce *= muViscosity / p->density() * p->mass(); + + p->setNonPressureForce(viscosityForce + gravityForce); +} + +inline void SPHSolver::calculatePressureForce(SPHParticle *p) { + glm::vec3 pressureForce(0); + float pDensity2 = p->density() * p->density(); + for (SPHParticle *n : *(p->neighbors())) { + pressureForce += n->mass() * ( + p->pressure() / pDensity2 + + n->pressure() / (n->density() * n->density())) * + kernelFunctions.computeSpikyGradient(p->position() - n->position()); + } + pressureForce *= -1 * p->mass(); + p->setPressureForce(pressureForce); +} + +inline void SPHSolver::enforceBounds(SPHParticle *p) { + glm::ivec3 violations; + if (!fluidContainer->intersects(p->position(), &violations)) { + glm::vec3 position = p->position(); + p->reverseVelocity(violations, 0.1f); // Bounce + + // TODO: Generalize + glm::vec3 scaleVec = fluidContainer->transform.scale(); + glm::vec3 minBounds = -0.5f * scaleVec; + glm::vec3 maxBounds = 0.5f * scaleVec; + float r = 0.f; + + if (violations.x < 0) position.x = minBounds.x; + else if (violations.x > 0) position.x = maxBounds.x; + + if (violations.y < 0) position.y = minBounds.y; + else if (violations.y > 0) position.y = maxBounds.y; + + if (violations.z < 0) position.z = minBounds.z; + else if (violations.z > 0) position.z = maxBounds.z; + + p->setPosition(position); + } +} + +inline void SPHSolver::prepNeighborSearch() { + if (nSearchType == NeighborSearchType::UniformGrid) { + nSearch->clear(); + for (SPHParticle &p : _particles) { + nSearch->addParticle(&p); + } + } else if (nSearchType == NeighborSearchType::IndexSortedUniformGrid || + nSearchType == NeighborSearchType::ZIndexSortedUniformGrid) { + IndexSortedUniformGridNeighborSearch *isugSearch = + static_cast(nSearch); + isugSearch->isuGrid->resetAndFillCells(true); + } else if ( + nSearchType == NeighborSearchType::IndexSortedUniformGridWithInsertion || + nSearchType == NeighborSearchType::ZIndexSortedUniformGridWithInsertion) { + IndexSortedUniformGridNeighborSearch *isugSearch = + static_cast(nSearch); + isugSearch->isuGrid->resetAndFillCells(false); // Use insertion sort + } + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_DEBUG + else { + std::cout << "DEBUG: No neighbor search prep being done" << std::endl; + } + #endif +} + +inline void SPHSolver::runNeighborSearch() { + iter_all_sphparticles_start + nSearch->findNeighbors(&p); + iter_all_sphparticles_end +} + +inline void SPHSolver::visualizeParticle(SPHParticle *p) { + switch (visualizationType) { + case FluidVisualizationType::Density: { + float densityNormalized = (p->density() - dRestDensity) / + visualizationMaxDensityDifference; + if (densityNormalized > 1.f) densityNormalized = 1.f; + if (densityNormalized < -1.f) densityNormalized = -1.f; + if (densityNormalized >= 0) { + p->color = (1.f - densityNormalized) * glm::vec3(1.f) + + densityNormalized * glm::vec3(0, 1, 0); + } else { + densityNormalized *= -1; + p->color = (1.f - densityNormalized) * glm::vec3(1.f) + + densityNormalized * glm::vec3(1, 0, 0); + } + break; + } + case FluidVisualizationType::Index: { + float indexNormalized = static_cast(p->index) / + static_cast(largestIndex); + p->color = (1.f - indexNormalized) * glm::vec3(1.f) + + indexNormalized * glm::vec3(1, 0, 0); + break; + } + case FluidVisualizationType::Pressure: { + float pressureNormalized = p->pressure() / visualizationMaxPressure; + if (pressureNormalized > 1.f) pressureNormalized = 1.f; + if (pressureNormalized < -1.f) pressureNormalized = -1.f; + if (pressureNormalized >= 0) { + p->color = (1.f - pressureNormalized) * glm::vec3(1.f) + + pressureNormalized * glm::vec3(0, 1, 0); + } else { + pressureNormalized *= -1; + p->color = (1.f - pressureNormalized) * glm::vec3(1.f) + + pressureNormalized * glm::vec3(1, 0, 0); + } + break; + } + case FluidVisualizationType::Velocity: { + float velocityNormalized = glm::length(p->velocity()) / + visualizationMaxVelocity; + if (velocityNormalized > 1.f) velocityNormalized = 1.f; + p->color = (1.f - velocityNormalized) * visualizationVelocityColor + + velocityNormalized * glm::vec3(1.f); + break; + } + case FluidVisualizationType::VelocityDir: { + p->color = glm::normalize(glm::abs(p->velocity())); + break; + } + default: + break; + } +} + +#endif // MFLUIDSOLVER_SPHSOLVER_INLINE_HPP_ diff --git a/src/fluidSolver/sphSolver/sphUniformGrid.cpp b/src/fluidSolver/sphSolver/sphUniformGrid.cpp new file mode 100644 index 00000000..fc24933b --- /dev/null +++ b/src/fluidSolver/sphSolver/sphUniformGrid.cpp @@ -0,0 +1,174 @@ +// Copyright 2016 Robert Zhou +// +// sphUniformGrid.cpp +// MFluidSolver + +#include "sphUniformGrid.hpp" + +#include +#include +#include + +#if MFluidSolver_USE_ASSERTS +#include +#endif +#if MFluidSolver_USE_OPENVDB +#include +#endif + +SPHUniformGrid::SPHUniformGrid( + const glm::vec3 &minBounds, const glm::vec3 &maxBounds, float cellSize) + : SPHGrid(minBounds, maxBounds, cellSize) { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + std:: cout << "INFO: Created grid with cell bounds: (" << + cellBounds.x << ", " << + cellBounds.y << ", " << cellBounds.z << ")" << std::endl; + #endif + + data.resize(numCells); +} + +SPHUniformGrid::~SPHUniformGrid() { +} + +void SPHUniformGrid::getNeighbors(SPHParticle *p) { + glm::ivec3 pC = getGridCoordinates(p->position()); + std::vector *neighbors = p->neighbors(); + neighbors->clear(); + + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_TRACE + std::cout << "TRACE: Getting neighbors for cell (" << + pC.x << ", " << pC.y << ", " << pC.z << ")" << std::endl; + #endif + + uint32_t index; + std::vector *vec; + + // Search neighboring grid cells + glm::ivec3 coords; + for (int i = -1; i < 2; ++i) { + for (int j = -1; j < 2; ++j) { + for (int k = -1; k < 2; ++k) { + coords = glm::ivec3(pC.x + i, pC.y + j, pC.z + k); + // Note: we check kernel radius and self in neighbor search + if (coords.x >= 0 && coords.y >= 0 && coords.z >= 0 && + coords.x < cellBounds.x && + coords.y < cellBounds.y && + coords.z < cellBounds.z) { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_TRACE + std::cout << "TRACE: Accessing cell (" << + coords.x << ", " << + coords.y << ", " << + coords.z << ")" << std::endl; + #endif + + index = getIndex(coords); + vec = &(data[index]); + + for (unsigned int l = 0; l < vec->size(); ++l) { + neighbors->push_back((*vec)[l]); + } // end for l + + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_TRACE + std::cout << "TRACE: Found " << + vec->size() << " neighbor(s)" << std::endl; + #endif + } + } // end for k + } // end for j + } // end for i +} + +void SPHUniformGrid::addParticle(SPHParticle *p) { + glm::ivec3 pC = getGridCoordinates(p->position()); + unsigned int pI = getIndex(pC); + p->index = pI; + data[pI].push_back(p); +} + +void SPHUniformGrid::updateParticle(SPHParticle *p) { + glm::ivec3 oldCoords = getGridCoordinates(p->oldPosition()); + glm::ivec3 newCoords = getGridCoordinates(p->position()); + if (oldCoords != newCoords) { + unsigned int oldIndex = getIndex(oldCoords); + unsigned int newIndex = getIndex(newCoords); + + std::vector *oldCell = &(data[oldIndex]); + auto oldIt = std::find(std::begin(*oldCell), std::end(*oldCell), p); + if (oldIt != std::end(*oldCell)) { + oldCell->erase(oldIt); + p->index = newIndex; + data[newIndex].push_back(p); + } else { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_ERROR + std::cerr << "ERROR: Could not find particle at (" << + p->position().x << ", " << + p->position().y << ", " << + p->position().z << ") based on old position (" << + p->oldPosition().x << ", " + << p->oldPosition().y << ", " + << p->oldPosition().z << "), index " << oldIndex << std::endl; + #endif + } + } +} + +void SPHUniformGrid::clear() { + for (unsigned int i = 0; i < numCells; ++i) { + data[i].clear(); + } +} + +void SPHUniformGrid::printDiagnostics() { + std::cout << "Grid bounds: (" << + cellBounds.x << ", " << + cellBounds.y << ", " << + cellBounds.z << ")" << std::endl; + std::cout << "Grid size: " << numCells << std::endl; + + for (unsigned int i = 0; i < cellBounds.x; ++i) { + for (unsigned int j = 0; j < cellBounds.y; ++j) { + for (unsigned int k = 0; k < cellBounds.z; ++k) { + std::printf( + "(%d, %d, %d) - %d\n", i, j, k, data[getIndex(i, j, k)].size()); + } + } + } +} + +#if MFluidSolver_USE_OPENVDB +void SPHUniformGrid::exportVDB( + const std::string &file, const std::string &gridName) { + // http://www.openvdb.org/documentation/doxygen/codeExamples.html + openvdb::FloatGrid::Ptr grid = + openvdb::FloatGrid::create(/*background value=*/0); + openvdb::FloatGrid::Accessor accessor = grid->getAccessor(); + openvdb::Coord coord; + + for (unsigned int i = 0; i < cellBounds.x; ++i) { + for (unsigned int j = 0; j < cellBounds.y; ++j) { + for (unsigned int k = 0; k < cellBounds.z; ++k) { + unsigned int numParticles = data[getIndex(i, j, k)].size(); + if (numParticles > 0) { + coord.reset(i, j, k); + accessor.setValue(coord, numParticles); + } + } + } + } + + grid->setName(gridName); + openvdb::io::File fileIO(file); + + // Add the grid pointer to a container. + openvdb::GridPtrVec grids; + grids.push_back(grid); + + fileIO.write(grids); + fileIO.close(); + + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + std::cout << "INFO: Exported grid VDB to: " << file << std::endl; + #endif +} +#endif diff --git a/src/fluidSolver/sphSolver/sphUniformGrid.hpp b/src/fluidSolver/sphSolver/sphUniformGrid.hpp new file mode 100644 index 00000000..9fcef15f --- /dev/null +++ b/src/fluidSolver/sphSolver/sphUniformGrid.hpp @@ -0,0 +1,35 @@ +// Copyright 2016 Robert Zhou +// +// sphUniformGrid.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_FLUIDSOLVER_SPHUNIFORMGRID_HPP_ +#define MFLUIDSOLVER_FLUIDSOLVER_SPHUNIFORMGRID_HPP_ + +#include +#include + +#include "sphGrid.hpp" + +// Grid type XZY (X elements together, then Z, then Y) +class SPHUniformGrid : public SPHGrid { + public: + SPHUniformGrid(const glm::vec3 &minBounds, const glm::vec3 &maxBounds, + float cellSize); + ~SPHUniformGrid(); + + virtual void getNeighbors(SPHParticle *p); + virtual void addParticle(SPHParticle *p); + virtual void updateParticle(SPHParticle *p); + virtual void clear(); + virtual void printDiagnostics(); + + #if MFluidSolver_USE_OPENVDB + virtual void exportVDB(const std::string &file, const std::string &gridName); + #endif + + private: + std::vector> data; +}; + +#endif // MFLUIDSOLVER_FLUIDSOLVER_SPHUNIFORMGRID_HPP_ diff --git a/src/fluidSolver/sphSolver/sphZIndexSortedUniformGrid.cpp b/src/fluidSolver/sphSolver/sphZIndexSortedUniformGrid.cpp new file mode 100644 index 00000000..ebe7f692 --- /dev/null +++ b/src/fluidSolver/sphSolver/sphZIndexSortedUniformGrid.cpp @@ -0,0 +1,24 @@ +// Copyright 2016 Robert Zhou +// +// sphZIndexSortedUniformGrid.cpp +// MFluidSolver + +#include + +#include "sphZIndexSortedUniformGrid.hpp" + +SPHZIndexSortedUniformGrid::SPHZIndexSortedUniformGrid( + const glm::vec3 &minBounds, const glm::vec3 &maxBounds, + float cellSize, std::vector *master) + : SPHIndexSortedUniformGrid(minBounds, maxBounds, cellSize, master) { + numCells = zCurve.initWithMax(cellBounds); + cells.resize(numCells); +} + +SPHZIndexSortedUniformGrid::~SPHZIndexSortedUniformGrid() { +} + +uint32_t SPHZIndexSortedUniformGrid::getIndex( + unsigned int x, unsigned int y, unsigned int z) { + return zCurve.getIndex(x, y, z); +} diff --git a/src/fluidSolver/sphSolver/sphZIndexSortedUniformGrid.hpp b/src/fluidSolver/sphSolver/sphZIndexSortedUniformGrid.hpp new file mode 100644 index 00000000..d4764106 --- /dev/null +++ b/src/fluidSolver/sphSolver/sphZIndexSortedUniformGrid.hpp @@ -0,0 +1,29 @@ +// Copyright 2016 Robert Zhou +// +// sphZIndexSortedUniformGrid.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_FLUIDSOLVER_SPHZINDEXSORTEDUNIFORMGRID_HPP_ +#define MFLUIDSOLVER_FLUIDSOLVER_SPHZINDEXSORTEDUNIFORMGRID_HPP_ + +#include + +#include "sphIndexSortedUniformGrid.hpp" +#include "fluidSolver/zCurve.hpp" + +// Grid type XZY (X elements together, then Z, then Y) +class SPHZIndexSortedUniformGrid : public SPHIndexSortedUniformGrid { + public: + SPHZIndexSortedUniformGrid( + const glm::vec3 &minBounds, const glm::vec3 &maxBounds, + float cellSize, std::vector *master); + ~SPHZIndexSortedUniformGrid(); + + virtual uint32_t getIndex( + unsigned int x, unsigned int y, unsigned int z); + + private: + ZCurve zCurve; +}; + +#endif // MFLUIDSOLVER_FLUIDSOLVER_SPHZINDEXSORTEDUNIFORMGRID_HPP_ diff --git a/src/fluidSolver/zCurve.cpp b/src/fluidSolver/zCurve.cpp new file mode 100644 index 00000000..e2e8cd5b --- /dev/null +++ b/src/fluidSolver/zCurve.cpp @@ -0,0 +1,121 @@ +// Copyright 2016 Robert Zhou +// +// zCurve.cpp +// MFluidSolver + +#include "zCurve.hpp" + +#include +#include + +#if MFluidSolver_USE_ASSERTS +#include +#endif + +#define LONG_COORD_LIMIT 2047 + +ZCurve::ZCurve() { +} + +uint32_t ZCurve::initWithMax(const glm::ivec3 &cellBounds) { + // Check sane bounds + if (cellBounds.x < 0 || cellBounds.y < 0 || cellBounds.z < 0) { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_FATAL + std::cerr << "FATAL: ZCurve bounds (" << + cellBounds.x << ", " << + cellBounds.y << ", " << + cellBounds.z << ") is negative" << std::endl; + #endif + throw ZCurveNegativeBoundsException(); + } + + if (cellBounds.x > LONG_COORD_LIMIT || + cellBounds.y > LONG_COORD_LIMIT || + cellBounds.z > LONG_COORD_LIMIT) { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_FATAL + std::cerr << "FATAL: ZCurve bounds (" << + cellBounds.x << ", " << + cellBounds.y << ", " << + cellBounds.z << ") exceeds limit " << + LONG_COORD_LIMIT << std::endl; + #endif + throw ZCurveTooLargeException(); + } + + // Cache Z-index if asked + #if MFluidSolver_ZCURVE_CACHING + maxIndex = calculateIndex(cellBounds); + #else + maxIndex = getIndex(cellBounds); + #endif + + #if MFluidSolver_ZCURVE_CACHING + for (unsigned int i = 0; i < cellBounds.x; ++i) { + cache.push_back(std::vector>()); + for (unsigned int j = 0; j < cellBounds.y; ++j) { + cache[i].push_back(std::vector()); + for (unsigned int k = 0; k < cellBounds.z; ++k) { + cache[i][j].push_back(calculateIndex(i, j, k)); + } + } + } + #endif + + return maxIndex; +} + +#if MFluidSolver_ZCURVE_CACHING +uint32_t ZCurve::getIndex(const glm::ivec3 &p) { + #if MFluidSolver_USE_ASSERTS + assert(p.x >= 0 && p.y >= 0 && p.z >= 0); + assert(p.x < cache.size() && + p.y < cache[0].size() && + p.z < cache[0][0].size()); + #endif + return cache[p.x][p.y][p.z]; +} + +uint32_t ZCurve::getIndex(unsigned int i, unsigned int j, unsigned int k) { + #if MFluidSolver_USE_ASSERTS + assert(i <= LONG_COORD_LIMIT && + j <= LONG_COORD_LIMIT && + k <= LONG_COORD_LIMIT); + #endif + + return getIndex(glm::ivec3(i, j, k)); +} + +uint32_t ZCurve::calculateIndex(const glm::ivec3 &p) { + #if MFluidSolver_USE_ASSERTS + assert(p.x >= 0 && p.y >= 0 && p.z >= 0); + #endif + + return calculateIndex(p.x, p.y, p.z); +} +#else +uint32_t ZCurve::getIndex(const glm::ivec3 &p) { + #if MFluidSolver_USE_ASSERTS + assert(p.x >= 0 && p.y >= 0 && p.z >= 0); + #endif + + return getIndex(p.x, p.y, p.z); +} +#endif + +#if MFluidSolver_ZCURVE_CACHING +uint32_t ZCurve::calculateIndex( + unsigned int i, unsigned int j, unsigned int k) { +#else +uint32_t ZCurve::getIndex(unsigned int i, unsigned int j, unsigned int k) { +#endif + return (splitBits(k) << 2) | (splitBits(j) << 1) | splitBits(i); +} + +// http://stackoverflow.com/questions/1024754/how-to-compute-a-3d-morton-number-interleave-the-bits-of-3-ints +uint32_t ZCurve::splitBits(uint32_t i) { + i = (i | (i << 16)) & 0x030000FF; + i = (i | (i << 8)) & 0x0300F00F; + i = (i | (i << 4)) & 0x030C30C3; + i = (i | (i << 2)) & 0x09249249; + return i; +} diff --git a/src/fluidSolver/zCurve.hpp b/src/fluidSolver/zCurve.hpp new file mode 100644 index 00000000..9e61ae46 --- /dev/null +++ b/src/fluidSolver/zCurve.hpp @@ -0,0 +1,61 @@ +// Copyright 2016 Robert Zhou +// +// zCurve.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_FLUIDSOLVER_ZCURVE_HPP_ +#define MFLUIDSOLVER_FLUIDSOLVER_ZCURVE_HPP_ + +#include +#include + +#include "MFluidSolverConfig.hpp" + +#if MFluidSolver_ZCURVE_CACHING +#include +#endif + +#if MFluidSolver_ZCURVE_CACHING +struct simpleIvec3Comparator { + bool operator()(const glm::ivec3 &a, const glm::ivec3 &b) const { + if (a.x != b.x) return a.x < b.x; + if (a.y != b.y) return a.y < b.y; + return a.z < b.z; + } +}; +#endif + +struct ZCurveNegativeBoundsException : std::exception { + const char *what() const noexcept { + return "ZCurve cannot be constructed with negative cell bounds.\n"; + }; +}; + +struct ZCurveTooLargeException : std::exception { + const char *what() const noexcept { + return "ZCurve cannot use such high indices"; + }; +}; + +class ZCurve { + public: + ZCurve(); + uint32_t initWithMax(const glm::ivec3 &cellBounds); + uint32_t getIndex(const glm::ivec3 &p); + uint32_t getIndex(unsigned int i, unsigned int j, unsigned int k); + uint32_t splitBits(uint32_t i); + + #if MFluidSolver_ZCURVE_CACHING + uint32_t calculateIndex(const glm::ivec3 &p); + uint32_t calculateIndex(unsigned int i, unsigned int j, unsigned int k); + #endif + + private: + uint32_t maxIndex; + + #if MFluidSolver_ZCURVE_CACHING + std::vector>> cache; + #endif +}; + +#endif // MFLUIDSOLVER_FLUIDSOLVER_ZCURVE_HPP_ diff --git a/src/geom/cube.cpp b/src/geom/cube.cpp new file mode 100644 index 00000000..67f3db6c --- /dev/null +++ b/src/geom/cube.cpp @@ -0,0 +1,171 @@ +// Copyright 2016 Robert Zhou +// +// cube.cpp +// MFluidSolver + +#include "cube.hpp" + +#include +#include + +#define GLM_FORCE_RADIANS +#include + +Cube::Cube(const glm::vec3 &color) : _color(color) { + const unsigned int IDX_COUNT = 24; + const unsigned int VERT_COUNT = 8; + + // Vertices + glm::vec3 vert_pos[VERT_COUNT]; + vert_pos[0] = glm::vec3(0.5f, 0.5f, 0.5f); + vert_pos[1] = glm::vec3(-0.5f, 0.5f, 0.5f); + vert_pos[2] = glm::vec3(0.5f, -0.5f, 0.5f); + vert_pos[3] = glm::vec3(-0.5f, -0.5f, 0.5f); + vert_pos[4] = glm::vec3(0.5f, 0.5f, -0.5f); + vert_pos[5] = glm::vec3(-0.5f, 0.5f, -0.5f); + vert_pos[6] = glm::vec3(0.5f, -0.5f, -0.5f); + vert_pos[7] = glm::vec3(-0.5f, -0.5f, -0.5f); + + // Indices + GLuint idx[IDX_COUNT]; + idx[0] = 0; idx[1] = 1; + idx[2] = 0; idx[3] = 2; + idx[4] = 0; idx[5] = 4; + idx[6] = 1; idx[7] = 3; + idx[8] = 1; idx[9] = 5; + idx[10] = 2; idx[11] = 3; + idx[12] = 2; idx[13] = 6; + idx[14] = 3; idx[15] = 7; + idx[16] = 4; idx[17] = 5; + idx[18] = 4; idx[19] = 6; + idx[20] = 5; idx[21] = 7; + idx[22] = 6; idx[23] = 7; + + // Color + glm::vec3 vert_col[VERT_COUNT]; + for (unsigned int i = 0; i < VERT_COUNT; ++i) { + vert_col[i] = _color; + } + + _idxCount = IDX_COUNT; + + // Bind + glGenBuffers(1, &vertexIndexArrBufferID); + glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, vertexIndexArrBufferID); + glBufferData(GL_ELEMENT_ARRAY_BUFFER, IDX_COUNT * sizeof(GLuint), idx, + GL_STATIC_DRAW); + + glGenBuffers(1, &vertexPositionArrBufferID); + glBindBuffer(GL_ARRAY_BUFFER, vertexPositionArrBufferID); + glBufferData(GL_ARRAY_BUFFER, VERT_COUNT * sizeof(glm::vec3), vert_pos, + GL_STATIC_DRAW); + + glGenBuffers(1, &vertexColorArrBufferID); + glBindBuffer(GL_ARRAY_BUFFER, vertexColorArrBufferID); + glBufferData(GL_ARRAY_BUFFER, VERT_COUNT * sizeof(glm::vec3), vert_col, + GL_STATIC_DRAW); +} + +GLenum Cube::drawMode() { + return GL_LINES; +} + +bool Cube::intersects(const glm::vec3 &point) const { + glm::vec3 localPos = glm::vec3(transform.invT() * glm::vec4(point, 1.f)); + if (localPos.x >= -0.5f && localPos.x <= 0.5f && + localPos.y >= -0.5f && localPos.y <= 0.5f && + localPos.z >= -0.5f && localPos.z <= 0.5f) + return true; + return false; +} + +void Cube::spawnParticlesInVolume( + FluidSolver *solver, ParticleSpawnMethod spawnMethod) const { + // Store old particle count for number of particles we added and max particle check + unsigned int oldParticleCount = solver->numParticles(); + if (oldParticleCount == solver->maxParticles()) { + std::cout << "INFO: Reached max number of particles (" << + solver->maxParticles() << ")!" << std::endl; + return; + } + + unsigned int particlesLeft = solver->maxParticles() - oldParticleCount; + glm::vec3 minBound, maxBound; + getBoundsByTransformedMinMax( + glm::vec3(-0.5f), glm::vec3(0.5f), &minBound, &maxBound); + + float particleSeparation = solver->particleSeparation(); + unsigned int count = 0; + + // Calculate padding due to particle separation + glm::vec3 padding = maxBound - minBound; + glm::ivec3 numParticlesOnAxis = (glm::ivec3) (padding / particleSeparation); + padding -= (glm::vec3)numParticlesOnAxis * particleSeparation; + minBound += padding; + + // Create particles + if (spawnMethod == ParticleSpawnMethod::Uniform) { + for (float i = minBound.x; i <= maxBound.x; i += particleSeparation) { + for (float j = minBound.y; j <= maxBound.y; j += particleSeparation) { + for (float k = minBound.z; k <= maxBound.z; k += particleSeparation) { + solver->addParticleAt(glm::vec3(i, j, k)); + ++count; + if (count > particlesLeft) break; + } + if (count > particlesLeft) break; + } + if (count > particlesLeft) break; + } + } else if (spawnMethod == ParticleSpawnMethod::Jittered) { + std::random_device rd; + std::mt19937 gen(rd()); + std::uniform_real_distribution unifDist(0, particleSeparation); + const float halfSeparation = particleSeparation * 0.5f; + + for (float i = minBound.x; i <= maxBound.x; i += particleSeparation) { + for (float j = minBound.y; j <= maxBound.y; j += particleSeparation) { + for (float k = minBound.z; k <= maxBound.z; k += particleSeparation) { + const float xOffset = unifDist(gen) - halfSeparation; + const float yOffset = unifDist(gen) - halfSeparation; + const float zOffset = unifDist(gen) - halfSeparation; + const float factor = 0.27f; + solver->addParticleAt(glm::vec3( + i + xOffset * factor, j + yOffset * factor, k + zOffset * factor)); + ++count; + if (count > particlesLeft) break; + } + if (count > particlesLeft) break; + } + if (count > particlesLeft) break; + } + } + + if (solver->numParticles() == solver->maxParticles()) { + std::cout << "INFO: Reached max number of particles (" << + solver->maxParticles() << ")!" << std::endl; + } + + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + std::cout << "INFO: Cube (" << name << ") seeded " << + (solver->numParticles() - oldParticleCount) << " particles" << std::endl; + #endif +} + +bool Cube::intersects(const glm::vec3 &point, glm::ivec3 *violations) const { + glm::vec3 localPos = glm::vec3(transform.invT() * glm::vec4(point, 1.f)); + + // Check each axis locally + violations->x = 0; + if (localPos.x <= -0.5f) violations->x = -1; + else if (localPos.x > 0.5f) violations->x = 1; + + violations->y = 0; + if (localPos.y <= -0.5f) violations->y = -1; + else if (localPos.y > 0.5f) violations->y = 1; + + violations->z = 0; + if (localPos.z <= -0.5f) violations->z = -1; + else if (localPos.z > 0.5f) violations->z = 1; + + return violations->x == 0 && violations->y == 0 && violations->z == 0; +} diff --git a/src/geom/cube.hpp b/src/geom/cube.hpp new file mode 100644 index 00000000..9227b064 --- /dev/null +++ b/src/geom/cube.hpp @@ -0,0 +1,32 @@ +// Copyright 2016 Robert Zhou +// +// cube.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_GEOM_CUBE_HPP_ +#define MFLUIDSOLVER_GEOM_CUBE_HPP_ + +#include "MFluidSolverConfig.hpp" + +#include "geom.hpp" + +class Cube : public Geometry { + public: + Cube() + : Cube(glm::vec3( + MFluidSolver_DEFAULT_GEOMETRY_COLOR_R, + MFluidSolver_DEFAULT_GEOMETRY_COLOR_G, + MFluidSolver_DEFAULT_GEOMETRY_COLOR_B)) {} + explicit Cube(const glm::vec3 &color); + + virtual GLenum drawMode(); + virtual bool intersects(const glm::vec3 &point) const; + virtual bool intersects(const glm::vec3 &point, glm::ivec3 *violations) const; + virtual void spawnParticlesInVolume( + FluidSolver *solver, ParticleSpawnMethod spawnMethod) const; + + private: + glm::vec3 _color; +}; + +#endif // MFLUIDSOLVER_GEOM_CUBE_HPP_ diff --git a/src/geom/drawable.cpp b/src/geom/drawable.cpp new file mode 100644 index 00000000..b953d8f6 --- /dev/null +++ b/src/geom/drawable.cpp @@ -0,0 +1,45 @@ +// Copyright 2016 Robert Zhou +// +// drawable.cpp +// MFluidSolver + +#include "drawable.hpp" + +Drawable::Drawable() + : vertexIndexArrBufferID(-1), + vertexColorArrBufferID(-1), + vertexPositionArrBufferID(-1) { +} + +Drawable::~Drawable() { + if (vertexIndexArrBufferID != -1) glDeleteBuffers(1, &vertexIndexArrBufferID); + if (vertexColorArrBufferID != -1) glDeleteBuffers(1, &vertexColorArrBufferID); + if (vertexPositionArrBufferID != -1) + glDeleteBuffers(1, &vertexPositionArrBufferID); +} + +bool Drawable::bindIndexBuffer() { + if (vertexIndexArrBufferID == -1) return false; + glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, vertexIndexArrBufferID); + return true; +} + +bool Drawable::bindColorBuffer() { + if (vertexColorArrBufferID == -1) return false; + glBindBuffer(GL_ARRAY_BUFFER, vertexColorArrBufferID); + return true; +} + +bool Drawable::bindPositionBuffer() { + if (vertexPositionArrBufferID == -1) return false; + glBindBuffer(GL_ARRAY_BUFFER, vertexPositionArrBufferID); + return true; +} + +GLenum Drawable::drawMode() { + return GL_TRIANGLES; +} + +GLenum Drawable::idxCount() { + return _idxCount; +} diff --git a/src/geom/drawable.hpp b/src/geom/drawable.hpp new file mode 100644 index 00000000..432d3572 --- /dev/null +++ b/src/geom/drawable.hpp @@ -0,0 +1,28 @@ +// Copyright 2016 Robert Zhou +// +// drawable.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_GEOM_DRAWABLE_HPP_ +#define MFLUIDSOLVER_GEOM_DRAWABLE_HPP_ + +#include + +class Drawable { + public: + Drawable(); + ~Drawable(); + bool bindIndexBuffer(); + bool bindColorBuffer(); + bool bindPositionBuffer(); + virtual GLenum drawMode(); + GLenum idxCount(); + + protected: + GLuint vertexIndexArrBufferID; + GLuint vertexColorArrBufferID; + GLuint vertexPositionArrBufferID; + int _idxCount; +}; + +#endif // MFLUIDSOLVER_GEOM_DRAWABLE_HPP_ diff --git a/src/geom/geom.cpp b/src/geom/geom.cpp index c4438fdc..d88c1ff1 100644 --- a/src/geom/geom.cpp +++ b/src/geom/geom.cpp @@ -1,5 +1,46 @@ +// Copyright 2016 Robert Zhou // // geom.cpp -// Thanda +// MFluidSolver + +#include +#include #include "geom.hpp" + +void Geometry::getBoundsByTransformedMinMax( + const glm::vec3 &min, const glm::vec3 &max, + glm::vec3 *outMin, glm::vec3 *outMax) const { + // Transforms min and max vertices, + // then gets min and max after they're transformed + std::vector points; + points.push_back(glm::vec3(min.x, min.y, min.z)); + points.push_back(glm::vec3(min.x, min.y, max.z)); + points.push_back(glm::vec3(min.x, max.y, min.z)); + points.push_back(glm::vec3(min.x, max.y, max.z)); + points.push_back(glm::vec3(max.x, min.y, min.z)); + points.push_back(glm::vec3(max.x, min.y, max.z)); + points.push_back(glm::vec3(max.x, max.y, min.z)); + points.push_back(glm::vec3(max.x, max.y, max.z)); + + for (size_t i = 0; i < points.size(); ++i) { + points[i] = glm::vec3(transform.T() * glm::vec4(points[i], 1.0)); + } + + for (size_t i = 0; i < 3; ++i) { + (*outMin)[i] = std::numeric_limits::infinity(); + (*outMax)[i] = -std::numeric_limits::infinity(); + } + + for (std::vector::iterator it = points.begin(); + it != points.end(); ++it) { + for (int i = 0; i < 3; ++i) { + if ((*it)[i] < (*outMin)[i]) { + (*outMin)[i] = (*it)[i]; + } + if ((*it)[i] > (*outMax)[i]) { + (*outMax)[i] = (*it)[i]; + } + } + } +} diff --git a/src/geom/geom.hpp b/src/geom/geom.hpp index 91e658f7..380017e3 100644 --- a/src/geom/geom.hpp +++ b/src/geom/geom.hpp @@ -1,9 +1,38 @@ +// Copyright 2016 Robert Zhou // // geom.hpp -// Thanda +// MFluidSolver -#ifndef geom_hpp -#define geom_hpp +#ifndef MFLUIDSOLVER_GEOM_HPP_ +#define MFLUIDSOLVER_GEOM_HPP_ +class Geometry; +class FluidSolver; -#endif /* geom_hpp */ +#include +#include + +#include + +#include "fluidSolver/fluidSolver.hpp" +#include "drawable.hpp" +#include "transform.hpp" + +enum ParticleSpawnMethod {Uniform, Jittered, PoissonDisk}; + +class Geometry : public Drawable { + public: + void getBoundsByTransformedMinMax( + const glm::vec3 &min, const glm::vec3 &max, + glm::vec3 *outMin, glm::vec3 *outMax) const; + virtual bool intersects(const glm::vec3 &position) const = 0; + virtual bool intersects(const glm::vec3 &position, + glm::ivec3 *violations) const = 0; + virtual void spawnParticlesInVolume(FluidSolver *solver, + ParticleSpawnMethod spawnMethod) const = 0; + + std::string name; + Transform transform; +}; + +#endif // MFLUIDSOLVER_GEOM_HPP_ diff --git a/src/geom/transform.cpp b/src/geom/transform.cpp new file mode 100644 index 00000000..3d1c976a --- /dev/null +++ b/src/geom/transform.cpp @@ -0,0 +1,66 @@ +// Copyright 2016 Robert Zhou +// +// transform.cpp +// MFluidSolver + +#include "transform.hpp" + +#define GLM_FORCE_RADIANS +#include + +Transform::Transform(const glm::vec3 &translation, + const glm::vec3 &rotation, + const glm::vec3 &scale) + : _translation(translation), + _rotation(rotation), + _scale(scale) { + calculateMatrices(); +} + +void Transform::calculateMatrices() { + worldTransform = glm::translate(glm::mat4(1.0f), _translation) * + glm::rotate(glm::mat4(1.0f), glm::radians(_rotation.x), glm::vec3(1, 0, 0)) * + glm::rotate(glm::mat4(1.0f), glm::radians(_rotation.y), glm::vec3(0, 1, 0)) * + glm::rotate(glm::mat4(1.0f), glm::radians(_rotation.z), glm::vec3(0, 0, 1)) * + glm::scale(glm::mat4(1.0f), _scale); + invWorldTransform = glm::inverse(worldTransform); +} + +void Transform::setTransform( + const glm::vec3 &translation, + const glm::vec3 &rotation, + const glm::vec3 &scale) { + _translation = translation; + _rotation = rotation; + _scale = scale; + calculateMatrices(); +} + +void Transform::setTranslation(const glm::vec3 &translation) { + _translation = translation; + calculateMatrices(); +} +void Transform::setRotation(const glm::vec3 &rotation) { + _rotation = rotation; + calculateMatrices(); +} +void Transform::setScale(const glm::vec3 &scale) { + _scale = scale; + calculateMatrices(); +} + +const glm::mat4 &Transform::T() const { + return worldTransform; +} +const glm::mat4 &Transform::invT() const { + return invWorldTransform; +} +const glm::vec3 &Transform::position() const { + return _translation; +} +const glm::vec3 &Transform::rotation() const { + return _rotation; +} +const glm::vec3 &Transform::scale() const { + return _scale; +} diff --git a/src/geom/transform.hpp b/src/geom/transform.hpp new file mode 100644 index 00000000..dc4d0675 --- /dev/null +++ b/src/geom/transform.hpp @@ -0,0 +1,42 @@ +// Copyright 2016 Robert Zhou +// +// transform.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_GEOM_TRANSFORM_HPP_ +#define MFLUIDSOLVER_GEOM_TRANSFORM_HPP_ + +#include + +class Transform { + public: + Transform() : Transform(glm::vec3(0), glm::vec3(0), glm::vec3(0)) {} + Transform(const glm::vec3 &translation, + const glm::vec3 &rotation, + const glm::vec3 &scale); + + void setTransform(const glm::vec3 &translation, + const glm::vec3 &rotation, + const glm::vec3 &scale); + void setTranslation(const glm::vec3 &translation); + void setRotation(const glm::vec3 &rotation); + void setScale(const glm::vec3 &scale); + + const glm::mat4 &T() const; + const glm::mat4 &invT() const; + const glm::vec3 &position() const; + const glm::vec3 &rotation() const; + const glm::vec3 &scale() const; + + private: + void calculateMatrices(); + + glm::vec3 _translation; + glm::vec3 _rotation; + glm::vec3 _scale; + + glm::mat4 worldTransform; + glm::mat4 invWorldTransform; +}; + +#endif // MFLUIDSOLVER_GEOM_TRANSFORM_HPP_ diff --git a/src/main.cpp b/src/main.cpp index 2a00d1e4..e2126546 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -1,11 +1,170 @@ -#include "main.hpp" +// Copyright 2016 Robert Zhou +// +// main.cpp +// MFluidSolver -using namespace std; +// OpenGL Includes +#include +#include +// Misc C Includes +#include +#include +#include -int main() -{ +#include "MFluidSolverConfig.hpp" - return 0; -} +// Optional C Includes +#if MFluidSolver_USE_OPENVDB +#include +#endif +#if MFluidSolver_USE_TBB +#include +#endif + +// Misc C++ Includes +#include + +// Project Includes +#include "viewer/input.hpp" +#include "viewer/particleShaderProgram.hpp" + +#define MFluidSolver_MAIN_CATCH_EXCEPTIONS 0 + +int main() { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + // Print version info + std::cout << "INFO: Version " << + MFluidSolver_VERSION_MAJOR << "." << + MFluidSolver_VERSION_MINOR << std::endl; + #endif + + // Initialize GLFW + if (!glfwInit()) { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_ERROR + std::cerr << "FATAL: Failed to initiaize GLFW" << std::endl; + #endif + + getchar(); // Wait for key before quit + return -1; + } + + // Initialize non-OpenGL + std::srand(std::time(NULL)); + #if MFluidSolver_USE_OPENVDB + openvdb::initialize(); + #endif + #if MFluidSolver_USE_TBB + tbb::task_scheduler_init init(tbb::task_scheduler_init::automatic); + #endif + + // Set JSON file location + std::string configJSON = MFluidSolver_DEFAULT_CONFIG_FILE; + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + std::cout << "INFO: Loading config file: " << configJSON << std::endl; + #endif + + // Open JSON file + Json::Reader reader; + Json::Value root; + std::ifstream sceneStream(configJSON, std::ifstream::binary); + bool success = reader.parse(sceneStream, root, false); + if (!success) { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_ERROR + std::cerr << "ERROR: Failed to parse config file " << + configJSON.c_str() << std::endl; + #endif + } + + // Parse JSON + std::string sceneJSON = root.get( + "sceneJSON", MFluidSolver_DEFAULT_SCENE_FILE).asString(); + std::string wireVShader = root.get( + "wireVShader", MFluidSolver_DEFAULT_WIRE_VERT_FILE).asString(); + std::string wireFShader = root.get( + "wireFShader", MFluidSolver_DEFAULT_WIRE_FRAG_FILE).asString(); + std::string particleVShader = root.get( + "particleVShader", MFluidSolver_DEFAULT_PARTICLE_VERT_FILE).asString(); + std::string particleFShader = root.get( + "particleFShader", MFluidSolver_DEFAULT_PARTICLE_FRAG_FILE).asString(); + std::string particleTexture = root.get( + "particleTexture", MFluidSolver_DEFAULT_PARTICLE_TEX_FILE).asString(); + bool autoRender = root.get( + "autoRender", MFluidSolver_DEFAULT_AUTORENDER).asBool(); + unsigned int renderSkip = root.get( + "renderSkip", MFluidSolver_DEFAULT_RENDERSKIP).asInt(); + // Initialize project objects + Viewer viewer; + Input::viewer = &viewer; + try { + viewer.init(); + } catch (std::exception &e) { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_FATAL + std::cerr << "FATAL: " << e.what() << std::endl; + #endif + return -1; + } + + // Print OpenGL info + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + std::cout << "INFO: OpenGL Version " << glGetString(GL_VERSION) << std::endl; + #endif + + // Set a few settings/modes in OpenGL rendering + glEnable(GL_DEPTH_TEST); + glDepthFunc(GL_LESS); + glEnable(GL_LINE_SMOOTH); + glEnable(GL_POLYGON_SMOOTH); + glHint(GL_LINE_SMOOTH_HINT, GL_NICEST); + glHint(GL_POLYGON_SMOOTH_HINT, GL_NICEST); + + // Set the size with which points should be rendered + glPointSize(MFluidSolver_DEFAULT_POINT_SIZE); + + // Set background color + glClearColor(MFluidSolver_DEFAULT_BG_COLOR_R, + MFluidSolver_DEFAULT_BG_COLOR_G, + MFluidSolver_DEFAULT_BG_COLOR_B, + 0.0f); + + // Bind vertex array object + GLuint vaoID; + glGenVertexArrays(1, &vaoID); + glBindVertexArray(vaoID); + + // After vertex array object init + viewer.wireShader = new ShaderProgram(wireVShader, wireFShader); + viewer.particleShader = + new ParticleShaderProgram(&(viewer.scene.solver), + particleVShader, particleFShader, particleTexture); + viewer.scene.solver.loadConfig(configJSON); + viewer.configureScreenshot(autoRender, renderSkip); + + // Run! Catch and print exceptions as necessary + int returnCode = 0; + #if MFluidSolver_MAIN_CATCH_EXCEPTIONS + try { + #endif + viewer.scene.loadJSON(sceneJSON); + viewer.run(); + #if MFluidSolver_MAIN_CATCH_EXCEPTIONS + } catch (std::exception &e) { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_FATAL + std::cerr << "FATAL: " << e.what() << std::endl; + #endif + returnCode = -1; + } + #endif + + // Performance + #if MFluidSolver_RECORD_PERFORMANCE + viewer.scene.solver.endSimulation(); + viewer.scene.solver.printPerformanceStats(); + #endif + + // Cleanup + glDeleteVertexArrays(1, &vaoID); + + return returnCode; +} diff --git a/src/main.hpp b/src/main.hpp deleted file mode 100644 index fe1f64d5..00000000 --- a/src/main.hpp +++ /dev/null @@ -1,9 +0,0 @@ -// -// main.hpp -// Thanda - -#ifndef main_hpp -#define main_hpp - -#endif /* main_hpp */ - diff --git a/src/scene/camera.cpp b/src/scene/camera.cpp new file mode 100644 index 00000000..81653732 --- /dev/null +++ b/src/scene/camera.cpp @@ -0,0 +1,149 @@ +// Copyright 2016 Robert Zhou +// +// camera.cpp +// MFluidSolver + +#include "camera.hpp" + +#include + +#define GLM_FORCE_RADIANS +#include +#include +#include + +Camera::Camera(unsigned int width, unsigned int height, + const glm::vec3 &eye, const glm::vec3 &ref, + const glm::vec3 &worldUp) + : nearClip(0.1f), farClip(1000.f), + _width(width), _height(height), + _eye(eye), _ref(ref), _worldUp(worldUp), + _fovy(45.f), arcballZoom(1.f), + panSpeed(0.01f), zoomSpeed(0.005f) { + arcballRotationMat = glm::mat4(1.f); + recomputeAttributes(); + recomputeLocalAxes(); + recomputeEyeAndRef(); +} + +glm::mat4 Camera::getViewProjection() { + glm::mat4 perpMat = + glm::perspective(glm::radians(_fovy), _aspect, nearClip, farClip); + glm::mat4 lookMat = + glm::lookAt( + glm::vec3(arcballEye), glm::vec3(arcballRef), glm::vec3(arcballUp)); + return perpMat * lookMat; +} + +void Camera::setEyeRef(const glm::vec3 &eye, const glm::vec3 &ref) { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + std::cout << "INFO: Set camera to eye (" << + eye.x << ", " << + eye.y << ", " << + eye.z << ") and ref (" << + ref.x << ", " << + ref.y << ", " << + ref.z << ")" << std::endl; + #endif + + _eye = eye; + _ref = ref; + arcballZoom = 1.f; + arcballPan = glm::vec4(0); + arcballRotationMat = glm::mat4(1.f); + + recomputeAttributes(); + recomputeLocalAxes(); + recomputeEyeAndRef(); +} + +void Camera::recomputeAttributes() { + // Vectors + _look = glm::normalize(_ref - _eye); + _right = glm::normalize(glm::cross(_look, _worldUp)); + _up = glm::normalize(glm::cross(_right, _look)); + + // Other + _aspect = static_cast(_width) / static_cast(_height); + + // Arcball + centerX = _width / 2.f; + centerY = _height / 2.f; + arcballRadius = _width * 0.30; // Arbitrary radius +} + +void Camera::recomputeLocalAxes() { + arcballLook = glm::normalize(arcballRotationMat * glm::vec4(_look, 0)); + arcballUp = glm::normalize(arcballRotationMat * glm::vec4(_up, 0)); + arcballRight = glm::normalize(arcballRotationMat * glm::vec4(_right, 0)); +} + +void Camera::recomputeEyeAndRef() { + arcballEye = arcballRotationMat * glm::vec4(_eye, 0); + arcballEye += arcballPan; + arcballRef = glm::vec4(_ref, 1.f) + arcballPan; + + // Apply zoom + glm::vec4 difference = arcballEye - arcballRef; + difference *= arcballZoom; + arcballEye = arcballRef + difference; +} + +void Camera::arcball(const glm::dvec2 &p1, const glm::dvec2 &p2) { + glm::vec3 pt1 = computeSpherePoint(p1); + glm::vec3 pt2 = computeSpherePoint(p2); + + glm::vec3 ptCross = glm::cross(pt1, pt2); + glm::quat dq = glm::quat(glm::dot(pt1, pt2), ptCross.x, ptCross.y, ptCross.z); + + arcballRotationMat *= glm::mat4_cast(dq * glm::quat(1.f, 0, 0, 0)); + + recomputeLocalAxes(); + recomputeEyeAndRef(); +} + +glm::vec3 Camera::computeSpherePoint(const glm::dvec2 &p) const { + glm::vec3 p2(0); + p2.x = (p.x - centerX) / arcballRadius; + p2.y = (p.y - centerY) / arcballRadius; + + float r = (p2.x * p2.x) + (p2.y * p2.y); + if (r > 1) { + p2 = glm::normalize(p2); + } else { + p2.x *= -1; // Reverse x position + p2.z = sqrt(1 - r); + } + + return p2; +} + +void Camera::pan(const glm::dvec2 &p1, const glm::dvec2 &p2) { + float dx = p2.x - p1.x; + float dy = p2.y - p1.y; + + arcballPan -= dx * panSpeed * arcballRight; + arcballPan += dy * panSpeed * arcballUp; + + recomputeEyeAndRef(); +} + +void Camera::zoom(double delta) { + // A typical mouse has 15 * 8 delta (eigth-degrees) per scroll wheel tick + arcballZoom += delta * zoomSpeed; + + // Prevent zooming past ref + if (arcballZoom < 0) { + arcballZoom = 0; + } + + recomputeEyeAndRef(); +} + +const glm::vec3 Camera::right() const { + return glm::vec3(arcballRight); +} + +const glm::vec3 Camera::up() const { + return glm::vec3(arcballUp); +} diff --git a/src/scene/camera.hpp b/src/scene/camera.hpp new file mode 100644 index 00000000..1ca0b5dd --- /dev/null +++ b/src/scene/camera.hpp @@ -0,0 +1,81 @@ +// Copyright 2016 Robert Zhou +// +// camera.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_SCENE_CAMERA_HPP_ +#define MFLUIDSOLVER_SCENE_CAMERA_HPP_ + +#include + +#include "MFluidSolverConfig.hpp" + +class Camera { + public: + Camera() : Camera(1024, 768) {} + Camera(unsigned int width, unsigned int height) : Camera( + width, height, + glm::vec3(0, 0, 12), glm::vec3(0, 0, 0), glm::vec3(0, 1, 0)) {} + Camera(unsigned int width, unsigned int height, + const glm::vec3 &eye, const glm::vec3 &ref, const glm::vec3 &worldUp); + + // Matrices + glm::mat4 getViewProjection(); + + void setEyeRef(const glm::vec3 &eye, const glm::vec3 &ref); + void recomputeAttributes(); + void recomputeLocalAxes(); + void recomputeEyeAndRef(); + + // Arcball + // Rotation interface for the mouse + void arcball(const glm::dvec2 &p1, const glm::dvec2 &p2); + // Panning interface for the mouse + void pan(const glm::dvec2 &p1, const glm::dvec2 &p2); + // Zooming interface for the scroll wheel + void zoom(double delta); + + // Properties + float nearClip, + farClip; + + const glm::vec3 right() const; + const glm::vec3 up() const; + + private: + // Dimensions + unsigned int _width, + _height; + + // Orientation vars + glm::vec3 _eye, + _ref, + _look, + _up, + _right, + _worldUp, + _V, + _H; + + // Arcball + float centerX, + centerY, + arcballRadius, + panSpeed, + zoomSpeed; + glm::mat4 arcballRotationMat; + glm::vec4 arcballEye, + arcballRef, + arcballLook, + arcballUp, + arcballRight; + glm::vec4 arcballPan; + float arcballZoom; + glm::vec3 computeSpherePoint(const glm::dvec2 &p) const; + + // Other properties + float _fovy, + _aspect; +}; + +#endif // MFLUIDSOLVER_SCENE_CAMERA_HPP_ diff --git a/src/scene/scene.cpp b/src/scene/scene.cpp index 32300ca9..4d718d7d 100644 --- a/src/scene/scene.cpp +++ b/src/scene/scene.cpp @@ -1,5 +1,139 @@ +// Copyright 2016 Robert Zhou // // scene.cpp -// Thanda +// MFluidSolver #include "scene.hpp" + +#include + +#include +#include +#include + +#include "geom/cube.hpp" +#include "utils.hpp" + +Scene::Scene() + : spawnMethod(ParticleSpawnMethod::Jittered) { +} + +Scene::~Scene() { + for (Geometry *g : objects) { + delete g; + } +} + +void Scene::loadJSON(const std::string &file) { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + std::cout << "INFO: Loading scene file: " << file << std::endl; + #endif + + // Read JSON file + Json::Reader reader; + Json::Value root; + std::ifstream sceneStream(file, std::ifstream::binary); + + bool success = reader.parse(sceneStream, root, false); + if (!success) { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_FATAL + std::cerr << "FATAL: Failed to parse scene file " << std::endl; + #endif + + throw InvalidSceneException(); + } + + // Load container info + glm::vec3 containerDim; + containerDim.x = root["containerDim"].get( + "scaleX", MFluidSolver_DEFAULT_SCENE_CONTAINER_SCALEX).asFloat(); + containerDim.y = root["containerDim"].get( + "scaleY", MFluidSolver_DEFAULT_SCENE_CONTAINER_SCALEY).asFloat(); + containerDim.z = root["containerDim"].get( + "scaleZ", MFluidSolver_DEFAULT_SCENE_CONTAINER_SCALEZ).asFloat(); + + // Load source info + glm::vec3 particleDim; + particleDim.x = root["particleDim"].get( + "scaleX", MFluidSolver_DEFAULT_SCENE_SOURCE_SCALEX).asFloat(); + particleDim.y = root["particleDim"].get( + "scaleY", MFluidSolver_DEFAULT_SCENE_SOURCE_SCALEY).asFloat(); + particleDim.z = root["particleDim"].get( + "scaleZ", MFluidSolver_DEFAULT_SCENE_SOURCE_SCALEZ).asFloat(); + + glm::vec3 particleConPos; + particleConPos.x = root["particleDim"].get( + "posX", MFluidSolver_DEFAULT_SCENE_SOURCE_POSX).asFloat(); + particleConPos.y = root["particleDim"].get( + "posY", MFluidSolver_DEFAULT_SCENE_SOURCE_POSY).asFloat(); + particleConPos.z = root["particleDim"].get( + "posZ", MFluidSolver_DEFAULT_SCENE_SOURCE_POSZ).asFloat(); + + // Get spawn method info + std::string spawningMethodString = root.get( + "spawnMethod", MFluidSolver_DEFAULT_SPAWNMETHODSTRING).asString(); + MUtils::toLowerInplace(&spawningMethodString); + if (spawningMethodString == "jittered") { + spawnMethod = ParticleSpawnMethod::Jittered; + } else if (spawningMethodString == "poissondisk") { + spawnMethod = ParticleSpawnMethod::PoissonDisk; + } else if (spawningMethodString == "uniform") { + spawnMethod = ParticleSpawnMethod::Uniform; + } + + // Load camera info + glm::vec3 cameraEye; + cameraEye.x = root["camera"].get( + "eyeX", MFluidSolver_DEFAULT_SCENE_CAMERA_EYEX).asFloat(); + cameraEye.y = root["camera"].get( + "eyeY", MFluidSolver_DEFAULT_SCENE_CAMERA_EYEY).asFloat(); + cameraEye.z = root["camera"].get( + "eyeZ", MFluidSolver_DEFAULT_SCENE_CAMERA_EYEZ).asFloat(); + + glm::vec3 cameraRef; + cameraRef.x = root["camera"].get( + "refX", MFluidSolver_DEFAULT_SCENE_CAMERA_REFX).asFloat(); + cameraRef.y = root["camera"].get( + "refY", MFluidSolver_DEFAULT_SCENE_CAMERA_REFY).asFloat(); + cameraRef.z = root["camera"].get( + "refZ", MFluidSolver_DEFAULT_SCENE_CAMERA_REFZ).asFloat(); + + // Set attributes + camera.setEyeRef(cameraEye, cameraRef); + + float particleSeparation = root.get( + "particleSeparation", MFluidSolver_DEFAULT_PARTICLE_SEPARATION).asFloat(); + solver.setParticleSeparation(particleSeparation); + + int maxParticles = root.get( + "maxParticles", MFluidSolver_DEFAULT_MAX_PARTICLES).asInt(); + solver.setMaxParticles(maxParticles); + + solver.init(containerDim * -0.5f, containerDim * 0.5f); + + // Create geometry + solver.fluidContainer = new Cube(glm::vec3(0)); + solver.fluidContainer->name = "Fluid Container"; + solver.fluidSource = new Cube(glm::vec3(0)); + solver.fluidSource->name = "Fluid Source"; + objects.push_back(solver.fluidContainer); + objects.push_back(solver.fluidSource); + + // Change geometry scale + solver.fluidContainer->transform.setScale(containerDim); + solver.fluidSource->transform.setTransform( + particleConPos, glm::vec3(0), particleDim); + + seedScene(); + + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + std::cout << "INFO: Particle count: " << + solver.numParticles() << " / " << solver.maxParticles() << std::endl; + #endif + + solver.sceneLoaded(); +} + +void Scene::seedScene() { + solver.fluidSource->spawnParticlesInVolume(&solver, spawnMethod); +} diff --git a/src/scene/scene.hpp b/src/scene/scene.hpp index 84022a43..3dcdd777 100644 --- a/src/scene/scene.hpp +++ b/src/scene/scene.hpp @@ -1,3 +1,37 @@ +// Copyright 2016 Robert Zhou // // scene.hpp -// Thanda \ No newline at end of file +// MFluidSolver + +#ifndef MFLUIDSOLVER_SCENE_HPP_ +#define MFLUIDSOLVER_SCENE_HPP_ + +#include +#include +#include + +#include "MFluidSolverConfig.hpp" + +#include "fluidSolver/sphSolver/iiSphSolver.hpp" +#include "camera.hpp" + +struct InvalidSceneException : std::exception { + const char *what() const noexcept {return "Invalid scene.\n";}; +}; + +class Scene { + public: + Scene(); + ~Scene(); + void loadJSON(const std::string &file); + void seedScene(); + + IISPHSolver solver; + std::vector objects; + Camera camera; + + private: + ParticleSpawnMethod spawnMethod; +}; + +#endif // MFLUIDSOLVER_SCENE_HPP_ diff --git a/src/scene/scene.json b/src/scene/scene.json deleted file mode 100644 index 789b4a42..00000000 --- a/src/scene/scene.json +++ /dev/null @@ -1,13 +0,0 @@ -{ - "containerDim" : { - "scaleX" : 5.0, - "scaleY" : 5.0, - "scaleZ" : 5.0 - }, - "particleDim" : { - "boundX" : 3.7, - "boundY" : 3.7, - "boundZ" : 3.7 - }, - "particleSeparation" : 0.1 -} \ No newline at end of file diff --git a/src/utils.cpp b/src/utils.cpp new file mode 100644 index 00000000..c04b695a --- /dev/null +++ b/src/utils.cpp @@ -0,0 +1,45 @@ +// Copyright 2016 Robert Zhou +// +// utils.cpp +// MFluidSolver + +#include "utils.hpp" + +#include +#include +#include +#include + +namespace MUtils { + std::string toHMS(double seconds) { + std::ostringstream strstm; + + int minutes = seconds / 60; + int hours = minutes / 60; + int iseconds = static_cast(seconds - minutes * 60); + minutes = minutes % 60; + + double dseconds = seconds - static_cast(seconds); + std::string truncateDSeconds = std::to_string(dseconds); + truncateDSeconds = truncateDSeconds.substr(1); + + strstm << zeroPad(hours, 1) << ":" << + zeroPad(minutes, 2) << ":" << + zeroPad(iseconds, 2) << + truncateDSeconds; + + return strstm.str(); + } + + void toLowerInplace(std::string *input) { + for (unsigned int i = 0; i < input->length(); ++i) { + (*input)[i] = std::tolower((*input)[i]); + } + } + + std::string zeroPad(int number, unsigned int digits) { + std::ostringstream ss; + ss << std::setw(digits) << std::setfill('0') << number; + return ss.str(); + } +} // namespace MUtils diff --git a/src/utils.hpp b/src/utils.hpp new file mode 100644 index 00000000..e5d007d3 --- /dev/null +++ b/src/utils.hpp @@ -0,0 +1,45 @@ +// Copyright 2016 Robert Zhou +// +// utils.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_UTILS_HPP_ +#define MFLUIDSOLVER_UTILS_HPP_ + +#include +#include + +namespace MUtils { + // Float approximate-equality comparison + template inline bool fequal(T a, T b, T epsilon = 0.0001) { + if (a == b) { + // Shortcut + return true; + } + + const T diff = std::abs(a - b); + if (a * b == 0) { + // a or b or both are zero; relative error is not meaningful here + return diff < (epsilon * epsilon); + } + + return diff / (std::abs(a) + std::abs(b)) < epsilon; + } + + // Based on pseudocode from https://en.wikipedia.org/wiki/Insertion_sort + template + inline void insertionSort(RandomIt first, RandomIt last, Compare comp) { + for (RandomIt i = first + 1; i < last; ++i) { + for (RandomIt j = i; j > first && comp(*j, *(j - 1)); --j) { + std::iter_swap(j, j - 1); + } + } + } + + // Non-template functions + std::string toHMS(double seconds); + void toLowerInplace(std::string *input); + std::string zeroPad(int number, unsigned int digits); +} // namespace MUtils + +#endif // MFLUIDSOLVER_UTILS_HPP_ diff --git a/src/viewer/input.cpp b/src/viewer/input.cpp new file mode 100644 index 00000000..ddae3fba --- /dev/null +++ b/src/viewer/input.cpp @@ -0,0 +1,59 @@ +// Copyright 2016 Robert Zhou +// +// input.cpp +// MFluidSolver + +#include "input.hpp" + +Viewer *Input::viewer = nullptr; + +void Input::computeArcballScrollCb( + GLFWwindow* window, double xoffset, double yoffset) { + if (viewer == nullptr) return; + viewer->scene.camera.zoom(yoffset); +} + +void Input::checkKeys( + GLFWwindow* window, int key, int scancode, int action, int mods) { + if (viewer == nullptr) return; + if (action == GLFW_RELEASE) { + switch (key) { + case GLFW_KEY_ESCAPE: { + viewer->stop(); + break; + } + case GLFW_KEY_E: { + #if MFluidSolver_USE_OPENVDB + viewer->scene.solver.exportVDB(); + #endif + break; + } + case GLFW_KEY_N: { + viewer->scene.solver.visualizeRandomParticlesNeighbors(); + break; + } + case GLFW_KEY_P: { + viewer->togglePause(); + break; + } + case GLFW_KEY_R: { + viewer->scene.seedScene(); + break; + } + case GLFW_KEY_S: { + viewer->screenshot(true); + break; + } + case GLFW_KEY_W: { + #if MFluidSolver_USE_PARTIO + viewer->scene.solver.exportBgeo(); + #endif + break; + } + case GLFW_KEY_RIGHT: { + viewer->scene.solver.updateStep(); + break; + } + } // end switch + } // end if RELEASE +} diff --git a/src/viewer/input.hpp b/src/viewer/input.hpp new file mode 100644 index 00000000..6e488790 --- /dev/null +++ b/src/viewer/input.hpp @@ -0,0 +1,22 @@ +// Copyright 2016 Robert Zhou +// +// input.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_VIEWER_INPUT_HPP_ +#define MFLUIDSOLVER_VIEWER_INPUT_HPP_ + +#include "MFluidSolverConfig.hpp" + +#include "viewer.hpp" + +class Input { + public: + static Viewer *viewer; + static void computeArcballScrollCb( + GLFWwindow *window, double xoffset, double yoffset); + static void checkKeys( + GLFWwindow *window, int key, int scancode, int action, int mods); +}; + +#endif // MFLUIDSOLVER_VIEWER_INPUT_HPP_ diff --git a/src/viewer/particleShaderProgram.cpp b/src/viewer/particleShaderProgram.cpp new file mode 100644 index 00000000..53389713 --- /dev/null +++ b/src/viewer/particleShaderProgram.cpp @@ -0,0 +1,191 @@ +// Copyright 2016 Robert Zhou +// +// particleShaderProgram.cpp +// MFluidSolver + +#include "particleShaderProgram.hpp" + +#include +#include +#include + +ParticleShaderProgram::ParticleShaderProgram(SPHSolver *solver, + const std::string &vertexShader, const std::string &fragmentShader, + const std::string &billboardDDS) + : ShaderProgram(vertexShader, fragmentShader), solver(solver), + aBillboardVertexArrID(-1), billboardVertexArrBufferID(-1), + uCameraRightVecID(-1), uCameraUpVecID(-1) { + // Get vertex buffer IDs + aBillboardVertexArrID = glGetAttribLocation(programID, "avs_Billboard"); + + // Get uniform IDs + uBillboardTextureSamplerID = + glGetUniformLocation(programID, "u_BillboardTextureSampler"); + uCameraRightVecID = glGetUniformLocation(programID, "u_CameraRight"); + uCameraUpVecID = glGetUniformLocation(programID, "u_CameraUp"); + uParticleSizeFloatID = glGetUniformLocation(programID, "u_ParticleSize"); + + // Load default particle size + setParticleSize(0.1f); + + // Load billboard texture + billboardTextureID = ShaderProgram::loadDDS(billboardDDS); + + // Fill billboard vertex buffer + glm::vec3 billboardData[] = { + glm::vec3(-0.5f, -0.5f, 0), + glm::vec3(0.5f, -0.5f, 0), + glm::vec3(-0.5f, 0.5f, 0), + glm::vec3(0.5f, 0.5f, 0) + }; + glGenBuffers(1, &billboardVertexArrBufferID); + glBindBuffer(GL_ARRAY_BUFFER, billboardVertexArrBufferID); + glBufferData(GL_ARRAY_BUFFER, 4 * sizeof(glm::vec3), billboardData, + GL_STATIC_DRAW); + + // Pre-allocate particle buffers + glGenBuffers(1, &particleColorArrBufferID); + glBindBuffer(GL_ARRAY_BUFFER, particleColorArrBufferID); + glBufferData(GL_ARRAY_BUFFER, solver->maxParticles() * sizeof(glm::vec3), + NULL, GL_STREAM_DRAW); + + glGenBuffers(1, &particlePositionArrBufferID); + glBindBuffer(GL_ARRAY_BUFFER, particlePositionArrBufferID); + glBufferData(GL_ARRAY_BUFFER, solver->maxParticles() * sizeof(glm::vec3), + NULL, GL_STREAM_DRAW); + + // Allocate particle space on RAM + particleColorArray = new glm::vec3[solver->maxParticles()]; + particlePositionArray = new glm::vec3[solver->maxParticles()]; + + // Print debug info, particularly info about unbounded variables + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_DEBUG + std::cout << "DEBUG:SHADER: Program " << programID << + " is type particles" << std::endl; + #endif + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_WARN + if (aBillboardVertexArrID == -1) { + std::cout << "WARN: avs_Billboard is not bound for program " << + programID << std::endl; + } + if (uBillboardTextureSamplerID == -1) { + std::cout << "WARN: u_BillboardTextureSampler is not bound for program " << + programID << std::endl; + } + if (uCameraRightVecID == -1) { + std::cout << "WARN: u_CameraRight is not bound for program " << + programID << std::endl; + } + if (uCameraUpVecID == -1) { + std::cout << "WARN: u_CameraUp is not bound for program " << + programID << std::endl; + } + if (uParticleSizeFloatID == -1) { + std::cout << "WARN: u_ParticleSize is not bound for program " << + programID << std::endl; + } + #endif +} + +ParticleShaderProgram::~ParticleShaderProgram() { + if (billboardVertexArrBufferID != -1) + glDeleteBuffers(1, &billboardVertexArrBufferID); + if (particleColorArrBufferID != -1) + glDeleteBuffers(1, &particleColorArrBufferID); + if (particlePositionArrBufferID != -1) + glDeleteBuffers(1, &particlePositionArrBufferID); + + delete particleColorArray; + delete particlePositionArray; +} + +void ParticleShaderProgram::setCameraVectors( + const glm::vec3 &right, const glm::vec3 &up) { + glUseProgram(programID); + + if (uCameraRightVecID != -1) { + // http://glm.g-truc.net/0.9.2/api/a00001.html + glUniform3f(uCameraRightVecID, right.x, right.y, right.z); + } + + if (uCameraUpVecID != -1) { + // http://glm.g-truc.net/0.9.2/api/a00001.html + glUniform3f(uCameraUpVecID, up.x, up.y, up.z); + } +} + +void ParticleShaderProgram::setParticleSize(float size) { + glUseProgram(programID); + if (uParticleSizeFloatID != -1) { + glUniform1f(uParticleSizeFloatID, size); + } +} + +void ParticleShaderProgram::draw() { + glEnable(GL_BLEND); + glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); + + glUseProgram(programID); + + // Copy data from particles to RAM buffers + std::vector &particles = solver->particles(); + for (unsigned int i = 0; i < particles.size(); ++i) { + particleColorArray[i] = glm::vec3(particles[i].color); + particlePositionArray[i] = particles[i].position(); + } + + // Set billboard texture to texture unit 0 + if (billboardTextureID != -1) { + glActiveTexture(GL_TEXTURE0); + glBindTexture(GL_TEXTURE_2D, billboardTextureID); + glUniform1i(uBillboardTextureSamplerID, 0); + } + + // Update buffers + if (particleColorArrBufferID != -1) { + glBindBuffer(GL_ARRAY_BUFFER, particleColorArrBufferID); + glBufferData(GL_ARRAY_BUFFER, // Buffer orphaning for performance + solver->maxParticles() * sizeof(glm::vec3), NULL, GL_STREAM_DRAW); + glBufferSubData(GL_ARRAY_BUFFER, 0, + solver->numParticles() * sizeof(glm::vec3), particleColorArray); + } + if (particlePositionArrBufferID != -1) { + glBindBuffer(GL_ARRAY_BUFFER, particlePositionArrBufferID); + glBufferData(GL_ARRAY_BUFFER, // Buffer orphaning for performance + solver->maxParticles() * sizeof(glm::vec3), NULL, GL_STREAM_DRAW); + glBufferSubData(GL_ARRAY_BUFFER, 0, + solver->numParticles() * sizeof(glm::vec3), particlePositionArray); + } + + // Insert data to attribute variables + if (aBillboardVertexArrID != -1 && billboardVertexArrBufferID != -1) { + glEnableVertexAttribArray(aBillboardVertexArrID); + glBindBuffer(GL_ARRAY_BUFFER, billboardVertexArrBufferID); + glVertexAttribPointer(aBillboardVertexArrID, 3, GL_FLOAT, GL_FALSE, 0, NULL); + // glVertexAttribDivisor(aBillboardVertexArrID, 0); + } + if (aVertexColorArrID != -1 && particleColorArrBufferID != -1) { + glEnableVertexAttribArray(aVertexColorArrID); + glBindBuffer(GL_ARRAY_BUFFER, particleColorArrBufferID); + glVertexAttribPointer(aVertexColorArrID, 3, GL_FLOAT, GL_FALSE, 0, NULL); + // glVertexAttribDivisor(aVertexColorArrID, 1); + } + if (aVertexPositionArrID != -1 && particlePositionArrBufferID != -1) { + glEnableVertexAttribArray(aVertexPositionArrID); + glBindBuffer(GL_ARRAY_BUFFER, particlePositionArrBufferID); + glVertexAttribPointer(aVertexPositionArrID, 3, GL_FLOAT, GL_FALSE, 0, NULL); + // glVertexAttribDivisor(aVertexPositionArrID, 1); + } + + // Draw billboards on all them particles! + // glDrawArraysInstanced(GL_TRIANGLE_STRIP, 0, 4, solver->numParticles()); + glDrawArrays(GL_POINTS, 0, solver->numParticles()); + + // Disable attributes + if (aBillboardVertexArrID != -1) + glDisableVertexAttribArray(aBillboardVertexArrID); + if (aVertexColorArrID != -1) + glDisableVertexAttribArray(aVertexColorArrID); + if (aVertexPositionArrID != -1) + glDisableVertexAttribArray(aVertexPositionArrID); +} diff --git a/src/viewer/particleShaderProgram.hpp b/src/viewer/particleShaderProgram.hpp new file mode 100644 index 00000000..7838559b --- /dev/null +++ b/src/viewer/particleShaderProgram.hpp @@ -0,0 +1,48 @@ +// Copyright 2016 Robert Zhou +// +// particleShaderProgram.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_VIEWER_PARTICLESHADERPROGRAM_HPP_ +#define MFLUIDSOLVER_VIEWER_PARTICLESHADERPROGRAM_HPP_ + +#include +#include +#include + +#include "MFluidSolverConfig.hpp" + +#include "fluidSolver/sphSolver/sphSolver.hpp" +#include "shaderProgram.hpp" + +class ParticleShaderProgram : public ShaderProgram { + public: + ParticleShaderProgram(SPHSolver *solver, + const std::string &vertexShader, const std::string &fragmentShader, + const std::string &billboardDDS); + ~ParticleShaderProgram(); + void draw(); + void setCameraVectors(const glm::vec3 &right, const glm::vec3 &up); + void setParticleSize(float size); + + glm::vec3 *particleColorArray; + glm::vec3 *particlePositionArray; + + private: + SPHSolver *solver; + GLuint billboardTextureID; + + // Variable IDs + GLuint uBillboardTextureSamplerID; + GLuint aBillboardVertexArrID; + GLuint uCameraRightVecID; + GLuint uCameraUpVecID; + GLuint uParticleSizeFloatID; + + // Buffer IDs + GLuint billboardVertexArrBufferID; + GLuint particleColorArrBufferID; + GLuint particlePositionArrBufferID; +}; + +#endif // MFLUIDSOLVER_VIEWER_PARTICLESHADERPROGRAM_HPP_ diff --git a/src/viewer/shaderProgram.cpp b/src/viewer/shaderProgram.cpp new file mode 100644 index 00000000..6b5dbb44 --- /dev/null +++ b/src/viewer/shaderProgram.cpp @@ -0,0 +1,283 @@ +// Copyright 2016 Robert Zhou +// +// shaderProgram.cpp +// MFluidSolver + +#include "shaderProgram.hpp" + +#include +#include +#include +#include +#include + +ShaderProgram::ShaderProgram( + const std::string &vertexShader, const std::string &fragmentShader) + : _vertexShader(vertexShader), _fragmentShader(fragmentShader), + programID(-1), + vertexShaderID(-1), fragmentShaderID(-1), + aVertexColorArrID(-1), aVertexPositionArrID(-1), + uModelMatID(-1), uViewProjectionMatID(-1) { + compile(); + + // Get vertex buffer IDs + aVertexColorArrID = glGetAttribLocation(programID, "avs_Color"); + aVertexPositionArrID = glGetAttribLocation(programID, "avs_Position"); + + // Get uniform IDs + uModelMatID = glGetUniformLocation(programID, "u_Model"); + uViewProjectionMatID = glGetUniformLocation(programID, "u_ViewProjection"); + + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_DEBUG + std::cout << "DEBUG:SHADER: Created program ID " << programID << std::endl; + #endif + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_WARN + if (aVertexColorArrID == -1) { + std::cout << "WARN: avs_Color is not bound for program " << + programID << std::endl; + } + if (aVertexPositionArrID == -1) { + std::cout << "WARN: avs_Position is not bound for program " << + programID << std::endl; + } + if (uModelMatID == -1) { + std::cout << "WARN: u_Model is not bound for program " << + programID << std::endl; + } + if (uViewProjectionMatID == -1) { + std::cout << "WARN: u_ViewProjection is not bound for program " << + programID << std::endl; + } + #endif +} + +ShaderProgram::~ShaderProgram() { + glDeleteProgram(programID); +} + +void ShaderProgram::draw(Drawable *d) { + glUseProgram(programID); + + // Insert data to attribute variables + if (aVertexColorArrID != -1 && d->bindColorBuffer()) { + glEnableVertexAttribArray(aVertexColorArrID); + glVertexAttribPointer(aVertexColorArrID, 3, GL_FLOAT, GL_FALSE, 0, NULL); + } + if (aVertexPositionArrID != -1 && d->bindPositionBuffer()) { + glEnableVertexAttribArray(aVertexPositionArrID); + glVertexAttribPointer(aVertexPositionArrID, 3, GL_FLOAT, GL_FALSE, 0, NULL); + } + + // Draw (indices) + d->bindIndexBuffer(); + glDrawElements(d->drawMode(), d->idxCount(), GL_UNSIGNED_INT, NULL); + + // Disable attributes + if (aVertexColorArrID != -1) + glDisableVertexAttribArray(aVertexColorArrID); + if (aVertexPositionArrID != -1) + glDisableVertexAttribArray(aVertexPositionArrID); +} + +void ShaderProgram::setModelMat(const glm::mat4 &modelMat) { + glUseProgram(programID); + + if (uModelMatID != -1) { + // http://glm.g-truc.net/0.9.2/api/a00001.html + glUniformMatrix4fv(uModelMatID, 1, GL_FALSE, &modelMat[0][0]); + } +} + +void ShaderProgram::setViewProjectionMat(const glm::mat4 &viewProjectionMat) { + glUseProgram(programID); + + if (uViewProjectionMatID != -1) { + // http://glm.g-truc.net/0.9.2/api/a00001.html + glUniformMatrix4fv(uViewProjectionMatID, 1, GL_FALSE, + &viewProjectionMat[0][0]); + } +} + +// Source: https://github.com/opengl-tutorials/ogl/blob/master/common/texture.cpp +#define FOURCC_DXT1 0x31545844 // Equivalent to "DXT1" in ASCII +#define FOURCC_DXT3 0x33545844 // Equivalent to "DXT3" in ASCII +#define FOURCC_DXT5 0x35545844 // Equivalent to "DXT5" in ASCII +GLuint ShaderProgram::loadDDS(const std::string &file) { + unsigned char header[124]; + FILE *fp; + + // Try to open the file + fp = fopen(file.c_str(), "rb"); + if (fp == NULL) { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_ERROR + std::cerr << "ERROR: " << file << " could not be opened" << std::endl; + #endif + + getchar(); + return 0; + } + + // Verify the type of file + char filecode[4]; + fread(filecode, 1, 4, fp); + if (std::strncmp(filecode, "DDS ", 4) != 0) { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_ERROR + std::cerr << "ERROR: " << file << " is not a DDS file" << std::endl; + #endif + + getchar(); + fclose(fp); + return 0; + } + + // Get the surface desc + fread(&header, 124, 1, fp); + + unsigned int height = *(unsigned int *)&(header[8 ]); + unsigned int width = *(unsigned int *)&(header[12]); + unsigned int linearSize = *(unsigned int *)&(header[16]); + unsigned int mipMapCount = *(unsigned int *)&(header[24]); + unsigned int fourCC = *(unsigned int *)&(header[80]); + + unsigned char *buffer; + unsigned int bufsize; + // How big is it going to be including all mipmaps? + bufsize = mipMapCount > 1 ? linearSize * 2 : linearSize; + buffer = (unsigned char *) malloc(bufsize * sizeof(unsigned char)); + fread(buffer, 1, bufsize, fp); + // Close the file pointer + fclose(fp); + + unsigned int components = (fourCC == FOURCC_DXT1) ? 3 : 4; + unsigned int format; + switch (fourCC) { + case FOURCC_DXT1: { + format = GL_COMPRESSED_RGBA_S3TC_DXT1_EXT; + break; + } + case FOURCC_DXT3: { + format = GL_COMPRESSED_RGBA_S3TC_DXT3_EXT; + break; + } + case FOURCC_DXT5: { + format = GL_COMPRESSED_RGBA_S3TC_DXT5_EXT; + break; + } + default: { + free(buffer); + return 0; + } + } + + // Create one OpenGL texture + GLuint textureID; + glGenTextures(1, &textureID); + + // "Bind" the newly created texture + // all future texture functions will modify this texture + glBindTexture(GL_TEXTURE_2D, textureID); + glPixelStorei(GL_UNPACK_ALIGNMENT, 1); + + unsigned int blockSize = + (format == GL_COMPRESSED_RGBA_S3TC_DXT1_EXT) ? 8 : 16; + unsigned int offset = 0; + + // Load the mipmaps + for (unsigned int level = 0; + level < mipMapCount && (width || height); ++level) { + unsigned int size = ((width + 3) / 4) * ((height + 3) / 4) * blockSize; + glCompressedTexImage2D(GL_TEXTURE_2D, level, format, width, height, + 0, size, buffer + offset); + + offset += size; + width /= 2; + height /= 2; + + // Deal with Non-Power-Of-Two textures. + // This code is not included in the webpage to reduce clutter. + if (width < 1) width = 1; + if (height < 1) height = 1; + } + + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + std::cout << "INFO: Loaded texture " << file << std::endl; + #endif + + free(buffer); + return textureID; +} + +void ShaderProgram::compile() { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + std::cout << "INFO: Loading shader program (VS:" << _vertexShader << + ", FS:" << _fragmentShader << ")" << std::endl; + #endif + + vertexShaderID = glCreateShader(GL_VERTEX_SHADER); + fragmentShaderID = glCreateShader(GL_FRAGMENT_SHADER); + + // Read shader files to string + std::ifstream vsStream(_vertexShader); + std::string vsText((std::istreambuf_iterator(vsStream)), + std::istreambuf_iterator()); + std::ifstream fsStream(_fragmentShader); + std::string fsText((std::istreambuf_iterator(fsStream)), + std::istreambuf_iterator()); + + // Compile + GLint compileResult = GL_FALSE; + int compileLogLen = 0; + + char const *vsTextPtr = vsText.c_str(); + glShaderSource(vertexShaderID, 1, &(vsTextPtr), NULL); + glCompileShader(vertexShaderID); + glGetShaderiv(vertexShaderID, GL_COMPILE_STATUS, &compileResult); + glGetShaderiv(vertexShaderID, GL_INFO_LOG_LENGTH, &compileLogLen); + if (compileLogLen > 1) { + std::vector vsErrMessage(compileLogLen + 1); + glGetShaderInfoLog(vertexShaderID, compileLogLen, NULL, &vsErrMessage[0]); + + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_ERROR + std::cerr << "ERROR: Error compiling vertex shader " << _vertexShader << + ": " << &vsErrMessage[0] << std::endl; + #endif + } + + char const *fsTextPtr = fsText.c_str(); + glShaderSource(fragmentShaderID, 1, &(fsTextPtr), NULL); + glCompileShader(fragmentShaderID); + glGetShaderiv(fragmentShaderID, GL_COMPILE_STATUS, &compileResult); + glGetShaderiv(fragmentShaderID, GL_INFO_LOG_LENGTH, &compileLogLen); + if (compileLogLen > 1) { + std::vector fsErrMessage(compileLogLen + 1); + glGetShaderInfoLog(fragmentShaderID, compileLogLen, NULL, &fsErrMessage[0]); + + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_ERROR + std::cerr << "ERROR: Error compiling fragment shader " << _fragmentShader << + ": " << &fsErrMessage[0] << std::endl; + #endif + } + + // Link the two programs + programID = glCreateProgram(); + glAttachShader(programID, vertexShaderID); + glAttachShader(programID, fragmentShaderID); + glLinkProgram(programID); + glGetProgramiv(programID, GL_LINK_STATUS, &compileResult); + glGetProgramiv(programID, GL_INFO_LOG_LENGTH, &compileLogLen); + if (compileLogLen > 1) { + std::vector programErrMessage(compileLogLen + 1); + glGetProgramInfoLog(programID, compileLogLen, NULL, &programErrMessage[0]); + + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_ERROR + std::cerr << "ERROR: Error linking program (VS:" << _vertexShader << + ", FS:" << _fragmentShader << "): " << &programErrMessage[0] << std::endl; + #endif + } + + // Free space: http://gamedev.stackexchange.com/questions/47910/after-a-succesful-gllinkprogram-should-i-delete-detach-my-shaders + glDetachShader(programID, vertexShaderID); + glDetachShader(programID, fragmentShaderID); + glDeleteShader(vertexShaderID); + glDeleteShader(fragmentShaderID); +} diff --git a/src/viewer/shaderProgram.hpp b/src/viewer/shaderProgram.hpp new file mode 100644 index 00000000..d11d8726 --- /dev/null +++ b/src/viewer/shaderProgram.hpp @@ -0,0 +1,45 @@ +// Copyright 2016 Robert Zhou +// +// shaderProgram.hpp +// MFluidSolver + +#ifndef MFLUIDSOLVER_VIEWER_SHADERPROGRAM_HPP_ +#define MFLUIDSOLVER_VIEWER_SHADERPROGRAM_HPP_ + +#include +#include +#include +#include + +#include "MFluidSolverConfig.hpp" + +#include "geom/geom.hpp" + +class ShaderProgram { + public: + ShaderProgram(const std::string &vertexShader, + const std::string &fragmentShader); + ~ShaderProgram(); + void draw(Drawable *d); + void setModelMat(const glm::mat4 &modelMat); + void setViewProjectionMat(const glm::mat4 &viewProjectionMat); + + static GLuint loadDDS(const std::string &file); + + // Required for derived classes + void compile(); + + GLuint programID; + GLuint vertexShaderID; + GLuint fragmentShaderID; + GLuint aVertexColorArrID; + GLuint aVertexPositionArrID; + GLuint uViewProjectionMatID; + + std::string _vertexShader, _fragmentShader; + + private: + GLuint uModelMatID; +}; + +#endif // MFLUIDSOLVER_VIEWER_SHADERPROGRAM_HPP_ diff --git a/src/viewer/viewer.cpp b/src/viewer/viewer.cpp index 60d1868a..a6574824 100644 --- a/src/viewer/viewer.cpp +++ b/src/viewer/viewer.cpp @@ -1,5 +1,384 @@ +// Copyright 2016 Robert Zhou // // viewer.cpp -// Thanda +// MFluidSolver #include "viewer.hpp" + +#include +#include +#include +#include +#include +#include + +#include + +#if MFluidSolver_USE_TBB +#include +#include +#endif + +#include "input.hpp" + +#if MFluidSolver_OPENGL_DEBUG +void APIENTRY openglDebugCallbackFunction(GLenum source, GLenum type, GLuint id, + GLenum severity, GLsizei length, const GLchar* message, void* userParam) { + std::string sType; + std::string sSeverity; + std::ostream *outStream = &std::cout; + + switch (type) { + case GL_DEBUG_TYPE_ERROR: { + sType = "ERROR"; + outStream = &std::cerr; + break; + } + case GL_DEBUG_TYPE_DEPRECATED_BEHAVIOR: { + sType = "DEPRECATED"; + break; + } + case GL_DEBUG_TYPE_UNDEFINED_BEHAVIOR: { + sType = "UNDEFINED"; + break; + } + case GL_DEBUG_TYPE_PORTABILITY: { + sType = "PORTABILITY"; + break; + } + case GL_DEBUG_TYPE_PERFORMANCE: { + sType = "PERFORMANCE"; + break; + } + case GL_DEBUG_TYPE_OTHER: { + sType = "OTHER"; + break; + } + } + switch (severity) { + case GL_DEBUG_SEVERITY_LOW: { + sSeverity = "LOW"; + #if MFluidSolver_LOG_LEVEL > MFluidSolver_LOG_WARN + return; + #endif + break; + } + case GL_DEBUG_SEVERITY_MEDIUM: { + sSeverity = "MEDIUM"; + #if MFluidSolver_LOG_LEVEL > MFluidSolver_LOG_ERROR + return; + #endif + break; + } + case GL_DEBUG_SEVERITY_HIGH: { + sSeverity = "HIGH"; + #if MFluidSolver_LOG_LEVEL > MFluidSolver_LOG_ERROR + return; + #endif + break; + } + default: { + sSeverity = "OTHER"; + #if MFluidSolver_LOG_LEVEL > MFluidSolver_LOG_INFO + return; + #endif + } + } + *outStream << "GL:" << sType << + "(" << sSeverity << "): " << message << std::endl; +} +#endif + +Viewer::Viewer() + : wireShader(nullptr), particleShader(nullptr), + oldLeftState(GLFW_RELEASE), oldRightState(GLFW_RELEASE), + paused(true), shouldStop(false), + autoRender(MFluidSolver_DEFAULT_AUTORENDER), + renderSkip(MFluidSolver_DEFAULT_RENDERSKIP) { +} + +Viewer::~Viewer() { + delete wireShader; +} + +void Viewer::init(int width, int height) { + _width = width; + _height = height; + + // http://www.opengl-tutorial.org/beginners-tutorials/tutorial-1-opening-a-window/ + glfwWindowHint(GLFW_SAMPLES, 4); // 4x antialiasing + glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 3); + glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 3); + // To make MacOS happy; should not be needed + glfwWindowHint(GLFW_OPENGL_FORWARD_COMPAT, GL_TRUE); + // We don't want the old OpenGL + glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE); + + #if MFluidSolver_OPENGL_DEBUG + glfwWindowHint(GLFW_OPENGL_DEBUG_CONTEXT, GL_TRUE); + #endif + + window = glfwCreateWindow(width, height, "MFluidSolver", NULL, NULL); + if (window == NULL) { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_FATAL + std::cerr << "FATAL: Failed to open GLFW window" << std::endl; + #endif + + getchar(); // Wait for key before quit + glfwTerminate(); + throw GLFWWindowInitException(); + } + + glfwMakeContextCurrent(window); // Initialize GLEW + + glewExperimental = GL_TRUE; + if (glewInit() != GLEW_OK) { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_FATAL + std::cerr << "FATAL: Failed to initialize GLEW" << std::endl; + #endif + + getchar(); // Wait for key before quit + glfwTerminate(); + throw GLEWInitException(); + } + + // Sticky keys and mouse buttons + glfwSetInputMode(window, GLFW_STICKY_MOUSE_BUTTONS, GL_TRUE); + + // Key callback + glfwSetKeyCallback(window, Input::checkKeys); + + // Mouse callbacks + glfwSetScrollCallback(window, Input::computeArcballScrollCb); + + #if MFluidSolver_OPENGL_DEBUG + // Set debug callback + if (glDebugMessageCallback) { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + std::cout << "DEBUG: Registering OpenGL debug callback" << std:: endl; + #endif + + glEnable(GL_DEBUG_OUTPUT_SYNCHRONOUS); + glDebugMessageCallback(openglDebugCallbackFunction, nullptr); + GLuint unusedIds = 0; + glDebugMessageControl(GL_DONT_CARE, GL_DONT_CARE, GL_DEBUG_SEVERITY_HIGH, + 0, &unusedIds, true); + } else { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + std::cout << "DEBUG: OpenGL debug callback not available" << std::endl; + #endif + } + #endif + + // Init screenshot array + screenshotArray = std::vector(3 * width * height); + screenshotArrayFlipped = std::vector(3 * width * height); + + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + if (paused) std::cout << "INFO: We are currently paused!" << std::endl; + if (!paused) std::cout << "INFO: We are currently playing!" << std::endl; + #endif + + // Create necessary output folders + #if MFluidSolver_PARTICLE_STATS_FILES + boost::filesystem::path particleStatsDir("particlestats"); + boost::filesystem::create_directory(particleStatsDir); + if (!boost::filesystem::is_directory(particleStatsDir)) { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_ERROR + std::cout << "ERROR: particlestats folder could not be created" << + std::endl; + #endif + } + #endif +} + +void Viewer::run() { + double oldTime = glfwGetTime(); + double currTime, deltaT; + GLenum glErrorCode; + const GLubyte *glErrorString; + + do { + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); + + // Get elapsed time + currTime = glfwGetTime(); + deltaT = currTime - oldTime; + oldTime = currTime; + + // Update and render particles + if (!paused) { + scene.solver.update(deltaT); + if (scene.solver.hasEndedSimulation()) { + togglePause(); + } + } + + // Convenience method for stopping before blowing up + #if MFluidSolver_PARTICLE_STATS + if (scene.solver.numFlyaways > 0) { + paused = true; + } + #endif + + // Render particles + particleShader->setViewProjectionMat(scene.camera.getViewProjection()); + particleShader->setCameraVectors(scene.camera.right(), scene.camera.up()); + particleShader->draw(); + + // Render boxes + wireShader->setViewProjectionMat(scene.camera.getViewProjection()); + // for (Geometry *g : scene.objects) { + Geometry *g = scene.solver.fluidContainer; + wireShader->setModelMat(g->transform.T()); + wireShader->draw(g); + // } + + // Parse mouse button input + leftState = glfwGetMouseButton(window, GLFW_MOUSE_BUTTON_LEFT); + rightState = glfwGetMouseButton(window, GLFW_MOUSE_BUTTON_RIGHT); + if (leftState == GLFW_PRESS) { + if (oldLeftState == GLFW_RELEASE) { + glfwGetCursorPos(window, &oldPos.x, &oldPos.y); + } else { + glfwGetCursorPos(window, &newPos.x, &newPos.y); + scene.camera.arcball(oldPos, newPos); + oldPos = glm::vec2(newPos); + } + } + if (rightState == GLFW_PRESS) { + if (oldRightState == GLFW_RELEASE) { + glfwGetCursorPos(window, &oldPos.x, &oldPos.y); + } else { + glfwGetCursorPos(window, &newPos.x, &newPos.y); + scene.camera.pan(oldPos, newPos); + oldPos = glm::vec2(newPos); + } + } + oldLeftState = leftState; + oldRightState = rightState; + + // Swap buffers + glfwSwapBuffers(window); + + // Render + if (autoRender) { + if (renderSkip > 1 && scene.solver.updateNumber() % renderSkip == 0) { + screenshot(false); + } + } + + glfwPollEvents(); + + #if MFluidSolver_OPENGL_DEBUG == 0 + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_ERROR + // Check for errors + if ((glErrorCode = glGetError()) != GL_NO_ERROR) { + glErrorString = gluErrorString(glErrorCode); + std::cerr << "GL:ERROR: " << glErrorString << std::endl; + } + #endif + #endif + } while (!shouldStop && + glfwWindowShouldClose(window) == 0); +} + +inline unsigned int Viewer::screenshotIndex(unsigned int x, unsigned int y) { + return y * _width + x; +} + +void Viewer::screenshot(bool manual) { + // Capture + glReadPixels(0, 0, _width, _height, + GL_RGB, GL_UNSIGNED_BYTE, screenshotArray.data()); + + // Flip + std::string outString = + MUtils::zeroPad(scene.solver.updateNumber(), 6) + ".tga"; + if (manual) { + outString = "screenshot_" + outString; + } else { + outString = "render/render_" + outString; + } + #if MFluidSolver_USE_TBB + tbb::parallel_for(tbb::blocked_range2d(0, _width, 0, _height), + [&](const tbb::blocked_range2d &r) { + for (unsigned int x = r.rows().begin(); x != r.rows().end(); ++x) { + for (unsigned int y = r.cols().begin(); y != r.cols().end(); ++y) { + #else + for (unsigned int x = 0; x < _width; ++x) { + for (unsigned int y = 0; y < _height / 2; ++y) { + #endif + for (unsigned int i = 0; i < 3; ++i) { + screenshotArrayFlipped[3 * screenshotIndex(x, y) + i] = + screenshotArray[3 * screenshotIndex(x, _height - 1 - y) + i]; + screenshotArrayFlipped[3 * screenshotIndex(x, _height - 1 - y) + i] = + screenshotArray[3 * screenshotIndex(x, y) + i]; + } + #if MFluidSolver_USE_TBB + } + } + } + ); + #else + } + } + #endif + + // Write + int result = stbi_write_tga(outString.c_str(), + _width, _height, 3, screenshotArrayFlipped.data()); + if (manual) { + if (result != 0) { // Success + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + std::cout << "INFO: Wrote rendered image to: " << + outString << std::endl; + #endif + } else { // Failure + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_ERROR + std::cout << "ERROR: Failed to write rendered image to: " << + outString << std::endl; + #endif + } + } +} + +void Viewer::configureScreenshot(bool render, unsigned int skip) { + autoRender = render; + renderSkip = skip; + + boost::filesystem::path renderDir("render"); + boost::filesystem::create_directory(renderDir); + if (!boost::filesystem::is_directory(renderDir)) { + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_ERROR + std::cout << + "ERROR: Failed to create render folder. Disabling autorender" << + std::endl; + #endif + autoRender = false; + } + + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + if (autoRender) { + if (renderSkip > 1) { + std::cout << "INFO: Rendering 1 frame for every " << + renderSkip << " to file!" << std::endl; + } else { + std::cout << "INFO: Rendering every frame to file!" << std::endl; + } + } + #endif +} + +void Viewer::stop() { + shouldStop = true; +} + +void Viewer::togglePause() { + paused = !paused; + + #if MFluidSolver_LOG_LEVEL <= MFluidSolver_LOG_INFO + if (paused) std::cout << "INFO: Paused!" << std::endl; + if (!paused) std::cout << "INFO: Unpaused!" << std::endl; + #endif +} diff --git a/src/viewer/viewer.hpp b/src/viewer/viewer.hpp index abb78a17..0085b9da 100644 --- a/src/viewer/viewer.hpp +++ b/src/viewer/viewer.hpp @@ -1,9 +1,63 @@ +// Copyright 2016 Robert Zhou // // viewer.hpp -// Thanda +// MFluidSolver -#ifndef viewer_hpp -#define viewer_hpp +#ifndef MFLUIDSOLVER_VIEWER_HPP_ +#define MFLUIDSOLVER_VIEWER_HPP_ +#include +#include +#include -#endif /* viewer_hpp */ +#include "MFluidSolverConfig.hpp" + +#include "scene/scene.hpp" +#include "particleShaderProgram.hpp" +#include "utils.hpp" + +struct GLFWWindowInitException : std::exception { + const char *what() const noexcept { + return "Failed to initialize GLFW window.\n"; + }; +}; +struct GLEWInitException : std::exception { + const char *what() const noexcept { + return "Failed to initialize GLEW.\n"; + }; +}; + +class Viewer { + public: + Viewer(); + ~Viewer(); + void init(int width = MFluidSolver_DEFAULT_WINDOW_WIDTH, + int height = MFluidSolver_DEFAULT_WINDOW_HEIGHT); + void run(); + void screenshot(bool manual = false); + void stop(); + void togglePause(); + + void configureScreenshot(bool render, unsigned int skip); + + GLFWwindow *window; + Scene scene; + ShaderProgram *wireShader; + ParticleShaderProgram *particleShader; + + private: + int _width, _height; + + int leftState, rightState; + int oldLeftState, oldRightState; + glm::dvec2 oldPos, newPos; + bool paused, shouldStop; + + bool autoRender; + unsigned int renderSkip; + std::vector screenshotArray; + std::vector screenshotArrayFlipped; + inline unsigned int screenshotIndex(unsigned int x, unsigned int y); +}; + +#endif // MFLUIDSOLVER_VIEWER_HPP_ diff --git a/texture/particle.dds b/texture/particle.dds new file mode 100644 index 00000000..7dd08563 Binary files /dev/null and b/texture/particle.dds differ diff --git a/tools/particle_attribute_plotter.py b/tools/particle_attribute_plotter.py new file mode 100644 index 00000000..6b74facf --- /dev/null +++ b/tools/particle_attribute_plotter.py @@ -0,0 +1,28 @@ +#!/usr/bin/env python +import numpy as np +import matplotlib.pyplot as plt +import sys + +spread = np.random.rand(50) * 100 +center = np.ones(25) * 50 +flier_high = np.random.rand(10) * 100 + 100 +flier_low = np.random.rand(10) * -100 +data = np.concatenate((spread, center, flier_high, flier_low), 0) + +pressure_data = np.loadtxt("../build/particlestats/psr-" + sys.argv[1] + ".txt") +pressure_force_data = np.loadtxt("../build/particlestats/pfm-" + sys.argv[1] + ".txt") +velocity_data = np.loadtxt("../build/particlestats/vel-" + sys.argv[1] + ".txt") + +fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(18, 12)) + +axes[0].boxplot(pressure_data) +axes[0].set_title('Pressure') + +axes[1].boxplot(pressure_force_data) +axes[1].set_title('Pressure Force') + +axes[2].boxplot(velocity_data) +axes[2].set_title('Velocity') + +plt.tight_layout() +plt.show() \ No newline at end of file diff --git a/unittest/CMakeLists.txt b/unittest/CMakeLists.txt new file mode 100644 index 00000000..4f0a162a --- /dev/null +++ b/unittest/CMakeLists.txt @@ -0,0 +1,23 @@ +######################## +# MFluidSolver # +######################## + ##### Tests ##### + +# Boost Testing Framework +find_package(Boost COMPONENTS unit_test_framework REQUIRED) +include_directories(${Boost_INCLUDE_DIRS}) +add_definitions(-DBOOST_TEST_DYN_LINK) + +set(MTEST_SRC + main.cpp # All cpp test files are included in main +) + +set(MTEST_LIBS + MKernelFunctions + MNeighborSearch + MZCurve + ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY} +) + +add_executable(MTest ${MTEST_SRC}) +target_link_libraries(MTest ${MTEST_LIBS}) \ No newline at end of file diff --git a/unittest/main.cpp b/unittest/main.cpp new file mode 100644 index 00000000..bb5ad086 --- /dev/null +++ b/unittest/main.cpp @@ -0,0 +1,11 @@ +// Copyright 2016 Robert Zhou +// +// main.cpp +// MFluidSolver + +#define BOOST_NORETURN +#define BOOST_TEST_MODULE "MFluidSolver C++ Unit Tests" + +#include "testKernelFunctions.cpp" +#include "testNeighborSearch.cpp" +#include "testZCurve.cpp" diff --git a/unittest/testKernelFunctions.cpp b/unittest/testKernelFunctions.cpp new file mode 100644 index 00000000..1e19c0a3 --- /dev/null +++ b/unittest/testKernelFunctions.cpp @@ -0,0 +1,135 @@ +// Copyright 2016 Robert Zhou +// +// testKernelFunctions.cpp +// MFluidSolver + +#include + +#include "fluidSolver/sphSolver/kernelFunctions.hpp" + +BOOST_AUTO_TEST_SUITE(KernelFunctionsTests) + +BOOST_AUTO_TEST_CASE(KernelFunctions_Poly6_InsideKernelTest) { + KernelFunctions kernelFunctions; + double out; + + kernelFunctions.setKernelRadius(0.1f); + + out = kernelFunctions.computePoly6(glm::vec3(0.05f, 0.05f, 0.05f)); + BOOST_CHECK_CLOSE(out, 24.4794, 0.01); + + out = kernelFunctions.computePoly6(glm::vec3(0.03f, 0.04f, 0.05f)); + BOOST_CHECK_CLOSE(out, 195.835, 0.01); + + out = kernelFunctions.computePoly6(glm::vec3(0.003f, 0.002f, 0.001f)); + BOOST_CHECK_CLOSE(out, 1560.11, 0.01); +} + +/*BOOST_AUTO_TEST_CASE(KernelFunctions_Poly6_OutsideKernelTest) { + KernelFunctions kernelFunctions; + double out; + + kernelFunctions.setKernelRadius(0.1f); + + out = kernelFunctions.computePoly6(glm::vec3(0.1f, 0.1f, 0.1f)); + BOOST_CHECK_CLOSE(out, 0, 0.001); +}*/ + +BOOST_AUTO_TEST_CASE(KernelFunctions_Spiky_InsideKernelTest) { + KernelFunctions kernelFunctions; + double out; + + kernelFunctions.setKernelRadius(0.1f); + + out = kernelFunctions.computeSpiky(glm::vec3(0.05f, 0.05f, 0.05f)); + BOOST_CHECK_CLOSE(out, 11.4819, 0.01); + + out = kernelFunctions.computeSpiky(glm::vec3(0.03f, 0.04f, 0.05f)); + BOOST_CHECK_CLOSE(out, 119.969, 0.01); + + out = kernelFunctions.computeSpiky(glm::vec3(0.003f, 0.002f, 0.001f)); + BOOST_CHECK_CLOSE(out, 4258.50, 0.01); +} + +/*BOOST_AUTO_TEST_CASE(KernelFunctions_Spiky_OutsideKernelTest) { + KernelFunctions kernelFunctions; + double out; + + kernelFunctions.setKernelRadius(0.1f); + + out = kernelFunctions.computeSpiky(glm::vec3(0.1f, 0.1f, 0.1f)); + BOOST_CHECK_CLOSE(out, 0, 0.001); +}*/ + +BOOST_AUTO_TEST_CASE(KernelFunctions_SpikyGradient_InsideKernelTest) { + KernelFunctions kernelFunctions; + glm::vec3 out; + + kernelFunctions.setKernelRadius(0.1f); + + out = kernelFunctions.computeSpikyGradient(glm::vec3(0.05f, 0.05f, 0.05f)); + BOOST_CHECK_CLOSE(out.x, -1484.40, 0.01); + out = kernelFunctions.computeSpikyGradient(glm::vec3(0.05f, 0.05f, 0.05f)); + BOOST_CHECK_CLOSE(out.y, -1484.40, 0.01); + out = kernelFunctions.computeSpikyGradient(glm::vec3(0.05f, 0.05f, 0.05f)); + BOOST_CHECK_CLOSE(out.z, -1484.40, 0.01); + + out = kernelFunctions.computeSpikyGradient(glm::vec3(0.03f, 0.04f, 0.05f)); + BOOST_CHECK_CLOSE(out.x, -5213.35, 0.01); + out = kernelFunctions.computeSpikyGradient(glm::vec3(0.03f, 0.04f, 0.05f)); + BOOST_CHECK_CLOSE(out.y, -6951.13, 0.01); + out = kernelFunctions.computeSpikyGradient(glm::vec3(0.03f, 0.04f, 0.05f)); + BOOST_CHECK_CLOSE(out.z, -8688.91, 0.01); + + out = kernelFunctions.computeSpikyGradient(glm::vec3(0.003f, 0.002f, 0.001f)); + BOOST_CHECK_CLOSE(out.x, -106413.397, 0.01); + out = kernelFunctions.computeSpikyGradient(glm::vec3(0.003f, 0.002f, 0.001f)); + BOOST_CHECK_CLOSE(out.y, -70942.265, 0.01); + out = kernelFunctions.computeSpikyGradient(glm::vec3(0.003f, 0.002f, 0.001f)); + BOOST_CHECK_CLOSE(out.z, -35471.133, 0.01); +} + +BOOST_AUTO_TEST_CASE(KernelFunctions_Viscous_InsideKernelTest) { + KernelFunctions kernelFunctions; + double out; + + kernelFunctions.setKernelRadius(0.1f); + + out = kernelFunctions.computeViscous(glm::vec3(0.05f, 0.05f, 0.05f)); + BOOST_CHECK_CLOSE(out, 6.18494, 0.01); + + out = kernelFunctions.computeViscous(glm::vec3(0.03f, 0.04f, 0.05f)); + BOOST_CHECK_CLOSE(out, 72.4076, 0.01); + + out = kernelFunctions.computeViscous(glm::vec3(0.003f, 0.002f, 0.001f)); + BOOST_CHECK_CLOSE(out, 29517.8941, 0.01); +} + +/*BOOST_AUTO_TEST_CASE(KernelFunctions_Viscous_OutsideKernelTest) { + KernelFunctions kernelFunctions; + double out; + + kernelFunctions.setKernelRadius(0.1f); + + out = kernelFunctions.computeViscous(glm::vec3(0.1f, 0.1f, 0.1f)); + BOOST_CHECK_CLOSE(out, 0, 0.001); +}*/ + +BOOST_AUTO_TEST_CASE(KernelFunctions_ViscousLaplacian_Test) { + KernelFunctions kernelFunctions; + double out; + + kernelFunctions.setKernelRadius(0.1f); + + out = kernelFunctions.computeViscousLaplacian(glm::vec3(0.05f, 0.05f, 0.05f)); + BOOST_CHECK_CLOSE(out, 191905.05, 1); + + out = kernelFunctions.computeViscousLaplacian(glm::vec3(0.03f, 0.04f, 0.05f)); + BOOST_CHECK_CLOSE(out, 419538.32, 1); + + out = + kernelFunctions.computeViscousLaplacian(glm::vec3(0.003f, 0.002f, 0.001f)); + BOOST_CHECK_CLOSE(out, 1378799, 1); +} + +BOOST_AUTO_TEST_SUITE_END() diff --git a/unittest/testNeighborSearch.cpp b/unittest/testNeighborSearch.cpp new file mode 100644 index 00000000..3b0302a1 --- /dev/null +++ b/unittest/testNeighborSearch.cpp @@ -0,0 +1,702 @@ +// Copyright 2016 Robert Zhou +// +// testNeighborSearch.cpp +// MFluidSolver + +#include +#include + +#include + +#include "fluidSolver/sphSolver/neighborSearch.hpp" + +BOOST_AUTO_TEST_SUITE(NeighborSearchTests) + +/*******************************************/ +/****************** Naive ******************/ +/*******************************************/ + +BOOST_AUTO_TEST_CASE(NeighborSearch_Naive_BasicSuccess) { + NaiveNeighborSearch nSearch(1.f); + + SPHParticle a(glm::vec3(1.f)); + SPHParticle b(glm::vec3(1.f, 1.5f, 1.5f)); + + nSearch.addParticle(&a); + nSearch.addParticle(&b); + + // Check a's neighbors + nSearch.findNeighbors(&a); + + BOOST_CHECK_EQUAL(1, a.neighbors()->size()); + if (a.neighbors()->size() > 0) + BOOST_CHECK_EQUAL(&b, a.neighbors()->at(0)); + + // Check b's neighbors + nSearch.findNeighbors(&b); + + BOOST_CHECK_EQUAL(1, b.neighbors()->size()); + if (b.neighbors()->size() > 0) + BOOST_CHECK_EQUAL(&a, b.neighbors()->at(0)); +} + +BOOST_AUTO_TEST_CASE(NeighborSearch_Naive_BasicFail) { + NaiveNeighborSearch nSearch(1.f); + + SPHParticle a(glm::vec3(1.f)); + SPHParticle b(glm::vec3(3.f)); + + nSearch.addParticle(&a); + nSearch.addParticle(&b); + + // Check a's neighbors + std::vector neighbors; + nSearch.findNeighbors(&a); + + BOOST_CHECK_EQUAL(0, a.neighbors()->size()); + + // Check b's neighbors + neighbors.clear(); + nSearch.findNeighbors(&b); + + BOOST_CHECK_EQUAL(0, b.neighbors()->size()); +} + +BOOST_AUTO_TEST_CASE(NeighborSearch_Naive_StressTest) { + NaiveNeighborSearch nSearch(0.2f); + std::vector particles; + + // Init particles + glm::vec3 minBound(0); + glm::vec3 maxBound(2.5); + float particleSeparation = 0.1f; + minBound += particleSeparation / 2.f; + maxBound -= particleSeparation / 2.f; + for (float i = minBound.x; i <= maxBound.x; i += particleSeparation) { + for (float j = minBound.y; j <= maxBound.y; j += particleSeparation) { + for (float k = minBound.z; k <= maxBound.z; k += particleSeparation) { + SPHParticle *p = new SPHParticle(glm::vec3(i, j, k)); + particles.push_back(p); + nSearch.addParticle(p); + } + } + } + + std::clock_t startTime = std::clock(); + for (SPHParticle *p : particles) { + nSearch.findNeighbors(p); + } + clock_t endTime = clock(); + + double timeInSeconds = + (endTime - startTime) / static_cast(CLOCKS_PER_SEC); + std::cout << "INFO: Naive neighbor search took " << + timeInSeconds << " seconds for " << + particles.size() << " particles\n" << std::endl; +} + +/*************************************************/ +/****************** UniformGrid ******************/ +/*************************************************/ + +BOOST_AUTO_TEST_CASE(NeighborSearch_UniformGrid_SameCell) { + UniformGridNeighborSearch nSearch(0.5f, glm::vec3(0), glm::vec3(1), 0.5f); + + SPHParticle a(glm::vec3(0.25f)); + SPHParticle b(glm::vec3(0.2f)); + + nSearch.addParticle(&a); + nSearch.addParticle(&b); + + // Check a's neighbors + nSearch.findNeighbors(&a); + + BOOST_CHECK_EQUAL(1, a.neighbors()->size()); + if (a.neighbors()->size() > 0) + BOOST_CHECK_EQUAL(&b, a.neighbors()->at(0)); + + // Check b's neighbors + nSearch.findNeighbors(&b); + + BOOST_CHECK_EQUAL(1, b.neighbors()->size()); + if (b.neighbors()->size() > 0) + BOOST_CHECK_EQUAL(&a, b.neighbors()->at(0)); +} + +BOOST_AUTO_TEST_CASE(NeighborSearch_UniformGrid_DifferentCellStillNeighborsX) { + UniformGridNeighborSearch nSearch(0.5f, glm::vec3(0), glm::vec3(1), 0.5f); + + SPHParticle a(glm::vec3(0.25f)); + SPHParticle b(glm::vec3(0.6f, 0.25f, 0.25f)); + + nSearch.addParticle(&a); + nSearch.addParticle(&b); + + // Check a's neighbors + nSearch.findNeighbors(&a); + + BOOST_CHECK_EQUAL(1, a.neighbors()->size()); + if (a.neighbors()->size() > 0) + BOOST_CHECK_EQUAL(&b, a.neighbors()->at(0)); + + // Check b's neighbors + nSearch.findNeighbors(&b); + + BOOST_CHECK_EQUAL(1, b.neighbors()->size()); + if (b.neighbors()->size() > 0) + BOOST_CHECK_EQUAL(&a, b.neighbors()->at(0)); +} + +BOOST_AUTO_TEST_CASE(NeighborSearch_UniformGrid_DifferentCellStillNeighborsY) { + UniformGridNeighborSearch nSearch(0.5f, glm::vec3(0), glm::vec3(1), 0.5f); + + SPHParticle a(glm::vec3(0.25f)); + SPHParticle b(glm::vec3(0.25f, 0.6f, 0.25f)); + + nSearch.addParticle(&a); + nSearch.addParticle(&b); + + // Check a's neighbors + nSearch.findNeighbors(&a); + + BOOST_CHECK_EQUAL(1, a.neighbors()->size()); + if (a.neighbors()->size() > 0) + BOOST_CHECK_EQUAL(&b, a.neighbors()->at(0)); + + // Check b's neighbors + nSearch.findNeighbors(&b); + + BOOST_CHECK_EQUAL(1, b.neighbors()->size()); + if (b.neighbors()->size() > 0) + BOOST_CHECK_EQUAL(&a, b.neighbors()->at(0)); +} + +BOOST_AUTO_TEST_CASE(NeighborSearch_UniformGrid_DifferentCellStillNeighborsZ) { + UniformGridNeighborSearch nSearch(0.5f, glm::vec3(0), glm::vec3(1), 0.5f); + + SPHParticle a(glm::vec3(0.25f)); + SPHParticle b(glm::vec3(0.25f, 0.25f, 0.6f)); + + nSearch.addParticle(&a); + nSearch.addParticle(&b); + + // Check a's neighbors + nSearch.findNeighbors(&a); + + BOOST_CHECK_EQUAL(1, a.neighbors()->size()); + if (a.neighbors()->size() > 0) + BOOST_CHECK_EQUAL(&b, a.neighbors()->at(0)); + + // Check b's neighbors + nSearch.findNeighbors(&b); + + BOOST_CHECK_EQUAL(1, b.neighbors()->size()); + if (b.neighbors()->size() > 0) + BOOST_CHECK_EQUAL(&a, b.neighbors()->at(0)); +} + +BOOST_AUTO_TEST_CASE(NeighborSearch_UniformGrid_DifferentCellNotNeighborsX) { + UniformGridNeighborSearch nSearch(0.2f, glm::vec3(0), glm::vec3(1), 0.2f); + + SPHParticle a(glm::vec3(0.1f)); + SPHParticle b(glm::vec3(0.9f, 0.1f, 0.1f)); + + nSearch.addParticle(&a); + nSearch.addParticle(&b); + + // Check a's neighbors + nSearch.findNeighbors(&a); + + BOOST_CHECK_EQUAL(0, a.neighbors()->size()); + + // Check b's neighbors + nSearch.findNeighbors(&b); + + BOOST_CHECK_EQUAL(0, b.neighbors()->size()); +} + +BOOST_AUTO_TEST_CASE(NeighborSearch_UniformGrid_DifferentCellNotNeighborsY) { + UniformGridNeighborSearch nSearch(0.2f, glm::vec3(0), glm::vec3(1), 0.2f); + + SPHParticle a(glm::vec3(0.1f)); + SPHParticle b(glm::vec3(0.1f, 0.9f, 0.1f)); + + nSearch.addParticle(&a); + nSearch.addParticle(&b); + + // Check a's neighbors + nSearch.findNeighbors(&a); + + BOOST_CHECK_EQUAL(0, a.neighbors()->size()); + + // Check b's neighbors + nSearch.findNeighbors(&b); + + BOOST_CHECK_EQUAL(0, b.neighbors()->size()); +} + +BOOST_AUTO_TEST_CASE(NeighborSearch_UniformGrid_DifferentCellNotNeighborsZ) { + UniformGridNeighborSearch nSearch(0.2f, glm::vec3(0), glm::vec3(1), 0.2f); + + SPHParticle a(glm::vec3(0.1f)); + SPHParticle b(glm::vec3(0.1f, 0.1f, 0.9f)); + + nSearch.addParticle(&a); + nSearch.addParticle(&b); + + // Check a's neighbors + nSearch.findNeighbors(&a); + + BOOST_CHECK_EQUAL(0, a.neighbors()->size()); + + // Check b's neighbors + nSearch.findNeighbors(&b); + + BOOST_CHECK_EQUAL(0, b.neighbors()->size()); +} + +BOOST_AUTO_TEST_CASE(NeighborSearch_UniformGrid_DifferentCellStillNeighborsPositiveOffsetX) { + UniformGridNeighborSearch nSearch(0.5f, glm::vec3(2), glm::vec3(3), 0.5f); + + SPHParticle a(glm::vec3(2.25f)); + SPHParticle b(glm::vec3(2.6f, 2.25f, 2.25f)); + + nSearch.addParticle(&a); + nSearch.addParticle(&b); + + // Check a's neighbors + nSearch.findNeighbors(&a); + + BOOST_CHECK_EQUAL(1, a.neighbors()->size()); + if (a.neighbors()->size() > 0) + BOOST_CHECK_EQUAL(&b, a.neighbors()->at(0)); + + // Check b's neighbors + nSearch.findNeighbors(&b); + + BOOST_CHECK_EQUAL(1, b.neighbors()->size()); + if (b.neighbors()->size() > 0) + BOOST_CHECK_EQUAL(&a, b.neighbors()->at(0)); +} + +BOOST_AUTO_TEST_CASE(NeighborSearch_UniformGrid_DifferentCellStillNeighborsNegativeOffsetX) { + UniformGridNeighborSearch nSearch(0.5f, glm::vec3(-3), glm::vec3(-2), 0.5f); + + SPHParticle a(glm::vec3(-2.75f)); + SPHParticle b(glm::vec3(-2.4f, -2.75f, -2.75f)); + + nSearch.addParticle(&a); + nSearch.addParticle(&b); + + // Check a's neighbors + nSearch.findNeighbors(&a); + + BOOST_CHECK_EQUAL(1, a.neighbors()->size()); + if (a.neighbors()->size() > 0) + BOOST_CHECK_EQUAL(&b, a.neighbors()->at(0)); + + // Check b's neighbors + nSearch.findNeighbors(&b); + + BOOST_CHECK_EQUAL(1, b.neighbors()->size()); + if (b.neighbors()->size() > 0) + BOOST_CHECK_EQUAL(&a, b.neighbors()->at(0)); +} + +BOOST_AUTO_TEST_CASE(NeighborSearch_UniformGrid_StressTest) { + glm::vec3 minBound(0); + glm::vec3 maxBound(2.5); + + UniformGridNeighborSearch nSearch(0.2f, minBound, maxBound, 0.2f); + std::vector particles; + + // Init particles + float particleSeparation = 0.1f; + minBound += particleSeparation / 2.f; + maxBound -= particleSeparation / 2.f; + for (float i = minBound.x; i <= maxBound.x; i += particleSeparation) { + for (float j = minBound.y; j <= maxBound.y; j += particleSeparation) { + for (float k = minBound.z; k <= maxBound.z; k += particleSeparation) { + SPHParticle *p = new SPHParticle(glm::vec3(i, j, k)); + particles.push_back(p); + nSearch.addParticle(p); + } + } + } + + std::clock_t startTime = std::clock(); + for (SPHParticle *p : particles) { + nSearch.findNeighbors(p); + } + clock_t endTime = clock(); + + double timeInSeconds = + (endTime - startTime) / static_cast(CLOCKS_PER_SEC); + std:: cout << "INFO: Uniform grid neighbor search took" << + timeInSeconds << " seconds for " << + particles.size() << " particles\n" << std::endl; +} + +/************************************************************/ +/****************** IndexSortedUniformGrid ******************/ +/************************************************************/ + +BOOST_AUTO_TEST_CASE(NeighborSearch_IndexSortedUniformGrid_SameCell) { + std::vector pList; + pList.push_back(SPHParticle(glm::vec3(0.25f))); + pList.push_back(SPHParticle(glm::vec3(0.2f))); + + SPHParticle *a = &(pList.at(0)); + SPHParticle *b = &(pList.at(1)); + + IndexSortedUniformGridNeighborSearch nSearch( + 0.5f, glm::vec3(0), glm::vec3(1), 0.5f, &pList, false); + nSearch.isuGrid->resetAndFillCells(); + + // Check a's neighbors + nSearch.findNeighbors(a); + + BOOST_CHECK_EQUAL(1, a->neighbors()->size()); + if (a->neighbors()->size() > 0) + BOOST_CHECK_EQUAL(b, a->neighbors()->at(0)); + + // Check b's neighbors + nSearch.findNeighbors(b); + + BOOST_CHECK_EQUAL(1, b->neighbors()->size()); + if (b->neighbors()->size() > 0) + BOOST_CHECK_EQUAL(a, b->neighbors()->at(0)); +} + +BOOST_AUTO_TEST_CASE(NeighborSearch_IndexSortedUniformGrid_DifferentCellStillNeighborsX) { + std::vector pList; + pList.push_back(SPHParticle(glm::vec3(0.25f))); + pList.push_back(SPHParticle(glm::vec3(0.6f, 0.25f, 0.25f))); + + SPHParticle *a = &(pList.at(0)); + SPHParticle *b = &(pList.at(1)); + + IndexSortedUniformGridNeighborSearch nSearch( + 0.5f, glm::vec3(0), glm::vec3(1), 0.5f, &pList, false); + nSearch.isuGrid->resetAndFillCells(); + + // Check a's neighbors + nSearch.findNeighbors(a); + + BOOST_CHECK_EQUAL(1, a->neighbors()->size()); + if (a->neighbors()->size() > 0) + BOOST_CHECK_EQUAL(b, a->neighbors()->at(0)); + + // Check b's neighbors + nSearch.findNeighbors(b); + + BOOST_CHECK_EQUAL(1, b->neighbors()->size()); + if (b->neighbors()->size() > 0) + BOOST_CHECK_EQUAL(a, b->neighbors()->at(0)); +} + +BOOST_AUTO_TEST_CASE(NeighborSearch_IndexSortedUniformGrid_DifferentCellStillNeighborsY) { + std::vector pList; + pList.push_back(SPHParticle(glm::vec3(0.25f))); + pList.push_back(SPHParticle(glm::vec3(0.26f, 0.6f, 0.25f))); + + SPHParticle *a = &(pList.at(0)); + SPHParticle *b = &(pList.at(1)); + + IndexSortedUniformGridNeighborSearch nSearch( + 0.5f, glm::vec3(0), glm::vec3(1), 0.5f, &pList, false); + nSearch.isuGrid->resetAndFillCells(); + + // Check a's neighbors + nSearch.findNeighbors(a); + + BOOST_CHECK_EQUAL(1, a->neighbors()->size()); + if (a->neighbors()->size() > 0) + BOOST_CHECK_EQUAL(b, a->neighbors()->at(0)); + + // Check b's neighbors + nSearch.findNeighbors(b); + + BOOST_CHECK_EQUAL(1, b->neighbors()->size()); + if (b->neighbors()->size() > 0) + BOOST_CHECK_EQUAL(a, b->neighbors()->at(0)); +} + +BOOST_AUTO_TEST_CASE(NeighborSearch_IndexSortedUniformGrid_DifferentCellStillNeighborsZ) { + std::vector pList; + pList.push_back(SPHParticle(glm::vec3(0.25f))); + pList.push_back(SPHParticle(glm::vec3(0.26f, 0.25f, 0.6f))); + + SPHParticle *a = &(pList.at(0)); + SPHParticle *b = &(pList.at(1)); + + IndexSortedUniformGridNeighborSearch nSearch( + 0.5f, glm::vec3(0), glm::vec3(1), 0.5f, &pList, false); + nSearch.isuGrid->resetAndFillCells(); + + // Check a's neighbors + nSearch.findNeighbors(a); + + BOOST_CHECK_EQUAL(1, a->neighbors()->size()); + if (a->neighbors()->size() > 0) + BOOST_CHECK_EQUAL(b, a->neighbors()->at(0)); + + // Check b's neighbors + nSearch.findNeighbors(b); + + BOOST_CHECK_EQUAL(1, b->neighbors()->size()); + if (b->neighbors()->size() > 0) + BOOST_CHECK_EQUAL(a, b->neighbors()->at(0)); +} + +BOOST_AUTO_TEST_CASE(NeighborSearch_IndexSortedUniformGrid_DifferentCellNotNeighborsX) { + std::vector pList; + pList.push_back(SPHParticle(glm::vec3(0.1f))); + pList.push_back(SPHParticle(glm::vec3(0.9f, 0.1f, 0.1f))); + + SPHParticle *a = &(pList.at(0)); + SPHParticle *b = &(pList.at(1)); + + IndexSortedUniformGridNeighborSearch nSearch( + 0.2f, glm::vec3(0), glm::vec3(1), 0.2f, &pList, false); + nSearch.isuGrid->resetAndFillCells(); + + // Check a's neighbors + nSearch.findNeighbors(a); + + BOOST_CHECK_EQUAL(0, a->neighbors()->size()); + + // Check b's neighbors + nSearch.findNeighbors(b); + + BOOST_CHECK_EQUAL(0, b->neighbors()->size()); +} + +BOOST_AUTO_TEST_CASE(NeighborSearch_IndexSortedUniformGrid_DifferentCellNotNeighborsY) { + std::vector pList; + pList.push_back(SPHParticle(glm::vec3(0.1f))); + pList.push_back(SPHParticle(glm::vec3(0.1f, 0.9f, 0.1f))); + + SPHParticle *a = &(pList.at(0)); + SPHParticle *b = &(pList.at(1)); + + IndexSortedUniformGridNeighborSearch nSearch( + 0.2f, glm::vec3(0), glm::vec3(1), 0.2f, &pList, false); + nSearch.isuGrid->resetAndFillCells(); + + // Check a's neighbors + nSearch.findNeighbors(a); + + BOOST_CHECK_EQUAL(0, a->neighbors()->size()); + + // Check b's neighbors + nSearch.findNeighbors(b); + + BOOST_CHECK_EQUAL(0, b->neighbors()->size()); +} + +BOOST_AUTO_TEST_CASE(NeighborSearch_IndexSortedUniformGrid_DifferentCellNotNeighborsZ) { + std::vector pList; + pList.push_back(SPHParticle(glm::vec3(0.1f))); + pList.push_back(SPHParticle(glm::vec3(0.1f, 0.1f, 0.9f))); + + SPHParticle *a = &(pList.at(0)); + SPHParticle *b = &(pList.at(1)); + + IndexSortedUniformGridNeighborSearch nSearch( + 0.2f, glm::vec3(0), glm::vec3(1), 0.2f, &pList, false); + nSearch.isuGrid->resetAndFillCells(); + + // Check a's neighbors + nSearch.findNeighbors(a); + + BOOST_CHECK_EQUAL(0, a->neighbors()->size()); + + // Check b's neighbors + nSearch.findNeighbors(b); + + BOOST_CHECK_EQUAL(0, b->neighbors()->size()); +} + +/*************************************************************/ +/****************** ZIndexSortedUniformGrid ******************/ +/*************************************************************/ + +BOOST_AUTO_TEST_CASE(NeighborSearch_ZIndexSortedUniformGrid_SameCell) { + std::vector pList; + pList.push_back(SPHParticle(glm::vec3(0.25f))); + pList.push_back(SPHParticle(glm::vec3(0.2f))); + + SPHParticle *a = &(pList.at(0)); + SPHParticle *b = &(pList.at(1)); + + IndexSortedUniformGridNeighborSearch nSearch( + 0.5f, glm::vec3(0), glm::vec3(1), 0.5f, &pList, true); + nSearch.isuGrid->resetAndFillCells(); + + // Check a's neighbors + nSearch.findNeighbors(a); + + BOOST_CHECK_EQUAL(1, a->neighbors()->size()); + if (a->neighbors()->size() > 0) + BOOST_CHECK_EQUAL(b, a->neighbors()->at(0)); + + // Check b's neighbors + nSearch.findNeighbors(b); + + BOOST_CHECK_EQUAL(1, b->neighbors()->size()); + if (b->neighbors()->size() > 0) + BOOST_CHECK_EQUAL(a, b->neighbors()->at(0)); +} + +BOOST_AUTO_TEST_CASE(NeighborSearch_ZIndexSortedUniformGrid_DifferentCellStillNeighborsX) { + std::vector pList; + pList.push_back(SPHParticle(glm::vec3(0.25f))); + pList.push_back(SPHParticle(glm::vec3(0.6f, 0.25f, 0.25f))); + + SPHParticle *a = &(pList.at(0)); + SPHParticle *b = &(pList.at(1)); + + IndexSortedUniformGridNeighborSearch nSearch( + 0.5f, glm::vec3(0), glm::vec3(1), 0.5f, &pList, true); + nSearch.isuGrid->resetAndFillCells(); + + // Check a's neighbors + nSearch.findNeighbors(a); + + BOOST_CHECK_EQUAL(1, a->neighbors()->size()); + if (a->neighbors()->size() > 0) + BOOST_CHECK_EQUAL(b, a->neighbors()->at(0)); + + // Check b's neighbors + nSearch.findNeighbors(b); + + BOOST_CHECK_EQUAL(1, b->neighbors()->size()); + if (b->neighbors()->size() > 0) + BOOST_CHECK_EQUAL(a, b->neighbors()->at(0)); +} + +BOOST_AUTO_TEST_CASE(NeighborSearch_ZIndexSortedUniformGrid_DifferentCellStillNeighborsY) { + std::vector pList; + pList.push_back(SPHParticle(glm::vec3(0.25f))); + pList.push_back(SPHParticle(glm::vec3(0.26f, 0.6f, 0.25f))); + + SPHParticle *a = &(pList.at(0)); + SPHParticle *b = &(pList.at(1)); + + IndexSortedUniformGridNeighborSearch nSearch( + 0.5f, glm::vec3(0), glm::vec3(1), 0.5f, &pList, true); + nSearch.isuGrid->resetAndFillCells(); + + // Check a's neighbors + nSearch.findNeighbors(a); + + BOOST_CHECK_EQUAL(1, a->neighbors()->size()); + if (a->neighbors()->size() > 0) + BOOST_CHECK_EQUAL(b, a->neighbors()->at(0)); + + // Check b's neighbors + nSearch.findNeighbors(b); + + BOOST_CHECK_EQUAL(1, b->neighbors()->size()); + if (b->neighbors()->size() > 0) + BOOST_CHECK_EQUAL(a, b->neighbors()->at(0)); +} + +BOOST_AUTO_TEST_CASE(NeighborSearch_ZIndexSortedUniformGrid_DifferentCellStillNeighborsZ) { + std::vector pList; + pList.push_back(SPHParticle(glm::vec3(0.25f))); + pList.push_back(SPHParticle(glm::vec3(0.26f, 0.25f, 0.6f))); + + SPHParticle *a = &(pList.at(0)); + SPHParticle *b = &(pList.at(1)); + + IndexSortedUniformGridNeighborSearch nSearch( + 0.5f, glm::vec3(0), glm::vec3(1), 0.5f, &pList, true); + nSearch.isuGrid->resetAndFillCells(); + + // Check a's neighbors + nSearch.findNeighbors(a); + + BOOST_CHECK_EQUAL(1, a->neighbors()->size()); + if (a->neighbors()->size() > 0) + BOOST_CHECK_EQUAL(b, a->neighbors()->at(0)); + + // Check b's neighbors + nSearch.findNeighbors(b); + + BOOST_CHECK_EQUAL(1, b->neighbors()->size()); + if (b->neighbors()->size() > 0) + BOOST_CHECK_EQUAL(a, b->neighbors()->at(0)); +} + +BOOST_AUTO_TEST_CASE(NeighborSearch_ZIndexSortedUniformGrid_DifferentCellNotNeighborsX) { + std::vector pList; + pList.push_back(SPHParticle(glm::vec3(0.1f))); + pList.push_back(SPHParticle(glm::vec3(0.9f, 0.1f, 0.1f))); + + SPHParticle *a = &(pList.at(0)); + SPHParticle *b = &(pList.at(1)); + + IndexSortedUniformGridNeighborSearch nSearch( + 0.2f, glm::vec3(0), glm::vec3(1), 0.2f, &pList, true); + nSearch.isuGrid->resetAndFillCells(); + + // Check a's neighbors + nSearch.findNeighbors(a); + + BOOST_CHECK_EQUAL(0, a->neighbors()->size()); + + // Check b's neighbors + nSearch.findNeighbors(b); + + BOOST_CHECK_EQUAL(0, b->neighbors()->size()); +} + +BOOST_AUTO_TEST_CASE(NeighborSearch_ZIndexSortedUniformGrid_DifferentCellNotNeighborsY) { + std::vector pList; + pList.push_back(SPHParticle(glm::vec3(0.1f))); + pList.push_back(SPHParticle(glm::vec3(0.1f, 0.9f, 0.1f))); + + SPHParticle *a = &(pList.at(0)); + SPHParticle *b = &(pList.at(1)); + + IndexSortedUniformGridNeighborSearch nSearch( + 0.2f, glm::vec3(0), glm::vec3(1), 0.2f, &pList, true); + nSearch.isuGrid->resetAndFillCells(); + + // Check a's neighbors + nSearch.findNeighbors(a); + + BOOST_CHECK_EQUAL(0, a->neighbors()->size()); + + // Check b's neighbors + nSearch.findNeighbors(b); + + BOOST_CHECK_EQUAL(0, b->neighbors()->size()); +} + +BOOST_AUTO_TEST_CASE(NeighborSearch_ZIndexSortedUniformGrid_DifferentCellNotNeighborsZ) { + std::vector pList; + pList.push_back(SPHParticle(glm::vec3(0.1f))); + pList.push_back(SPHParticle(glm::vec3(0.1f, 0.1f, 0.9f))); + + SPHParticle *a = &(pList.at(0)); + SPHParticle *b = &(pList.at(1)); + + IndexSortedUniformGridNeighborSearch nSearch( + 0.2f, glm::vec3(0), glm::vec3(1), 0.2f, &pList, true); + nSearch.isuGrid->resetAndFillCells(); + + // Check a's neighbors + nSearch.findNeighbors(a); + + BOOST_CHECK_EQUAL(0, a->neighbors()->size()); + + // Check b's neighbors + nSearch.findNeighbors(b); + + BOOST_CHECK_EQUAL(0, b->neighbors()->size()); +} + +BOOST_AUTO_TEST_SUITE_END() diff --git a/unittest/testZCurve.cpp b/unittest/testZCurve.cpp new file mode 100644 index 00000000..425ec48a --- /dev/null +++ b/unittest/testZCurve.cpp @@ -0,0 +1,58 @@ +// +// testZCurve.cpp +// MFluidSolver + +#include + +#include "fluidSolver/zCurve.hpp" + +BOOST_AUTO_TEST_SUITE(ZCurveTests) + +BOOST_AUTO_TEST_CASE(ZCurve_NegativeXFail) { + ZCurve z; + BOOST_TEST_MESSAGE("Expecting negative bounds exception"); + BOOST_CHECK_THROW(z.initWithMax(glm::ivec3(-1, 10, 10)), + ZCurveNegativeBoundsException); +} + +BOOST_AUTO_TEST_CASE(ZCurve_NegativeYFail) { + ZCurve z; + BOOST_TEST_MESSAGE("Expecting negative bounds exception"); + BOOST_CHECK_THROW(z.initWithMax(glm::ivec3(10, -1, 10)), + ZCurveNegativeBoundsException); +} + +BOOST_AUTO_TEST_CASE(ZCurve_NegativeZFail) { + ZCurve z; + BOOST_TEST_MESSAGE("Expecting negative bounds exception"); + BOOST_CHECK_THROW(z.initWithMax(glm::ivec3(10, 10, -1)), + ZCurveNegativeBoundsException); +} + +BOOST_AUTO_TEST_CASE(ZCurve_LargeXFail) { + ZCurve z; + BOOST_TEST_MESSAGE("Expecting exceeded bounds exception"); + BOOST_CHECK_THROW(z.initWithMax(glm::ivec3(2050, 10, 10)), + ZCurveTooLargeException); +} + +BOOST_AUTO_TEST_CASE(ZCurve_LargeYFail) { + ZCurve z; + BOOST_TEST_MESSAGE("Expecting exceeded bounds exception"); + BOOST_CHECK_THROW(z.initWithMax(glm::ivec3(10, 2050, 10)), + ZCurveTooLargeException); +} + +BOOST_AUTO_TEST_CASE(ZCurve_LargeZFail) { + ZCurve z; + BOOST_TEST_MESSAGE("Expecting exceeded bounds exception"); + BOOST_CHECK_THROW(z.initWithMax(glm::ivec3(10, 10, 2050)), + ZCurveTooLargeException); +} + +BOOST_AUTO_TEST_CASE(ZCurve_Small1) { + ZCurve z; + BOOST_CHECK_EQUAL(z.initWithMax(glm::ivec3(1, 2, 3)), 53); +} + +BOOST_AUTO_TEST_SUITE_END()