more NH fixes

Author: Ahdhn

apps/NeoHookean/barrier_energy.h

@@ -99,7 +99,6 @@ template <typename ProblemT,
           typename T = typename VAttrT::Type>
 void ceiling_barrier_energy(ProblemT&      problem,
                             VAttrT&        contact_area,
-                            const VAttrT&  x,
                             const T        h,  // time_step
                             const vec3<T>& ground_n,
                             const vec3<T>& ground_o,
@@ -119,10 +118,10 @@ void ceiling_barrier_energy(ProblemT&      problem,
 
         const VertexHandle c1 = iter[1];
 
-        //???
-        const Eigen::Vector3<ActiveT> xi = iter_val<ActiveT, 3>(id, iter, x, 0);
+        const Eigen::Vector3<ActiveT> xi =
+            iter_val<ActiveT, 3>(id, iter, obj, 0);
 
-        const Eigen::Vector3<T> x_dbc = x.template to_eigen<3>(c1);
+        const Eigen::Vector3<T> x_dbc = obj.template to_eigen<3>(c1);
 
 
         // ceiling
@@ -138,15 +137,6 @@ void ceiling_barrier_energy(ProblemT&      problem,
                 return ActiveT(std::numeric_limits<PassiveT>::max());
             }
 
-            if constexpr (is_scalar_v<ActiveT>) {
-                printf("\n Ceil AC: T= %d, d= %f, s= %f",
-                       threadIdx.x,
-                       d.val(),
-                       s.val());
-            } else {
-                printf("\n Ceil PA: T= %d, d= %f, s= %f", threadIdx.x, d, s);
-            }
-
             E = h_sq * contact_area(c0) * dhat * T(0.5) * kappa * (s - 1) *
                 log(s);
         }

apps/NeoHookean/neo_hookean.cu

@@ -53,7 +53,7 @@ void neo_hookean(RXMeshStatic& rx, T dx)
 
     // TODO the limits and velocity should be different for different Dirichlet
     // nodes
-    const vec3<T> v_dbc_vel(0, -0.9, 0);        // Dirichlet node velocity
+    const vec3<T> v_dbc_vel(0, -1.0, 0);        // Dirichlet node velocity
     const vec3<T> v_dbc_limit(0, -0.7, 0);      // Dirichlet node limit position
     const vec3<T> ground_o(0.0f, -1.0f, 0.0f);  // a point on the slope
     const vec3<T> ground_n =
@@ -166,7 +166,7 @@ void neo_hookean(RXMeshStatic& rx, T dx)
                          kappa);
 
     ceiling_barrier_energy(
-        problem, contact_area, x, time_step, ground_n, ground_o, dhat, kappa);
+        problem, contact_area, time_step, ground_n, ground_o, dhat, kappa);
 
 
     T line_search_init_step = 0;
@@ -182,8 +182,7 @@ void neo_hookean(RXMeshStatic& rx, T dx)
     //                ground_n);
 
     // add neo hooken energy
-    neo_hookean_energy(
-        problem, x, is_dbc, volume, inv_b, mu_lame, time_step, lam);
+    neo_hookean_energy(problem, is_dbc, volume, inv_b, mu_lame, time_step, lam);
 
 
     int steps = 0;
@@ -273,9 +272,7 @@ void neo_hookean(RXMeshStatic& rx, T dx)
             line_search_init_step = std::min(nh_step, bar_step);
 
             // TODO: line search should pass the step to the friction energy
-            if (!newton_solver.line_search(line_search_init_step, 0.5, 5)) {
-                break;
-            }
+            newton_solver.line_search(line_search_init_step, 0.5);
 
             // evaluate energy
             add_contact(rx, problem.vv_pairs, v_dbc[0], is_dbc, x, dhat);
@@ -309,7 +306,7 @@ void neo_hookean(RXMeshStatic& rx, T dx)
             iter++;
         }
 
-        RXMESH_INFO("===================");
+        RXMESH_INFO("\n===================\n");
 
         //  update velocity
         rx.for_each_vertex(

apps/NeoHookean/neo_hookean_energy.h

@@ -5,13 +5,11 @@
 using namespace rxmesh;
 
 template <typename ProblemT,
-          typename VAttrT,
           typename VAttrI,
           typename VAttrF,
           typename VAttrFM,
           typename T>
 void neo_hookean_energy(ProblemT&      problem,
-                        const VAttrT&  x,
                         const VAttrI&  is_dbc,
                         const VAttrF&  volume,
                         const VAttrFM& inv_b,
@@ -34,9 +32,9 @@ void neo_hookean_energy(ProblemT&      problem,
                 return ActiveT();
             }
 
-            Eigen::Vector3<ActiveT> x0 = iter_val<ActiveT, 3>(fh, iter, x, 0);
-            Eigen::Vector3<ActiveT> x1 = iter_val<ActiveT, 3>(fh, iter, x, 1);
-            Eigen::Vector3<ActiveT> x2 = iter_val<ActiveT, 3>(fh, iter, x, 2);
+            Eigen::Vector3<ActiveT> x0 = iter_val<ActiveT, 3>(fh, iter, obj, 0);
+            Eigen::Vector3<ActiveT> x1 = iter_val<ActiveT, 3>(fh, iter, obj, 1);
+            Eigen::Vector3<ActiveT> x2 = iter_val<ActiveT, 3>(fh, iter, obj, 2);
 
             Eigen::Vector3<ActiveT> e0 = x2 - x0;
             Eigen::Vector3<ActiveT> e1 = x1 - x0;