From 03e9156cb23068144dd26c59891370e45f7429c3 Mon Sep 17 00:00:00 2001
From: adany292 <adany292@student.liu.se>
Date: Sun, 8 Dec 2024 18:47:51 +0100
Subject: [PATCH] Rewrite/reeefactored code to an application for eeasieer
 maintnance, next refactor cuda code as well

---
 CMakeLists.txt                   |  43 +--
 src/Camera.cpp                   | 145 ++++----
 src/FluidRenderer.cpp            | 189 ++++++++++
 src/FluidSimulation.cpp          |  66 ++++
 src/FluidSimulationApp.cpp       |  75 ++++
 src/GUI_fluid_sim.cpp            |   4 +-
 src/MainGUI.cpp                  | 176 +++++++++
 src/include/Camera.h             |  56 +--
 src/include/FluidRenderer.h      |  54 +++
 src/include/FluidSimulation.h    |  32 ++
 src/include/FluidSimulationApp.h |  36 ++
 src/include/PCISPH_sim.h         |   2 +-
 src/kernels_cuda.cu              | 606 -------------------------------
 src/settings/constants.h         |   2 +
 14 files changed, 749 insertions(+), 737 deletions(-)
 create mode 100644 src/FluidRenderer.cpp
 create mode 100644 src/FluidSimulation.cpp
 create mode 100644 src/FluidSimulationApp.cpp
 create mode 100644 src/MainGUI.cpp
 create mode 100644 src/include/FluidRenderer.h
 create mode 100644 src/include/FluidSimulation.h
 create mode 100644 src/include/FluidSimulationApp.h
 delete mode 100644 src/kernels_cuda.cu

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 1063d17..f6daabb 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -9,20 +9,8 @@ enable_language(CUDA)
 # ------------------------------
 find_package(OpenGL REQUIRED)
 find_package(GLEW REQUIRED)
-find_package(glm REQUIRED)
 find_package(glfw3 3.3 REQUIRED)
 
-# If the above 'find_package(glfw3 3.3 REQUIRED)' does not work on your system:
-# - You may need pkg-config or a custom Findglfw3.cmake.
-# Example:
-# find_package(PkgConfig REQUIRED)
-# pkg_check_modules(GLFW3 REQUIRED glfw3)
-# and then link using ${GLFW3_LIBRARIES} and include using ${GLFW3_INCLUDE_DIRS}.
-
-# CUDA is handled by enable_language(CUDA), assuming it finds a suitable CUDA compiler.
-# If you need finer control:
-# find_package(CUDA REQUIRED)
-
 # ------------------------------
 # 2. Include ImGui
 # ------------------------------
@@ -61,6 +49,9 @@ set(COMMON_SOURCES
     src/spawn_particles.cpp
     src/Camera.cpp
     src/ImGuiManager.cpp
+    src/FluidSimulationApp.cpp 
+    src/FluidRenderer.cpp
+    src/FluidSimulation.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/common/GL_utilities.c
     # Add other common source files if needed
 )
