avoiding duplicates pairs in hessian

Author: Ahdhn

apps/NeoHookean/barrier_energy.h

@@ -55,17 +55,19 @@ void floor_barrier_energy(ProblemT&      problem,
         });
 }
 
-template <typename VAttrT, typename VAttrB, typename T = typename VAttrT::Type>
+template <typename VAttrT,
+          typename VAttrB,
+          typename PairT,
+          typename T = typename VAttrT::Type>
 void add_contact(RXMeshStatic&      rx,
-                 CandidatePairsVV&  contact_pairs,
+                 PairT&             contact_pairs,
                  const VertexHandle dbc_vertex,
+                 const VertexHandle dbc_vertex1,
+                 const VertexHandle dbc_vertex2,
                  const VAttrB&      is_dbc,
                  const VAttrT&      x,
                  const T            dhat)
 {
-    if (contact_pairs.num_pairs() > 0) {
-        return;
-    }
     contact_pairs.reset();
 
     rx.for_each_vertex(DEVICE, [=] __device__(const VertexHandle& vh) mutable {
@@ -79,12 +81,22 @@ void add_contact(RXMeshStatic&      rx,
         T d = (xi - x_dbc).dot(normal);
 
         if (d < dhat) {
+            // if (vh.local_id() == 20 || vh.local_id() == 21 || vh.local_id()
+            // == 22 || vh.local_id() == 23 || vh.local_id() == 24) {
             bool inserted = contact_pairs.insert(vh, dbc_vertex);
             assert(inserted);
 
-            printf("\n addec contact pair between %d, %d",
+            inserted = contact_pairs.insert(vh, dbc_vertex1);
+            assert(inserted);
+
+            inserted = contact_pairs.insert(vh, dbc_vertex2);
+            assert(inserted);
+
+            printf("\n addec contact pair between %d, (%d, %d, %d)",
                    vh.local_id(),
-                   dbc_vertex.local_id());
+                   dbc_vertex.local_id(),
+                   dbc_vertex1.local_id(),
+                   dbc_vertex2.local_id());
         }
     });
 

apps/NeoHookean/draw.h

@@ -4,11 +4,13 @@
 
 using namespace rxmesh;
 
-template <typename FunT, typename VAttrT>
+template <typename FunT, typename VAttrT, typename DirT>
 void draw(RXMeshStatic& rx,
           VAttrT&       x,
           VAttrT&       velocity,
           FunT&         step_forward,
+          DirT&         dir_mat,
+          DirT&         grad_mat,
           int&          step)
 {
 #if USE_POLYSCOPE
@@ -18,13 +20,25 @@ void draw(RXMeshStatic& rx,
 
     bool is_running = false;
 
+    auto dir  = *rx.add_vertex_attribute<float>("Dir", 3);
+    auto grad = *rx.add_vertex_attribute<float>("Grad", 3);
+
     auto ps_callback = [&]() mutable {
         auto step_and_update = [&]() {
             step_forward();
             x.move(DEVICE, HOST);
             velocity.move(DEVICE, HOST);
             auto vel = rx.get_polyscope_mesh()->addVertexVectorQuantity(
                 "Velocity", velocity);
+
+            dir_mat.move(DEVICE, HOST);
+            dir.from_matrix(&dir_mat);
+            rx.get_polyscope_mesh()->addVertexVectorQuantity("Direction", dir);
+
+            grad_mat.move(DEVICE, HOST);
+            grad.from_matrix(&grad_mat);
+            rx.get_polyscope_mesh()->addVertexVectorQuantity("Grad", grad);
+
             rx.get_polyscope_mesh()->updateVertexPositions(x);
         };
         if (ImGui::Button("Step")) {

apps/NeoHookean/neo_hookean.cu

@@ -168,6 +168,11 @@ void neo_hookean(RXMeshStatic& rx, T dx)
     ceiling_barrier_energy(
         problem, contact_area, time_step, ground_n, ground_o, dhat, kappa);
 
+    DenseMatrix<T, Eigen::RowMajor> dir(
+        rx, problem.grad.rows(), problem.grad.cols());
+
+    DenseMatrix<T, Eigen::RowMajor> grad(
+        rx, problem.grad.rows(), problem.grad.cols());
 
     T line_search_init_step = 0;
 
@@ -193,6 +198,11 @@ void neo_hookean(RXMeshStatic& rx, T dx)
     timer.add("LinearSolver");
     timer.add("Diff");
 
+    //add_contact(
+    //    rx, problem.vv_pairs, v_dbc[0], v_dbc[1], v_dbc[2], is_dbc, x, dhat);
+    //problem.update_hessian();
+    //printf("\n");
+
     auto step_forward = [&]() {
         // x_tilde = x + v*h
         timer.start("Step");
@@ -221,14 +231,22 @@ void neo_hookean(RXMeshStatic& rx, T dx)
             rx, is_dbc, x, v_dbc_vel, v_dbc_limit, time_step, dbc_target);
 
         // evaluate energy
-        add_contact(rx, problem.vv_pairs, v_dbc[0], is_dbc, x, dhat);
+        add_contact(rx,
+                    problem.vv_pairs,
+                    v_dbc[0],
+                    v_dbc[1],
+                    v_dbc[2],
+                    is_dbc,
+                    x,
+                    dhat);
         problem.update_hessian();
         problem.eval_terms();
         {
             CUDA_ERROR(cudaGetLastError());
             CUDA_ERROR(cudaDeviceSynchronize());
         }
 
+        grad.copy_from(problem.grad, DEVICE, DEVICE);
         // DBC satisfied
         check_dbc_satisfied(
             rx, is_dbc_satisfied, x, is_dbc, dbc_target, time_step, tol);
@@ -244,18 +262,25 @@ void neo_hookean(RXMeshStatic& rx, T dx)
         // get newton direction
         newton_solver.compute_direction();
 
+        dir.copy_from(newton_solver.dir, DEVICE, DEVICE);
         // residual is abs_max(newton_dir)/ h
         T residual = newton_solver.dir.abs_max() / time_step;
 
         T f = problem.get_current_loss();
         RXMESH_INFO(
-            "*******Step: {}, Energy: {}, Residual: {}", steps, f, residual);
+            "*******Step: {}, Energy: {}, Residual: {},  DBC_satisfied= {}",
+            steps,
+            f,
+            residual,
+            num_satisfied);
 
         int iter = 0;
         while (residual > tol || num_satisfied != num_dbc_vertices) {
 
             if (residual <= tol && num_satisfied != num_dbc_vertices) {
                 dbc_stiff.multiply(T(2));
+                problem.eval_terms();
+                newton_solver.compute_direction();
             }
 
             T nh_step = neo_hookean_step_size(rx, x, newton_solver.dir, alpha);
@@ -272,10 +297,22 @@ 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
-            newton_solver.line_search(line_search_init_step, 0.5);
+            bool ls_success =
+                newton_solver.line_search(line_search_init_step, 0.5, 64, 0.0);
+
+            if (!ls_success) {
+                RXMESH_WARN("Line search failed!");
+            }
 
             // evaluate energy
-            add_contact(rx, problem.vv_pairs, v_dbc[0], is_dbc, x, dhat);
+            add_contact(rx,
+                        problem.vv_pairs,
+                        v_dbc[0],
+                        v_dbc[1],
+                        v_dbc[2],
+                        is_dbc,
+                        x,
+                        dhat);
             problem.update_hessian();
             problem.eval_terms();
             {
@@ -301,9 +338,20 @@ void neo_hookean(RXMeshStatic& rx, T dx)
             // residual is abs_max(newton_dir)/ h
             residual = newton_solver.dir.abs_max() / time_step;
 
-            RXMESH_INFO("  Subsetp: {}, F: {}, R: {}", iter, f, residual);
+            RXMESH_INFO(
+                "  Subsetp: {}, F: {}, R: {}, line_search_init_step={}, "
+                "DBC_satisfied= {}",
+                iter,
+                f,
+                residual,
+                line_search_init_step,
+                num_satisfied);
 
             iter++;
+
+            if (iter > 10) {
+                break;
+            }
         }
 
         RXMESH_INFO("\n===================\n");
@@ -325,7 +373,7 @@ void neo_hookean(RXMeshStatic& rx, T dx)
     };
 
 #if USE_POLYSCOPE
-    draw(rx, x, velocity, step_forward, steps);
+    draw(rx, x, velocity, step_forward, dir, grad, steps);
 #else
     while (steps < 5) {
         step_forward();

include/rxmesh/attribute.h

@@ -190,9 +190,10 @@ class Attribute : public AttributeBase
      * rows represent the number of mesh elements of this attribute and number
      * of columns is the number of attributes
      */
-    std::shared_ptr<DenseMatrix<T>> to_matrix() const
+    template <int Order = Eigen::ColMajor>
+    std::shared_ptr<DenseMatrix<T, Order>> to_matrix() const
     {
-        std::shared_ptr<DenseMatrix<T>> mat =
+        std::shared_ptr<DenseMatrix<T, Order>> mat =
             std::make_shared<DenseMatrix<T>>(*m_rxmesh, rows(), cols());
 
         if constexpr (std::is_same_v<HandleT, VertexHandle>) {
@@ -230,7 +231,8 @@ class Attribute : public AttributeBase
      * attribute on the host side
      * @param mat
      */
-    void from_matrix(DenseMatrix<T>* mat)
+    template <int Order>
+    void from_matrix(DenseMatrix<T, Order>* mat)
     {
         assert(mat->rows() == rows());
         assert(mat->cols() == cols());

include/rxmesh/diff/candidate_pairs.h

@@ -7,9 +7,11 @@ namespace rxmesh {
 /**
  * @brief Storing a pair of candidate contact pairs using their handles
  */
-template <typename HandleT0, typename HandleT1, typename IndexT>
+template <typename HandleT0, typename HandleT1, typename HessMatT>
 struct CandidatePairs
 {
+    using IndexT = int;
+
     template <typename T, int V, typename O, bool W>
     friend struct DiffScalarProblem;
 
@@ -36,11 +38,12 @@ struct CandidatePairs
      * contribute to -- but this is already allocated in the hessian matrix
      */
     __host__ CandidatePairs(int            max_capacity,
-                            int            variable_dim,
+                            HessMatT&      hess,
                             const Context& ctx)
-        : m_variable_dim(variable_dim),
+        : m_hess(hess),
+          m_variable_dim(hess.K_),
           m_pairs_id(DenseMatrix<IndexT, Eigen::ColMajor>(
-              max_capacity * variable_dim * variable_dim * 2,
+              max_capacity * m_variable_dim * m_variable_dim * 2,
               2)),
           m_pairs_handle(DenseMatrix<PairT, Eigen::ColMajor>(max_capacity, 1)),
           m_current_num_pairs(DenseMatrix<int>(1, 1)),
@@ -53,15 +56,23 @@ struct CandidatePairs
 
     /**
      * @brief Insert a candidate pair and return true if it succeeded. Otherwise
-     * return false, i.e., in case of exceeding the max capacity
+     * return false, i.e., in case of exceeding the max capacity.
+     * Here, we insert two things 1) the handles in a buffer to iterate over
+     * later 2) the corresponding new indices of the local Hessian block in the
+     * sparse Hessian matrix that these two (possibly) contacting pair
+     * generates.
+     * We always insert new handles (so long as the buffer is not overflown).
+     * For 2) we check if the block is already allocated. If not, we add the
+     * indices to another buffer which is later used in the update_hessian in
+     * the Problem
      */
     __device__ __host__ bool insert(const HandleT0& c0, const HandleT1& c1)
     {
 
         // first check if we have a space
         /// if we have a space to insert new pairs
 
-        auto add_candidate = [&](int id) {
+        auto add_candidate_with_indices = [&](int id) {
             // add the handles
             m_pairs_handle(id).first  = c0;
             m_pairs_handle(id).second = c1;
@@ -86,12 +97,24 @@ struct CandidatePairs
             }
         };
 
+
+        auto add_candidate = [&](int id) {
+            // add the handles
+            m_pairs_handle(id).first  = c0;
+            m_pairs_handle(id).second = c1;
+        };
+
+
 #ifdef __CUDA_ARCH__
         int id = ::atomicAdd(m_current_num_pairs.data(DEVICE), 1);
-        if (id < m_pairs_handle.rows()) {
-            ::atomicAdd(m_current_num_index.data(DEVICE),
-                        m_variable_dim * m_variable_dim * 2);
-            add_candidate(id);
+        if (id < m_pairs_handle.rows()) {            
+            if (m_hess.is_non_zero(c0, c1)) {
+                add_candidate(id);
+            } else {
+                ::atomicAdd(m_current_num_index.data(DEVICE),
+                            m_variable_dim * m_variable_dim * 2);
+                add_candidate_with_indices(id);
+            }
             return true;
         } else {
             ::atomicAdd(m_current_num_pairs.data(DEVICE), -1);
@@ -100,9 +123,13 @@ struct CandidatePairs
 #else
         if (m_current_num_pairs(0) < m_pairs_handle.rows()) {
             int id = m_current_num_pairs(0);
-            m_current_num_pairs(0)++;
-            m_current_num_index(0) += m_variable_dim * m_variable_dim * 2;
-            add_candidate(id);
+            m_current_num_pairs(0)++;            
+            if (m_hess.is_non_zero(c0, c1)) {
+                add_candidate(id);
+            } else {
+                m_current_num_index(0) += m_variable_dim * m_variable_dim * 2;
+                add_candidate_with_indices(id);
+            }
             return true;
         } else {
             return false;
@@ -168,7 +195,6 @@ struct CandidatePairs
         m_current_num_index(0) = 0;
 
 #ifndef __CUDA_ARCH__
-
         m_current_num_pairs.move(HOST, DEVICE);
         m_current_num_index.move(HOST, DEVICE);
 #endif
@@ -194,12 +220,16 @@ struct CandidatePairs
     // track the number of pairs (not the number of indices)
     DenseMatrix<int> m_current_num_pairs;
     DenseMatrix<int> m_current_num_index;
-    int              m_variable_dim;
+
+    // TODO remove this and use m_hess.k_
+    int m_variable_dim;
+
+    HessMatT m_hess;
 
     Context m_context;
 };
 
-
-using CandidatePairsVV = CandidatePairs<VertexHandle, VertexHandle, int>;
+template <typename HessMatT>
+using CandidatePairsVV = CandidatePairs<VertexHandle, VertexHandle, HessMatT>;
 
 }  // namespace rxmesh