7

Author: Roushelfy

simulators/7_self_contact/src/BarrierEnergy.cu

@@ -160,7 +160,7 @@ const DeviceBuffer<T> &BarrierEnergy<T, dim>::grad()
 								   for (int j = 0; j < dim; j++)
 								   {
 									   T grad =device_contact_area(i) * dhat * (kappa / 2 * (log(s) / dhat + (s - 1) / d)) * device_n_ceil(j);
-									   device_grad(i * dim + j) += grad;
+									   atomicAdd(&device_grad(i * dim + j), grad);
 									   atomicAdd(&device_grad((N-1) * dim + j), -grad);
 								   }
 							   } })
@@ -381,11 +381,14 @@ T BarrierEnergy<T, dim>::init_step_size(const DeviceBuffer<T> &p)
 				   if (bbox_overlap(p, e0, e1, dp, de0, de1, current_alpha))
 				   {
 					   T toc = narrow_phase_CCD(p, e0, e1, dp, de0, de1, current_alpha);
+					   printf("toc: %f\n", toc);
 					   device_alpha1(i) = min(device_alpha1(i), toc);
 				   }
 				   } })
 		.wait();
-	return min(min_vector(device_alpha1), current_alpha);
+	T a = min(min_vector(device_alpha1), current_alpha);
+	printf("alpha: %f\n", a);
+	return a;
 }
 template class BarrierEnergy<float, 2>;
 template class BarrierEnergy<float, 3>;

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

@@ -34,7 +34,10 @@ __device__ __host__ T narrow_phase_CCD(Eigen::Matrix<T, 2, 1> p, Eigen::Matrix<T
     {
         return toc_upperbound;
     }
-
+    printf("maxDispMag: %f\n", maxDispMag);
+    printf("dp: %f, %f\n", dp[0], dp[1]);
+    printf("de0: %f, %f\n", de0[0], de0[1]);
+    printf("de1: %f, %f\n", de1[0], de1[1]);
     T eta = 0.1; // calculate the toc that first brings the distance to 0.1x the current distance
     T dist2_cur = PointEdgeDistanceVal(p, e0, e1);
     T dist_cur = std::sqrt(dist2_cur);
@@ -51,6 +54,7 @@ __device__ __host__ T narrow_phase_CCD(Eigen::Matrix<T, 2, 1> p, Eigen::Matrix<T
         e1 += tocLowerBound * de1;
         dist2_cur = PointEdgeDistanceVal(p, e0, e1);
         dist_cur = std::sqrt(dist2_cur);
+        printf("dist_cur: %f\n", dist_cur);
         if (toc != 0 && dist_cur < gap)
         {
             break;

simulators/7_self_contact/src/main.cpp

@@ -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 = 3, h = 0.01, side_len = 0.45, tol = 0.01, mu = 0.11;
+		   k = 4e4, initial_stretch = 1, n_seg = 1, 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();