@@ -68,30 +59,24 @@ set(COMMON_SOURCES
 add_library(common_lib STATIC ${COMMON_SOURCES})
 
 target_include_directories(common_lib PUBLIC
-    ${CMAKE_CURRENT_SOURCE_DIR}/src           # For src/*.h if needed
-    ${CMAKE_CURRENT_SOURCE_DIR}/src/settings           # For src/*.h if needed
-    /usr/local/cuda/include                   # CUDA include directory
+    ${CMAKE_CURRENT_SOURCE_DIR}/src
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/settings
+    /usr/local/cuda/include
     ${OPENGL_INCLUDE_DIRS}
     ${GLEW_INCLUDE_DIRS}
-    ${CMAKE_CURRENT_SOURCE_DIR}/common        # Common directory (GL_utilities, etc.)
-    ${CMAKE_CURRENT_SOURCE_DIR}/common/Linux/ # Linux-specific if needed
+    ${CMAKE_CURRENT_SOURCE_DIR}/common
+    ${CMAKE_CURRENT_SOURCE_DIR}/common/Linux/
     ${IMGUI_DIR}
     ${IMGUI_DIR}/backends
-    # Add other include directories as needed, e.g. ${CMAKE_CURRENT_SOURCE_DIR}/include
 )
 
-# Link necessary libraries to the common library
-# CUDA_LIBRARIES might be implicitly used if using newer CUDA integrations.
-# If not defined, you can manually specify -lcuda or rely on language CUDA integration.
 target_link_libraries(common_lib PUBLIC
     OpenGL::GL
     ${GLEW_LIBRARIES}
     glfw
     GLU
-    glm::glm
 )
 
-# Optionally, if you need CUDA separable compilation:
 set_target_properties(common_lib PROPERTIES CUDA_SEPARABLE_COMPILATION ON)
 
 # ------------------------------
@@ -106,6 +91,10 @@ target_link_libraries(cuda_sph_main PRIVATE common_lib imgui)
 add_executable(screen_spaced src/screen_spaced_main.cpp)
 target_link_libraries(screen_spaced PRIVATE common_lib imgui)
 
+# New executable for MainGUI
+add_executable(fluidSimulator src/MainGUI.cpp)
+target_link_libraries(fluidSimulator PRIVATE common_lib imgui)
+
 # ------------------------------
 # 5. Set Target Properties
 # ------------------------------
@@ -119,7 +108,7 @@ set_target_properties(screen_spaced PROPERTIES
     CXX_STANDARD_REQUIRED YES
 )
 
-# If needed, set CUDA architectures or flags:
-# For example:
-# set_property(TARGET common_lib PROPERTY CUDA_ARCHITECTURES 52)
-# This depends on your GPU and CUDA version.
+set_target_properties(fluidSimulator PROPERTIES
+    CXX_STANDARD 17
+    CXX_STANDARD_REQUIRED YES
+)
diff --git a/src/Camera.cpp b/src/Camera.cpp
index eab52c4..f44b438 100644
--- a/src/Camera.cpp
+++ b/src/Camera.cpp
@@ -1,88 +1,85 @@
-// Camera.cpp
 #include "include/Camera.h"
 #include <cstdio>
 #include <cstdlib>
 
-static Camera* g_camera = nullptr;
+namespace FluidSimulation {
 
-Camera::Camera(vec3 position, vec3 target, vec3 up, float fov_, float aspect_, float zNear_, float zFar_)
-    : campos(position), forward(VectorSub(target, position)), upVector(up),
-      fov(fov_), aspectRatio(aspect_), zNear(zNear_), zFar(zFar_), rotationMatrix(IdentityMatrix()) 
-{
-    float norm = sqrtf(forward.x*forward.x + forward.y*forward.y + forward.z*forward.z);
-    forward = ScalarMult(forward, 1.0f / norm);
-}
+    static Camera* g_camera = nullptr;
 
-// Processes keyboard input to control the camera
-void Camera::processKey(unsigned char key)
-{
-    switch (key)
+    Camera::Camera(vec3 position, vec3 target, vec3 up, float fov_, float aspect_, float zNear_, float zFar_)
+        : campos(position), forward(VectorSub(target, position)), upVector(up),
+          fov(fov_), aspectRatio(aspect_), zNear(zNear_), zFar(zFar_), rotationMatrix(IdentityMatrix()) 
     {
-        case 'w':
-            campos += ScalarMult(forward, 0.5f);
-            break;
-        case 's':
-            campos -= ScalarMult(forward, 0.5f);
-            break;
-        case 'a':
-            campos -= ScalarMult(CrossProduct(forward, upVector), 0.5f);
-            break;
-        case 'd':
-            campos += ScalarMult(CrossProduct(forward, upVector), 0.5f);
-            break;
-        case 'j':
-            rotateSide(0.03f);
-            break;
-        case 'l':
-            rotateSide(-0.03f);
-            break;
-        case 'i':
-            // Rotate upwards
-            rotateUp(0.05f);
-            break;
-        case 'k':
-            // Rotate downwards
-            rotateUp(-0.05f);
-            break;
-        case 32: // SPACEBAR
-            campos += ScalarMult(upVector, 1.0f);
-            break;
-        case 'c':
-            campos -= ScalarMult(upVector, 1.0f);
-            break;
-        case 'v':
-            // Toggle boundary display
-            toggleDisplayBorder();
-            break;
-        case 0x1b: // ESC
-            exit(0);
+        float norm = sqrtf(forward.x*forward.x + forward.y*forward.y + forward.z*forward.z);
+        forward = ScalarMult(forward, 1.0f / norm);
     }
-}
 
-void Camera::reshape(int width, int height) {
-    glViewport(0, 0, width, height);
-    aspectRatio = (float)width / (float)height;
-}
+    void Camera::processKey(unsigned char key) {
+        switch (key) {
+            case 'w':
+                campos += ScalarMult(forward, 0.5f);
+                break;
+            case 's':
+                campos -= ScalarMult(forward, 0.5f);
+                break;
+            case 'a':
+                campos -= ScalarMult(CrossProduct(forward, upVector), 0.5f);
+                break;
+            case 'd':
+                campos += ScalarMult(CrossProduct(forward, upVector), 0.5f);
+                break;
+            case 'j':
+                rotateSide(0.03f);
+                break;
+            case 'l':
+                rotateSide(-0.03f);
+                break;
+            case 'i':
+                rotateUp(0.05f);
+                break;
+            case 'k':
+                rotateUp(-0.05f);
+                break;
+            case 32: // SPACEBAR
+                campos += ScalarMult(upVector, 1.0f);
+                break;
+            case 'c':
+                campos -= ScalarMult(upVector, 1.0f);
+                break;
+            case 'v':
+                toggleDisplayBorder();
+                break;
+            case 0x1b: // ESC
+                exit(0);
+        }
+    }
 
-mat4 Camera::getViewMatrix() const {
-    return lookAt(campos, VectorAdd(campos, forward), upVector);
-}
+    void Camera::reshape(int width, int height) {
+        glViewport(0, 0, width, height);
+        aspectRatio = (float)width / (float)height;
+    }
 
-mat4 Camera::getProjectionMatrix() const {
-    return perspective(fov, aspectRatio, zNear, zFar);
-}
+    mat4 Camera::getViewMatrix() const {
+        return lookAt(campos, VectorAdd(campos, forward), upVector);
+    }
 
-void Camera::rotateSide(float angle) {
-    mat3 rot = mat4tomat3(Ry(angle));
-    forward = MultMat3Vec3(rot, forward);
-}
+    mat4 Camera::getProjectionMatrix() const {
+        return perspective(fov, aspectRatio, zNear, zFar);
+    }
 
-void Camera::rotateUp(float angle) {
-    vec3 side = CrossProduct(forward, upVector);
-    mat3 rot = mat4tomat3(ArbRotate(side, angle));
-    forward = MultMat3Vec3(rot, forward);
-}
+    void Camera::rotateSide(float angle) {
+        mat3 rot = mat4tomat3(Ry(angle));
+        forward = MultMat3Vec3(rot, forward);
+    }
+
+    void Camera::rotateUp(float angle) {
+        vec3 side = CrossProduct(forward, upVector);
+        mat3 rot = mat4tomat3(ArbRotate(side, angle));
+        forward = MultMat3Vec3(rot, forward);
+    }
+
+    void initCamera(Camera& camera) {
+        g_camera = &camera;
+    }
 
-void initCamera(Camera& camera) {
-    g_camera = &camera;
-}
+} // namespace FluidSimulation
\ No newline at end of file
diff --git a/src/FluidRenderer.cpp b/src/FluidRenderer.cpp
new file mode 100644
index 0000000..55cd677
--- /dev/null
+++ b/src/FluidRenderer.cpp
@@ -0,0 +1,189 @@
+#include "include/FluidRenderer.h"
+#include "GL_utilities.h" // For loadShaders
+#include <iostream>
+
+#define MAIN
+#include "VectorUtils4.h"
+#undef MAIN
+
+namespace FluidSimulation {
+
+    FluidRenderer::FluidRenderer(int width, int height, Camera* camera)
+        : m_fbWidth(width), m_fbHeight(height), m_camera(camera) {
+
+            init();
+    }
+
+    FluidRenderer::~FluidRenderer() {
+        glDeleteTextures(1, &m_fboTexture);
+        glDeleteRenderbuffers(1, &m_depthRBO);
+        glDeleteFramebuffers(1, &m_fbo);
+        std::cout << "FluidRenderer destroyed." << std::endl;
+    }
+
+    void FluidRenderer::init() {
+        initGLEW();
+        initOpenGLSettings();
+        initShaders();
+        initShaderUniforms();
+        initFBOs();
+    }
+
+    void FluidRenderer::initGLEW() {
+        glewExperimental = GL_TRUE;
+        GLenum err = glewInit();
+        if (err != GLEW_OK) {
+            fprintf(stderr, "Error initializing GLEW: %s\n", glewGetErrorString(err));
+            exit(1);
+        }
+        std::cout << "GLEW initialized successfully." << std::endl;
+    }
+
+    void FluidRenderer::initOpenGLSettings() {
+        glEnable(GL_PROGRAM_POINT_SIZE);
+        glEnable(GL_DEPTH_TEST);
+        glDisable(GL_CULL_FACE);
+
+        // Create and bind a VAO - mandatory for opeeeenGl 3.3+ ??
+        GLuint vao;
+        glGenVertexArrays(1, &vao);
+        glBindVertexArray(vao);
+
+        std::cout << "OpenGL settings initialized." << std::endl;
+    }
+
+    void FluidRenderer::initShaders() {
+        m_shaderProgram = loadShaders("../shaders/particle.vert", "../shaders/particle.frag");
+        glUseProgram(m_shaderProgram);
+        std::cout << "Shaders initialized successfully." << std::endl;
+    }
+
+    void FluidRenderer::initShaderUniforms() {
+        m_positionAttribLoc = glGetAttribLocation(m_shaderProgram, "aPosition");
+        if (m_positionAttribLoc == -1) {
+            fprintf(stderr, "Could not find attribute 'aPosition'\n");
+        } else {
+            printf("Attribute 'aPosition' location: %d\n", m_positionAttribLoc);
+        }
+
+        m_modelMatrixLoc = glGetUniformLocation(m_shaderProgram, "modelMatrix");
+        m_viewMatrixLoc = glGetUniformLocation(m_shaderProgram, "viewMatrix");
+        m_projectionMatrixLoc = glGetUniformLocation(m_shaderProgram, "projectionMatrix");
+        m_pointRadiusLoc = glGetUniformLocation(m_shaderProgram, "pointRadius");
+        m_pointScaleLoc = glGetUniformLocation(m_shaderProgram, "pointScale");
+        m_uColorLoc = glGetUniformLocation(m_shaderProgram, "uColor");
+
+        // Should be moved to a setting and be controllable with GUI
+        if (m_pointRadiusLoc != -1) glUniform1f(m_pointRadiusLoc, 20.0f); // Example radius
+        if (m_pointScaleLoc != -1) glUniform1f(m_pointScaleLoc, 10.0f);   // Example scale
+
+        mat4 modelMatrix = IdentityMatrix();
+        glUniformMatrix4fv(m_modelMatrixLoc, 1, GL_TRUE, modelMatrix.m);
+
+        std::cout << "Shader uniforms initialized." << std::endl;
+    }
+
+    void FluidRenderer::initFBOs() {
+        glGenFramebuffers(1, &m_fbo);
+        glBindFramebuffer(GL_FRAMEBUFFER, m_fbo);
+
+        // Create texture for FBO
+        glGenTextures(1, &m_fboTexture);
+        glBindTexture(GL_TEXTURE_2D, m_fboTexture);
+        glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA8, m_fbWidth, m_fbHeight, 0, GL_RGBA, GL_UNSIGNED_BYTE, NULL);
+        glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
+        glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, GL_TEXTURE_2D, m_fboTexture, 0);
+
+        // Create depth renderbuffer
+        glGenRenderbuffers(1, &m_depthRBO);
+        glBindRenderbuffer(GL_RENDERBUFFER, m_depthRBO);
+        glRenderbufferStorage(GL_RENDERBUFFER, GL_DEPTH_COMPONENT24, m_fbWidth, m_fbHeight);
+        glFramebufferRenderbuffer(GL_FRAMEBUFFER, GL_DEPTH_ATTACHMENT, GL_RENDERBUFFER, m_depthRBO);
+
+        if (glCheckFramebufferStatus(GL_FRAMEBUFFER) != GL_FRAMEBUFFER_COMPLETE) {
+            fprintf(stderr, "FBO not complete: 0x%X\n", glCheckFramebufferStatus(GL_FRAMEBUFFER));
+        }
+
+        glBindFramebuffer(GL_FRAMEBUFFER, 0); // Unbind FBO
+        std::cout << "FBO initialized successfully." << std::endl;
+    }
+
+    void FluidRenderer::resizeFBO(int newWidth, int newHeight) {
+        if (newWidth <= 0 || newHeight <= 0) return;
+
+        m_fbWidth = newWidth;
+        m_fbHeight = newHeight;
+
+        // Bind the existing FBO
+        glBindFramebuffer(GL_FRAMEBUFFER, m_fbo);
+
+        // Resize the color attachment texture
+        glBindTexture(GL_TEXTURE_2D, m_fboTexture);
+        glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA8, m_fbWidth, m_fbHeight, 0, GL_RGBA, GL_UNSIGNED_BYTE, NULL);
+
+        // Resize the depth renderbuffer
+        glBindRenderbuffer(GL_RENDERBUFFER, m_depthRBO);
+        glRenderbufferStorage(GL_RENDERBUFFER, GL_DEPTH_COMPONENT24, m_fbWidth, m_fbHeight);
+
+        // Check framebuffer completeness
+        if (glCheckFramebufferStatus(GL_FRAMEBUFFER) != GL_FRAMEBUFFER_COMPLETE) {
+            fprintf(stderr, "FBO not complete after resize: 0x%X\n", glCheckFramebufferStatus(GL_FRAMEBUFFER));
+        }
+
+        glBindFramebuffer(GL_FRAMEBUFFER, 0); // Unbind FBO
+        //std::cout << "FBO resized to " << newWidth << "x" << newHeight << "." << std::endl;
+    }
+
+
+    void FluidRenderer::setVBOs(GLuint particleVBO, GLuint boundaryVBO) {
+        m_particleVBO = particleVBO;
+        m_boundaryVBO = boundaryVBO;
+    }
+
+    void FluidRenderer::renderParticles(size_t particleCount) {
+        glUniform3f(m_uColorLoc, 0.29f, 0.573f, 1.0f);
+        glBindBuffer(GL_ARRAY_BUFFER, m_particleVBO);
+        glEnableVertexAttribArray(m_positionAttribLoc);
+        glVertexAttribPointer(m_positionAttribLoc, 3, GL_FLOAT, GL_FALSE, sizeof(float3), (void*)0);
+        glDrawArrays(GL_POINTS, 0, particleCount); // Use actual particle count
+        glDisableVertexAttribArray(m_positionAttribLoc);
+    }
+
+    void FluidRenderer::renderBoundaries(size_t boundaryCount) {
+        glUniform3f(m_uColorLoc, 0.29f, 1.0f, 0.57f);
+        glBindBuffer(GL_ARRAY_BUFFER, m_boundaryVBO);
+        glEnableVertexAttribArray(m_positionAttribLoc);
+        glVertexAttribPointer(m_positionAttribLoc, 3, GL_FLOAT, GL_FALSE, sizeof(float3), (void*)0);
+        glDrawArrays(GL_POINTS, 0, boundaryCount); // Use actual boundary particle count
+        glDisableVertexAttribArray(m_positionAttribLoc);
+    }
+
+    void FluidRenderer::renderFrame(size_t particleCount, size_t boundaryCount) {
+        glClearColor(0.9f, 0.9f, 0.9f, 1.0f);
+
+        glBindFramebuffer(GL_FRAMEBUFFER, m_fbo);
+        glViewport(0, 0, m_fbWidth, m_fbHeight);
+        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+        glDepthMask(GL_TRUE);
+
+        glUseProgram(m_shaderProgram);
+
+        // Update shader with camera matrices
+        mat4 projectionMatrix = m_camera->getProjectionMatrix();
+        glUniformMatrix4fv(m_projectionMatrixLoc, 1, GL_TRUE, projectionMatrix.m);
+
+        mat4 viewMatrix = m_camera->getViewMatrix();
+        glUniformMatrix4fv(m_viewMatrixLoc, 1, GL_TRUE, viewMatrix.m);
+
+        // Render particles
+        renderParticles(particleCount);
+
+        // Render boundaries if enabled
+        if (m_camera->shouldDisplayBorder()) {
+            renderBoundaries(boundaryCount);
+        }
+
+        glBindFramebuffer(GL_FRAMEBUFFER, 0);
+    }
+
+}
diff --git a/src/FluidSimulation.cpp b/src/FluidSimulation.cpp
new file mode 100644
index 0000000..8f58b8f
--- /dev/null
+++ b/src/FluidSimulation.cpp
@@ -0,0 +1,66 @@
+#include "include/FluidSimulation.h"
+#include "include/spawn_particles.h" // For createBoundaryParticles, initParticles
+#include "include/PCISPH_sim.h"      // For CUDA simulation functions
+#include <iostream>
+
+namespace FluidSimulation {
+
+    FluidSimulation::FluidSimulation(int numParticles)
+        : m_numParticles(numParticles), m_numBoundaryParticles(0),
+          m_vbo_particles(0), m_vbo_boundary(0), m_timeStep(0) {
+    }
+
+    FluidSimulation::~FluidSimulation() {
+        unregisterVBOFromCUDA();
+        unregisterBoundaryVBOFromCUDA();
+        glDeleteBuffers(1, &m_vbo_particles);
+        glDeleteBuffers(1, &m_vbo_boundary);
+    }
+
+    void FluidSimulation::initSimulation() {
+        // Upload constants to the GPU
+        initializeConstants();
+
+        // Initialize particles
+        std::vector<float3> h_positions(m_numParticles);
+        std::vector<float3> h_velocities(m_numParticles);
+        float gridSpacing = 0.25f;
+
+        initParticles(h_positions, h_velocities, m_numParticles, gridSpacing);
+
+        // Initialize particle device arrays
+        initializeFluidDeviceArrays(h_positions.data(), h_velocities.data(), m_numParticles);
+
+        // Create particle VBO
+        glGenBuffers(1, &m_vbo_particles);
+        glBindBuffer(GL_ARRAY_BUFFER, m_vbo_particles);
+        glBufferData(GL_ARRAY_BUFFER, h_positions.size() * sizeof(float3), h_positions.data(), GL_DYNAMIC_DRAW);
+        glBindBuffer(GL_ARRAY_BUFFER, 0);
+        registerVBOWithCUDA(m_vbo_particles);
+
+        // Initialize boundary
+        std::vector<float3> boundary_positions;
+        std::vector<float3> boundary_velocities;
+        int boundaryLayers = 3;
+        float spacing = 0.25f;
+
+        createBoundaryParticles(boundary_positions, boundary_velocities, boundaryLayers, spacing);
+        m_numBoundaryParticles = boundary_positions.size();
+
+        initializeBoundryDeviceArrays(boundary_positions.data(), boundary_velocities.data(), m_numBoundaryParticles);
+
+        // Create boundary VBO
+        glGenBuffers(1, &m_vbo_boundary);
+        glBindBuffer(GL_ARRAY_BUFFER, m_vbo_boundary);
+        glBufferData(GL_ARRAY_BUFFER, boundary_positions.size() * sizeof(float3), boundary_positions.data(), GL_DYNAMIC_DRAW);
+        glBindBuffer(GL_ARRAY_BUFFER, 0);
+        registerBoundaryVBOWithCUDA(m_vbo_boundary);
+    }
+
+    void FluidSimulation::updateSimulation() {
+        for (int i = 0; i < 4; ++i) {
+            updateParticles(0.001);
+        }
+    }
+
+}
diff --git a/src/FluidSimulationApp.cpp b/src/FluidSimulationApp.cpp
new file mode 100644
index 0000000..6371edd
--- /dev/null
+++ b/src/FluidSimulationApp.cpp
@@ -0,0 +1,75 @@
+
+#include "include/FluidSimulationApp.h"
+#include <iostream>
+#include "imgui.h"
+#include "settings/constants.h"
+
+namespace FluidSimulation {
+
+    // Constructor
+    FluidSimulatorApp::FluidSimulatorApp(appSettings* appSettings)
+        : m_appsettings(appSettings) {
+        // Initialize FluidRenderer with the settings
+        m_camera = new Camera(
+            vec3(5.0f, 15.0f, 20.0f),
+            vec3(0.0f, 4.0f, 0.0f),
+            vec3(0.0f, 1.0f, 0.0f),
+            45.0f,
+            appSettings->width / (float)appSettings->height,
+            0.1f,
+            500.0f
+        );
+
+        m_fluidRenderer = new FluidRenderer(appSettings->width, appSettings->height, m_camera);
+        
+        m_fluidSimulation = new FluidSimulation(NUM_PARTICLES);
+        m_fluidSimulation->initSimulation();
+
+        m_fluidRenderer->setVBOs(m_fluidSimulation->getParticleVBO(), m_fluidSimulation->getBoundaryVBO());
+
+        //m_fluidRenderer->init();
+        std::cout << "FluidSimulatorApp created with settings: "
+                  << appSettings->width << "x" << appSettings->height << std::endl;
+    }
+
+    // Destructor
+    FluidSimulatorApp::~FluidSimulatorApp() {
+        delete m_fluidRenderer;
+        std::cout << "FluidSimulatorApp destroyed!" << std::endl;
+    }
+
+    // Render function
+    void FluidSimulatorApp::Render() {
+        // Render fluid simulation to FBO
+
+        m_fluidSimulation->updateSimulation();
+        m_fluidRenderer->renderFrame(
+            m_fluidSimulation->getParticleCount(),
+            m_fluidSimulation->getBoundaryParticleCount()
+        );
+
+        /*
+        // Render ImGui controls
+        ImGui::Begin("Fluid Simulation Controls");
+        ImGui::Text("Adjust simulation parameters here.");
+        // Add sliders, buttons, etc., to control the simulation
+        ImGui::End();
+        */
+        
+        // Display FBO texture in ImGui window
+        ImGui::Begin("Simulation");
+
+        ImVec2 availableSize = ImGui::GetContentRegionAvail();
+        //ImGui::Image(()(intptr_)m_fluidRenderer->getFBOTexture(),
+        //             ImVec2(m_appsettings->width, m_appsettings->height));
+        ImGui::Image((ImTextureID)(uintptr_t)m_fluidRenderer->getFBOTexture(), availableSize, ImVec2(0, 1), ImVec2(1, 0)); // Flip Y-axis
+
+        ImGui::End();
+    }
+
+    // Factory function
+    FluidSimulatorApp* creatFluidSimulatorApp(appSettings* settings) {
+        return new FluidSimulatorApp(settings);
+    }
+
+} // namespace FluidSimulation
diff --git a/src/GUI_fluid_sim.cpp b/src/GUI_fluid_sim.cpp
index a4a4d5a..5c33f1e 100644
--- a/src/GUI_fluid_sim.cpp
+++ b/src/GUI_fluid_sim.cpp
@@ -141,7 +141,7 @@ void display() {
     glfwGetFramebufferSize(glfwGetCurrentContext(), &winWidth, &winHeight);
     glViewport(0, 0, winWidth, winHeight);
 
-    printFrameStats();
+    //printFrameStats();
 }
 
 void initializeBoundaryVBO(const std::vector<float3>& h_boundary_positions) {
@@ -353,6 +353,8 @@ void GLFWFramebufferSizeCallback(GLFWwindow* window, int width, int height)
     cam->reshape(width, height);
 }
 
