bug fix

Author: Roushelfy

simulators/1_mass_spring/include/SparseMatrix.h

@@ -1,6 +1,5 @@
 #pragma once
 #include <Eigen/Dense>
-#include "SparseMatrix.h"
 #include <muda/cub/device/device_reduce.h>
 #include <muda/container.h>
 #include <muda/muda.h>

simulators/1_mass_spring/include/uti.h

@@ -137,6 +137,7 @@ const DeviceTripletMatrix<T, 1> &MassSpringEnergy<T, dim>::hess()
 	auto device_hess_row_idx = device_hess.row_indices();
 	auto device_hess_col_idx = device_hess.col_indices();
 	auto device_hess_val = device_hess.values();
+	device_hess_val.fill(0);
 	ParallelFor(256).apply(N, [device_x = device_x.cviewer(), device_e = device_e.cviewer(), device_l2 = device_l2.cviewer(), device_k = device_k.cviewer(), device_hess_val = device_hess_val.viewer(), device_hess_row_idx = device_hess_row_idx.viewer(), device_hess_col_idx = device_hess_col_idx.viewer(), N] __device__(int i) mutable
 						   {
 		int idx[2] = {device_e(2 * i), device_e(2 * i + 1)}; // First node index
@@ -162,7 +163,7 @@ const DeviceTripletMatrix<T, 1> &MassSpringEnergy<T, dim>::hess()
 				int indStart = i * 4*dim*dim + (ni * 2 + nj) * dim*dim;
 				for (int d1 = 0; d1 < dim; d1++)
 					for (int d2 = 0; d2 < dim; d2++){
-						device_hess_row_idx(indStart + d1 * dim + d2)= idx[ni]*dim + d1;
+						device_hess_row_idx(indStart + d1 * dim + d2)= idx[ni] * dim + d1;
 						device_hess_col_idx(indStart + d1 * dim + d2)= idx[nj] * dim + d2;
 						device_hess_val(indStart + d1 * dim + d2) = H_local(ni * dim + d1, nj * dim + d2);
 					}

simulators/1_mass_spring/src/MassSpringEnergy.cu

@@ -31,7 +31,8 @@ void __device__ make_PSD(const Eigen::Matrix<T, Size, Size> &hess, Eigen::Matrix
     Eigen::Matrix<T, Size, Size> lamDiag;
     lamDiag.setZero();
     for (int i = 0; i < Size; i++)
-        lamDiag(i, i) = lam(i);
+        if (lam(i) > 0)
+            lamDiag(i, i) = lam(i);
 
     Eigen::Matrix<T, Size, Size> VT = V.transpose();
 

simulators/1_mass_spring/src/SparseMatrix.cu

@@ -2,8 +2,7 @@
 
 int main()
 {
-	float rho = 1000, k = 1e5, initial_stretch = 1.4, n_seg = 10, h = 0.004, 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);
-	MassSpringSimulator<float, 2> simulator(rho, side_len, initial_stretch, k, h, tol, n_seg);
+	double rho = 1000, k = 1e5, initial_stretch = 1.4, n_seg = 1, h = 0.004, side_len = 1, tol = 0.01;
+	MassSpringSimulator<double, 2> simulator(rho, side_len, initial_stretch, k, h, tol, n_seg);
 	simulator.run();
 }

simulators/1_mass_spring/src/device_uti.cu

@@ -79,6 +79,7 @@ void MassSpringSimulator<T, dim>::run()
     assert(dim == 2);
     bool running = true;
     auto &window = pimpl_->window;
+    int time_step = 0;
     while (running)
     {
         sf::Event event;
@@ -89,7 +90,7 @@ void MassSpringSimulator<T, dim>::run()
         }
 
         pimpl_->draw(); // Draw the current state
-
+        std::cout << "Time step " << time_step++ << "\n";
         // Update the simulation state
         pimpl_->step_forward();
 
@@ -113,6 +114,7 @@ void MassSpringSimulator<T, dim>::Impl::step_forward()
     // std::cout << "Initial residual " << residual << "\n";
     while (residual > tol)
     {
+        std::cout << "Iteration " << iter << " residual " << residual << "E_last" << E_last << "\n";
         // Line search
         T alpha = 1;
         DeviceBuffer<T> x0 = x;
@@ -122,9 +124,8 @@ void MassSpringSimulator<T, dim>::Impl::step_forward()
             alpha /= 2;
             update_x(add_vector<T>(x0, p, 1.0, alpha));
         }
-        // std::cout << "step size = " << alpha << "\n";
+        std::cout << "step size = " << alpha << "\n";
         E_last = IP_val();
-        // std::cout << "Iteration " << iter << " residual " << residual << "E_last" << E_last << "\n";
         p = search_direction();
         residual = max_vector(p) / h;
         iter += 1;

simulators/1_mass_spring/src/main.cpp

@@ -1,6 +1,5 @@
 #pragma once
 #include <Eigen/Dense>
-#include "SparseMatrix.h"
 #include <muda/cub/device/device_reduce.h>
 #include <muda/container.h>
 #include <muda/muda.h>