Per-contact interaction adjoint support

src/smith/physics/contact/contact_data.cpp

@@ -155,9 +155,10 @@ std::unique_ptr<mfem::BlockOperator> ContactData::mergedJacobian() const
   }
 
   for (size_t i{0}; i < interactions_.size(); ++i) {
-    // this is the BlockOperator for one of the contact interactions
-    auto interaction_J = interactions_[i].jacobian();
+    // this is the BlockOperator for one of the contact interactions, post-processed for Smith's assembly conventions
+    auto interaction_J = interactions_[i].jacobianContribution();
     interaction_J->owns_blocks = false;  // we'll manage the ownership of the blocks on our own...
+
     // add the contact interaction's contribution to df_(contact)/dx (the 0, 0 block)
     if (!interaction_J->IsZeroBlock(0, 0)) {
       SLIC_ERROR_ROOT_IF(!dynamic_cast<mfem::HypreParMatrix*>(&interaction_J->GetBlock(0, 0)),
@@ -171,42 +172,16 @@ std::unique_ptr<mfem::BlockOperator> ContactData::mergedJacobian() const
         delete &interaction_J->GetBlock(0, 0);
       }
     }
-    // add the contact interaction's (other) contribution to df_(contact)/dx (for penalty) or to df_(contact)/dp and
-    // dg/dx (for Lagrange multipliers)
+
+    // add the contact interaction's contribution to df_(contact)/dp and dg/dx (for Lagrange multipliers)
     if (!interaction_J->IsZeroBlock(1, 0) && !interaction_J->IsZeroBlock(0, 1)) {
       auto dgdu = dynamic_cast<mfem::HypreParMatrix*>(&interaction_J->GetBlock(1, 0));
       auto dfdp = dynamic_cast<mfem::HypreParMatrix*>(&interaction_J->GetBlock(0, 1));
       SLIC_ERROR_ROOT_IF(!dgdu, "Only HypreParMatrix constraint matrix blocks are currently supported.");
       SLIC_ERROR_ROOT_IF(!dfdp, "Only HypreParMatrix constraint matrix blocks are currently supported.");
-      // zero out rows and cols not in the active set
-      auto inactive_dofs = interactions_[i].inactiveDofs();
-      dgdu->EliminateRows(inactive_dofs);
-      auto dfdp_elim = std::unique_ptr<mfem::HypreParMatrix>(dfdp->EliminateCols(inactive_dofs));
-      if (interactions_[i].getContactOptions().enforcement == ContactEnforcement::Penalty) {
-        // compute contribution to df_(contact)/dx (the 0, 0 block) for penalty
-        std::unique_ptr<mfem::HypreParMatrix> BTB(mfem::ParMult(dfdp, dgdu, true));
-        delete &interaction_J->GetBlock(1, 0);
-        delete &interaction_J->GetBlock(0, 1);
-        if (block_J->IsZeroBlock(0, 0)) {
-          mfem::Vector penalty(reference_nodes_->ParFESpace()->GetTrueVSize());
-          penalty = interactions_[i].getContactOptions().penalty;
-          BTB->ScaleRows(penalty);
-          block_J->SetBlock(0, 0, BTB.release());
-        } else {
-          block_J->SetBlock(0, 0,
-                            mfem::Add(1.0, static_cast<mfem::HypreParMatrix&>(block_J->GetBlock(0, 0)),
-                                      interactions_[i].getContactOptions().penalty, *BTB));
-        }
-      } else  // enforcement == ContactEnforcement::LagrangeMultiplier
-      {
-        // compute contribution to off-diagonal blocks for Lagrange multiplier
-        dgdu_blocks(static_cast<int>(i), 0) = dgdu;
-        dfdp_blocks(0, static_cast<int>(i)) = dfdp;
-      }
-      if (!interaction_J->IsZeroBlock(1, 1)) {
-        // we track our own active set, so get rid of the tribol inactive dof block
-        delete &interaction_J->GetBlock(1, 1);
-      }
+
+      dgdu_blocks(static_cast<int>(i), 0) = dgdu;
+      dfdp_blocks(0, static_cast<int>(i)) = dfdp;
     }
   }
   if (haveLagrangeMultipliers()) {

src/smith/physics/contact/contact_interaction.cpp

@@ -107,6 +107,68 @@ std::unique_ptr<mfem::BlockOperator> ContactInteraction::jacobian() const
   return tribol::getMfemBlockJacobian(getInteractionId());
 }
 