+
+// Just KEEP IMGUI setup in this file. Move simulation/aplication to anothere file.
 int main(int argc, char** argv) {
     // Initialize GLFW
     if (!glfwInit()) {
diff --git a/src/MainGUI.cpp b/src/MainGUI.cpp
new file mode 100644
index 0000000..58d3a0f
--- /dev/null
+++ b/src/MainGUI.cpp
@@ -0,0 +1,176 @@
+#include <GL/glew.h>
+#include <GLFW/glfw3.h>   // Use GLFW instead of GLUT
+
+#include "imgui.h"
+#include "imgui_impl_glfw.h"     // Use the GLFW backend
+#include "imgui_impl_opengl3.h"  // Use the OpenGL3 backend
+
+#include <cstdlib> // For exit()
+#include <cstdio>
+
+#include "include/FluidSimulationApp.h"
+
+void GLFWKeyCallback(GLFWwindow* window, int key, int scancode, int action, int mods) {
+    // Retrieve FluidSimulatorApp from window user pointer
+    FluidSimulation::FluidSimulatorApp* app = static_cast<FluidSimulation::FluidSimulatorApp*>(glfwGetWindowUserPointer(window));
+    if (!app) return;
+
+    // Retrieve the camera from the app
+    FluidSimulation::Camera* cam = app->getCamera();
+    if (!cam) return;
+
+    // Map GLFW keys to ASCII keys and process them
+    if (action == GLFW_PRESS || action == GLFW_REPEAT) {
+        unsigned char ascii_key;
+        switch (key) {
+            case GLFW_KEY_W: ascii_key = 'w'; break;
+            case GLFW_KEY_S: ascii_key = 's'; break;
+            case GLFW_KEY_A: ascii_key = 'a'; break;
+            case GLFW_KEY_D: ascii_key = 'd'; break;
+            case GLFW_KEY_J: ascii_key = 'j'; break;
+            case GLFW_KEY_L: ascii_key = 'l'; break;
+            case GLFW_KEY_I: ascii_key = 'i'; break;
+            case GLFW_KEY_K: ascii_key = 'k'; break;
+            case GLFW_KEY_SPACE: ascii_key = ' '; break;
+            case GLFW_KEY_C: ascii_key = 'c'; break;
+            case GLFW_KEY_V: ascii_key = 'v'; break;
+            case GLFW_KEY_ESCAPE: ascii_key = 0x1b; break;
+            default: return; // Key not handled
+        }
+
+        cam->processKey(ascii_key);
+    }
+}
+
+void GLFWFramebufferSizeCallback(GLFWwindow* window, int width, int height) {
+    // Retrieve FluidSimulatorApp from window user pointer
+    FluidSimulation::FluidSimulatorApp* app = static_cast<FluidSimulation::FluidSimulatorApp*>(glfwGetWindowUserPointer(window));
+    if (!app) return;
+
+    // Update camera aspect ratio
+    FluidSimulation::Camera* cam = app->getCamera();
+    if (cam) cam->reshape(width, height);
+
+    // Resize FBO
+    FluidSimulation::FluidRenderer* renderer = app->getRenderer(); // Add this getter
+    if (renderer) renderer->resizeFBO(width, height);
+}
+
+
+int main(int argc, char** argv) {
+    // Initialize GLFW
+    if (!glfwInit()) {
+        fprintf(stderr, "Failed to initialize GLFW\n");
+        return -1;
+    }
+
+    // Set OpenGL version (3.3 Core)
+    glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 3);
+    glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 3);
+    glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE);
+
+    // Create GLFW window
+    GLFWwindow* window = glfwCreateWindow(1280, 720, "Simulation with ImGui Docking & FBO", NULL, NULL);
+    if (!window) {
+        fprintf(stderr, "Failed to create GLFW window\n");
+        glfwTerminate();
+        return -1;
+    }
+
+    glfwMakeContextCurrent(window);
+    glfwSwapInterval(1); // Enable V-sync
+
+    // Set initial viewport
+    int initialWidth = 1280, initialHeight = 720;
+    glViewport(0, 0, initialWidth, initialHeight);
+
+
+    // Initialize FBO with initial window size (1280x720)
+
+    // Instantiate ImGuiManager
+    IMGUI_CHECKVERSION();
+    ImGui::CreateContext();
+    ImGuiIO& io = ImGui::GetIO(); (void)io;
+
+    // Enable Docking
+    io.ConfigFlags |= ImGuiConfigFlags_DockingEnable;
+
+    // Optional: Enable Multi-Viewport / Platform Windows
+    // io.ConfigFlags |= ImGuiConfigFlags_ViewportsEnable;
+
+    // Setup Dear ImGui style
+    ImGui::StyleColorsDark();
+    // ImGui::StyleColorsClassic();
+
+    // Setup Platform/Renderer backends
+    ImGui_ImplGlfw_InitForOpenGL(window, true);
+    ImGui_ImplOpenGL3_Init("#version 330");
+
+    // Set GLFW callbacks
+    //glfwSetWindowUserPointer(window, cameraPtr);
+    //glfwSetKeyCallback(window, GLFWKeyCallback);
+    //glfwSetFramebufferSizeCallback(window, GLFWFramebufferSizeCallback);
+    FluidSimulation::appSettings appSettings = {1280, 720};
+    FluidSimulation::FluidSimulatorApp* fluidSimulator = FluidSimulation::creatFluidSimulatorApp(&appSettings);
+
+    // Set GLFW callbacks
+    glfwSetWindowUserPointer(window, fluidSimulator);
+    glfwSetKeyCallback(window, GLFWKeyCallback);
+    glfwSetFramebufferSizeCallback(window, GLFWFramebufferSizeCallback);
+
+    // Main loop
+    while (!glfwWindowShouldClose(window)) {
+        // Poll events
+        glfwPollEvents();
+
+        // Start ImGui frame
+        ImGui_ImplOpenGL3_NewFrame();
+        ImGui_ImplGlfw_NewFrame();
+        ImGui::NewFrame();
+
+        // Obtain ImGuiIO reference
+        //ImGuiIO& io = imguiManager.getIO();
+
+        // Create a full-screen dockspace
+        ImGui::SetNextWindowPos(ImVec2(0, 0));
+        ImGui::SetNextWindowSize(io.DisplaySize);
+        ImGuiWindowFlags dockspace_flags = ImGuiWindowFlags_NoTitleBar | ImGuiWindowFlags_NoCollapse |
+                                           ImGuiWindowFlags_NoResize | ImGuiWindowFlags_NoMove |
+                                           ImGuiWindowFlags_NoBringToFrontOnFocus | ImGuiWindowFlags_NoNavFocus;
+        
+        ImGui::Begin("MainDockSpaceWindow", NULL, dockspace_flags);
+        ImGuiID dockspace_id = ImGui::GetID("MainDockSpace");
+        ImGui::DockSpace(dockspace_id, ImVec2(0.0f, 0.0f), ImGuiDockNodeFlags_None);
+        ImGui::End();
+
+        // Render Applicatiion Here
+        fluidSimulator->Render();
+
+        // Move control to app laters
+        ImGui::Begin("Controls");
+
+        ImGui::Text("Fluid Simulation Parameters");
+
+        ImGui::End();
+        // Render ImGui
+        //imguiManager.render();
+        ImGui::Render();
+        ImGui_ImplOpenGL3_RenderDrawData(ImGui::GetDrawData());
+
+        // Swap buffers
+        glfwSwapBuffers(window);
+    }
+
+    // Cleanup ImGui and OpenGL resources is handled by ImGuiManager destructor
+
+    // Cleanup your own resources
+    //cleanup();
+    //glDeleteFramebuffers(1, &fbo);
+    //glDeleteTextures(1, &fboTexture);
+    //glDeleteRenderbuffers(1, &depthRBO);
+    delete fluidSimulator;
+    // Terminate GLFW
+    glfwDestroyWindow(window);
+    glfwTerminate();
+    return 0;
+}
\ No newline at end of file
diff --git a/src/include/Camera.h b/src/include/Camera.h
index 8c96803..f1d9a94 100644
--- a/src/include/Camera.h
+++ b/src/include/Camera.h
@@ -1,35 +1,35 @@
-// Camera.h
-#ifndef CAMERA_H
-#define CAMERA_H
+#pragma once
 
 #include "VectorUtils4.h"
 
