diff --git a/CMakeLists.txt b/CMakeLists.txt index cf54de96e..5739ec1f2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -41,6 +41,7 @@ option(WITH_BOX2D_EXAMPLE "Build Box2D integration example" OFF) option(WITH_BULLET_EXAMPLE "Build Bullet integration example (requires the BulletIntegration library)" OFF) cmake_dependent_option(WITH_CUBEMAP_EXAMPLE "Build CubeMap example (requires some JPEG importer plugin)" OFF "NOT MAGNUM_TARGET_GLES" OFF) option(WITH_DART_EXAMPLE "Build DART integration example (requires the DartIntegration library)" OFF) +option(WITH_FLUIDSIMULATION2D_EXAMPLE "Build 2D Fluid Simulation example (requires the ImGui integration)" OFF) option(WITH_FLUIDSIMULATION3D_EXAMPLE "Build 3D Fluid Simulation example (requires the ImGui integration)" OFF) option(WITH_IMGUI_EXAMPLE "Build ImGui example" OFF) option(WITH_LEAPMOTION_EXAMPLE "Build LeapMotion example" OFF) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 7a98cc148..7d7165605 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -62,6 +62,10 @@ if(WITH_DART_EXAMPLE) add_subdirectory(dart) endif() +if(WITH_FLUIDSIMULATION2D_EXAMPLE) + add_subdirectory(fluidsimulation2d) +endif() + if(WITH_FLUIDSIMULATION3D_EXAMPLE) add_subdirectory(fluidsimulation3d) endif() diff --git a/src/fluidsimulation2d/CMakeLists.txt b/src/fluidsimulation2d/CMakeLists.txt new file mode 100644 index 000000000..fe2d3a629 --- /dev/null +++ b/src/fluidsimulation2d/CMakeLists.txt @@ -0,0 +1,87 @@ +# +# This file is part of Magnum. +# +# Original authors — credit is appreciated but not required: +# +# 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 — +# Vladimír Vondruš +# 2019 — Nghia Truong +# +# This is free and unencumbered software released into the public domain. +# +# Anyone is free to copy, modify, publish, use, compile, sell, or distribute +# this software, either in source code form or as a compiled binary, for any +# purpose, commercial or non-commercial, and by any means. +# +# In jurisdictions that recognize copyright laws, the author or authors of +# this software dedicate any and all copyright interest in the software to +# the public domain. We make this dedication for the benefit of the public +# at large and to the detriment of our heirs and successors. We intend this +# dedication to be an overt act of relinquishment in perpetuity of all +# present and future rights to this software under copyright law. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL +# THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER +# IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN +# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. +# + +cmake_minimum_required(VERSION 3.4) +project(magnum-fluidsimulation2d) + +# Add module path in case this is project root +if(PROJECT_SOURCE_DIR STREQUAL CMAKE_SOURCE_DIR) + set(CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/../../modules/" ${CMAKE_MODULE_PATH}) +endif() + +find_package(Corrade REQUIRED Main) +find_package(Magnum REQUIRED + GL + MeshTools + Primitives + SceneGraph + Shaders + Sdl2Application) +find_package(MagnumIntegration REQUIRED ImGui) + +set_directory_properties(PROPERTIES CORRADE_USE_PEDANTIC_FLAGS ON) + +corrade_add_resource(FluidSimulation2D_RESOURCES resources.conf) + +add_executable(magnum-fluidsimulation2d WIN32 + FluidSimulation2DExample.cpp + DataStructures/Array2X.h + DataStructures/MathHelpers.h + DataStructures/PCGSolver.h + DataStructures/SDFObject.h + DataStructures/SparseMatrix.h + DrawableObjects/FlatShadeObject2D.h + DrawableObjects/ParticleGroup2D.h + DrawableObjects/ParticleGroup2D.cpp + DrawableObjects/WireframeObject2D.h + FluidSolver/SolverData.h + FluidSolver/ApicSolver2D.h + FluidSolver/ApicSolver2D.cpp + Shaders/ParticleSphereShader2D.h + Shaders/ParticleSphereShader2D.cpp + ${FluidSimulation2D_RESOURCES}) +target_link_libraries(magnum-fluidsimulation2d PRIVATE + Corrade::Main + Magnum::Application + Magnum::GL + Magnum::Magnum + Magnum::MeshTools + Magnum::Primitives + Magnum::SceneGraph + Magnum::Shaders + MagnumIntegration::ImGui) +target_include_directories(magnum-fluidsimulation2d PRIVATE + ${CMAKE_CURRENT_SOURCE_DIR} + ${CMAKE_CURRENT_BINARY_DIR}) + +install(TARGETS magnum-fluidsimulation2d DESTINATION ${MAGNUM_BINARY_INSTALL_DIR}) + +# Make the executable a default target to build & run in Visual Studio +set_property(DIRECTORY ${PROJECT_SOURCE_DIR} PROPERTY VS_STARTUP_PROJECT magnum-fluidsimulation2d) diff --git a/src/fluidsimulation2d/DataStructures/Array2X.h b/src/fluidsimulation2d/DataStructures/Array2X.h new file mode 100644 index 000000000..81aeda264 --- /dev/null +++ b/src/fluidsimulation2d/DataStructures/Array2X.h @@ -0,0 +1,185 @@ +#ifndef Magnum_Examples_FluidSimulation2D_Array2X_h +#define Magnum_Examples_FluidSimulation2D_Array2X_h +/* + This file is part of Magnum. + + Original authors — credit is appreciated but not required: + + 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 — + Vladimír Vondruš + 2019 — Nghia Truong + + This is free and unencumbered software released into the public domain. + + Anyone is free to copy, modify, publish, use, compile, sell, or distribute + this software, either in source code form or as a compiled binary, for any + purpose, commercial or non-commercial, and by any means. + + In jurisdictions that recognize copyright laws, the author or authors of + this software dedicate any and all copyright interest in the software to + the public domain. We make this dedication for the benefit of the public + at large and to the detriment of our heirs and successors. We intend this + dedication to be an overt act of relinquishment in perpetuity of all + present and future rights to this software under copyright law. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER + IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN + CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include +#include +#include +#include + +#include "MathHelpers.h" + +namespace Magnum { namespace Examples { +template +class Array2X { +public: + Array2X() = default; + Array2X(std::size_t nx, std::size_t ny) : _size{nx, ny}, _data(nx * ny) {} + Array2X(std::size_t nx, std::size_t ny, const T& value) : _size{nx, ny}, _data(nx * ny, value) {} + + Array2X& operator=(const Array2X& other) { + if(other.size_x() != size_x() || other.size_y() != size_y()) { + Fatal{} << "Copy array with different size!"; + } + std::memcpy(data(), other.data(), count() * sizeof(T)); + return *this; + } + + /******************************* Accessors *******************************/ + template + const T& operator()(IntType i, IntType j) const { + assert(i >= 0 + && j >= 0 + && static_cast(i) < _size[0] + && static_cast(j) < _size[1]); + return _data[i + _size[0] * j]; + } + + template + T& operator()(IntType i, IntType j) { + assert(i >= 0 + && j >= 0 + && static_cast(i) < _size[0] + && static_cast(j) < _size[1]); + return _data[i + _size[0] * j]; + } + + template + const T& operator()(const Math::Vector2& coord) const { return (*this)(coord[0], coord[1]); } + + template + T& operator()(const Math::Vector2& coord) { return (*this)(coord[0], coord[1]); } + + const T* data() const { return _data.data(); } + T* data() { return _data.data(); } + + std::size_t size_x() const { return _size[0]; } + std::size_t size_y() const { return _size[1]; } + std::size_t count() const { return _data.size(); } + + /******************************* Modifiers *******************************/ + void assign(const T& value) { _data.assign(_data.size(), value); } + void set_zero() { _data.assign(_data.size(), T(0)); } + + template + void resize(IntType nx, IntType ny) { + _size[0] = static_cast(nx); + _size[1] = static_cast(ny); + _data.resize(_size[0] * _size[1]); + } + + template + void resize(IntType nx, IntType ny, const T& value) { + _size = { static_cast(nx), static_cast(ny) }; + _data.resize(_size[0] * _size[1], value); + } + + void swapContent(Array2X& other) { + /* Only allow to swap content of array having the same sizes */ + if(other.size_x() != size_x() || other.size_y() != size_y()) { + Fatal{} << "Swap content of arrays having different sizes!"; + } + _data.swap(other._data); + } + + /*************************** Data manimupations ***************************/ + template + void loop1D(Function&& func) const { + for(std::size_t i = 0, iend = count(); i < iend; ++i) { + func(i); + } + } + + template + void loop2D(Function&& func) const { + for(std::size_t j = 0; j < size_y(); ++j) { + for(std::size_t i = 0; i < size_x(); ++i) { + func(i, j); + } + } + } + + T interpolateValue(const Math::Vector2& point) const { + Int i, j; + T fx, fy; + barycentric(point[0], i, fx, 0, Int(size_x())); + barycentric(point[1], j, fy, 0, Int(size_y())); + T v00 = (*this)(i, j); + T v10 = (*this)(i + 1, j); + T v01 = (*this)(i, j + 1); + T v11 = (*this)(i + 1, j + 1); + return bilerp(v00, v10, v01, v11, fx, fy); + } + + Math::Vector2 affineInterpolateValue(const Math::Vector2& point) const { + Int i, j; + T fx, fy; + barycentric(point[0], i, fx, 0, Int(size_x())); + barycentric(point[1], j, fy, 0, Int(size_y())); + T v00 = (*this)(i, j); + T v10 = (*this)(i + 1, j); + T v01 = (*this)(i, j + 1); + T v11 = (*this)(i + 1, j + 1); + return bilerpGradient(v00, v10, v01, v11, fx, fy); + } + + Math::Vector2 interpolateGradient(const Math::Vector2& point) const { + Int i, j; + T fx, fy; + barycentric(point[0], i, fx, 0, Int(size_x())); + barycentric(point[1], j, fy, 0, Int(size_y())); + T v00 = (*this)(i, j); + T v01 = (*this)(i, j + 1); + T v10 = (*this)(i + 1, j); + T v11 = (*this)(i + 1, j + 1); + T ddy0 = (v01 - v00); + T ddy1 = (v11 - v10); + T ddx0 = (v10 - v00); + T ddx1 = (v11 - v01); + + Math::Vector2 grad( + Math::lerp(ddx0, ddx1, fy), + Math::lerp(ddy0, ddy1, fx) + ); + const auto magSqr = grad.dot(); + if(magSqr > T(1e-20)) { + grad /= std::sqrt(magSqr); + } + return grad; + } + +private: + std::size_t _size[2] { 0, 0 }; + std::vector _data; +}; +} } + +#endif diff --git a/src/fluidsimulation2d/DataStructures/MathHelpers.h b/src/fluidsimulation2d/DataStructures/MathHelpers.h new file mode 100644 index 000000000..557031a9a --- /dev/null +++ b/src/fluidsimulation2d/DataStructures/MathHelpers.h @@ -0,0 +1,104 @@ +#ifndef Magnum_Examples_FluidSimulation2D_MathHelpers_h +#define Magnum_Examples_FluidSimulation2D_MathHelpers_h +/* + This file is part of Magnum. + + Original authors — credit is appreciated but not required: + + 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 — + Vladimír Vondruš + 2019 — Nghia Truong + + This is free and unencumbered software released into the public domain. + + Anyone is free to copy, modify, publish, use, compile, sell, or distribute + this software, either in source code form or as a compiled binary, for any + purpose, commercial or non-commercial, and by any means. + + In jurisdictions that recognize copyright laws, the author or authors of + this software dedicate any and all copyright interest in the software to + the public domain. We make this dedication for the benefit of the public + at large and to the detriment of our heirs and successors. We intend this + dedication to be an overt act of relinquishment in perpetuity of all + present and future rights to this software under copyright law. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER + IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN + CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include +#include +#include +#include + +namespace Magnum { namespace Examples { +template +inline T fractionInside(T phi_left, T phi_right) { + if(phi_left < 0 && phi_right < 0) { + return T(1); + } + if(phi_left < 0 && phi_right >= 0) { + return phi_left / (phi_left - phi_right); + } + if(phi_left >= 0 && phi_right < 0) { + return phi_right / (phi_right - phi_left); + } else { + return T(0); + } +} + +template +inline void barycentric(T x, int& i, T& f, int i_low, int i_high) { + T s = std::floor(x); + i = static_cast(s); + if(i < i_low) { + i = i_low; + f = 0; + } else if(i > i_high - 2) { + i = i_high - 2; + f = 1; + } else { + f = T(x - s); + } +} + +template +inline T bilerp(const T& v00, const T& v10, + const T& v01, const T& v11, + T fx, T fy) { + return Math::lerp(Math::lerp(v00, v10, fx), + Math::lerp(v01, v11, fx), + fy); +} + +template +inline Math::Vector2 bilerpGradient(const T& v00, const T& v10, + const T& v01, const T& v11, + T fx, T fy) { + const Math::Vector2 f00(fy - T(1), fx - T(1)); + const Math::Vector2 f10(T(1) - fy, -fx); + const Math::Vector2 f01(-fy, T(1) - fx); + const Math::Vector2 f11(fy, fx); + return f00 * v00 + f10 * v10 + f01 * v01 + f11 * v11; +} + +template +inline T smoothKernel(T r2, T h2) { + const auto t = T(1) - r2 / h2; + const auto t_exp3 = t * t * t; + return Math::max(t_exp3, T(0)); +} + +template +inline T linearKernel(const Math::Vector2& d, T hInv) { + const auto tx = T(1) - std::abs(d.x() * hInv); + const auto ty = T(1) - std::abs(d.y() * hInv); + return Math::max(tx * ty, T(0)); +} +} } + +#endif diff --git a/src/fluidsimulation2d/DataStructures/PCGSolver.h b/src/fluidsimulation2d/DataStructures/PCGSolver.h new file mode 100644 index 000000000..ad6ae7d18 --- /dev/null +++ b/src/fluidsimulation2d/DataStructures/PCGSolver.h @@ -0,0 +1,276 @@ +#ifndef Magnum_Examples_FluidSimulation2D_PCGSolver_h +#define Magnum_Examples_FluidSimulation2D_PCGSolver_h +/* + This file is part of Magnum. + + Original authors — credit is appreciated but not required: + + 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 — + Vladimír Vondruš + 2019 — Nghia Truong + + This is free and unencumbered software released into the public domain. + + Anyone is free to copy, modify, publish, use, compile, sell, or distribute + this software, either in source code form or as a compiled binary, for any + purpose, commercial or non-commercial, and by any means. + + In jurisdictions that recognize copyright laws, the author or authors of + this software dedicate any and all copyright interest in the software to + the public domain. We make this dedication for the benefit of the public + at large and to the detriment of our heirs and successors. We intend this + dedication to be an overt act of relinquishment in perpetuity of all + present and future rights to this software under copyright law. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER + IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN + CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include + +#include "SparseMatrix.h" + +namespace Magnum { namespace Examples { +template +class PCGSolver { +public: + PCGSolver(T toleranceFactor_ = T(1e-10), uint32_t maxIterations_ = 1000) : + _toleranceFactor{toleranceFactor_}, _maxIterations{maxIterations_} {} + + bool solve(SparseMatrix& matrix, const std::vector& rhs, std::vector& result) { + uint32_t rows = matrix.size; + if(rows == 0) { + return false; + } + if(_m.size() != rows) { + _m.resize(rows); + _s.resize(rows); + _z.resize(rows); + _r.resize(rows); + } + + _r = rhs; + _lastResidual = maxAbs(_r); + if(!(_lastResidual > 0)) { + _lastIterationCount = 0; + return true; + } + + formPreconditioner(matrix); + applyPreconditioner(_r, _z); + T rho = dotProduct(_z, _r); + if(!(rho > 0) || rho != rho) { + _lastIterationCount = 0; + return false; + } + + _s = _z; + matrix.compressData(); /* must prepare compact data for Matrix * Vector operation */ + + const T tolerance = _toleranceFactor * _lastResidual; + uint32_t iter { 0 }; + for(; iter < _maxIterations; ++iter) { + matrix.multiply(_s, _z); + const T alpha = rho / dotProduct(_s, _z); + addScaled(alpha, _s, result); + addScaled(-alpha, _z, _r); + _lastResidual = maxAbs(_r); + if(_lastResidual < tolerance) { + _lastIterationCount = iter + 1; + return true; + } + applyPreconditioner(_r, _z); + const T rho_new = dotProduct(_z, _r); + const T beta = rho_new / rho; + addScaled(beta, _s, _z); + _s.swap(_z); + rho = rho_new; + } + + /* Failed to converge */ + _lastIterationCount = _maxIterations; + return false; + } + + /* API to query last solve */ + uint32_t lastIterationCount() const { return _lastIterationCount; } + T lastResidual() const { return _lastResidual; } + +private: + void formPreconditioner(const SparseMatrix& matrix) { + static constexpr T s_modification_parameter = T(0.97); + + const auto size = matrix.size; + _precond.reset(size); + + /* Copy data from the matrix */ + for(uint32_t i = 0; i < size; ++i) { + _precond.colStartIdx[i] = static_cast(_precond.colIndices.size()); + const auto& indices = matrix.rowIndices[i]; + const auto& values = matrix.rowValues[i]; + for(uint32_t j = 0, jend = indices.size(); j < jend; ++j) { + if(indices[j] > i) { + _precond.colIndices.push_back(indices[j]); + _precond.colValues.push_back(values[j]); + } else if(indices[j] == i) { + _precond.invdiag[i] = values[j]; + } + } + } + _precond.colStartIdx[size] = static_cast(_precond.colIndices.size()); + + for(uint32_t k = 0; k < size; ++k) { + auto invdiag = _precond.invdiag[k]; + if(invdiag == T(0)) { + continue; /* null row/column */ + } + invdiag = T(1) / std::sqrt(invdiag); + _precond.invdiag[k] = invdiag; + + const auto pStart = _precond.colStartIdx[k]; + const auto pEnd = _precond.colStartIdx[k + 1]; + for(uint32_t p = pStart; p < pEnd; ++p) { + _precond.colValues[p] *= invdiag; + } + + /* Process the lower elements of column k */ + for(uint32_t p = pStart; p < pEnd; ++p) { + const auto j = _precond.colIndices[p]; + const auto multiplier = _precond.colValues[p]; + T missing = 0; + uint32_t a = pStart; + + uint32_t b = 0; + while(a < pEnd && _precond.colIndices[a] < j) { + while(b < matrix.rowIndices[j].size()) { + if(matrix.rowIndices[j][b] < _precond.colIndices[a]) { + ++b; + } else if(matrix.rowIndices[j][b] == _precond.colIndices[a]) { + break; + } else { + missing += _precond.colValues[a]; + break; + } + } + ++a; + } + + invdiag = _precond.invdiag[j]; + if(a < pEnd && _precond.colIndices[a] == j) { + invdiag -= multiplier * _precond.colValues[a]; + } + ++a; + + b = _precond.colStartIdx[j]; + const auto jEnd = _precond.colStartIdx[j + 1]; + while(a < pEnd && b < jEnd) { + if(_precond.colIndices[b] < _precond.colIndices[a]) { + ++b; + } else if(_precond.colIndices[b] == _precond.colIndices[a]) { + _precond.colValues[b] -= multiplier * _precond.colValues[a]; + ++a; + ++b; + } else { + missing += _precond.colValues[a++]; + } + } + + while(a < pEnd) { + missing += _precond.colValues[a++]; + } + + _precond.invdiag[j] = invdiag - s_modification_parameter * multiplier * missing; + } + } + } + + void applyPreconditioner(const std::vector& rhs, std::vector& result) const { + const uint32_t rows = _precond.rows; + + /* Solve L * result = rhs */ + result = rhs; + for(uint32_t i = 0; i < rows; ++i) { + const auto colStart = _precond.colStartIdx[i]; + const auto colEnd = _precond.colStartIdx[i + 1]; + const auto tmp = result[i] * _precond.invdiag[i]; + result[i] = tmp; + for(uint32_t j = colStart; j < colEnd; ++j) { + result[_precond.colIndices[j]] -= _precond.colValues[j] * tmp; + } + } + + /* solve L^T * result = result */ + uint32_t i = rows; + do { + --i; + const auto colStart = _precond.colStartIdx[i]; + const auto colEnd = _precond.colStartIdx[i + 1]; + auto tmp = result[i]; + for(uint32_t j = colStart; j < colEnd; ++j) { + tmp -= _precond.colValues[j] * result[_precond.colIndices[j]]; + } + tmp *= _precond.invdiag[i]; + result[i] = tmp; + } while (i != 0); + } + + T dotProduct(const std::vector& x, const std::vector& y) const { + T sum = 0; + for(std::size_t i = 0; i < x.size(); ++i) { + sum += x[i] * y[i]; + } + return sum; + } + + T maxAbs(const std::vector& x) const { + T maxVal = 0; + for(std::size_t i = 0; i < x.size(); ++i) { + const auto absVal = std::abs(x[i]); + if(absVal > maxVal) { + maxVal = absVal; + } + } + return maxVal; + } + + void addScaled(T alpha, const std::vector& x, std::vector& y) const { + for(std::size_t i = 0; i < x.size(); ++i) { + y[i] += alpha * x[i]; + } + } + + /* Solver parameters */ + const uint32_t _maxIterations; + const T _toleranceFactor; + + /* Preconditioner - triangular matrix */ + struct { + uint32_t rows; + std::vector invdiag; /* inverse of diagonal elements */ + std::vector colValues; /* values below the diagonal, listed column by column */ + std::vector colIndices; /* a list of all row indices, for each column in turn */ + std::vector colStartIdx; /* where each column begins in rowindex (plus an extra entry at the end, of #nonzeros) */ + + void reset(uint32_t rows_) { + rows = rows_; + invdiag.assign(rows, T(0)); /* important: must set zero */ + colValues.resize(0); + colIndices.resize(0); + colStartIdx.resize(rows + 1); + } + } _precond; + + /* Solver temporary variables */ + std::vector _m, _z, _s, _r; + + /* Status of last solve */ + uint32_t _lastIterationCount { 0 }; + T _lastResidual { 0 }; +}; +} } + +#endif diff --git a/src/fluidsimulation2d/DataStructures/SDFObject.h b/src/fluidsimulation2d/DataStructures/SDFObject.h new file mode 100644 index 000000000..1414308f8 --- /dev/null +++ b/src/fluidsimulation2d/DataStructures/SDFObject.h @@ -0,0 +1,114 @@ +#ifndef Magnum_Examples_FluidSimulation2D_SDFObject_h +#define Magnum_Examples_FluidSimulation2D_SDFObject_h +/* + This file is part of Magnum. + + Original authors — credit is appreciated but not required: + + 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 — + Vladimír Vondruš + 2019 — Nghia Truong + + This is free and unencumbered software released into the public domain. + + Anyone is free to copy, modify, publish, use, compile, sell, or distribute + this software, either in source code form or as a compiled binary, for any + purpose, commercial or non-commercial, and by any means. + + In jurisdictions that recognize copyright laws, the author or authors of + this software dedicate any and all copyright interest in the software to + the public domain. We make this dedication for the benefit of the public + at large and to the detriment of our heirs and successors. We intend this + dedication to be an overt act of relinquishment in perpetuity of all + present and future rights to this software under copyright law. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER + IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN + CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include +#include +#include +#include +#include + +namespace Magnum { namespace Examples { +struct SDFObject { + enum class ObjectType { + /* Basic primitives */ + Circle = 0, + Box, + + /* Boolean operations */ + Intersection, + Subtraction, + Union + }; + + SDFObject() : center{0, 0}, radii{1}, type{ObjectType::Circle}, negativeInside{true} {} + + SDFObject(const Vector2& center_, Float radius_, ObjectType type_, bool negativeInside_ = true) : + SDFObject(center_, Vector2{ radius_, 0 }, type_, negativeInside_) {} + + SDFObject(const Vector2& center_, const Vector2& radii_, ObjectType type_, bool negativeInside_ = true) : + center{center_}, radii{radii_}, type{type_}, negativeInside{negativeInside_} { + if(type != ObjectType::Circle + && type != ObjectType::Box) { + Fatal{} << "Invalid object type"; + } + } + + SDFObject(SDFObject* obj1_, SDFObject* obj2_, ObjectType type_) : + type{type_}, obj1{obj1_}, obj2{obj2_} { + if(type != ObjectType::Intersection + && type != ObjectType::Subtraction + && type != ObjectType::Union) { + Fatal{} << "Invalid boolean operation"; + } + } + + Float signedDistance(const Vector2& pos) const { + switch(type) { + case ObjectType::Circle: { + const auto dist = (pos - center).length() - radii[0]; + return negativeInside ? dist : -dist; + } + case ObjectType::Box: { + const auto dx = std::abs(pos[0] - center[0]) - radii[0]; + const auto dy = std::abs(pos[1] - center[1]) - radii[1]; + Float dist { 0 }; + if(dx < 0 && dy < 0) { + dist = Math::max(dx, dy); + } else { + const Float dax = Math::max(dx, 0.0f); + const Float day = Math::max(dy, 0.0f); + dist = std::sqrt(dax * dax + day * day); + } + return negativeInside ? dist : -dist; + } + + case ObjectType::Intersection: + return Math::max(obj1->signedDistance(pos), obj2->signedDistance(pos)); + case ObjectType::Subtraction: + return Math::max(obj1->signedDistance(pos), -obj2->signedDistance(pos)); + case ObjectType::Union: + return Math::min(obj1->signedDistance(pos), obj2->signedDistance(pos)); + } + return 0; + } + + Vector2 center; + Vector2 radii; + ObjectType type; + bool negativeInside { true }; + + Containers::Pointer obj1 { nullptr }; + Containers::Pointer obj2 { nullptr }; +}; +} } + +#endif diff --git a/src/fluidsimulation2d/DataStructures/SparseMatrix.h b/src/fluidsimulation2d/DataStructures/SparseMatrix.h new file mode 100644 index 000000000..9643bc6b2 --- /dev/null +++ b/src/fluidsimulation2d/DataStructures/SparseMatrix.h @@ -0,0 +1,120 @@ +#ifndef Magnum_Examples_FluidSimulation2D_SparseMatrix_h +#define Magnum_Examples_FluidSimulation2D_SparseMatrix_h +/* + This file is part of Magnum. + + Original authors — credit is appreciated but not required: + + 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 — + Vladimír Vondruš + 2019 — Nghia Truong + + This is free and unencumbered software released into the public domain. + + Anyone is free to copy, modify, publish, use, compile, sell, or distribute + this software, either in source code form or as a compiled binary, for any + purpose, commercial or non-commercial, and by any means. + + In jurisdictions that recognize copyright laws, the author or authors of + this software dedicate any and all copyright interest in the software to + the public domain. We make this dedication for the benefit of the public + at large and to the detriment of our heirs and successors. We intend this + dedication to be an overt act of relinquishment in perpetuity of all + present and future rights to this software under copyright law. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER + IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN + CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include +#include + +namespace Magnum { namespace Examples { +template +struct SparseMatrix { + explicit SparseMatrix(uint32_t size_ = 0) { resize(size_); } + + void resize(std::size_t size_) { + size = size_; + rowIndices.resize(size); + rowValues.resize(size); + rowStartIdx.resize(size + 1); + } + + void clear() { + for(auto& indices : rowIndices) { + indices.resize(0); + } + for(auto& vals : rowValues) { + vals.resize(0); + } + } + + template + void addToElement(IntType i_, IntType j_, U inc_val_) { + const auto i = static_cast(i_); + const auto j = static_cast(j_); + const auto val = static_cast(inc_val_); + auto& currentRowIndices = rowIndices[i]; + + auto iter = std::lower_bound(currentRowIndices.begin(), currentRowIndices.end(), j); + if((iter != currentRowIndices.end()) && (*iter == j)) { + const auto k = std::distance(currentRowIndices.begin(), iter); + rowValues[i][k] += val; + } else { + auto insert_sorted_vector = + [](std::vector& vec, uint32_t item) { + return vec.insert(std::upper_bound(vec.begin(), vec.end(), item), item); + }; + iter = insert_sorted_vector(currentRowIndices, j); + const auto k = std::distance(currentRowIndices.begin(), iter); + rowValues[i].insert(rowValues[i].begin() + k, val); + } + } + + void compressData() { + rowStartIdx[0] = 0; + for(uint32_t i = 0; i < size; ++i) { + rowStartIdx[i + 1] = rowStartIdx[i] + static_cast(rowIndices[i].size()); + } + + auto copyData = [](auto& dst, const auto& src) { + dst.resize(0); + for(const auto& vec : src) { + dst.insert(dst.end(), vec.begin(), vec.end()); + } + }; + + copyData(compactIndices, rowIndices); + copyData(compactValues, rowValues); + } + + void multiply(const std::vector& x, std::vector& result) const { + result.resize(size); + for(uint32_t i = 0; i < size; ++i) { + T tmp = 0; + for(uint32_t j = rowStartIdx[i], jend = rowStartIdx[i + 1]; j < jend; ++j) { + tmp += compactValues[j] * x[compactIndices[j]]; + } + result[i] = tmp; + } + } + + uint32_t size; + + /* Sparse matrix data */ + std::vector> rowIndices; + std::vector> rowValues; + + /* Compact data */ + std::vector compactValues; + std::vector compactIndices; + std::vector rowStartIdx; +}; +} } + +#endif diff --git a/src/fluidsimulation2d/DrawableObjects/FlatShadeObject2D.h b/src/fluidsimulation2d/DrawableObjects/FlatShadeObject2D.h new file mode 100644 index 000000000..b8d5d8b4b --- /dev/null +++ b/src/fluidsimulation2d/DrawableObjects/FlatShadeObject2D.h @@ -0,0 +1,67 @@ +#ifndef Magnum_Examples_DrawableObjects_FlatShadeObject2D_h +#define Magnum_Examples_DrawableObjects_FlatShadeObject2D_h +/* + This file is part of Magnum. + + Original authors — credit is appreciated but not required: + + 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 — + Vladimír Vondruš + 2019 — Nghia Truong + + This is free and unencumbered software released into the public domain. + + Anyone is free to copy, modify, publish, use, compile, sell, or distribute + this software, either in source code form or as a compiled binary, for any + purpose, commercial or non-commercial, and by any means. + + In jurisdictions that recognize copyright laws, the author or authors of + this software dedicate any and all copyright interest in the software to + the public domain. We make this dedication for the benefit of the public + at large and to the detriment of our heirs and successors. We intend this + dedication to be an overt act of relinquishment in perpetuity of all + present and future rights to this software under copyright law. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER + IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN + CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include +#include +#include +#include +#include +#include + +namespace Magnum { namespace Examples { +using Object2D = SceneGraph::Object; + +class FlatShadeObject2D : public SceneGraph::Drawable2D { +public: + explicit FlatShadeObject2D(Object2D& object, Shaders::Flat2D& shader, const Color3& color, GL::Mesh& mesh, SceneGraph::DrawableGroup2D* const drawables) : + SceneGraph::Drawable2D{object, drawables}, _shader(shader), _color(color), _mesh(mesh) {} + + void draw(const Matrix3& transformation, SceneGraph::Camera2D& camera) override { + if(_bEnabled) { + _shader.setColor(_color) + .setTransformationProjectionMatrix(camera.projectionMatrix() * transformation); + _mesh.draw(_shader); + } + } + + FlatShadeObject2D& setColor(const Color3& color) { _color = color; return *this; } + FlatShadeObject2D& setEnabled(bool bEnabled) { _bEnabled = bEnabled; return *this; } + +private: + Shaders::Flat2D& _shader; + Color3 _color; + GL::Mesh& _mesh; + bool _bEnabled { true }; +}; +} } + +#endif diff --git a/src/fluidsimulation2d/DrawableObjects/ParticleGroup2D.cpp b/src/fluidsimulation2d/DrawableObjects/ParticleGroup2D.cpp new file mode 100644 index 000000000..8072dcf3a --- /dev/null +++ b/src/fluidsimulation2d/DrawableObjects/ParticleGroup2D.cpp @@ -0,0 +1,77 @@ +/* + This file is part of Magnum. + + Original authors — credit is appreciated but not required: + + 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 — + Vladimír Vondruš + 2019 — Nghia Truong + + This is free and unencumbered software released into the public domain. + + Anyone is free to copy, modify, publish, use, compile, sell, or distribute + this software, either in source code form or as a compiled binary, for any + purpose, commercial or non-commercial, and by any means. + + In jurisdictions that recognize copyright laws, the author or authors of + this software dedicate any and all copyright interest in the software to + the public domain. We make this dedication for the benefit of the public + at large and to the detriment of our heirs and successors. We intend this + dedication to be an overt act of relinquishment in perpetuity of all + present and future rights to this software under copyright law. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER + IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN + CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include "DrawableObjects/ParticleGroup2D.h" + +#include +#include +#include +#include +#include +#include +#include + +namespace Magnum { namespace Examples { +using namespace Math::Literals; + +ParticleGroup2D::ParticleGroup2D(const std::vector& points, float particleRadius) : + _points{points}, + _particleRadius{particleRadius}, + _meshParticles{GL::MeshPrimitive::Points} { + _meshParticles.addVertexBuffer(_bufferParticles, 0, Shaders::Generic2D::Position{}); + _particleShader.reset(new ParticleSphereShader2D); +} + +ParticleGroup2D& ParticleGroup2D::draw(Containers::Pointer& camera, Int screenHeight, Int projectionHeight) { + if(_points.empty()) { return *this; } + + if(_dirty) { + Containers::ArrayView data(reinterpret_cast(&_points[0]), _points.size() * 2); + _bufferParticles.setData(data); + _meshParticles.setCount(static_cast(_points.size())); + _dirty = false; + } + + (*_particleShader) + /* particle data */ + .setNumParticles(static_cast(_points.size())) + .setParticleRadius(_particleRadius) + /* sphere render data */ + .setColorMode(_colorMode) + .setColor(_color) + /* view/prj matrices and size */ + .setViewProjectionMatrix(camera->projectionMatrix() * camera->cameraMatrix()) + .setScreenHeight(screenHeight) + .setDomainHeight(projectionHeight); + + _meshParticles.draw(*_particleShader); + return *this; +} +} } diff --git a/src/fluidsimulation2d/DrawableObjects/ParticleGroup2D.h b/src/fluidsimulation2d/DrawableObjects/ParticleGroup2D.h new file mode 100644 index 000000000..0c427a6e6 --- /dev/null +++ b/src/fluidsimulation2d/DrawableObjects/ParticleGroup2D.h @@ -0,0 +1,75 @@ +#ifndef Magnum_Examples_FluidSimulation2D_Particle2DGroup_h +#define Magnum_Examples_FluidSimulation2D_Particle2DGroup_h +/* + This file is part of Magnum. + + Original authors — credit is appreciated but not required: + + 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 — + Vladimír Vondruš + 2019 — Nghia Truong + + This is free and unencumbered software released into the public domain. + + Anyone is free to copy, modify, publish, use, compile, sell, or distribute + this software, either in source code form or as a compiled binary, for any + purpose, commercial or non-commercial, and by any means. + + In jurisdictions that recognize copyright laws, the author or authors of + this software dedicate any and all copyright interest in the software to + the public domain. We make this dedication for the benefit of the public + at large and to the detriment of our heirs and successors. We intend this + dedication to be an overt act of relinquishment in perpetuity of all + present and future rights to this software under copyright law. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER + IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN + CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include +#include +#include +#include +#include +#include + +#include "Shaders/ParticleSphereShader2D.h" + +namespace Magnum { namespace Examples { +class ParticleGroup2D { +public: + explicit ParticleGroup2D(const std::vector& points, float particleRadius); + + ParticleGroup2D& draw(Containers::Pointer& camera, Int screenHeight, Int projectionHeight); + + bool isDirty() const { return _dirty; } + ParticleGroup2D& setDirty() { _dirty = true; return *this; } + + Float particleRadius() const { return _particleRadius; } + ParticleGroup2D& setParticleRadius(Float radius) { _particleRadius = radius; return *this; } + + ParticleSphereShader2D::ColorMode colorMode() const { return _colorMode; } + ParticleGroup2D& setColorMode(ParticleSphereShader2D::ColorMode colorMode) { _colorMode = colorMode; return *this; } + + const Color3& color() const { return _color; } + ParticleGroup2D& setColor(const Color3& color) { _color = color; return *this; } + +private: + const std::vector& _points; + bool _dirty { false }; + + Float _particleRadius { 1.0f }; + ParticleSphereShader2D::ColorMode _colorMode { ParticleSphereShader2D::ColorMode::RampColorById }; + Color3 _color{ 0.1f }; + + GL::Buffer _bufferParticles; + GL::Mesh _meshParticles; + Containers::Pointer _particleShader; +}; +} } + +#endif diff --git a/src/fluidsimulation2d/DrawableObjects/WireframeObject2D.h b/src/fluidsimulation2d/DrawableObjects/WireframeObject2D.h new file mode 100644 index 000000000..1362bfbf1 --- /dev/null +++ b/src/fluidsimulation2d/DrawableObjects/WireframeObject2D.h @@ -0,0 +1,65 @@ +#ifndef Magnum_Examples_DrawableObjects_WireframeObject2D_h +#define Magnum_Examples_DrawableObjects_WireframeObject2D_h +/* + This file is part of Magnum. + + Original authors — credit is appreciated but not required: + + 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 — + Vladimír Vondruš + 2019 — Nghia Truong + + This is free and unencumbered software released into the public domain. + + Anyone is free to copy, modify, publish, use, compile, sell, or distribute + this software, either in source code form or as a compiled binary, for any + purpose, commercial or non-commercial, and by any means. + + In jurisdictions that recognize copyright laws, the author or authors of + this software dedicate any and all copyright interest in the software to + the public domain. We make this dedication for the benefit of the public + at large and to the detriment of our heirs and successors. We intend this + dedication to be an overt act of relinquishment in perpetuity of all + present and future rights to this software under copyright law. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER + IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN + CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include +#include +#include +#include +#include + +#include "DrawableObjects/FlatShadeObject2D.h" + +namespace Magnum { namespace Examples { +using Scene2D = SceneGraph::Scene; + +class WireframeObject2D { +public: + explicit WireframeObject2D(Scene2D* const scene, SceneGraph::DrawableGroup2D* const drawableGroup, GL::Mesh&& mesh) : + _mesh(std::move(mesh)) { + _obj2D.reset(new Object2D{ scene }); + _flatShader = Shaders::Flat2D{}; + _drawableObj.reset(new FlatShadeObject2D{ *_obj2D, _flatShader, Color3{ 1.0f }, _mesh, drawableGroup }); + } + + WireframeObject2D& setColor(const Color3& color) { _drawableObj->setColor(color); return *this; } + WireframeObject2D& setTransformation(const Matrix3& matrix) { _obj2D->setTransformation(matrix); return *this; } + WireframeObject2D& setEnabled(bool bEnabled) { _drawableObj->setEnabled(bEnabled); return *this; } + +protected: + GL::Mesh _mesh{ NoCreate }; + Shaders::Flat2D _flatShader{ NoCreate }; + Containers::Pointer _obj2D; + Containers::Pointer _drawableObj; +}; +} } + +#endif diff --git a/src/fluidsimulation2d/FluidSimulation2DExample.cpp b/src/fluidsimulation2d/FluidSimulation2DExample.cpp new file mode 100644 index 000000000..996bbb8be --- /dev/null +++ b/src/fluidsimulation2d/FluidSimulation2DExample.cpp @@ -0,0 +1,480 @@ +/* + This file is part of Magnum. + + Original authors — credit is appreciated but not required: + + 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 — + Vladimír Vondruš + 2019 — Nghia Truong + + This is free and unencumbered software released into the public domain. + + Anyone is free to copy, modify, publish, use, compile, sell, or distribute + this software, either in source code form or as a compiled binary, for any + purpose, commercial or non-commercial, and by any means. + + In jurisdictions that recognize copyright laws, the author or authors of + this software dedicate any and all copyright interest in the software to + the public domain. We make this dedication for the benefit of the public + at large and to the detriment of our heirs and successors. We intend this + dedication to be an overt act of relinquishment in perpetuity of all + present and future rights to this software under copyright law. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER + IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN + CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "DrawableObjects/ParticleGroup2D.h" +#include "DrawableObjects/WireframeObject2D.h" +#include "FluidSolver/ApicSolver2D.h" + +namespace Magnum { namespace Examples { +class FluidSimulation2DExample : public Platform::Application { +public: + explicit FluidSimulation2DExample(const Arguments& arguments); + virtual ~FluidSimulation2DExample() = default; + +protected: + void viewportEvent(ViewportEvent& event) override; + void keyPressEvent(KeyEvent& event) override; + void keyReleaseEvent(KeyEvent& event) override; + void mousePressEvent(MouseEvent& event) override; + void mouseReleaseEvent(MouseEvent& event) override; + void mouseMoveEvent(MouseMoveEvent& event) override; + void mouseScrollEvent(MouseScrollEvent& event) override; + void textInputEvent(TextInputEvent& event) override; + void drawEvent() override; + + /* Fluid simulation helper functions */ + void resetSimulation(); + Vector2 windowPos2WorldPos(const Vector2i& winPos); + + /* Window control */ + void showMenu(); + bool _showMenu = true; + ImGuiIntegration::Context _imGuiContext{ NoCreate }; + + /* Scene and drawable group must be constructed before camera and other + drawble objects */ + Containers::Pointer _scene; + Containers::Pointer _drawableGroup; + + /* Camera helpers */ + Containers::Pointer _objCamera; + Containers::Pointer _camera; + + /* Fluid simulation system */ + Containers::Pointer _fluidSolver; + Containers::Pointer _drawableParticles; + Containers::Pointer _drawableBoundary; + Float _speed { 2.0f }; + Float _evolvedTime { 0 }; + Int _numEmission { 0 }; + bool _bAutoEmitParticles { true }; + bool _pausedSimulation { false }; + + /* Mouse-Fluid interaction */ + Containers::Pointer _drawablePointer; + Timeline _timeline; + Vector2 _lastMousePressedWorldPos; + Float _mouseInteractionRadius { 5 }; + Float _mouseInteractionMagnitude { 5 }; + bool _bMouseInteraction { true }; +}; + +namespace { +static constexpr Float GridCellLength = 1.0f; /* length of 1 grid cell */ +static constexpr Int NumGridCellX = 100; /* number of cells in x dimension */ +static constexpr Int NumGridCellY = 100; /* number of cells in y dimension */ +static constexpr Vector2 GridStart { -50, -50 }; /* lower corner of the grid */ +static constexpr Int radiusCircleBoundary = 45; /* radius of the boundary circle */ + +/* Viewport will display this window */ +static constexpr Float ProjectionScale = 1.05f; +static constexpr Int DomainDisplayW = Int(NumGridCellX * GridCellLength * ProjectionScale); +static constexpr Int DomainDisplayH = Int(NumGridCellY * GridCellLength * ProjectionScale); + +Vector2 gridCenter() { + return Vector2(NumGridCellX, NumGridCellY) * GridCellLength * 0.5f + GridStart; +} +} + +FluidSimulation2DExample::FluidSimulation2DExample(const Arguments& arguments) : Platform::Application{arguments, NoCreate} { + /* Setup window */ + { + const Vector2 dpiScaling = this->dpiScaling({}); + Configuration conf; + conf.setTitle("Magnum 2D Fluid Simulation Example") + .setSize(conf.size(), dpiScaling) + .setWindowFlags(Configuration::WindowFlag::Resizable); + GLConfiguration glConf; + glConf.setSampleCount(dpiScaling.max() < 2.0f ? 8 : 2); + if(!tryCreate(conf, glConf)) { + create(conf, glConf.setSampleCount(0)); + } + } + + /* Setup ImGui, load a better font */ + { + ImGui::CreateContext(); + ImGui::StyleColorsDark(); + + ImFontConfig fontConfig; + fontConfig.FontDataOwnedByAtlas = false; + const Vector2 size = Vector2{ windowSize() } / dpiScaling(); + Utility::Resource rs{ "data" }; + Containers::ArrayView font = rs.getRaw("SourceSansPro-Regular.ttf"); + ImGui::GetIO().Fonts->AddFontFromMemoryTTF( + const_cast(font.data()), Int(font.size()), 16.0f * framebufferSize().x() / size.x(), &fontConfig); + + _imGuiContext = ImGuiIntegration::Context(*ImGui::GetCurrentContext(), + Vector2{ windowSize() } / dpiScaling(), windowSize(), framebufferSize()); + + /* Setup proper blending to be used by ImGui */ + GL::Renderer::setBlendFunction( + GL::Renderer::BlendFunction::SourceAlpha, + GL::Renderer::BlendFunction::OneMinusSourceAlpha); + } + + /* Setup scene objects and camera */ + { + /* Setup scene objects */ + _scene.reset(new Scene2D{}); + _drawableGroup.reset(new SceneGraph::DrawableGroup2D{}); + + /* Configure camera */ + _objCamera.reset(new Object2D{ _scene.get() }); + _objCamera->setTransformation(Matrix3::translation(gridCenter())); + + _camera.reset(new SceneGraph::Camera2D{ *_objCamera }); + _camera->setAspectRatioPolicy(SceneGraph::AspectRatioPolicy::Extend) + .setProjectionMatrix(Matrix3::projection(Vector2{ DomainDisplayW, DomainDisplayH })) + .setViewport(GL::defaultFramebuffer.viewport().size()); + } + + /* Setup fluid solver */ + { + SceneObjects* sceneObjs = new SceneObjects; + sceneObjs->emitterT0 = SDFObject(gridCenter() + Vector2(10, 10), 30, SDFObject::ObjectType::Circle); + sceneObjs->emitter = SDFObject(gridCenter() + Vector2(15, 20), 15, SDFObject::ObjectType::Circle); + sceneObjs->boundary = SDFObject(gridCenter(), radiusCircleBoundary, SDFObject::ObjectType::Circle, false); + _fluidSolver.reset(new ApicSolver2D{ GridStart, GridCellLength, NumGridCellX, NumGridCellY, sceneObjs }); + + /* Drawable particles */ + _drawableParticles.reset(new ParticleGroup2D{ _fluidSolver->particlePositions(), _fluidSolver->particleRadius() }); + _drawableParticles->setColor(Color3(85.0f / 255, 200.0f / 255, 245.0f / 255)); + + /* Drawable boundary*/ + _drawableBoundary.reset(new WireframeObject2D(_scene.get(), _drawableGroup.get(), + MeshTools::compile(Primitives::circle2DWireframe(128)))); + _drawableBoundary->setTransformation(Matrix3::scaling(Vector2{ radiusCircleBoundary + _fluidSolver->particleRadius() })); + _drawableBoundary->setColor(Color3(1, 1, 1)); + + /* Visualize mouse pointer for mouse-fluid interaction */ + _drawablePointer.reset(new WireframeObject2D(_scene.get(), _drawableGroup.get(), + MeshTools::compile(Primitives::circle2DWireframe(32)))); + _drawablePointer->setColor(Color3(0, 1, 0)); + _drawablePointer->setEnabled(false); + } + + /* Enable depth test, render particles as sprites */ + GL::Renderer::enable(GL::Renderer::Feature::DepthTest); + GL::Renderer::enable(GL::Renderer::Feature::ProgramPointSize); + + /* Start the timer, loop at 60 Hz max */ + setSwapInterval(1); + setMinimalLoopPeriod(16); + _timeline.start(); +} + +void FluidSimulation2DExample::drawEvent() { + GL::defaultFramebuffer.clear(GL::FramebufferClear::Color | GL::FramebufferClear::Depth); + _imGuiContext.newFrame(); + + /* Enable text input, if needed */ + if(ImGui::GetIO().WantTextInput && !isTextInputActive()) { + startTextInput(); + } else if(!ImGui::GetIO().WantTextInput && isTextInputActive()) { + stopTextInput(); + } + + /* Draw objects */ + { + /* Trigger drawable object to update the particles to the GPU */ + _drawableParticles->setDirty(); + _drawableParticles->draw(_camera, GL::defaultFramebuffer.viewport().size().y(), DomainDisplayH); + + /* Draw other objects (boundary mesh, pointer mesh) */ + _camera->draw(*_drawableGroup); + } + + if(!_pausedSimulation) { + static constexpr Float frameTime = 1.0f / 60.0f; + if(_evolvedTime > 1.0f) { /* pause for a while before starting simulation */ + _fluidSolver->advanceFrame(frameTime * _speed); + } + _evolvedTime += frameTime; + + /* Emit particles automatically */ + if(_bAutoEmitParticles && _evolvedTime > 10.0f) { + static Float lastTime { _evolvedTime }; + if(_evolvedTime - lastTime > 1.5f /* emit every 1.5 second */ + && _numEmission < 5) { /* emit 5 times */ + _fluidSolver->emitParticles(); + lastTime = _evolvedTime; + ++_numEmission; + } + } + } + + /* Menu for parameters */ + if(_showMenu) { showMenu(); } + + /* Update application cursor */ + _imGuiContext.updateApplicationCursor(*this); + + /* Render ImGui window */ + { + GL::Renderer::enable(GL::Renderer::Feature::Blending); + GL::Renderer::disable(GL::Renderer::Feature::FaceCulling); + GL::Renderer::disable(GL::Renderer::Feature::DepthTest); + GL::Renderer::enable(GL::Renderer::Feature::ScissorTest); + + _imGuiContext.drawFrame(); + + GL::Renderer::disable(GL::Renderer::Feature::ScissorTest); + GL::Renderer::enable(GL::Renderer::Feature::DepthTest); + GL::Renderer::enable(GL::Renderer::Feature::FaceCulling); + GL::Renderer::disable(GL::Renderer::Feature::Blending); + } + + swapBuffers(); + + /* Run next frame immediately */ + redraw(); +} + +void FluidSimulation2DExample::viewportEvent(ViewportEvent& event) { + /* Resize the main framebuffer */ + GL::defaultFramebuffer.setViewport({ {}, event.framebufferSize() }); + + /* Relayout ImGui */ + _imGuiContext.relayout(Vector2{ event.windowSize() } / event.dpiScaling(), event.windowSize(), event.framebufferSize()); + + /* Recompute the camera's projection matrix */ + _camera->setViewport(event.framebufferSize()); +} + +void FluidSimulation2DExample::keyPressEvent(Platform::Sdl2Application::KeyEvent& event) { + switch(event.key()) { + case KeyEvent::Key::E: + _fluidSolver->emitParticles(); + break; + case KeyEvent::Key::H: + _showMenu ^= true; + event.setAccepted(true); + break; + case KeyEvent::Key::R: + resetSimulation(); + event.setAccepted(true); + break; + case KeyEvent::Key::Space: + _pausedSimulation ^= true; + event.setAccepted(true); + break; + default: + if(_imGuiContext.handleKeyPressEvent(event)) { + event.setAccepted(true); + } + } +} + +void FluidSimulation2DExample::keyReleaseEvent(KeyEvent& event) { + if(_imGuiContext.handleKeyReleaseEvent(event)) { + event.setAccepted(true); + return; + } +} + +void FluidSimulation2DExample::mousePressEvent(MouseEvent& event) { + if(_imGuiContext.handleMousePressEvent(event)) { + event.setAccepted(true); + return; + } + _lastMousePressedWorldPos = windowPos2WorldPos(event.position()); + if(_bMouseInteraction) { + _timeline.nextFrame(); + _drawablePointer->setEnabled(true); + _drawablePointer->setTransformation(Matrix3::translation(_lastMousePressedWorldPos) * + Matrix3::scaling(Vector2{ _mouseInteractionRadius })); + event.setAccepted(); + } +} + +void FluidSimulation2DExample::mouseReleaseEvent(MouseEvent& event) { + if(_imGuiContext.handleMouseReleaseEvent(event)) { + event.setAccepted(true); + } + if(_bMouseInteraction) { + _drawablePointer->setEnabled(false); + event.setAccepted(); + } +} + +void FluidSimulation2DExample::mouseMoveEvent(MouseMoveEvent& event) { + if(_imGuiContext.handleMouseMoveEvent(event)) { + event.setAccepted(true); + return; + } + if(!event.buttons()) { return; } + + const Vector2 currentPos = windowPos2WorldPos(event.position()); + if(_bMouseInteraction) { + _timeline.nextFrame(); + const auto dt = _timeline.previousFrameDuration(); + if(dt > 1e-4f) { + _fluidSolver->addRepulsiveVelocity(_lastMousePressedWorldPos, currentPos, dt, + _mouseInteractionRadius, _mouseInteractionMagnitude * 0.01f); + _drawablePointer->setTransformation(Matrix3::translation(currentPos) * + Matrix3::scaling(Vector2{ _mouseInteractionRadius })); + event.setAccepted(); + } + _lastMousePressedWorldPos = currentPos; + } +} + +void FluidSimulation2DExample::mouseScrollEvent(MouseScrollEvent& event) { + const Float delta = event.offset().y(); + if(Math::abs(delta) < 1.0e-2f) { + return; + } + if(_imGuiContext.handleMouseScrollEvent(event)) { + /* Prevent scrolling the page */ + event.setAccepted(); + return; + } +} + +void FluidSimulation2DExample::textInputEvent(TextInputEvent& event) { + if(_imGuiContext.handleTextInputEvent(event)) { + event.setAccepted(true); + } +} + +void FluidSimulation2DExample::showMenu() { + ImGui::SetNextWindowPos({ 10.0f, 10.0f }, ImGuiCond_FirstUseEver); + ImGui::SetNextWindowBgAlpha(0.5f); + ImGui::Begin("Options", nullptr); + + /* General information */ + ImGui::Text("Hide/show menu: H"); + ImGui::Text("Num. particles: %d", Int(_fluidSolver->numParticles())); + ImGui::Text("Rendering: %3.2f FPS", Double(ImGui::GetIO().Framerate)); + ImGui::Spacing(); + + /* Rendering parameters */ + if(ImGui::TreeNode("Particle Rendering")) { + ImGui::PushID("Particle Rendering"); + { + constexpr const char* items[] = { "Uniform", "Ramp by ID" }; + static Int colorMode = 1; + ImGui::PushItemWidth(ImGui::GetWindowWidth() * 0.5f); + if(ImGui::Combo("Color Mode", &colorMode, items, 2)) { + _drawableParticles->setColorMode(ParticleSphereShader2D::ColorMode(colorMode)); + } + ImGui::PopItemWidth(); + if(colorMode == 0) { /* Uniform color */ + static Color3 color = _drawableParticles->color(); + if(ImGui::ColorEdit3("Color", color.data())) { + _drawableParticles->setColor(color); + } + } + } + ImGui::PopID(); + ImGui::TreePop(); + } + ImGui::Spacing(); + ImGui::Separator(); + + /* Simulation parameters */ + if(ImGui::TreeNodeEx("Simulation", ImGuiTreeNodeFlags_DefaultOpen)) { + ImGui::PushID("Simulation"); + ImGui::PushItemWidth(ImGui::GetWindowWidth() * 0.3f); + ImGui::InputFloat("Speed", &_speed); + ImGui::Checkbox("Auto emit particles 5 times", &_bAutoEmitParticles); + ImGui::PopItemWidth(); + ImGui::BeginGroup(); + ImGui::Checkbox("Mouse interaction", &_bMouseInteraction); + if(_bMouseInteraction) { + ImGui::PushItemWidth(ImGui::GetWindowWidth() * 0.5f); + ImGui::SliderFloat("Radius", &_mouseInteractionRadius, 1.0f, 10.0f); + ImGui::SliderFloat("Magnitude", &_mouseInteractionMagnitude, 1.0f, 10.0f); + ImGui::PopItemWidth(); + } + ImGui::EndGroup(); + ImGui::PopID(); + ImGui::TreePop(); + } + ImGui::Spacing(); + ImGui::Separator(); + + /* Reset */ + ImGui::Spacing(); + if(ImGui::Button("Emit Particles")) { + _fluidSolver->emitParticles(); + } + ImGui::SameLine(); + if(ImGui::Button(_pausedSimulation ? "Play Sim" : "Pause Sim")) { + _pausedSimulation ^= true; + } + ImGui::SameLine(); + if(ImGui::Button("Reset Sim")) { + resetSimulation(); + } + ImGui::End(); +} + +void FluidSimulation2DExample::resetSimulation() { + _fluidSolver->reset(); + _pausedSimulation = false; + _evolvedTime = 0; + _numEmission = 0; +} + +Vector2 FluidSimulation2DExample::windowPos2WorldPos(const Vector2i& winPos) { + /* Compute inverted model view projection matrix */ + const Matrix3 invViewProjMat = (_camera->projectionMatrix() * _camera->cameraMatrix()).inverted(); + + /* Compute the world coordinate from window coordinate */ + const Vector2i flippedPos = Vector2i(winPos.x(), GL::defaultFramebuffer.viewport().size().y() - winPos.y()); + const Vector2 ndcPos = Vector2(flippedPos) / Vector2(GL::defaultFramebuffer.viewport().size()) * Vector2(2) - Vector2(1); + const Vector3 worldPos = invViewProjMat * Vector3(ndcPos, 1); + return Vector2(worldPos.x() / worldPos.z(), worldPos.y() / worldPos.z()); +} +} } + +MAGNUM_APPLICATION_MAIN(Magnum::Examples::FluidSimulation2DExample) diff --git a/src/fluidsimulation2d/FluidSolver/ApicSolver2D.cpp b/src/fluidsimulation2d/FluidSolver/ApicSolver2D.cpp new file mode 100644 index 000000000..f5e7cd716 --- /dev/null +++ b/src/fluidsimulation2d/FluidSolver/ApicSolver2D.cpp @@ -0,0 +1,520 @@ +/* + This file is part of Magnum. + + Original authors — credit is appreciated but not required: + + 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 — + Vladimír Vondruš + 2019 — Nghia Truong + + This is free and unencumbered software released into the public domain. + + Anyone is free to copy, modify, publish, use, compile, sell, or distribute + this software, either in source code form or as a compiled binary, for any + purpose, commercial or non-commercial, and by any means. + + In jurisdictions that recognize copyright laws, the author or authors of + this software dedicate any and all copyright interest in the software to + the public domain. We make this dedication for the benefit of the public + at large and to the detriment of our heirs and successors. We intend this + dedication to be an overt act of relinquishment in perpetuity of all + present and future rights to this software under copyright law. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER + IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN + CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include "FluidSolver/ApicSolver2D.h" + +#include + +namespace Magnum { namespace Examples { +ApicSolver2D::ApicSolver2D(const Vector2& origin, Float cellSize, Int Ni, Int Nj, SceneObjects* sceneObjs) : + _objects{sceneObjs}, + _particles{cellSize}, + _grid{origin, cellSize, Ni, Nj} { + if(Ni < 1 || Nj < 1) { + Fatal{} << "Invalid grid resolution"; + } + if(!sceneObjs) { + Fatal{} << "Invalid scene object"; + } + + /* Initialize data */ + initBoundary(); + generateParticles(_objects->emitterT0, 0); +} + +/* This function should be called again every time the boundary changes */ +void ApicSolver2D::initBoundary() { + _grid.boundarySDF.loop2D( + [&](std::size_t i, std::size_t j) { + _grid.boundarySDF(i, j) = _objects->boundary.signedDistance(_grid.getWorldPos(i, j)); + }); + + /* Initialize the fluid cell weights from boundary signed distance field */ + _grid.u_weights.loop2D( + [&](std::size_t i, std::size_t j) { + _grid.u_weights(i, j) = Float(1) - fractionInside(_grid.boundarySDF(i, j + 1), _grid.boundarySDF(i, j)); + _grid.u_weights(i, j) = Math::clamp(_grid.u_weights(i, j), Float(0), Float(1)); + }); + _grid.v_weights.loop2D( + [&](std::size_t i, std::size_t j) { + _grid.v_weights(i, j) = Float(1) - fractionInside(_grid.boundarySDF(i + 1, j), _grid.boundarySDF(i, j)); + _grid.v_weights(i, j) = Math::clamp(_grid.v_weights(i, j), Float(0), Float(1)); + }); +} + +void ApicSolver2D::generateParticles(const SDFObject& sdfObj, Float initialVelocity_y) { + using Distribution = std::uniform_real_distribution; + const Float rndScale = _particles.particleRadius * 0.5f; + std::mt19937 gen((std::random_device{})()); + Distribution distr(-rndScale, rndScale); + + /* Generate new particles */ + std::vector newParticles; + _grid.fluidSDF.loop2D( + [&](std::size_t i, std::size_t j) { + const Vector2 cellCenter = _grid.getWorldPos(i + 0.5f, j + 0.5f); + for(Int k = 0; k < 2; ++k) { + const Vector2 ppos = cellCenter + Vector2(distr(gen), distr(gen)); + if(sdfObj.signedDistance(ppos) < 0) { + newParticles.push_back(ppos); + } + } + }); + + /* Insert into the system */ + _particles.addParticles(newParticles, initialVelocity_y); +} + +void ApicSolver2D::addRepulsiveVelocity(const Vector2& p0, const Vector2& p1, Float dt, Float radius, Float magnitude) { + Vector2 movingVel = p1 - p0; + const auto movingDist = movingVel.length(); + if(movingDist < _particles.particleRadius) { + return; + } + movingVel /= dt; + + Vector2 from, to; + for(std::size_t i = 0; i < 2; ++i) { + from[i] = Math::min(p0[i], p1[i]); + to[i] = Math::max(p0[i], p1[i]); + } + + const Int span = Int(radius / _grid.cellSize); + const Vector2i fromCell = _grid.getCellIdx(from) - Vector2i(span); + const Vector2i toCell = _grid.getCellIdx(to) + Vector2i(span); + const Vector2 p01 = p1 - p0; + const Float p01DistSqr = p01.dot(); + + const auto distToSegment = + [&](const Vector2& pos) { + const Float t = Math::max(0.0f, Math::min(1.0f, Math::dot(pos - p0, p01) / p01DistSqr)); + const Vector2 prj = p0 + t * p01; + return (pos - prj).length(); + }; + + for(Int j = fromCell.y(); j <= toCell.y(); ++j) { + for(Int i = fromCell.x(); i <= toCell.x(); ++i) { + if(!_grid.isValidCellIdx(i, j)) { + continue; + } + + const auto& particleIdxs = _grid.cellParticles(i, j); + for(auto p : particleIdxs) { + const auto dist = distToSegment(_particles.positions[p]); + const auto t = dist / radius; + if(t < 1.0f) { + const auto w = Math::lerp(0.0f, magnitude, t); + _particles.velocities[p] += movingVel * w; + } + } + } + } +} + +void ApicSolver2D::advanceFrame(Float frameDuration) { + Float frameTime = 0; + + while(frameTime < frameDuration) { + auto substep = timestepCFL(); + const auto remainingTime = frameDuration - frameTime; + if(frameTime + substep > frameDuration) { + substep = remainingTime; + } else if(frameTime + Float(1.5) * substep > frameDuration) { + substep = remainingTime * Float(0.5); + } + frameTime += substep; + + /* Advect particles */ + moveParticles(substep); + + /* Particles => grid */ + collectParticlesToCells(); + particleVelocity2Grid(); + + /* Update grid velocity */ + extrapolate(_grid.u, _grid.u_tmp, _grid.u_valid, _grid.u_old_valid); + extrapolate(_grid.v, _grid.v_tmp, _grid.v_valid, _grid.v_old_valid); + addGravity(substep); + computeFluidSDF(); + solvePressures(substep); + + /* Enforce boundary condition */ + constrainVelocity(); + + /* Grid => particles */ + relaxParticlePositions(substep); + gridVelocity2Particle(); + } +} + +Float ApicSolver2D::timestepCFL() const { + Float maxVel = 0; + _grid.u.loop1D([&](std::size_t i) { maxVel = Math::max(maxVel, std::abs(_grid.u.data()[i])); }); + _grid.v.loop1D([&](std::size_t i) { maxVel = Math::max(maxVel, std::abs(_grid.v.data()[i])); }); + return maxVel > 0 ? _grid.cellSize / maxVel * 3.0f : 1.0f; +} + +void ApicSolver2D::moveParticles(Float dt) { + _particles.loopAll( + [&](uint32_t p) { + const Vector2 newPos = _particles.positions[p] + _particles.velocities[p] * dt; + _particles.positions[p] = _grid.constrainBoundary(newPos); + }); +} + +void ApicSolver2D::collectParticlesToCells() { + _grid.cellParticles.loop1D([&](std::size_t i) { _grid.cellParticles.data()[i].resize(0); }); + _particles.loopAll([&](uint32_t p) { + const Vector2 ppos = _particles.positions[p]; + const Vector2i gridCoord = _grid.getValidCellIdx(ppos); + _grid.cellParticles(gridCoord).push_back(p); + }); +} + +void ApicSolver2D::particleVelocity2Grid() { + _grid.u.loop2D( + [&](std::size_t i, std::size_t j) { + Float sum_w = 0.0f; + Float sum_u = 0.0f; + const Vector2 nodePos = _grid.getWorldPos(i, j + 0.5f); + _grid.loopNeigborParticles(static_cast(i), static_cast(j), -1, 0, -1, 1, [&](uint32_t p) { + const Vector2 xpg = nodePos - _particles.positions[p]; + const auto w = linearKernel(xpg, _grid.invCellSize); + if(w > 0) { + sum_w += w; + sum_u += w * (_particles.velocities[p].x() + + Math::dot(_particles.affineMat[p][0], xpg)); + } + }); + _grid.u(i, j) = sum_w > 0 ? sum_u / sum_w : 0.0f; + _grid.u_valid(i, j) = sum_w > 0 ? 1 : 0; + }); + + _grid.v.loop2D( + [&](std::size_t i, std::size_t j) { + Float sum_w = 0.0; + Float sum_v = 0.0; + const Vector2 nodePos = _grid.getWorldPos(i + 0.5f, j); + _grid.loopNeigborParticles(static_cast(i), static_cast(j), -1, 1, -1, 0, [&](uint32_t p) { + const Vector2 xpg = nodePos - _particles.positions[p]; + const auto w = linearKernel(xpg, _grid.invCellSize); + if(w > 0) { + sum_w += w; + sum_v += w * (_particles.velocities[p].y() + + Math::dot(_particles.affineMat[p][1], xpg)); + } + }); + _grid.v(i, j) = sum_w > 0 ? sum_v / sum_w : 0.0f; + _grid.v_valid(i, j) = sum_w > 0 ? 1 : 0; + }); +} + +void ApicSolver2D::extrapolate(Array2X& grid, Array2X& tmp_grid, Array2X& valid, Array2X& old_valid) const { + tmp_grid = grid; + old_valid = valid; + + Array2X* pgrids[] = { &grid, &tmp_grid }; + Array2X* pvalids[] = { &valid, &old_valid }; + + for(int layers = 0; layers < 1; ++layers) { + auto pgrid_src = pgrids[layers & 1]; + auto pgrid_tgt = pgrids[!(layers & 1)]; + + auto pvalid_src = pvalids[layers & 1]; + auto pvalid_tgt = pvalids[!(layers & 1)]; + + grid.loop2D([&](std::size_t i, std::size_t j) { + if(i == 0 || i == grid.size_x() - 1 || + j == 0 || j == grid.size_y() - 1) { + return; + } + + Float sum = 0; + Int count = 0; + + if(!(*pvalid_src)(i, j)) { + const std::size_t rows[] = { i + 1, i - 1, i, i }; + const std::size_t cols[] = { j, j, j + 1, j - 1 }; + + for(std::size_t cell = 0; cell < 4; ++cell) { + if((*pvalid_src)(rows[cell], cols[cell])) { + sum += (*pgrid_src)(rows[cell], cols[cell]); + ++count; + } + } + + if(count > 0) { + (*pgrid_tgt)(i, j) = sum / static_cast(count); + (*pvalid_tgt)(i, j) = 1; + } + } + }); + + (*pgrid_src).swapContent(*pgrid_tgt); + (*pvalid_src).swapContent(*pvalid_tgt); + } +} + +void ApicSolver2D::addGravity(Float dt) { + _grid.v.loop2D( + [&](std::size_t i, std::size_t j) { + if(_grid.v_valid(i, j)) { + _grid.v(i, j) -= 9.81f * dt; /* gravity */ + } + }); +} + +void ApicSolver2D::computeFluidSDF() { + _grid.fluidSDF.assign(3 * _grid.cellSize); + + _particles.loopAll( + [&](uint32_t p) { + const Vector2 ppos = _particles.positions[p]; + const Vector2i gridPos = Vector2i(_grid.getGridPos(ppos) - Vector2(0.5)); + + for(int j = gridPos.y() - 2; j <= gridPos.y() + 2; ++j) { + for(int i = gridPos.x() - 2; i <= gridPos.x() + 2; ++i) { + if(!_grid.isValidCellIdx(i, j)) { + continue; + } + const Vector2 cellCenter = _grid.getWorldPos(i + 0.5f, j + 0.5f); + const Float sdfVal = (cellCenter - ppos).length() - _particles.particleRadius; + if(_grid.fluidSDF(i, j) > sdfVal) { + _grid.fluidSDF(i, j) = sdfVal; + } + } + } + }); + + _grid.fluidSDF.loop2D( + [&](std::size_t i, std::size_t j) { + const Vector2 cellCenter = _grid.getWorldPos(i + 0.5f, j + 0.5f); + const Float sdfVal = _objects->boundary.signedDistance(cellCenter); + if(_grid.fluidSDF(i, j) > sdfVal) { + _grid.fluidSDF(i, j) = sdfVal; + } + }); +} + +void ApicSolver2D::solvePressures(Float dt) { + const auto Ni = static_cast(_grid.Ni); + const auto Nj = static_cast(_grid.Nj); + const auto numCells = Ni * Nj; + + _pressureSolver.resize(numCells); + _pressureSolver.clear(); + + for(std::size_t j = 1; j < Nj - 1; ++j) { + for(std::size_t i = 1; i < Ni - 1; ++i) { + const auto row = i + Ni * j; + double rhsVal = 0; + const auto centerSDF = _grid.fluidSDF(i, j); + + if(centerSDF >= 0) { + _pressureSolver.rhs[row] = rhsVal; + continue; + } + + const Float cellsWeights[] = { + _grid.u_weights(i + 1, j), + _grid.u_weights(i, j), + _grid.v_weights(i, j + 1), + _grid.v_weights(i, j) + }; + const Float cellsSDF[] = { + _grid.fluidSDF(i + 1, j), + _grid.fluidSDF(i - 1, j), + _grid.fluidSDF(i, j + 1), + _grid.fluidSDF(i, j - 1) + }; + const Float cellsVel[] = { + -_grid.u(i + 1, j), /* minus velocity */ + _grid.u(i, j), + -_grid.v(i, j + 1), /* minus velocity */ + _grid.v(i, j) + }; + const std::size_t cols[] = { + row + 1, + row - 1, + row + Ni, + row - Ni + }; + + /* Fill-in matrix */ + for(std::size_t cell = 0; cell < 4; ++cell) { + rhsVal += static_cast(cellsWeights[cell] * cellsVel[cell]); + const auto term = cellsWeights[cell] * dt; + if(cellsSDF[cell] < 0) { + _pressureSolver.matrix.addToElement(row, row, term); + _pressureSolver.matrix.addToElement(row, cols[cell], -term); + } else { + const auto theta = Math::max(0.01f, fractionInside(centerSDF, cellsSDF[cell])); + _pressureSolver.matrix.addToElement(row, row, term / theta); + } + } + + /* Write rhs */ + _pressureSolver.rhs[row] = rhsVal; + } + } + + _pressureSolver.solve(); /* now solve the linear system for cells' pressure */ + + _grid.u.loop2D( + [&](std::size_t i, std::size_t j) { + /* Edges of the domain, or entirely in solid */ + if(i == 0 || i == _grid.u.size_x() - 1 || !(_grid.u_weights(i, j) > 0)) { + _grid.u(i, j) = 0; + return; + } + + const auto centerSDF = _grid.fluidSDF(i, j); + const auto leftSDF = _grid.fluidSDF(i - 1, j); + if(centerSDF < 0 || leftSDF < 0) { + Float theta = 1; + if(_grid.fluidSDF(i, j) >= 0 || leftSDF >= 0) { + theta = Math::max(0.01f, fractionInside(leftSDF, centerSDF)); + } + const auto row = i + j * Ni; + const auto pressure = static_cast(_pressureSolver.solution[row] - _pressureSolver.solution[row - 1]); + _grid.u(i, j) -= pressure * (dt / theta); + } + }); + + _grid.v.loop2D( + [&](std::size_t i, std::size_t j) { + /* Edges of the domain, or entirely in solid */ + if(j == 0 || j == _grid.v.size_y() - 1 || !(_grid.v_weights(i, j) > 0)) { + _grid.v(i, j) = 0; + return; + } + + const auto centerSDF = _grid.fluidSDF(i, j); + const auto bottomSDF = _grid.fluidSDF(i, j - 1); + if(centerSDF < 0 || bottomSDF < 0) { + Float theta = 1; + if(centerSDF >= 0 || bottomSDF >= 0) { + theta = Math::max(0.01f, fractionInside(bottomSDF, centerSDF)); + } + const auto row = i + j * Ni; + const auto pressure = static_cast(_pressureSolver.solution[row] - _pressureSolver.solution[row - Ni]); + _grid.v(i, j) -= pressure * (dt / theta); + } + }); +} + +void ApicSolver2D::constrainVelocity() { + _grid.u_tmp = _grid.u; + _grid.v_tmp = _grid.v; + + _grid.u.loop2D( + [&](std::size_t i, std::size_t j) { + if(_grid.u_weights(i, j) > 0) { /* not entirely in solid */ + return; + } + const Vector2 gridPos = Vector2(i, j + 0.5f); + const Vector2 normal = _grid.boundarySDF.interpolateGradient(gridPos); + Vector2 vel = _grid.velocityFromGridPos(gridPos); + Float perp_component = Math::dot(vel, normal); + vel -= perp_component * normal; + _grid.u_tmp(i, j) = vel[0]; + }); + + _grid.v.loop2D( + [&](std::size_t i, std::size_t j) { + if(_grid.v_weights(i, j) > 0) { /* not entirely in solid */ + return; + } + const Vector2 gridPos = Vector2(i + 0.5f, j); + const Vector2 normal = _grid.boundarySDF.interpolateGradient(gridPos); + Vector2 vel = _grid.velocityFromGridPos(gridPos); + Float perp_component = Math::dot(vel, normal); + vel -= perp_component * normal; + _grid.v_tmp(i, j) = vel[1]; + }); + + /* Only swap u_tmp and v_tmp after constraining both u and v */ + _grid.u.swapContent(_grid.u_tmp); + _grid.v.swapContent(_grid.v_tmp); +} + +void ApicSolver2D::relaxParticlePositions(Float dt) { + const Float restDist = _grid.cellSize / std::sqrt(2.0f) * 1.1f; + const Float restDistSqr = restDist * restDist; + const Float overlappedSqr = restDistSqr * 0.0001f; + const Float jitterMag = restDist / dt / 128.0f * 0.01f; + static constexpr Float stiffness = 5.0f; + + _particles.loopAll( + [&](uint32_t p) { + const Vector2 ppos = _particles.positions[p]; + const Vector2i gridCoord = _grid.getValidCellIdx(ppos); + Vector2 spring = Vector2(0); + _grid.loopNeigborParticles( + gridCoord.x(), gridCoord.y(), -1, 1, -1, 1, [&](uint32_t q) { + if(p == q) { + return; + } + const Vector2 xpq = ppos - _particles.positions[q]; + const auto distSqr = xpq.dot(); + const auto w = stiffness * smoothKernel(distSqr, restDistSqr); + if(distSqr > overlappedSqr) { + spring += xpq * (w / std::sqrt(distSqr) * restDist); + } else { + spring.x() += ((rand() & 255) - 128) * jitterMag; + spring.y() += ((rand() & 255) - 128) * jitterMag; + } + }); + + const Vector2 newPos = ppos + dt * spring; + _particles.tmp[p] = _grid.constrainBoundary(newPos); + }); + + _particles.positions.swap(_particles.tmp); +} + +void ApicSolver2D::gridVelocity2Particle() { + auto& u = _grid.u; + auto& v = _grid.v; + const auto dxInv = _grid.invCellSize; + + _particles.loopAll( + [&](uint32_t p) { + const Vector2 gridPos = _grid.getGridPos(_particles.positions[p]); + const Vector2 px = gridPos - Vector2(0, 0.5); + const Vector2 py = gridPos - Vector2(0.5, 0); + + _particles.velocities[p] = Vector2(u.interpolateValue(px), + v.interpolateValue(py)); + _particles.affineMat[p] = Matrix2x2(u.affineInterpolateValue(px) * dxInv, + v.affineInterpolateValue(py) * dxInv); + }); +} +} } diff --git a/src/fluidsimulation2d/FluidSolver/ApicSolver2D.h b/src/fluidsimulation2d/FluidSolver/ApicSolver2D.h new file mode 100644 index 000000000..87b174110 --- /dev/null +++ b/src/fluidsimulation2d/FluidSolver/ApicSolver2D.h @@ -0,0 +1,81 @@ +#ifndef Magnum_Examples_FluidSimulation2D_ApicSolver2D_h +#define Magnum_Examples_FluidSimulation2D_ApicSolver2D_h +/* + This file is part of Magnum. + + Original authors — credit is appreciated but not required: + + 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 — + Vladimír Vondruš + 2019 — Nghia Truong + + This is free and unencumbered software released into the public domain. + + Anyone is free to copy, modify, publish, use, compile, sell, or distribute + this software, either in source code form or as a compiled binary, for any + purpose, commercial or non-commercial, and by any means. + + In jurisdictions that recognize copyright laws, the author or authors of + this software dedicate any and all copyright interest in the software to + the public domain. We make this dedication for the benefit of the public + at large and to the detriment of our heirs and successors. We intend this + dedication to be an overt act of relinquishment in perpetuity of all + present and future rights to this software under copyright law. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER + IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN + CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include + +#include "FluidSolver/SolverData.h" + +namespace Magnum { namespace Examples { +/* + * 2D Affine Particle-in-Cell fluid solver + */ +class ApicSolver2D { +public: + ApicSolver2D(const Vector2& origin, Float cellSize, Int Ni, Int Nj, SceneObjects* sceneObjs); + + /* Manipulation */ + void reset() { _particles.reset(); _particles.addParticles(_particles.positionsT0, 0); } + void emitParticles() { generateParticles(_objects->emitter, 10); } + void addRepulsiveVelocity(const Vector2& p0, const Vector2& p1, Float dt, Float radius, Float magnitude); + void advanceFrame(Float frameDuration); + + /* Properties */ + UnsignedInt numParticles() const { return _particles.size(); } + Float particleRadius() const { return _particles.particleRadius; } + const std::vector& particlePositions() const { return _particles.positions; } + +private: + /* Initialization */ + void initBoundary(); + void generateParticles(const SDFObject& sdfObj, Float initialVelocity_y); + + /* Simulation */ + Float timestepCFL() const; + void moveParticles(Float dt); + void collectParticlesToCells(); + void particleVelocity2Grid(); + void extrapolate(Array2X& grid, Array2X& tmp_grid, Array2X& valid, Array2X& old_valid) const; + void addGravity(Float dt); + void computeFluidSDF(); + void solvePressures(Float dt); + void constrainVelocity(); + void relaxParticlePositions(Float dt); + void gridVelocity2Particle(); + + Containers::Pointer _objects; + ParticleData _particles; + GridData _grid; + LinearSystemSolver _pressureSolver; +}; +} } + +#endif diff --git a/src/fluidsimulation2d/FluidSolver/SolverData.h b/src/fluidsimulation2d/FluidSolver/SolverData.h new file mode 100644 index 000000000..32dcee1ba --- /dev/null +++ b/src/fluidsimulation2d/FluidSolver/SolverData.h @@ -0,0 +1,192 @@ +#ifndef Magnum_Examples_FluidSimulation2D_SolverData_h +#define Magnum_Examples_FluidSimulation2D_SolverData_h +/* + This file is part of Magnum. + + Original authors — credit is appreciated but not required: + + 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 — + Vladimír Vondruš + 2019 — Nghia Truong + + This is free and unencumbered software released into the public domain. + + Anyone is free to copy, modify, publish, use, compile, sell, or distribute + this software, either in source code form or as a compiled binary, for any + purpose, commercial or non-commercial, and by any means. + + In jurisdictions that recognize copyright laws, the author or authors of + this software dedicate any and all copyright interest in the software to + the public domain. We make this dedication for the benefit of the public + at large and to the detriment of our heirs and successors. We intend this + dedication to be an overt act of relinquishment in perpetuity of all + present and future rights to this software under copyright law. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER + IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN + CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include +#include + +#include "DataStructures/Array2X.h" +#include "DataStructures/SDFObject.h" +#include "DataStructures/PCGSolver.h" + +namespace Magnum { namespace Examples { +struct SceneObjects { + SDFObject emitterT0; /* emitter that is called once upon initialization */ + SDFObject emitter; /* emitter that is called by user when requested */ + SDFObject boundary; /* solid boundary */ +}; + +struct ParticleData { + explicit ParticleData(Float cellSize) : particleRadius{cellSize* 0.5f} {} + + uint32_t size() const { return static_cast(positions.size()); } + + void addParticles(const std::vector& newParticles, Float initialVelocity_y) { + if(positionsT0.size() == 0) { + positionsT0 = newParticles; + } + positions.insert(positions.end(), newParticles.begin(), newParticles.end()); + velocities.resize(size(), Vector2(0, -initialVelocity_y)); + affineMat.resize(size(), Matrix2x2(0)); + tmp.resize(size(), Vector2(0)); + } + + void reset() { + positions.resize(0); + velocities.resize(0); + affineMat.resize(0); + tmp.resize(0); + } + + template + void loopAll(Function&& func) const { + for(uint32_t p = 0, pend = size(); p < pend; ++p) { + func(p); + } + } + + const Float particleRadius; + std::vector positionsT0; + std::vector positions; + std::vector velocities; + std::vector affineMat; + std::vector tmp; +}; + +struct GridData { + GridData(const Vector2& origin_, Float cellSize_, Int Ni_, Int Nj_) : + origin{origin_}, Ni{Ni_}, Nj{Nj_}, + cellSize{cellSize_}, + invCellSize{Float(1.0) / cellSize} { + u.resize(Ni + 1, Nj); + v.resize(Ni, Nj + 1); + u_tmp.resize(Ni + 1, Nj); + v_tmp.resize(Ni, Nj + 1); + u_weights.resize(Ni + 1, Nj); + v_weights.resize(Ni, Nj + 1); + u_valid.resize(Ni + 1, Nj); + v_valid.resize(Ni, Nj + 1); + u_old_valid.resize(Ni + 1, Nj); + v_old_valid.resize(Ni, Nj + 1); + + fluidSDF.resize(Ni, Nj); + boundarySDF.resize(Ni + 1, Nj + 1); + cellParticles.resize(Ni, Nj); + } + + Vector2 getGridPos(const Vector2& worldPos) const { return (worldPos - origin) * invCellSize; } + Vector2 getWorldPos(Float grid_x, Float grid_y) const { return Vector2(grid_x, grid_y) * cellSize + origin; } + + bool isValidCellIdx(int x, int y) const { return x >= 0 && x < Ni && y >= 0 && y < Nj; } + Vector2i getCellIdx(const Vector2& worldPos) const { return Vector2i(getGridPos(worldPos)); } + Vector2i getValidCellIdx(const Vector2& worldPos) const { + auto tmp = getCellIdx(worldPos); + tmp.x() = Math::max(0, Math::min(Ni - 1, tmp.x())); + tmp.y() = Math::max(0, Math::min(Nj - 1, tmp.y())); + return tmp; + } + + Vector2 velocityFromGridPos(const Vector2& gridPos) const { + const Vector2 px = Vector2(gridPos[0], gridPos[1] - 0.5f); + const Vector2 py = Vector2(gridPos[0] - 0.5f, gridPos[1]); + return Vector2(u.interpolateValue(px), + v.interpolateValue(py)); + } + + Vector2 constrainBoundary(const Vector2& worldPos) const { + const Vector2 gridPos = getGridPos(worldPos); + const auto sdfVal = boundarySDF.interpolateValue(gridPos); + if(sdfVal < 0) { + const auto normal = boundarySDF.interpolateGradient(gridPos); + return worldPos - sdfVal * normal; + } else { + return worldPos; + } + } + + template + void loopNeigborParticles(int i, int j, int il, int ih, int jl, int jh, Function&& func) const { + for(int sj = j + jl; sj <= j + jh; ++sj) { + for(int si = i + il; si <= i + ih; ++si) { + if(si < 0 || si > Ni - 1 || sj < 0 || sj > Nj - 1) { continue; } + const auto& neighbors = cellParticles(si, sj); + for(auto p : neighbors) { + func(p); + } + } + } + } + + /* Grid spatial information */ + const Vector2 origin; + const Int Ni, Nj; + const Float cellSize; + const Float invCellSize; + + /* Nodes and cells' data */ + Array2X u, u_tmp, u_weights; + Array2X v, v_tmp, v_weights; + Array2X u_valid, v_valid, u_old_valid, v_old_valid; + Array2X boundarySDF; + Array2X fluidSDF; + + Array2X> cellParticles; +}; + +struct LinearSystemSolver { + void resize(std::size_t newSize) { + rhs.resize(newSize); + solution.resize(newSize); + matrix.resize(newSize); + } + + void clear() { + matrix.clear(); + solution.assign(solution.size(), 0); + } + + void solve() { + bool bSuccess = pcgSolver.solve(matrix, rhs, solution); + if(!bSuccess) { + Error{} << "Pressure solve failed!"; + } + } + + /* Use double for linear system (the solver converges slower if using float number) */ + using pcg_real = Double; + PCGSolver pcgSolver; + SparseMatrix matrix; + std::vector rhs; + std::vector solution; +}; +} } + +#endif diff --git a/src/fluidsimulation2d/Shaders/ParticleSphereShader2D.cpp b/src/fluidsimulation2d/Shaders/ParticleSphereShader2D.cpp new file mode 100644 index 000000000..250cc0cd1 --- /dev/null +++ b/src/fluidsimulation2d/Shaders/ParticleSphereShader2D.cpp @@ -0,0 +1,98 @@ +/* + This file is part of Magnum. + + Original authors — credit is appreciated but not required: + + 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 — + Vladimír Vondruš + 2019 — Nghia Truong + + This is free and unencumbered software released into the public domain. + + Anyone is free to copy, modify, publish, use, compile, sell, or distribute + this software, either in source code form or as a compiled binary, for any + purpose, commercial or non-commercial, and by any means. + + In jurisdictions that recognize copyright laws, the author or authors of + this software dedicate any and all copyright interest in the software to + the public domain. We make this dedication for the benefit of the public + at large and to the detriment of our heirs and successors. We intend this + dedication to be an overt act of relinquishment in perpetuity of all + present and future rights to this software under copyright law. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER + IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN + CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include "Shaders/ParticleSphereShader2D.h" + +#include +#include +#include +#include +#include +#include + +namespace Magnum { namespace Examples { +ParticleSphereShader2D::ParticleSphereShader2D() { + Utility::Resource rs("data"); + + GL::Shader vertShader{ GL::Version::GL330, GL::Shader::Type::Vertex }; + GL::Shader fragShader{ GL::Version::GL330, GL::Shader::Type::Fragment }; + vertShader.addSource(rs.get("ParticleSphereShader2D.vert")); + fragShader.addSource(rs.get("ParticleSphereShader2D.frag")); + + CORRADE_INTERNAL_ASSERT(GL::Shader::compile({ vertShader, fragShader })); + attachShaders({ vertShader, fragShader }); + CORRADE_INTERNAL_ASSERT(link()); + + _uNumParticles = uniformLocation("numParticles"); + _uParticleRadius = uniformLocation("particleRadius"); + + _uColorMode = uniformLocation("colorMode"); + _uColor = uniformLocation("uniformColor"); + + _uViewProjectionMatrix = uniformLocation("viewProjectionMatrix"); + _uScreenHeight = uniformLocation("screenHeight"); + _uDomainHeight = uniformLocation("domainHeight"); +} + +ParticleSphereShader2D& ParticleSphereShader2D::setNumParticles(Int numParticles) { + setUniform(_uNumParticles, numParticles); + return *this; +} + +ParticleSphereShader2D& ParticleSphereShader2D::setParticleRadius(Float radius) { + setUniform(_uParticleRadius, radius); + return *this; +} + +ParticleSphereShader2D& ParticleSphereShader2D::setColorMode(Int colorMode) { + setUniform(_uColorMode, colorMode); + return *this; +} + +ParticleSphereShader2D& ParticleSphereShader2D::setColor(const Color3& color) { + setUniform(_uColor, color); + return *this; +} + +ParticleSphereShader2D& ParticleSphereShader2D::setViewProjectionMatrix(const Matrix3& matrix) { + setUniform(_uViewProjectionMatrix, matrix); + return *this; +} + +ParticleSphereShader2D& ParticleSphereShader2D::setScreenHeight(Int height) { + setUniform(_uScreenHeight, height); + return *this; +} + +ParticleSphereShader2D& ParticleSphereShader2D::setDomainHeight(Int height) { + setUniform(_uDomainHeight, height); + return *this; +} +} } diff --git a/src/fluidsimulation2d/Shaders/ParticleSphereShader2D.frag b/src/fluidsimulation2d/Shaders/ParticleSphereShader2D.frag new file mode 100644 index 000000000..647da07c8 --- /dev/null +++ b/src/fluidsimulation2d/Shaders/ParticleSphereShader2D.frag @@ -0,0 +1,39 @@ +/* + This file is part of Magnum. + + Original authors — credit is appreciated but not required: + + 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 — + Vladimír Vondruš + 2019 — Nghia Truong + + This is free and unencumbered software released into the public domain. + + Anyone is free to copy, modify, publish, use, compile, sell, or distribute + this software, either in source code form or as a compiled binary, for any + purpose, commercial or non-commercial, and by any means. + + In jurisdictions that recognize copyright laws, the author or authors of + this software dedicate any and all copyright interest in the software to + the public domain. We make this dedication for the benefit of the public + at large and to the detriment of our heirs and successors. We intend this + dedication to be an overt act of relinquishment in perpetuity of all + present and future rights to this software under copyright law. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER + IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN + CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +flat in vec3 color; +layout(location = 0) out lowp vec4 fragmentColor; + +void main() { + vec2 point = gl_PointCoord.xy*vec2(2.0, -2.0) + vec2(-1.0, 1.0); + float mag = dot(point, point); + if(mag > 1.0) discard; /* outside the circle */ + fragmentColor = vec4(color, 1.0); +} diff --git a/src/fluidsimulation2d/Shaders/ParticleSphereShader2D.h b/src/fluidsimulation2d/Shaders/ParticleSphereShader2D.h new file mode 100644 index 000000000..ee251681c --- /dev/null +++ b/src/fluidsimulation2d/Shaders/ParticleSphereShader2D.h @@ -0,0 +1,67 @@ +#ifndef Magnum_Examples_FluidSimulation2D_ParticleSphereFlatShader_h +#define Magnum_Examples_FluidSimulation2D_ParticleSphereFlatShader_h +/* + This file is part of Magnum. + + Original authors — credit is appreciated but not required: + + 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 — + Vladimír Vondruš + 2019 — Nghia Truong + + This is free and unencumbered software released into the public domain. + + Anyone is free to copy, modify, publish, use, compile, sell, or distribute + this software, either in source code form or as a compiled binary, for any + purpose, commercial or non-commercial, and by any means. + + In jurisdictions that recognize copyright laws, the author or authors of + this software dedicate any and all copyright interest in the software to + the public domain. We make this dedication for the benefit of the public + at large and to the detriment of our heirs and successors. We intend this + dedication to be an overt act of relinquishment in perpetuity of all + present and future rights to this software under copyright law. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER + IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN + CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +#include + +namespace Magnum { namespace Examples { +class ParticleSphereShader2D : public GL::AbstractShaderProgram { +public: + enum ColorMode { + UniformDiffuseColor = 0, + RampColorById + }; + + explicit ParticleSphereShader2D(); + + ParticleSphereShader2D& setNumParticles(Int numParticles); + ParticleSphereShader2D& setParticleRadius(Float radius); + + ParticleSphereShader2D& setColorMode(Int colorMode); + ParticleSphereShader2D& setColor(const Color3& color); + + ParticleSphereShader2D& setViewport(const Vector2i& viewport); + ParticleSphereShader2D& setViewProjectionMatrix(const Matrix3& matrix); + ParticleSphereShader2D& setScreenHeight(Int height); + ParticleSphereShader2D& setDomainHeight(Int height); + +private: + Int _uNumParticles, + _uParticleRadius, + _uColorMode, + _uColor, + _uViewProjectionMatrix, + _uScreenHeight, + _uDomainHeight; +}; +} } + +#endif diff --git a/src/fluidsimulation2d/Shaders/ParticleSphereShader2D.vert b/src/fluidsimulation2d/Shaders/ParticleSphereShader2D.vert new file mode 100644 index 000000000..d6bc1cc29 --- /dev/null +++ b/src/fluidsimulation2d/Shaders/ParticleSphereShader2D.vert @@ -0,0 +1,70 @@ +/* + This file is part of Magnum. + + Original authors — credit is appreciated but not required: + + 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 — + Vladimír Vondruš + 2019 — Nghia Truong + + This is free and unencumbered software released into the public domain. + + Anyone is free to copy, modify, publish, use, compile, sell, or distribute + this software, either in source code form or as a compiled binary, for any + purpose, commercial or non-commercial, and by any means. + + In jurisdictions that recognize copyright laws, the author or authors of + this software dedicate any and all copyright interest in the software to + the public domain. We make this dedication for the benefit of the public + at large and to the detriment of our heirs and successors. We intend this + dedication to be an overt act of relinquishment in perpetuity of all + present and future rights to this software under copyright law. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL + THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER + IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN + CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. + */ + +uniform highp mat3 viewProjectionMatrix; +uniform int numParticles; +uniform int colorMode; +uniform int screenHeight; +uniform int domainHeight; +uniform float particleRadius; +uniform vec3 uniformColor; + + +layout(location = 0) in highp vec2 position; +flat out vec3 color; + +const vec3 colorRamp[] = vec3[] ( + vec3(1.0, 0.0, 0.0), + vec3(1.0, 0.5, 0.0), + vec3(1.0, 1.0, 0.0), + vec3(1.0, 0.0, 1.0), + vec3(0.0, 1.0, 0.0), + vec3(0.0, 1.0, 1.0), + vec3(0.0, 0.0, 1.0) +); + +vec3 generateVertexColor() { + if(colorMode == 1 ) { /* ramp color by particle id */ + float segmentSize = float(numParticles)/6.0f; + float segment = floor(float(gl_VertexID)/segmentSize); + float t = (float(gl_VertexID) - segmentSize*segment)/segmentSize; + vec3 startVal = colorRamp[int(segment)]; + vec3 endVal = colorRamp[int(segment) + 1]; + return mix(startVal, endVal, t); + } else { /* uniform diffuse color */ + return uniformColor; + } +} + +void main() { + color = generateVertexColor(); + gl_PointSize = particleRadius * float(screenHeight) / float(domainHeight); + gl_Position = mat4(viewProjectionMatrix) * vec4(position, 0, 1.0); +} diff --git a/src/fluidsimulation2d/SourceSansPro-Regular.ttf b/src/fluidsimulation2d/SourceSansPro-Regular.ttf new file mode 100644 index 000000000..b422bf432 Binary files /dev/null and b/src/fluidsimulation2d/SourceSansPro-Regular.ttf differ diff --git a/src/fluidsimulation2d/resources.conf b/src/fluidsimulation2d/resources.conf new file mode 100644 index 000000000..37053ec31 --- /dev/null +++ b/src/fluidsimulation2d/resources.conf @@ -0,0 +1,12 @@ +group=data + +[file] +filename=Shaders/ParticleSphereShader2D.vert +alias=ParticleSphereShader2D.vert + +[file] +filename=Shaders/ParticleSphereShader2D.frag +alias=ParticleSphereShader2D.frag + +[file] +filename=SourceSansPro-Regular.ttf