Refactors so that we can call contact with resets in between and get …

…the same answer.

src/serac/physics/contact/contact_data.cpp

@@ -41,6 +41,15 @@ void ContactData::addContactInteraction(int interaction_id, const std::set<int>&
   }
 }
 
+void ContactData::reset()
+{
+  for (auto& interaction : interactions_) {
+    FiniteElementState zero = interaction.pressure();
+    zero = 0.0;
+    interaction.setPressure(zero);
+  }
+}
+
 void ContactData::update(int cycle, double time, double& dt)
 {
   cycle_ = cycle;
@@ -240,6 +249,7 @@ void ContactData::residualFunction(const mfem::Vector& u_shape, const mfem::Vect
   mfem::Vector g_blk(r, disp_size, numPressureDofs());
 
   setDisplacements(u_shape, u_blk);
+
   // we need to call update first to update gaps
   update(cycle_, time_, dt_);
   // with updated gaps, we can update pressure for contact interactions with penalty enforcement

src/serac/physics/contact/contact_data.hpp

@@ -74,6 +74,12 @@ class ContactData {
    */
   void update(int cycle, double time, double& dt);
 
+  /**
+   * @brief Updates the positions, forces, and Jacobian contributions associated with contact
+   *
+   */
+  void reset();
+
   /**
    * @brief Get the contact constraint residual (i.e. nodal forces) from all contact interactions
    *

src/serac/physics/solid_mechanics_contact.hpp

@@ -110,6 +110,17 @@ class SolidMechanicsContact<order, dim, Parameters<parameter_space...>,
     duals_.push_back(&forces_);
   }
 
+  /// @overload
+  void resetStates(int cycle = 0, double time = 0.0) override
+  {
+    SolidMechanicsBase::resetStates(cycle, time);
+    forces_ = 0.0;
+    contact_.setDisplacements(shape_displacement_, displacement_);
+    contact_.reset();
+    double dt = 0.0;
+    contact_.update(cycle, time, dt);
+  }
+
   /// @brief Build the quasi-static operator corresponding to the total Lagrangian formulation
   std::unique_ptr<mfem_ext::StdFunctionOperator> buildQuasistaticOperator() override
   {
@@ -121,7 +132,7 @@ class SolidMechanicsContact<order, dim, Parameters<parameter_space...>,
       // NOTE this copy is required as the sundials solvers do not allow move assignments because of their memory
       // tracking strategy
       // See https://github.com/mfem/mfem/issues/3531
-      mfem::Vector r_blk(r, 0, displacement_.Size());
+      mfem::Vector r_blk(r, 0, displacement_.space().TrueVSize());
       r_blk = res;
 
       contact_.residualFunction(shape_displacement_, u, r);
@@ -178,6 +189,7 @@ class SolidMechanicsContact<order, dim, Parameters<parameter_space...>,
                                           *parameters_[parameter_indices].state...);
             J_             = assemble(drdu);
 
+            //contact_.setDisplacements(u_shape, u_blk);
             // get 11-block holding jacobian contributions
             auto block_J         = contact_.jacobianFunction(J_.release());
             block_J->owns_blocks = false;
@@ -296,13 +308,30 @@ class SolidMechanicsContact<order, dim, Parameters<parameter_space...>,
     }
 
     if (use_warm_start_) {
+
+      // Update the system residual
+      mfem::Vector augmented_residual(displacement_.space().TrueVSize() + contact_.numPressureDofs());
+      augmented_residual     = 0.0;
+      const mfem::Vector res = (*residual_)(time_ + dt, shape_displacement_, displacement_, acceleration_,
+                                            *parameters_[parameter_indices].state...);
+
+      contact_.setPressures(mfem::Vector(augmented_residual, displacement_.Size(), contact_.numPressureDofs()));
+      contact_.update(cycle_, time_, dt);
+      // TODO this copy is required as the sundials solvers do not allow move assignments because of their memory
+      // tracking strategy
+      // See https://github.com/mfem/mfem/issues/3531
+      mfem::Vector r_blk(augmented_residual, 0, displacement_.space().TrueVSize());
+      r_blk = res;
+
+      contact_.residualFunction(shape_displacement_, displacement_, augmented_residual);
+      r_blk.SetSubVector(bcs_.allEssentialTrueDofs(), 0.0);
+
       // use the most recently evaluated Jacobian
       auto [_, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_,
                                     *parameters_[parameter_indices].previous_state...);
       J_.reset();
       J_ = assemble(drdu);
 
-      contact_.update(cycle_, time_, dt);
       if (contact_.haveLagrangeMultipliers()) {
         J_offsets_    = mfem::Array<int>({0, displacement_.Size(), displacement_.Size() + contact_.numPressureDofs()});
         J_constraint_ = contact_.jacobianFunction(J_.release());
@@ -334,29 +363,16 @@ class SolidMechanicsContact<order, dim, Parameters<parameter_space...>,
         J_operator_ = J_.get();
       }
 
-      // Update the linearized Jacobian matrix
-      mfem::Vector augmented_residual(displacement_.space().TrueVSize() + contact_.numPressureDofs());
-      augmented_residual     = 0.0;
-      const mfem::Vector res = (*residual_)(time_ + dt, shape_displacement_, displacement_, acceleration_,
-                                            *parameters_[parameter_indices].state...);
-
-      // TODO this copy is required as the sundials solvers do not allow move assignments because of their memory
-      // tracking strategy
-      // See https://github.com/mfem/mfem/issues/3531
-      mfem::Vector r(augmented_residual, 0, displacement_.space().TrueVSize());
-      r = res;
-      r += contact_.forces();
-
       augmented_residual *= -1.0;
 
       mfem::Vector augmented_solution(displacement_.space().TrueVSize() + contact_.numPressureDofs());
       augmented_solution = 0.0;
       mfem::Vector du(augmented_solution, 0, displacement_.space().TrueVSize());
       du = du_;
-      mfem::EliminateBC(*J_, *J_e_, constrained_dofs, du, r);
+      mfem::EliminateBC(*J_, *J_e_, constrained_dofs, du, r_blk);
       for (int i = 0; i < constrained_dofs.Size(); i++) {
         int j = constrained_dofs[i];
-        r[j]  = du[j];
+        r_blk[j]  = du[j];
       }
 
       auto& lin_solver = nonlin_solver_->linearSolver();

src/serac/physics/tests/contact_solid_adjoint.cpp

@@ -63,7 +63,7 @@ std::unique_ptr<SolidMechT> createContactSolver(const NonlinearSolverOptions& no
 
   auto solid =
       std::make_unique<SolidMechT>(nonlinear_opts, solid_mechanics::direct_linear_options, dyn_opts,
-                                   physics_prefix + std::to_string(iter++), mesh_tag);
+                                   physics_prefix + std::to_string(iter++), mesh_tag, std::vector<std::string>{}, 0, 0.0, false, true);
   solid->setMaterial(mat);
 
   solid->setDisplacementBCs({2}, [](const mfem::Vector& /*X*/, double /*t*/, mfem::Vector& disp) { disp = 0.0; });
@@ -168,15 +168,14 @@ TEST_F(ContactSensitivityFixture, WhenShapeSensitivitiesCalledTwice_GetSameObjec
 {
   auto solid_solver               = createContactSolver(nonlinear_opts, dyn_opts, mat);
   auto [qoi1, shape_sensitivity1] = computeContactQoiSensitivities(*solid_solver, tsInfo);
+  auto [qoi2, shape_sensitivity2] = computeContactQoiSensitivities(*solid_solver, tsInfo);
 
+  EXPECT_EQ(qoi1, qoi2);
+
   solid_solver->resetStates();
   FiniteElementState derivative_direction(shape_sensitivity1.space(), "derivative_direction");
   fillDirection(derivative_direction);
 
-  auto [qoi2, shape_sensitivity2] = computeContactQoiSensitivities(*solid_solver, tsInfo);
-
-  EXPECT_EQ(qoi1, qoi2);
-
   double directional_deriv1 = innerProduct(derivative_direction, shape_sensitivity1);
   double directional_deriv2 = innerProduct(derivative_direction, shape_sensitivity2);
   EXPECT_EQ(directional_deriv1, directional_deriv2);