-class Camera {
-private:
-    vec3 campos;
-    vec3 forward;
-    vec3 upVector;
-    float fov;
-    float aspectRatio;
-    float zNear;
-    float zFar;
-    bool displayBorder = false;
-    mat4 rotationMatrix;
+namespace FluidSimulation {
 
-    void rotateSide(float angle);
-    void rotateUp(float angle);
+    class Camera {
+    private:
+        vec3 campos;
+        vec3 forward;
+        vec3 upVector;
+        float fov;
+        float aspectRatio;
+        float zNear;
+        float zFar;
+        bool displayBorder = false;
+        mat4 rotationMatrix;
 
-public:
-    Camera(vec3 position, vec3 target, vec3 up, float fov, float aspect, float zNear, float zFar);
-    void processKey(unsigned char key);
-    void reshape(int width, int height);
-    mat4 getViewMatrix() const;
-    mat4 getProjectionMatrix() const;
-    bool shouldDisplayBorder() const { return displayBorder; }
-    void toggleDisplayBorder() { displayBorder = !displayBorder; }
-    void setAspectRatio(float aspect) { aspectRatio = aspect; }
-};
+        void rotateSide(float angle);
+        void rotateUp(float angle);
 
-void initCamera(Camera& camera);
+    public:
+        Camera(vec3 position, vec3 target, vec3 up, float fov, float aspect, float zNear, float zFar);
+        void processKey(unsigned char key);
+        void reshape(int width, int height);
+        mat4 getViewMatrix() const;
+        mat4 getProjectionMatrix() const;
+        bool shouldDisplayBorder() const { return displayBorder; }
+        void toggleDisplayBorder() { displayBorder = !displayBorder; }
+        void setAspectRatio(float aspect) { aspectRatio = aspect; }
+    };
 
-#endif
+    void initCamera(Camera& camera);
+
+} // namespace FluidSimulation
diff --git a/src/include/FluidRenderer.h b/src/include/FluidRenderer.h
new file mode 100644
index 0000000..d62325a
--- /dev/null
+++ b/src/include/FluidRenderer.h
@@ -0,0 +1,54 @@
+#pragma once
+
+#include <GL/glew.h>
+#include <vector>
+#include <cuda_runtime.h>
+#include "include/Camera.h"
+
+namespace FluidSimulation {
+
+    class FluidRenderer {
+    private:
+        GLuint m_shaderProgram;
+
+        GLuint m_positionAttribLoc;
+        GLuint m_modelMatrixLoc;
+        GLuint m_viewMatrixLoc;
+        GLuint m_projectionMatrixLoc;
+        GLuint m_pointRadiusLoc;
+        GLuint m_pointScaleLoc;
+        GLuint m_uColorLoc;
+        
+        GLuint m_particleVBO;
+        GLuint m_boundaryVBO;
+
+        GLuint m_fbo;
+        GLuint m_fboTexture;
+        GLuint m_depthRBO;
+
+        int m_fbWidth;
+        int m_fbHeight;
+
+        Camera* m_camera;
+
+        void initGLEW();          // Initialize GLEW
+        void initOpenGLSettings(); // Set up OpenGL state
+        void initShaders();       // Load and configure shaders
+        void initFBOs();          // Set up framebuffers
+        void initShaderUniforms(); // Configure shader uniforms
+
+        void renderParticles(size_t particleCount);
+        void renderBoundaries(size_t boundaryCount);
+
+    public:
+        FluidRenderer(int width, int height, Camera* camera);
+        ~FluidRenderer();
+
+        void init();              // Consolidated initialization function
+        void resizeFBO(int newWidth, int newHeight);
+        void setVBOs(GLuint particleVBO, GLuint boundaryVBO);
+        void renderFrame(size_t particleCount, size_t boundaryCount);
+        GLuint getFBOTexture() const { return m_fboTexture; }
+    };
+
+}
diff --git a/src/include/FluidSimulation.h b/src/include/FluidSimulation.h
new file mode 100644
index 0000000..81cfe83
--- /dev/null
+++ b/src/include/FluidSimulation.h
@@ -0,0 +1,32 @@
+#pragma once
+
+#include <vector>
+#include <cuda_runtime.h>
+#include <GL/glew.h>
+
+namespace FluidSimulation {
+
+    class FluidSimulation {
+    private:
+        GLuint m_vbo_particles;
+        GLuint m_vbo_boundary;
+
+        int m_numParticles;
+        int m_numBoundaryParticles;
+
+        float m_timeStep;
+
+    public:
+        FluidSimulation(int numParticles);
+        ~FluidSimulation();
+
+        void initSimulation();  // Initialize simulation resources
+        void updateSimulation();  // Update particles using CUDA
+        GLuint getParticleVBO() const { return m_vbo_particles; }
+        GLuint getBoundaryVBO() const { return m_vbo_boundary; }
+
+        int getParticleCount() const { return m_numParticles; };
+        int getBoundaryParticleCount() const { return m_numBoundaryParticles; };
+    };
+
+}
diff --git a/src/include/FluidSimulationApp.h b/src/include/FluidSimulationApp.h
new file mode 100644
index 0000000..b7e083c
--- /dev/null
+++ b/src/include/FluidSimulationApp.h
@@ -0,0 +1,36 @@
+//#include "include/FluidConfig.h"
+#include "include/FluidRenderer.h"
+#include "include/FluidSimulation.h"
+#include "Camera.h"
+
+namespace FluidSimulation
+{
+    struct appSettings {
+        int width;
+        int height;
+    };
+
+    //FluidConfig fluidConfig;
+
+
+    class FluidSimulatorApp
+    {
+    private:
+        appSettings* m_appsettings;
+
+        Camera* m_camera;
+        FluidRenderer* m_fluidRenderer;
+        FluidSimulation* m_fluidSimulation;
+
+    public:
+        FluidSimulatorApp(appSettings* appSettings);
+        ~FluidSimulatorApp();
+        
+        void Render();
+
+        Camera* getCamera() const { return m_camera; }
+        FluidRenderer* getRenderer() const { return m_fluidRenderer; };
+    };
+    
+    FluidSimulatorApp* creatFluidSimulatorApp(appSettings* settings);
+}
\ No newline at end of file
diff --git a/src/include/PCISPH_sim.h b/src/include/PCISPH_sim.h
index 4ee42e6..96ea285 100644
--- a/src/include/PCISPH_sim.h
+++ b/src/include/PCISPH_sim.h
@@ -1,4 +1,4 @@
-// kernels_cuda.h
+// PCISPH_sim.h
 #ifndef KERNELS_CUDA_H
 #define KERNELS_CUDA_H
 