+std::unique_ptr<mfem::BlockOperator> ContactInteraction::jacobianContribution() const
+{
+  const int disp_size = current_coords_.ParFESpace()->GetTrueVSize();
+  const int pressure_size = numPressureDofs();
+
+  auto offsets = mfem::Array<int>({0, disp_size, disp_size + pressure_size});
+  auto out = std::make_unique<mfem::BlockOperator>(offsets);
+  out->owns_blocks = true;
+
+  auto interaction_J = jacobian();
+  interaction_J->owns_blocks = false;  // manage block ownership explicitly
+
+  std::unique_ptr<mfem::HypreParMatrix> dfdx;
+  if (!interaction_J->IsZeroBlock(0, 0)) {
+    auto* block_0_0 = dynamic_cast<mfem::HypreParMatrix*>(&interaction_J->GetBlock(0, 0));
+    SLIC_ERROR_ROOT_IF(!block_0_0, "Only HypreParMatrix constraint matrix blocks are currently supported.");
+    dfdx.reset(block_0_0);
+  }
+
+  if (!interaction_J->IsZeroBlock(1, 0) && !interaction_J->IsZeroBlock(0, 1)) {
+    auto* dgdu = dynamic_cast<mfem::HypreParMatrix*>(&interaction_J->GetBlock(1, 0));
+    auto* dfdp = dynamic_cast<mfem::HypreParMatrix*>(&interaction_J->GetBlock(0, 1));
+    SLIC_ERROR_ROOT_IF(!dgdu, "Only HypreParMatrix constraint matrix blocks are currently supported.");
+    SLIC_ERROR_ROOT_IF(!dfdp, "Only HypreParMatrix constraint matrix blocks are currently supported.");
+
+    // zero out rows and cols not in the active set
+    auto inactive_dofs = inactiveDofs();
+    dgdu->EliminateRows(inactive_dofs);
+    auto dfdp_elim = std::unique_ptr<mfem::HypreParMatrix>(dfdp->EliminateCols(inactive_dofs));
+
+    if (getContactOptions().enforcement == ContactEnforcement::Penalty) {
+      std::unique_ptr<mfem::HypreParMatrix> BTB(mfem::ParMult(dfdp, dgdu, true));
+      delete &interaction_J->GetBlock(1, 0);
+      delete &interaction_J->GetBlock(0, 1);
+
+      if (!dfdx) {
+        mfem::Vector penalty(disp_size);
+        penalty = getContactOptions().penalty;
+        BTB->ScaleRows(penalty);
+        dfdx = std::move(BTB);
+      } else {
+        dfdx.reset(mfem::Add(1.0, *dfdx, getContactOptions().penalty, *BTB));
+      }
+    } else  // ContactEnforcement::LagrangeMultiplier
+    {
+      out->SetBlock(1, 0, dgdu);
+      out->SetBlock(0, 1, dfdp);
+    }
+  }
+
+  if (dfdx) {
+    out->SetBlock(0, 0, dfdx.release());
+  }
+
+  // Smith builds the merged (1,1) inactive-dof identity itself; discard any Tribol-provided (1,1) block.
+  if (!interaction_J->IsZeroBlock(1, 1)) {
+    delete &interaction_J->GetBlock(1, 1);
+  }
+
+  return out;
+}
+
 int ContactInteraction::numPressureDofs() const
 {
   return getContactOptions().enforcement == ContactEnforcement::LagrangeMultiplier

src/smith/physics/contact/contact_interaction.hpp

@@ -109,6 +109,20 @@ class ContactInteraction {
    */
   std::unique_ptr<mfem::BlockOperator> jacobian() const;
 
+  /**
+   * @brief Get the Smith-processed Jacobian contribution for this interaction
+   *
+   * This method post-processes the Tribol Jacobian for Smith's assembly conventions:
+   * - For penalty enforcement, it folds the (df/dp, dg/dx) blocks into df/dx via:
+   *     df/dx += penalty * (df/dp)^T * dg/dx
+   *   and returns a BlockOperator with only the (0,0) block populated (pressure block has size 0).
+   * - For Lagrange multiplier enforcement, it returns (0,0), (1,0), and (0,1) blocks (with inactive dofs eliminated),
+   *   leaving construction of the merged (1,1) inactive-dof identity to ContactData.
+   *
+   * @return A new BlockOperator whose blocks are owned by the returned object.
+   */
+  std::unique_ptr<mfem::BlockOperator> jacobianContribution() const;
+
   /**
    * @brief Get the finite element space of the pressure DOFs
    *