Refactor ContactData::update()

src/smith/physics/contact/contact_data.cpp

@@ -74,20 +74,60 @@ void ContactData::reset()
   }
 }
 
-void ContactData::update(int cycle, double time, double& dt)
+void ContactData::updateGaps(int cycle, double time, double& dt,
+                             std::optional<std::reference_wrapper<const mfem::Vector>> u_shape,
+                             std::optional<std::reference_wrapper<const mfem::Vector>> u, bool eval_jacobian)
 {
   cycle_ = cycle;
   time_ = time;
   dt_ = dt;
+
+  if (u_shape && u) {
+    setDisplacements(u_shape->get(), u->get());
+  }
+
+  for (auto& interaction : interactions_) {
+    interaction.evalJacobian(eval_jacobian);
+  }
   // This updates the redecomposed surface mesh based on the current displacement, then transfers field quantities to
   // the updated mesh.
-  tribol::updateMfemParallelDecomposition();
-  // This function computes forces, gaps, and Jacobian contributions based on the current field quantities. Note the
-  // fields (with the exception of pressure) are stored on the redecomposed surface mesh until transferred by calling
-  // forces(), mergedGaps(), etc.
+  if (u_shape && u) {
+    tribol::updateMfemParallelDecomposition();
+  }
+  // This function computes gaps (and optionally geometric Jacobian blocks) based on the current mesh.
   tribol::update(cycle, time, dt);
 }
 
