7

Author: Roushelfy

simulators/7_self_contact/src/BarrierEnergy.cu

@@ -370,17 +370,19 @@ T BarrierEnergy<T, dim>::init_step_size(const DeviceBuffer<T> &p)
 			   {
 				   int xI = device_bp(i / Nbe);
 				   int eI0 = device_be(2*(i % Nbe)), eI1 = device_be(2*(i % Nbe) + 1);
-				   Eigen::Matrix<T, 2, 1> p, e0, e1,dp,de0,de1; 
+				   if (xI != eI0 && xI != eI1){
+					   Eigen::Matrix<T, 2, 1> p, e0, e1, dp, de0, de1; 
 				   p<<device_x(xI*dim),device_x(xI*dim+1);
 				   e0<<device_x(eI0*dim),device_x(eI0*dim+1);
 				   e1<<device_x(eI1*dim),device_x(eI1*dim+1);
 				   dp<<P(xI*dim),P(xI*dim+1);
 				   de0<<P(eI0*dim),P(eI0*dim+1);
 				   de1<<P(eI1*dim),P(eI1*dim+1); 
 				   if (bbox_overlap(p,e0,e1,dp,de0,de1,current_alpha)){
+					//printf("bbox_overlap at %d %d\n",xI,i);
 					 T toc=narrow_phase_CCD(p,e0,e1,dp,de0,de1,current_alpha);
 					 device_alpha1(i)=min(device_alpha1(i),toc);
-				   } })
+				   }} })
 		.wait();
 	return min(min_vector(device_alpha1), current_alpha);
 }

simulators/7_self_contact/src/NeoHookean_auto.cu

@@ -15,6 +15,7 @@ __device__ __host__ bool bbox_overlap(const Eigen::Matrix<T, 2, 1> &p, const Eig
     }
     else
     {
+        // printf("min_p: %f %f\n max_p: %f %f\n min_e: %f %f\n max_e: %f %f\n", min_p[0], min_p[1], max_p[0], max_p[1], min_e[0], min_e[1], max_e[0], max_e[1]);
         return true;
     }
 }

simulators/7_self_contact/src/distance/CCD.cu

@@ -6,7 +6,7 @@ int main()
 	double nu = 0.4, E = 1e5;
 	double Mu = E / (2 * (1 + nu)), Lam = E * nu / ((1 + nu) * (1 - 2 * nu));
 	double rho = 1000,
-		   k = 4e4, initial_stretch = 1, n_seg = 10, h = 0.01, side_len = 1, tol = 0.01, mu = 0.11;
+		   k = 4e4, initial_stretch = 1, n_seg = 3, h = 0.01, side_len = 0.45, tol = 0.01, mu = 0.11;
 	// 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);
 	InvFreeSimulator<double, 2> simulator(rho, side_len, initial_stretch, k, h, tol, mu, Mu, Lam, n_seg);
 	simulator.run();

simulators/7_self_contact/src/main.cpp

@@ -1,4 +1,5 @@
 #include "square_mesh.h"
+#include <stdio.h>
 template <typename T>
 void generate(T side_length, int n_seg, std::vector<T> &x, std::vector<int> &e)
 {
@@ -59,7 +60,7 @@ void find_boundary(const std::vector<int> &e, std::vector<int> &bp, std::vector<
 {
     // index all half-edges for fast query
     std::set<std::pair<int, int>> edge_set;
-    for (size_t i = 0; i < e.size()/3; i += 3)
+    for (size_t i = 0; i < e.size() / 3; i++)
     {
         for (int j = 0; j < 3; ++j)
         {