4_friction

Author: Roushelfy

simulators/2_dirichlet/src/main.cpp

@@ -2,7 +2,7 @@
 
 int main()
 {
-	float rho = 1000, k = 1e3, initial_stretch = 1, n_seg = 10, h = 0.02, side_len = 1, tol = 0.01;
+	float rho = 1000, k = 1e3, initial_stretch = 1, n_seg = 20, h = 0.02, side_len = 1, tol = 0.01;
 	// printf("Running mass-spring simulator with parameters: rho = %f, k = %f, initial_stretch = %f, n_seg = %d, h = %f, side_len = %f, tol = %f\n", rho, k, initial_stretch, n_seg, h, side_len, tol);
 	DirichletSimulator<float, 2> simulator(rho, side_len, initial_stretch, k, h, tol, n_seg);
 	simulator.run();

simulators/3_contact/CMakeLists.txt

@@ -0,0 +1,18 @@
+add_executable(3_contact)
+
+target_compile_options(3_contact PRIVATE -g)
+set_target_properties(3_contact PROPERTIES CUDA_SEPARABLE_COMPILATION ON)
+
+target_link_libraries(3_contact PRIVATE muda cusolver cublas cusparse )
+
+option(BUILD_SHARED_LIBS "Build shared libraries" OFF)
+
+include_directories(include)
+
+file(GLOB_RECURSE 3_contact_CU_SOURCE CONFIGURE_DEPENDS "src/*.cu")
+target_sources(3_contact PRIVATE ${3_contact_CU_SOURCE})
+
+file(GLOB_RECURSE 3_contact_CPP_SOURCE CONFIGURE_DEPENDS "src/*.cpp")
+target_sources(3_contact PRIVATE ${3_contact_CPP_SOURCE})
+
+target_link_libraries(3_contact PRIVATE sfml-graphics)

simulators/3_contact/include/BarrierEnergy.h

@@ -0,0 +1,29 @@
+#pragma once
+
+#include <memory>
+#include <Eigen/Dense>
+#include "uti.h"
+
+template <typename T, int dim>
+class BarrierEnergy
+{
+public:
+    BarrierEnergy(const std::vector<T> &x, const std::vector<T> &contact_area, T y_ground);
+    BarrierEnergy();
+    ~BarrierEnergy();
+    BarrierEnergy(BarrierEnergy &&rhs);
+    BarrierEnergy(const BarrierEnergy &rhs);
+    BarrierEnergy &operator=(BarrierEnergy &&rhs);
+
+    void update_x(const DeviceBuffer<T> &x);
+    T val();                                    // Calculate the value of the energy
+    const DeviceBuffer<T> &grad();              // Calculate the gradient of the energy
+    const DeviceTripletMatrix<T, 1> &hess();    // Calculate the Hessian matrix of the energy
+    T init_step_size(const DeviceBuffer<T> &p); // Calculate the initial step size for the line search
+
+private:
+    // The implementation details of the VecAdder class are placed in the implementation class declared here.
+    struct Impl;
+    // The private pointer to the implementation class Impl
+    std::unique_ptr<Impl> pimpl_;
+};

simulators/3_contact/include/GravityEnergy.h

@@ -0,0 +1,35 @@
+#pragma once
+
+#include <memory>
+#include <Eigen/Dense>
+#include "square_mesh.h"
+#include <muda/muda.h>
+#include <muda/container.h>
+#include <muda/ext/linear_system.h>
+
+using namespace muda;
+
+template <typename T, int dim>
+class GravityEnergy
+{
+public:
+    GravityEnergy(int N, T m);
+    GravityEnergy();
+    ~GravityEnergy();
+    GravityEnergy(GravityEnergy &&rhs);
+    GravityEnergy(const GravityEnergy &rhs);
+    GravityEnergy &operator=(GravityEnergy &&rhs);
+    GravityEnergy &operator=(const GravityEnergy &rhs);
+
+    void update_x(const DeviceBuffer<T> &x);
+    void update_x_tilde(const DeviceBuffer<T> &x_tilde);
+    void update_m(T m);
+    T val();                       // Calculate the value of the energy
+    const DeviceBuffer<T> &grad(); // Calculate the gradient of the energy
+
+private:
+    // The implementation details of the VecAdder class are placed in the implementation class declared here.
+    struct Impl;
+    // The private pointer to the implementation class Impl
+    std::unique_ptr<Impl> pimpl_;
+};

simulators/3_contact/include/InertialEnergy.h

@@ -0,0 +1,36 @@
+#pragma once
+
+#include <memory>
+#include <Eigen/Dense>
+#include "square_mesh.h"
+#include <muda/muda.h>
+#include <muda/container.h>
+#include <muda/ext/linear_system.h>
+
+using namespace muda;
+
+template <typename T, int dim>
+class InertialEnergy
+{
+public:
+    InertialEnergy(int N, T m);
+    InertialEnergy();
+    ~InertialEnergy();
+    InertialEnergy(InertialEnergy &&rhs);
+    InertialEnergy(const InertialEnergy &rhs);
+    InertialEnergy &operator=(InertialEnergy &&rhs);
+
+    void update_x(const DeviceBuffer<T> &x);
+    void generate_hess();
+    void update_x_tilde(const DeviceBuffer<T> &x_tilde);
+    void update_m(T m);
+    T val();                                 // Calculate the value of the energy
+    const DeviceBuffer<T> &grad();           // Calculate the gradient of the energy
+    const DeviceTripletMatrix<T, 1> &hess(); // Calculate the Hessian matrix of the energy
+
+private:
+    // The implementation details of the VecAdder class are placed in the implementation class declared here.
+    struct Impl;
+    // The private pointer to the implementation class Impl
+    std::unique_ptr<Impl> pimpl_;
+};

simulators/3_contact/include/MassSpringEnergy.h

@@ -0,0 +1,31 @@
+#pragma once
+
+#include <memory>
+#include <Eigen/Dense>
+#include "uti.h"
+
+template <typename T, int dim>
+class MassSpringEnergy
+{
+public:
+    MassSpringEnergy(const std::vector<T> &x, const std::vector<int> &e, const std::vector<T> &l2, const std::vector<T> &k);
+    MassSpringEnergy();
+    ~MassSpringEnergy();
+    MassSpringEnergy(MassSpringEnergy &&rhs);
+    MassSpringEnergy(const MassSpringEnergy &rhs);
+    MassSpringEnergy &operator=(MassSpringEnergy &&rhs);
+
+    void update_x(const DeviceBuffer<T> &x);
+    void update_e(const std::vector<int> &e);
+    void update_l2(const std::vector<T> &l2);
+    void update_k(const std::vector<T> &k);
+    T val();                                 // Calculate the value of the energy
+    const DeviceBuffer<T> &grad();           // Calculate the gradient of the energy
+    const DeviceTripletMatrix<T, 1> &hess(); // Calculate the Hessian matrix of the energy
+
+private:
+    // The implementation details of the VecAdder class are placed in the implementation class declared here.
+    struct Impl;
+    // The private pointer to the implementation class Impl
+    std::unique_ptr<Impl> pimpl_;
+};

simulators/3_contact/include/SparseMatrix.h

@@ -0,0 +1,35 @@
+#pragma once
+
+#include <vector>
+
+template <typename T>
+class SparseMatrix
+{
+public:
+    SparseMatrix(int size);
+    SparseMatrix();
+    ~SparseMatrix();
+    SparseMatrix(SparseMatrix<T> &&rhs);
+    SparseMatrix<T> &operator=(SparseMatrix<T> &&rhs);
+    SparseMatrix<T> &operator=(SparseMatrix<T> &rhs);
+    SparseMatrix<T> &operator*(const T &a);
+    SparseMatrix(const SparseMatrix<T> &rhs);
+    void set_value(int row, int col, T val, int loc);
+    void set_diagonal(T val);
+
+    const std::vector<int> &get_row_buffer() const;
+    const std::vector<int> &get_col_buffer() const;
+    const std::vector<T> &get_val_buffer() const;
+    std::vector<int> &set_row_buffer();
+    std::vector<int> &set_col_buffer();
+    std::vector<T> &set_val_buffer();
+    SparseMatrix<T> &combine(const SparseMatrix<T> &other);
+
+    int get_size() const;
+
+private:
+    int size;
+    std::vector<int> row_idx;
+    std::vector<int> col_idx;
+    std::vector<T> val;
+};

simulators/3_contact/include/device_uti.h

@@ -0,0 +1,10 @@
+#include <muda/muda.h>
+#include <muda/container.h>
+#include <muda/cub/device/device_reduce.h>
+#include <Eigen/Dense>
+// utility functions
+template <typename T>
+T devicesum(const muda::DeviceBuffer<T> &buffer);
+
+template <typename T, int Size>
+void __device__ make_PSD(const Eigen::Matrix<T, Size, Size> &hess, Eigen::Matrix<T, Size, Size> &PSD);

simulators/3_contact/include/simulator.h

@@ -0,0 +1,23 @@
+#pragma once
+
+#include <vector>
+#include <cmath>
+#include "square_mesh.h"
+#include <iostream>
+template <typename T, int dim>
+class ContactSimulator
+{
+public:
+    ContactSimulator();
+    ~ContactSimulator();
+    ContactSimulator(ContactSimulator &&rhs);
+    ContactSimulator &operator=(ContactSimulator &&rhs);
+    ContactSimulator(T rho, T side_len, T initial_stretch, T K, T h, T tol, int n_seg);
+    void run();
+
+private:
+    // The implementation details of the VecAdder class are placed in the implementation class declared here.
+    struct Impl;
+    // The private pointer to the implementation class Impl
+    std::unique_ptr<Impl> pimpl_;
+};

simulators/3_contact/include/square_mesh.h

@@ -0,0 +1,7 @@
+#pragma once
+#include <vector>
+#include <memory>
+
+// Function to generate mesh points and edges
+template <typename T>
+void generate(T side_length, int n_seg, std::vector<T> &x, std::vector<int> &e);