+void ContactData::update(int cycle, double time, double& dt,
+                         std::optional<std::reference_wrapper<const mfem::Vector>> u_shape,
+                         std::optional<std::reference_wrapper<const mfem::Vector>> u,
+                         std::optional<std::reference_wrapper<const mfem::Vector>> p)
+{
+  // First pass: update gaps if coordinates are provided
+  if (u_shape && u) {
+    updateGaps(cycle, time, dt, u_shape, u, false);
+  } else {
+    // Ensure internal timing is updated even if coordinates are not
+    cycle_ = cycle;
+    time_ = time;
+    dt_ = dt;
+  }
+
+  // second pass: update pressures and compute forces/Jacobians if p is provided
+  if (p) {
+    // with updated gaps, we can update pressure for contact interactions (active set detection and penalty)
+    setPressures(p->get());
+
+    for (auto& interaction : interactions_) {
+      interaction.evalJacobian(true);
+    }
+    // This second call is required to synchronize the updated pressures to Tribol's internal redecomposed surface mesh
+    // and to ensure Tribol's internal state is correctly reset for the second pass.
+    tribol::updateMfemParallelDecomposition();
+    tribol::update(cycle, time, dt);
+  }
+}
+
 FiniteElementDual ContactData::forces() const
 {
   FiniteElementDual f(*reference_nodes_->ParFESpace(), "contact force");
@@ -278,20 +318,7 @@ void ContactData::residualFunction(const mfem::Vector& u_shape, const mfem::Vect
   mfem::Vector r_blk(r, 0, disp_size);
   mfem::Vector g_blk(r, disp_size, numPressureDofs());
 
-  setDisplacements(u_shape, u_blk);
-
-  // we need to call update first to update gaps
-  for (auto& interaction : interactions_) {
-    interaction.evalJacobian(false);
-  }
-  update(cycle_, time_, dt_);
-  // with updated gaps, we can update pressure for contact interactions with penalty enforcement
-  setPressures(p_blk);
-  // call update again with the right pressures
-  for (auto& interaction : interactions_) {
-    interaction.evalJacobian(true);
-  }
-  update(cycle_, time_, dt_);
+  update(cycle_, time_, dt_, u_shape, u_blk, p_blk);
 
   r_blk += forces();
   // calling mergedGaps() with true will zero out gap on inactive dofs (so the residual converges and the linearized
@@ -456,7 +483,19 @@ void ContactData::addContactInteraction([[maybe_unused]] int interaction_id,
   SLIC_WARNING_ROOT("Smith built without Tribol support. No contact interaction will be added.");
 }
 
-void ContactData::update([[maybe_unused]] int cycle, [[maybe_unused]] double time, [[maybe_unused]] double& dt) {}
+void ContactData::updateGaps([[maybe_unused]] int cycle, [[maybe_unused]] double time, [[maybe_unused]] double& dt,
+                             [[maybe_unused]] std::optional<std::reference_wrapper<const mfem::Vector>> u_shape,
+                             [[maybe_unused]] std::optional<std::reference_wrapper<const mfem::Vector>> u,
+                             [[maybe_unused]] bool eval_jacobian)
+{
+}
+
+void ContactData::update([[maybe_unused]] int cycle, [[maybe_unused]] double time, [[maybe_unused]] double& dt,
+                         [[maybe_unused]] std::optional<std::reference_wrapper<const mfem::Vector>> u_shape,
+                         [[maybe_unused]] std::optional<std::reference_wrapper<const mfem::Vector>> u,
+                         [[maybe_unused]] std::optional<std::reference_wrapper<const mfem::Vector>> p)
+{
+}
 
 FiniteElementDual ContactData::forces() const
 {

src/smith/physics/contact/contact_data.hpp

@@ -15,6 +15,8 @@
 #include <memory>
 #include <set>
 #include <vector>
+#include <optional>
+#include <functional>
 
 #include "mfem.hpp"
 
@@ -69,14 +71,35 @@ class ContactData {
   void addContactInteraction(int interaction_id, const std::set<int>& bdry_attr_surf1,
                              const std::set<int>& bdry_attr_surf2, ContactOptions contact_opts);
 
+  /**
+   * @brief Updates only the gap contributions associated with contact
+   *
+   * @param cycle The current simulation cycle
+   * @param time The current time
+   * @param dt The timestep size to attempt
+   * @param u_shape Optional shape displacement vector
+   * @param u Optional current displacement dof values
+   * @param eval_jacobian Whether to also evaluate the Jacobian contributions (default false)
+   */
+  void updateGaps(int cycle, double time, double& dt,
+                  std::optional<std::reference_wrapper<const mfem::Vector>> u_shape = std::nullopt,
+                  std::optional<std::reference_wrapper<const mfem::Vector>> u = std::nullopt,
+                  bool eval_jacobian = false);
+
   /**
    * @brief Updates the positions, forces, and Jacobian contributions associated with contact
    *
    * @param cycle The current simulation cycle
    * @param time The current time
    * @param dt The timestep size to attempt
+   * @param u_shape Optional shape displacement vector
+   * @param u Optional current displacement dof values
+   * @param p Optional current pressure true dof values
    */
-  void update(int cycle, double time, double& dt);
+  void update(int cycle, double time, double& dt,
+              std::optional<std::reference_wrapper<const mfem::Vector>> u_shape = std::nullopt,
+              std::optional<std::reference_wrapper<const mfem::Vector>> u = std::nullopt,
+              std::optional<std::reference_wrapper<const mfem::Vector>> p = std::nullopt);
 
   /**
    * @brief Resets the contact pressures to zero
@@ -163,27 +186,6 @@ class ContactData {
    */
   std::unique_ptr<mfem::HypreParMatrix> contactSubspaceTransferOperator();
 
-  /**
-   * @brief Set the pressure field
-   *
-   * This sets Tribol's pressure degrees of freedom based on
-   *  1) the values in merged_pressure for Lagrange multiplier enforcement
-   *  2) the nodal gaps and penalty for penalty enforcement
-   *
-   * @note The nodal gaps must be up-to-date for penalty enforcement
-   *
-   * @param merged_pressures Current pressure true dof values in a merged mfem::Vector
-   */
-  void setPressures(const mfem::Vector& merged_pressures) const;
-
-  /**
-   * @brief Update the current coordinates based on the new displacement field
-   *
-   * @param u_shape Shape displacement vector
-   * @param u Current displacement dof values
-   */
-  void setDisplacements(const mfem::Vector& u_shape, const mfem::Vector& u);
-
   /**
    * @brief Have there been contact interactions added?
    *
@@ -223,6 +225,28 @@ class ContactData {
    */
   int numPressureDofs() const { return num_pressure_dofs_; };
 
+ protected:
+  /**
+   * @brief Set the pressure field
+   *
+   * This sets Tribol's pressure degrees of freedom based on
+   *  1) the values in merged_pressure for Lagrange multiplier enforcement
+   *  2) the nodal gaps and penalty for penalty enforcement
+   *
+   * @note The nodal gaps must be up-to-date for penalty enforcement
+   *
+   * @param merged_pressures Current pressure true dof values in a merged mfem::Vector
+   */
+  void setPressures(const mfem::Vector& merged_pressures) const;
+
+  /**
+   * @brief Update the current coordinates based on the new displacement field
+   *
+   * @param u_shape Shape displacement vector
+   * @param u Current displacement dof values
+   */
+  void setDisplacements(const mfem::Vector& u_shape, const mfem::Vector& u);
+
  private:
 #ifdef SMITH_USE_TRIBOL
   /**

src/smith/physics/contact_constraint.hpp

@@ -111,10 +111,9 @@ class ContactConstraint : public Constraint {
                         bool update_fields = true) const override
   {
     if (update_fields) {
-      contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]);
       // note: Tribol does not use cycle.
       int cycle = 0;
-      contact_.update(cycle, time, dt);
+      contact_.updateGaps(cycle, time, dt, *fields[ContactFields::SHAPE], *fields[ContactFields::DISP]);
       auto gaps_hpv = contact_.mergedGaps(false);
       // Note: this copy is needed to prevent the HypreParVector pointer from going out of scope.  see
       // https://github.com/mfem/mfem/issues/5029
@@ -140,19 +139,16 @@ class ContactConstraint : public Constraint {
                                                  [[maybe_unused]] bool fresh_derivative = true) const override
   {
     SLIC_ERROR_IF(direction != ContactFields::DISP, "requesting a non displacement-field derivative");
-    // if the fields are to be updated then we update the displacement field
-    if (update_fields) {
-      contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]);
-    }
     // if the field has been updated or we are requesting a fresh derivative
     // then re-compute the gap Jacobian
     // otherwise use previously cached Jacobian
     if (update_fields || fresh_derivative) {
-      for (auto& interaction : contact_.getContactInteractions()) {
-        interaction.evalJacobian(true);
-      }
       int cycle = 0;
-      contact_.update(cycle, time, dt);
+      if (update_fields) {
+        contact_.updateGaps(cycle, time, dt, *fields[ContactFields::SHAPE], *fields[ContactFields::DISP], true);
+      } else {
+        contact_.updateGaps(cycle, time, dt, std::nullopt, std::nullopt, true);
+      }
       J_contact_ = contact_.mergedJacobian();
     }
     // obtain (1, 0) block entry from the 2 x 2 block contact linear system
@@ -179,22 +175,12 @@ class ContactConstraint : public Constraint {
   {
     SLIC_ERROR_IF(direction != ContactFields::DISP, "requesting a non displacement-field derivative");
     int cycle = 0;
-    if (update_fields) {
-      contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]);
-      // first update gaps
-      for (auto& interaction : contact_.getContactInteractions()) {
-        interaction.evalJacobian(false);
-      }
-      contact_.update(cycle, time, dt);
-    }
     if (update_fields || fresh_derivative) {
-      // with updated gaps, then update pressure for contact interactions with penalty enforcement
-      contact_.setPressures(multipliers);
-      // call update again with the right pressures
-      for (auto& interaction : contact_.getContactInteractions()) {
-        interaction.evalJacobian(true);
+      if (update_fields) {
+        contact_.update(cycle, time, dt, *fields[ContactFields::SHAPE], *fields[ContactFields::DISP], multipliers);
+      } else {
+        contact_.update(cycle, time, dt, std::nullopt, std::nullopt, multipliers);
       }
-      contact_.update(cycle, time, dt);
     }
     return contact_.forces();
   };
@@ -220,22 +206,12 @@ class ContactConstraint : public Constraint {
     SLIC_ERROR_IF(direction != ContactFields::DISP, "requesting a non displacement-field derivative");
 
     int cycle = 0;
-    if (update_fields) {
-      contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]);
-      // first update gaps
-      for (auto& interaction : contact_.getContactInteractions()) {
-        interaction.evalJacobian(false);
-      }
-      contact_.update(cycle, time, dt);
-    }
     if (update_fields || fresh_derivative) {
-      // with updated gaps, we can update pressure for contact interactions with penalty enforcement
-      contact_.setPressures(multipliers);
-      // call update again with the right pressures
-      for (auto& interaction : contact_.getContactInteractions()) {
-        interaction.evalJacobian(true);
+      if (update_fields) {
+        contact_.update(cycle, time, dt, *fields[ContactFields::SHAPE], *fields[ContactFields::DISP], multipliers);
+      } else {
+        contact_.update(cycle, time, dt, std::nullopt, std::nullopt, multipliers);
       }
-      contact_.update(cycle, time, dt);
       J_contact_ = contact_.mergedJacobian();
     }
     // obtain (0, 0) block entry from the 2 x 2 block contact linear system
@@ -260,14 +236,13 @@ class ContactConstraint : public Constraint {
     SLIC_ERROR_IF(direction != ContactFields::DISP, "requesting a non displacement-field derivative");
 
     int cycle = 0;
-    if (update_fields) {
-      contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]);
-    }
     if (update_fields || fresh_derivative) {
-      for (auto& interaction : contact_.getContactInteractions()) {
-        interaction.evalJacobian(true);
+      mfem::Vector p = contact_.mergedPressures();
+      if (update_fields) {
+        contact_.update(cycle, time, dt, *fields[ContactFields::SHAPE], *fields[ContactFields::DISP], p);
+      } else {
+        contact_.update(cycle, time, dt, std::nullopt, std::nullopt, p);
       }
-      contact_.update(cycle, time, dt);
       J_contact_ = contact_.mergedJacobian();
     }
     // obtain (0, 1) block entry from the 2 x 2 block contact linear system

src/smith/physics/solid_mechanics_contact.hpp

@@ -115,10 +115,11 @@ class SolidMechanicsContact<order, dim, Parameters<parameter_space...>,
   {
     SolidMechanicsBase::resetStates(cycle, time);
     forces_ = 0.0;
-    contact_.setDisplacements(BasePhysics::shapeDisplacement(), displacement_);
     contact_.reset();
     double dt = 0.0;
-    contact_.update(cycle, time, dt);
+    mfem::Vector p(contact_.numPressureDofs());
+    p = 0.0;
+    contact_.update(cycle, time, dt, BasePhysics::shapeDisplacement(), displacement_, p);
   }
 
   /// @brief Build the quasi-static operator corresponding to the total Lagrangian formulation
@@ -236,7 +237,8 @@ class SolidMechanicsContact<order, dim, Parameters<parameter_space...>,
   void completeSetup() override
   {
     double dt = 0.0;
-    contact_.update(cycle_, time_, dt);
+    mfem::Vector p = pressure();
+    contact_.update(cycle_, time_, dt, BasePhysics::shapeDisplacement(), displacement_, p);
 
     SolidMechanicsBase::completeSetup();
   }
@@ -269,8 +271,6 @@ class SolidMechanicsContact<order, dim, Parameters<parameter_space...>,
     // solve the non-linear system resid = 0 and pressure * gap = 0
     nonlin_solver_->solve(augmented_solution);
     displacement_.Set(1.0, mfem::Vector(augmented_solution, 0, displacement_.Size()));
-    contact_.setPressures(mfem::Vector(augmented_solution, displacement_.Size(), contact_.numPressureDofs()));
-    contact_.update(cycle_, time_, dt);
     forces_.SetVector(contact_.forces(), 0);
   }
 
@@ -349,24 +349,29 @@ class SolidMechanicsContact<order, dim, Parameters<parameter_space...>,
       const mfem::Vector res = (*residual_)(time_ + dt, BasePhysics::shapeDisplacement(), displacement_, acceleration_,
                                             *parameters_[parameter_indices].state...);
 
-      contact_.setPressures(mfem::Vector(augmented_residual, displacement_.Size(), contact_.numPressureDofs()));
-      contact_.update(cycle_, time_, dt);
-      mfem::Vector r_blk(augmented_residual, 0, displacement_.space().TrueVSize());
-      r_blk = res;
-
       mfem::Vector augmented_solution(displacement_.space().TrueVSize() + contact_.numPressureDofs());
       augmented_solution = 0.0;
       mfem::Vector du(augmented_solution, 0, displacement_.space().TrueVSize());
       du = displacement_;
+      mfem::Vector p_blk(augmented_solution, displacement_.Size(), contact_.numPressureDofs());
+
+      // Perform a single update for the warm start evaluation.
+      // Note: we use time_ to match the previous Jacobian evaluation point.
+      contact_.update(cycle_, time_, dt, BasePhysics::shapeDisplacement(), displacement_, p_blk);
+
+      mfem::Vector r_blk(augmented_residual, 0, displacement_.space().TrueVSize());
+      r_blk = res;
+      r_blk += contact_.forces();
+
+      mfem::Vector g_blk(augmented_residual, displacement_.Size(), contact_.numPressureDofs());
+      g_blk.Set(1.0, contact_.mergedGaps(true));
 
-      contact_.residualFunction(BasePhysics::shapeDisplacement(), augmented_solution, augmented_residual);
       r_blk.SetSubVector(bcs_.allEssentialTrueDofs(), 0.0);
 
       // use the most recently evaluated Jacobian
       auto [_, drdu] = (*residual_)(time_, BasePhysics::shapeDisplacement(), differentiate_wrt(displacement_),
                                     acceleration_, *parameters_[parameter_indices].previous_state...);
 
-      contact_.update(cycle_, time_, dt);
       if (contact_.haveLagrangeMultipliers()) {
         J_offsets_ = mfem::Array<int>({0, displacement_.Size(), displacement_.Size() + contact_.numPressureDofs()});
         J_constraint_ = contact_.jacobianFunction(assemble(drdu));