diff --git a/src/kernels_cuda.cu b/src/kernels_cuda.cu
deleted file mode 100644
index d7dca02..0000000
--- a/src/kernels_cuda.cu
+++ /dev/null
@@ -1,606 +0,0 @@
-// kernels_cuda.cu
-#include "kernels_cuda.h"
-#include "constants.h"
-#include "include/float3_helper_math.h"
-#include <cuda_runtime.h>
-#include <cuda_gl_interop.h>
-#include <device_launch_parameters.h>
-#include <cmath>
-
-// Include Thrust for sorting
-#include <thrust/device_ptr.h>
-#include <thrust/sort.h>
-#include <thrust/sequence.h>
-
-// CUDA Graphics Resource (for OpenGL-CUDA interoperability)
-cudaGraphicsResource_t cuda_vbo_resource;
-
-// Constants
-__constant__ float PI;
-__constant__ float H;
-__constant__ float H2;
-__constant__ float POLY6_CONST;
-__constant__ float SPIKY_GRAD_CONST;
-__constant__ float VISC_LAP_CONST;
-
-__constant__ float PARTICLE_MASS = 1.0f;
-__constant__ float REST_DENSITY = 60.0f;
-__constant__ float GAS_CONSTANT;
-__constant__ float VISCOSITY;
-
-// Boundary constants
-__constant__ float BOUND_X_MIN;
-__constant__ float BOUND_X_MAX;
-__constant__ float BOUND_Y_MIN;
-__constant__ float BOUND_Y_MAX;
-__constant__ float BOUND_Z_MIN;
-__constant__ float BOUND_Z_MAX;
-
-// Device arrays
-float3* d_positions; // Positions mapped from OpenGL VBO
-float3* d_velocities;
-
-unsigned int* d_cellIndices;
-unsigned int* d_particleIndices;
-
-float3* d_positions_sorted;
-float3* d_velocities_sorted;
-
-unsigned int* d_cellStart;
-unsigned int* d_cellEnd;
-
-// Arrays for densities, pressures, accelerations
-float* d_densities;
-float* d_pressures;
-float3* d_accelerations;
-
-// Smoothing kernels
-__device__ float poly6Kernel(float r2)
-{
-    if (r2 >= 0.0f && r2 <= H2) {
-        float term = H2 - r2;
-        return POLY6_CONST * term * term * term;
-    } else {
-        return 0.0f;
-    }
-}
-
-__device__ float spikyKernelGrad(float r)
-{
-    if (r > 0.0f && r <= H) {
-        float term = H - r;
-        return SPIKY_GRAD_CONST * term * term;
-    } else {
-        return 0.0f;
-    }
-}
-
-__device__ float viscosityKernelLap(float r)
-{
-    if (r >= 0.0f && r <= H) {
-        return VISC_LAP_CONST * (H - r);
-    } else {
-        return 0.0f;
-    }
-}
-
-// Compute cell indices
-__global__ void computeCellIndices(
-    float3* positions,
-    unsigned int* cellIndices,
-    int numParticles,
-    float3 gridMin,
-    float3 gridCellSize,
-    int gridSizeX, int gridSizeY, int gridSizeZ
-)
-{
-    int idx = blockIdx.x * blockDim.x + threadIdx.x;
-
-    if (idx < numParticles) {
-        float3 pos = positions[idx];
-
-        int ix = (int)((pos.x - gridMin.x) / gridCellSize.x);
-        int iy = (int)((pos.y - gridMin.y) / gridCellSize.y);
-        int iz = (int)((pos.z - gridMin.z) / gridCellSize.z);
-
-        // Clamp indices to grid bounds
-        ix = max(0, min(ix, gridSizeX - 1));
-        iy = max(0, min(iy, gridSizeY - 1));
-        iz = max(0, min(iz, gridSizeZ - 1));
-
-        unsigned int cellIndex = ix + iy * gridSizeX + iz * gridSizeX * gridSizeY;
-
-        cellIndices[idx] = cellIndex;
-    }
-}
-
-// Reordering kernel
-__global__ void reorderData(
-    unsigned int* particleIndices,
-    float3* positions,
-    float3* velocities,
-    float3* positions_sorted,
-    float3* velocities_sorted,
-    int numParticles
-)
-{
-    int idx = blockIdx.x * blockDim.x + threadIdx.x;
-
-    if (idx < numParticles) {
-        unsigned int sortedIdx = particleIndices[idx];
-
-        positions_sorted[idx] = positions[sortedIdx];
-        velocities_sorted[idx] = velocities[sortedIdx];
-    }
-}
-
-// Find cell start and end indices
-// can bli reeplaceed med atomicAdd för varje ceell index??
-__global__ void findCellStartEnd(
-    unsigned int* cellIndices,
-    unsigned int* cellStart,
-    unsigned int* cellEnd,
-    int numParticles
-)
-{
-    int idx = blockIdx.x * blockDim.x + threadIdx.x;
-
-    if (idx < numParticles) {
-        unsigned int cellIndex = cellIndices[idx];
-
-        if (idx == 0) {
-            cellStart[cellIndex] = 0;
-        } else {
-            unsigned int prevCellIndex = cellIndices[idx - 1];
-            if (prevCellIndex != cellIndex) {
-                cellStart[cellIndex] = idx;
-                cellEnd[prevCellIndex] = idx;
-            }
-        }
-
-        if (idx == numParticles - 1) {
-            cellEnd[cellIndex] = numParticles;
-        }
-    }
-}
-
-// Compute density and pressure
-__global__ void computeDensityPressure(
-    float3* positions,
-    float* densities,
-    float* pressures,
-    unsigned int* cellIndices,
-    unsigned int* cellStart,
-    unsigned int* cellEnd,
-    int numParticles,
-    int gridSizeX, int gridSizeY, int gridSizeZ
-)
-{
-    int idx = blockIdx.x * blockDim.x + threadIdx.x;
-
-    if (idx < numParticles) {
-        float3 pos_i = positions[idx];
-        float density = 0.0f;
-
-        unsigned int cellIndex = cellIndices[idx];
-
-        // Compute cell coordinates
-        int cellX = cellIndex % gridSizeX;
-        int cellY = (cellIndex / gridSizeX) % gridSizeY;
-        int cellZ = cellIndex / (gridSizeX * gridSizeY);
-
-        // Iterate over neighboring cells
-        for (int dz = -1; dz <= 1; ++dz) {
-            int nz = cellZ + dz;
-            if (nz < 0 || nz >= gridSizeZ) continue;
-            for (int dy = -1; dy <= 1; ++dy) {
-                int ny = cellY + dy;
-                if (ny < 0 || ny >= gridSizeY) continue;
-                for (int dx = -1; dx <= 1; ++dx) {
-                    int nx = cellX + dx;
-                    if (nx < 0 || nx >= gridSizeX) continue;
-
-                    unsigned int neighborCellIndex = nx + ny * gridSizeX + nz * gridSizeX * gridSizeY;
-
-                    unsigned int startIdx = cellStart[neighborCellIndex];
-                    unsigned int endIdx = cellEnd[neighborCellIndex];
-
-                    if (startIdx != 0xFFFFFFFF) {
-                        for (unsigned int j = startIdx; j < endIdx; ++j) {
-                            float3 pos_j = positions[j];
-
-                            float3 r = pos_j - pos_i;
-                            float r2 = r.x * r.x + r.y * r.y + r.z * r.z;
-
-                            density += PARTICLE_MASS * poly6Kernel(r2);
-                        }
-                    }
-                }
-            }
-        }
-
-        densities[idx] = density;
-        //pressures[idx] = GAS_CONSTANT * fmaxf((density - REST_DENSITY), 0.0);
-        float dens = density/REST_DENSITY;
-        pressures[idx] = GAS_CONSTANT * ((dens*dens*dens*dens*dens*dens*dens)-1.0);
-        //pressures[idx] = GAS_CONSTANT * min((density - REST_DENSITY), 0.0);
-        //pressures[idx] = GAS_CONSTANT * (density - REST_DENSITY);
-    }
-}
-
-// Compute forces
-__global__ void computeForces(
-    float3* positions,
-    float3* velocities,
-    float* densities,
-    float* pressures,
-    float3* accelerations,
-    unsigned int* cellIndices,
-    unsigned int* cellStart,
-    unsigned int* cellEnd,
-    int numParticles,
-    int gridSizeX, int gridSizeY, int gridSizeZ
-)
-{
-    int idx = blockIdx.x * blockDim.x + threadIdx.x;
-
-    if (idx < numParticles) {
-        float3 pos_i = positions[idx];
-        float3 vel_i = velocities[idx];
-        float density_i = densities[idx];
-        float pressure_i = pressures[idx];
-
-        float3 fPressure = make_float3(0.0f, 0.0f, 0.0f);
-        float3 fViscosity = make_float3(0.0f, 0.0f, 0.0f);
-        float3 fGravity = make_float3(0.0f, -4.0f , 0.0f);
-
-        unsigned int cellIndex = cellIndices[idx];
-
-        // Compute cell coordinates
-        int cellX = cellIndex % gridSizeX;
-        int cellY = (cellIndex / gridSizeX) % gridSizeY;
-        int cellZ = cellIndex / (gridSizeX * gridSizeY);
-
-        // Iterate over neighboring cells
-        for (int dz = -1; dz <= 1; ++dz) {
-            int nz = cellZ + dz;
-            if (nz < 0 || nz >= gridSizeZ) continue;
-            for (int dy = -1; dy <= 1; ++dy) {
-                int ny = cellY + dy;
-                if (ny < 0 || ny >= gridSizeY) continue;
-                for (int dx = -1; dx <= 1; ++dx) {
-                    int nx = cellX + dx;
-                    if (nx < 0 || nx >= gridSizeX) continue;
-
-                    unsigned int neighborCellIndex = nx + ny * gridSizeX + nz * gridSizeX * gridSizeY;
-
-                    unsigned int startIdx = cellStart[neighborCellIndex];
-                    unsigned int endIdx = cellEnd[neighborCellIndex];
-
-                    if (startIdx != 0xFFFFFFFF) {
-                        for (unsigned int j = startIdx; j < endIdx; ++j) {
-                            if (j == idx) continue;
-
-                            float3 pos_j = positions[j];
-                            float3 r = (pos_i - pos_j);
-                            
-                            float distSq = r.x * r.x + r.y * r.y + r.z * r.z;
-
-                            if (distSq < H2 && distSq > 0.0001f) {
-                                float dist = sqrtf(distSq);
-
-                                // Pressure force
-                                float spikyGrad = spikyKernelGrad(dist);
-
-                                // Cache pressures[j]
-                                float density_j = densities[j];
-                                float pressure_j = pressures[j];
-
-                                // avreged pressure scaleed by inverse j density
-                                /*
-                                float pressureTerm = - PARTICLE_MASS * 0.5 * (pressure_i + pressure_j) / density_j * spikyGrad;
-                                fPressure += pressureTerm*(r/dist); // r/dist gives direction
-                                */
-
-                                // Symmetric pressure -- can cahnge to inverse dennsity in density/pressure calc to avoid devide by density here
-                                // Add particle_mass to pressure term according to equations !!! Added down below in aceleration calcualtion
-                                float pressureTerm = - (pressure_i / (density_i * density_i) + pressure_j / (density_j * density_j)) * spikyGrad;
-                                float3 pressureForce = pressureTerm * (r / dist); // Directional force
-                                fPressure += pressureForce;
-
-                                // Viscosity force
-                                float3 velDiff = velocities[j] - vel_i;
-                                float viscLap = viscosityKernelLap(dist);
-                                fViscosity  += VISCOSITY * velDiff * viscLap / density_j * PARTICLE_MASS;
-                            }
-                        }
-                    }
-                }
-            }
-        }
-
-        // Total acceleration
-        // According to paper v_i(t+dt) = v_(i+1) + (dt/m_i)*Fp_i
-        // While external and viscocisty force is calculated before But mass and nnot deeensity is used to deevidee force
-        // Deviding by mass, looks a lot better
-        accelerations[idx] = (fPressure * PARTICLE_MASS + fViscosity + fGravity) / PARTICLE_MASS;
-        //accelerations[idx] = (fPressure * PARTICLE_MASS + fViscosity + fGravity) / density_i;
-    }
-}
-
-// Integration
-__global__ void integrate(
-    float3* positions,
-    float3* velocities,
-    float3* accelerations,
-    float timeStep,
-    int numParticles
-)
-{
-    int idx = blockIdx.x * blockDim.x + threadIdx.x;
-
-    if (idx < numParticles) {
-        // Update velocity
-        float3 vel_i = velocities[idx] + (accelerations[idx] * timeStep);
-        velocities[idx] = vel_i;
-
-        // Update position
-        positions[idx] += vel_i * timeStep;
-    }
-}
-
-// Enforce boundary conditions
-__global__ void enforceBoundary(
-    float3* positions,
-    float3* velocities,
-    int numParticles,
-    float damping
-)
-{
-    int idx = blockIdx.x * blockDim.x + threadIdx.x;
-
-    if (idx < numParticles) {
-        // X-axis
-        if (positions[idx].x < BOUND_X_MIN) {
-            positions[idx].x = BOUND_X_MIN;
-            velocities[idx].x *= -damping;
-        }
-        if (positions[idx].x > BOUND_X_MAX) {
-            positions[idx].x = BOUND_X_MAX; // Should account for particle radius!!
-            velocities[idx].x *= -damping;
-        }
-
-        // Y-axis
-        if (positions[idx].y < BOUND_Y_MIN) {
-            positions[idx].y = BOUND_Y_MIN;
-            velocities[idx].y *= -damping;
-        }
-        if (positions[idx].y > BOUND_Y_MAX) {
-            positions[idx].y = BOUND_Y_MAX;
-            velocities[idx].y *= -damping;
-        }
-
-        // Z-axis
-        if (positions[idx].z < BOUND_Z_MIN) {
-            positions[idx].z = BOUND_Z_MIN;
-            velocities[idx].z *= -damping;
-        }
-        if (positions[idx].z > BOUND_Z_MAX) {
-            positions[idx].z = BOUND_Z_MAX;
-            velocities[idx].z *= -damping;
-        }
-    }
-}
-
-void updateParticles(float deltaTime)
-{
-    // Map the VBO resource to get a device pointer
-    size_t num_bytes;
-    cudaGraphicsMapResources(1, &cuda_vbo_resource, 0);
-    cudaGraphicsResourceGetMappedPointer((void**)&d_positions, &num_bytes, cuda_vbo_resource);
-
-    int threadsPerBlock = 256;
-    int blocksPerGrid = (NUM_PARTICLES + threadsPerBlock - 1) / threadsPerBlock;
-
-    // Compute grid size based on smoothing radius
-    float3 gridMin = make_float3(h_BOUND_X_MIN, h_BOUND_Y_MIN, h_BOUND_Z_MIN);
-    float3 gridMax = make_float3(h_BOUND_X_MAX, h_BOUND_Y_MAX, h_BOUND_Z_MAX);
-    float3 gridCellSize = make_float3(h_H, h_H, h_H);
-
-    int gridSizeX = (int)ceil((h_BOUND_X_MAX - h_BOUND_X_MIN) / h_H);
-    int gridSizeY = (int)ceil((h_BOUND_Y_MAX - h_BOUND_Y_MIN) / h_H);
-    int gridSizeZ = (int)ceil((h_BOUND_Z_MAX - h_BOUND_Z_MIN) / h_H);
-    int maxCells = gridSizeX * gridSizeY * gridSizeZ;
-
-    // Allocate or resize cellStart and cellEnd arrays if necessary
-    // (You may need to adjust your initialization code to handle dynamic sizes)
-
-    // Compute cell indices
-    computeCellIndices<<<blocksPerGrid, threadsPerBlock>>>(
-        d_positions, d_cellIndices, NUM_PARTICLES, gridMin, gridCellSize, gridSizeX, gridSizeY, gridSizeZ);
-    cudaDeviceSynchronize();
-
-    // Wrap device pointers with thrust iterators
-    thrust::device_ptr<unsigned int> dev_cellIndices_ptr(d_cellIndices);
-    thrust::device_ptr<unsigned int> dev_particleIndices_ptr(d_particleIndices);
-
-    // Initialize particle indices
-    thrust::sequence(dev_particleIndices_ptr, dev_particleIndices_ptr + NUM_PARTICLES);
-
-    // Sort particle indices based on cell indices
-    thrust::sort_by_key(dev_cellIndices_ptr, dev_cellIndices_ptr + NUM_PARTICLES, dev_particleIndices_ptr);
-
-    // Reorder positions and velocities
-    reorderData<<<blocksPerGrid, threadsPerBlock>>>(
-        d_particleIndices,
-        d_positions,
-        d_velocities,
-        d_positions_sorted,
-        d_velocities_sorted,
-        NUM_PARTICLES);
-    cudaDeviceSynchronize();
-
-    // Swap sorted data with original data
-    float3* temp_positions = d_positions;
-    d_positions = d_positions_sorted;
-    d_positions_sorted = temp_positions;
-
-    float3* temp_velocities = d_velocities;
-    d_velocities = d_velocities_sorted;
-    d_velocities_sorted = temp_velocities;
-
-    // Build cell start and end indices
-    // Initialize cellStart and cellEnd
-    cudaMemset(d_cellStart, 0xFF, maxCells * sizeof(unsigned int));
-    cudaMemset(d_cellEnd, 0xFF, maxCells * sizeof(unsigned int));
-
-    findCellStartEnd<<<blocksPerGrid, threadsPerBlock>>>(
-        d_cellIndices,
-        d_cellStart,
-        d_cellEnd,
-        NUM_PARTICLES);
-    cudaDeviceSynchronize();
-
-    // Compute density and pressure
-    computeDensityPressure<<<blocksPerGrid, threadsPerBlock>>>(
-        d_positions,
-        d_densities,
-        d_pressures,
-        d_cellIndices,
-        d_cellStart,
-        d_cellEnd,
-        NUM_PARTICLES,
-        gridSizeX, gridSizeY, gridSizeZ);
-    cudaDeviceSynchronize();
-
-    // Compute forces
-    computeForces<<<blocksPerGrid, threadsPerBlock>>>(
-        d_positions,
-        d_velocities,
-        d_densities,
-        d_pressures,
-        d_accelerations,
-        d_cellIndices,
-        d_cellStart,
-        d_cellEnd,
-        NUM_PARTICLES,
-        gridSizeX, gridSizeY, gridSizeZ);
-    cudaDeviceSynchronize();
-
-    // Integrate
-    float timeStep = 0.001f;
-    integrate<<<blocksPerGrid, threadsPerBlock>>>(
-        d_positions,
-        d_velocities,
-        d_accelerations,
-        timeStep,
-        //deltaTime,
-        NUM_PARTICLES);
-    cudaDeviceSynchronize();
-
-    // Enforce boundary conditions
-    float damping = 0.7f;
-    enforceBoundary<<<blocksPerGrid, threadsPerBlock>>>(
-        d_positions,
-        d_velocities,
-        NUM_PARTICLES,
-        damping);
-    cudaDeviceSynchronize();
-
-    // Unmap the VBO resource
-    cudaGraphicsUnmapResources(1, &cuda_vbo_resource, 0);
-}
-
-void initializeConstants()
-{
-    // Host constants
-    const float h_H2 = h_H * h_H;
-    const float h_POLY6_CONST = 315.0f / (64.0f * h_PI * powf(h_H, 9));
-    const float h_SPIKY_GRAD_CONST = -45.0f / (h_PI * powf(h_H, 6));
-    const float h_VISC_LAP_CONST = 45.0f / (h_PI * powf(h_H, 6));
-
-    // Copy constants to device
-    cudaMemcpyToSymbol(PI, &h_PI, sizeof(float));
-    cudaMemcpyToSymbol(H, &h_H, sizeof(float));
-    cudaMemcpyToSymbol(H2, &h_H2, sizeof(float));
-    cudaMemcpyToSymbol(POLY6_CONST, &h_POLY6_CONST, sizeof(float));
-    cudaMemcpyToSymbol(SPIKY_GRAD_CONST, &h_SPIKY_GRAD_CONST, sizeof(float));
-    cudaMemcpyToSymbol(VISC_LAP_CONST, &h_VISC_LAP_CONST, sizeof(float));
-
-    cudaMemcpyToSymbol(GAS_CONSTANT, &h_GAS_CONSTANT, sizeof(float));
-    cudaMemcpyToSymbol(VISCOSITY, &h_VISCOSITY, sizeof(float));
-
-    // Copy boundary constants to device
-    cudaMemcpyToSymbol(BOUND_X_MIN, &h_BOUND_X_MIN, sizeof(float));
-    cudaMemcpyToSymbol(BOUND_X_MAX, &h_BOUND_X_MAX, sizeof(float));
-    cudaMemcpyToSymbol(BOUND_Y_MIN, &h_BOUND_Y_MIN, sizeof(float));
-    cudaMemcpyToSymbol(BOUND_Y_MAX, &h_BOUND_Y_MAX, sizeof(float));
-    cudaMemcpyToSymbol(BOUND_Z_MIN, &h_BOUND_Z_MIN, sizeof(float));
-    cudaMemcpyToSymbol(BOUND_Z_MAX, &h_BOUND_Z_MAX, sizeof(float));
-}
-
-// Register VBO with CUDA
-void registerVBOWithCUDA(unsigned int vbo)
-{
-    cudaGraphicsGLRegisterBuffer(&cuda_vbo_resource, vbo, cudaGraphicsMapFlagsNone);
-}
-
-// Unregister VBO from CUDA
-void unregisterVBOFromCUDA()
-{
-    cudaGraphicsUnregisterResource(cuda_vbo_resource);
-}
-
-// Initialize device arrays
-void initializeDeviceArrays(
-    const float3* h_positions,
-    const float3* h_velocities,
-    int numParticles
-)
-{
-    // Allocate device memory for velocities
-    cudaMalloc(&d_velocities, numParticles * sizeof(float3));
-    cudaMemcpy(d_velocities, h_velocities, numParticles * sizeof(float3), cudaMemcpyHostToDevice);
-
-    // Allocate device memory for cell indices and particle indices
-    cudaMalloc(&d_cellIndices, numParticles * sizeof(unsigned int));
-    cudaMalloc(&d_particleIndices, numParticles * sizeof(unsigned int));
-
-    // Allocate device memory for sorted positions and velocities
-    cudaMalloc(&d_positions_sorted, numParticles * sizeof(float3));
-    cudaMalloc(&d_velocities_sorted, numParticles * sizeof(float3));
-
-    // Allocate cell start and end arrays
-    // Note: You need to ensure that 'maxCells' is available here or adjust the code to allocate after computing 'maxCells'
-    int gridSizeX = (int)ceil((h_BOUND_X_MAX - h_BOUND_X_MIN) / h_H);
-    int gridSizeY = (int)ceil((h_BOUND_Y_MAX - h_BOUND_Y_MIN) / h_H);
-    int gridSizeZ = (int)ceil((h_BOUND_Z_MAX - h_BOUND_Z_MIN) / h_H);
-    int maxCells = gridSizeX * gridSizeY * gridSizeZ;
-
-    cudaMalloc(&d_cellStart, maxCells * sizeof(unsigned int));
-    cudaMalloc(&d_cellEnd, maxCells * sizeof(unsigned int));
-
-    // Allocate device memory for densities, pressures, and accelerations
-    cudaMalloc(&d_densities, numParticles * sizeof(float));
-    cudaMalloc(&d_pressures, numParticles * sizeof(float));
-    cudaMalloc(&d_accelerations, numParticles * sizeof(float3));
-}
-
-// Cleanup device arrays
-void cleanupDeviceArrays()
-{
-    cudaFree(d_velocities);
-
-    cudaFree(d_cellIndices);
-    cudaFree(d_particleIndices);
-
-    cudaFree(d_positions_sorted);
-    cudaFree(d_velocities_sorted);
-
-    cudaFree(d_cellStart);
-    cudaFree(d_cellEnd);
-
-    cudaFree(d_densities);
-    cudaFree(d_pressures);
-    cudaFree(d_accelerations);
-}
\ No newline at end of file
diff --git a/src/settings/constants.h b/src/settings/constants.h
index 5f45e20..3faa9ca 100644
--- a/src/settings/constants.h
+++ b/src/settings/constants.h
@@ -18,4 +18,6 @@ const float h_REST_DENSITY = 3000.0f;
 const float h_PI = 3.14159265358979323846f;
 const float h_H = 0.35f; // Smoothing radius
 
+const int NUM_PARTICLES = 50000;
+
 #endif // CONSTANTS_H
-- 
GitLab