diff --git a/src/smith/physics/contact/contact_data.cpp b/src/smith/physics/contact/contact_data.cpp index ae8f10c45..60439c83f 100644 --- a/src/smith/physics/contact/contact_data.cpp +++ b/src/smith/physics/contact/contact_data.cpp @@ -155,9 +155,10 @@ std::unique_ptr 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(&interaction_J->GetBlock(0, 0)), @@ -171,42 +172,16 @@ std::unique_ptr 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(&interaction_J->GetBlock(1, 0)); auto dfdp = dynamic_cast(&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(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 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(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(i), 0) = dgdu; - dfdp_blocks(0, static_cast(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(i), 0) = dgdu; + dfdp_blocks(0, static_cast(i)) = dfdp; } } if (haveLagrangeMultipliers()) { diff --git a/src/smith/physics/contact/contact_interaction.cpp b/src/smith/physics/contact/contact_interaction.cpp index 855ac684f..4147519c6 100644 --- a/src/smith/physics/contact/contact_interaction.cpp +++ b/src/smith/physics/contact/contact_interaction.cpp @@ -107,6 +107,68 @@ std::unique_ptr ContactInteraction::jacobian() const return tribol::getMfemBlockJacobian(getInteractionId()); } +std::unique_ptr ContactInteraction::jacobianContribution() const +{ + const int disp_size = current_coords_.ParFESpace()->GetTrueVSize(); + const int pressure_size = numPressureDofs(); + + auto offsets = mfem::Array({0, disp_size, disp_size + pressure_size}); + auto out = std::make_unique(offsets); + out->owns_blocks = true; + + auto interaction_J = jacobian(); + interaction_J->owns_blocks = false; // manage block ownership explicitly + + std::unique_ptr dfdx; + if (!interaction_J->IsZeroBlock(0, 0)) { + auto* block_0_0 = dynamic_cast(&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(&interaction_J->GetBlock(1, 0)); + auto* dfdp = dynamic_cast(&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(dfdp->EliminateCols(inactive_dofs)); + + if (getContactOptions().enforcement == ContactEnforcement::Penalty) { + std::unique_ptr 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 diff --git a/src/smith/physics/contact/contact_interaction.hpp b/src/smith/physics/contact/contact_interaction.hpp index 6d179d235..4d54e78d4 100644 --- a/src/smith/physics/contact/contact_interaction.hpp +++ b/src/smith/physics/contact/contact_interaction.hpp @@ -109,6 +109,20 @@ class ContactInteraction { */ std::unique_ptr 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 jacobianContribution() const; + /** * @brief Get the finite element space of the pressure DOFs * diff --git a/src/smith/physics/solid_mechanics_contact.hpp b/src/smith/physics/solid_mechanics_contact.hpp index 713ad0e21..42844f052 100644 --- a/src/smith/physics/solid_mechanics_contact.hpp +++ b/src/smith/physics/solid_mechanics_contact.hpp @@ -12,6 +12,15 @@ #pragma once +#include +#include +#include +#include +#include +#include +#include +#include + #include "smith/physics/solid_mechanics.hpp" #include "smith/physics/contact/contact_data.hpp" @@ -213,9 +222,105 @@ class SolidMechanicsContact, { SLIC_ERROR_ROOT_IF(!is_quasistatic_, "Contact can only be applied to quasistatic problems."); SLIC_ERROR_ROOT_IF(order > 1, "Contact can only be applied to linear (order = 1) meshes."); + + const auto interaction_force_name = detail::addPrefix(name_, axom::fmt::format("contact_force_{}", interaction_id)); + + // Allow multiple calls with the same interaction_id; the new interaction overwrites the old one. + // Reuse previously-allocated objects so we don't invalidate pointers stored in BasePhysics vectors. + { + auto it = contact_interaction_forces_.find(interaction_id); + if (it == contact_interaction_forces_.end()) { + auto interaction_force_dual = + std::make_unique(StateManager::newDual(displacement_.space(), interaction_force_name)); + *interaction_force_dual = 0.0; + duals_.push_back(interaction_force_dual.get()); + contact_interaction_forces_.emplace(interaction_id, std::move(interaction_force_dual)); + } else { + *it->second = 0.0; + } + } + + { + const auto force_adjoint_bcs_name = + detail::addPrefix(name_, axom::fmt::format("contact_force_adjoint_bcs_{}", interaction_id)); + auto it = contact_interaction_force_adjoint_bcs_.find(interaction_id); + if (it == contact_interaction_force_adjoint_bcs_.end()) { + auto force_adjoint_bcs = std::make_unique(displacement_.space(), force_adjoint_bcs_name); + *force_adjoint_bcs = 0.0; + this->dual_adjoints_.push_back(force_adjoint_bcs.get()); + contact_interaction_force_adjoint_bcs_.emplace(interaction_id, std::move(force_adjoint_bcs)); + } else { + *it->second = 0.0; + } + } + + { + const auto sens_name = + detail::addPrefix(name_, axom::fmt::format("contact_shape_sensitivity_{}", interaction_id)); + auto it = contact_interaction_shape_sensitivities_.find(interaction_id); + if (it == contact_interaction_shape_sensitivities_.end()) { + auto sens = std::make_unique(BasePhysics::shapeDisplacementSensitivity().space(), sens_name); + *sens = 0.0; + contact_interaction_shape_sensitivities_.emplace(interaction_id, std::move(sens)); + } else { + *it->second = 0.0; + } + } + + { + auto insert_pos = std::lower_bound(contact_interaction_ids_sorted_.begin(), contact_interaction_ids_sorted_.end(), + interaction_id); + if (insert_pos == contact_interaction_ids_sorted_.end() || *insert_pos != interaction_id) { + contact_interaction_ids_sorted_.insert(insert_pos, interaction_id); + } + } + contact_.addContactInteraction(interaction_id, bdry_attr_surf1, bdry_attr_surf2, contact_opts); } + /// @overload + std::vector dualNames() const override + { + auto dual_names = SolidMechanicsBase::dualNames(); + dual_names.push_back("contact_forces"); + for (int interaction_id : contact_interaction_ids_sorted_) { + dual_names.push_back(axom::fmt::format("contact_force_{}", interaction_id)); + } + return dual_names; + } + + /// @overload + const FiniteElementDual& dual(const std::string& dual_name) const override + { + if (dual_name == "contact_forces" || dual_name == detail::addPrefix(this->name_, "contact_forces")) { + return forces_; + } + + const auto interaction_id = parseContactInteractionForceId(dual_name); + if (interaction_id.has_value()) { + auto it = contact_interaction_forces_.find(*interaction_id); + if (it != contact_interaction_forces_.end()) { + return *it->second; + } + } + + return SolidMechanicsBase::dual(dual_name); + } + + /// @overload + const FiniteElementState& dualAdjoint(const std::string& dual_name) const override + { + const auto interaction_id = parseContactInteractionForceId(dual_name); + if (interaction_id.has_value()) { + auto it = contact_interaction_force_adjoint_bcs_.find(*interaction_id); + if (it != contact_interaction_force_adjoint_bcs_.end()) { + return *it->second; + } + } + + return SolidMechanicsBase::dualAdjoint(dual_name); + } + /** * @brief create a contactSubspaceTransferOperator for AMGF */ @@ -248,7 +353,77 @@ class SolidMechanicsContact, */ mfem::HypreParVector pressure() const { return contact_.mergedPressures(); } +#ifdef SMITH_USE_TRIBOL + /** + * @brief Get a contact interaction by its interaction id + * + * @param interaction_id The unique identifier for the contact interaction + * @return Reference to the requested contact interaction + */ + const ContactInteraction& contactInteraction(int interaction_id) const + { + for (const auto& interaction : contact_.getContactInteractions()) { + if (interaction.getInteractionId() == interaction_id) { + return interaction; + } + } + SLIC_ERROR_ROOT(axom::fmt::format("No contact interaction found with interaction_id={}", interaction_id)); + return contact_.getContactInteractions().front(); + } +#endif + + /// @overload + void setDualAdjointBcs(std::unordered_map bcs) override + { + SLIC_ERROR_ROOT_IF(bcs.size() == 0, "Adjoint load container size must be greater than 0 in SolidMechanicsContact."); + + auto reaction_adjoint_load = bcs.find("reactions"); + if (reaction_adjoint_load != bcs.end()) { + SolidMechanicsBase::setDualAdjointBcs({{"reactions", reaction_adjoint_load->second}}); + } + + for (const auto& [name, bc] : bcs) { + if (name == "reactions") { + continue; + } + + const auto interaction_id = parseContactInteractionForceId(name); + SLIC_ERROR_ROOT_IF(!interaction_id.has_value(), + axom::fmt::format("Unknown dual adjoint BC '{}' for SolidMechanicsContact.", name)); + + auto it = contact_interaction_force_adjoint_bcs_.find(*interaction_id); + SLIC_ERROR_ROOT_IF( + it == contact_interaction_force_adjoint_bcs_.end(), + axom::fmt::format("No contact force adjoint BC registered for interaction_id={}", *interaction_id)); + + *it->second = bc; + } + } + protected: + /// @brief Converts a dual name into an interaction id (if it exists) + static std::optional parseContactInteractionForceId(std::string_view dual_name) + { + constexpr std::string_view prefix = "contact_force_"; + + // Accept both the bare name and the module-prefixed name, e.g. "solid_contact_force_0". + const auto idx = dual_name.rfind(prefix); + if (idx == std::string_view::npos) { + return std::nullopt; + } + + // This code converts everything after the prefix to a candidate id + const std::string_view id_str = dual_name.substr(idx + prefix.size()); + int interaction_id = -1; + const auto* begin = id_str.data(); + const auto* end = id_str.data() + id_str.size(); + auto [ptr, ec] = std::from_chars(begin, end, interaction_id); + if (ec != std::errc{} || ptr != end) { + return std::nullopt; + } + return interaction_id; + } + /// @brief Solve the Quasi-static Newton system void quasiStaticSolve(double dt) override { @@ -272,6 +447,15 @@ class SolidMechanicsContact, contact_.setPressures(mfem::Vector(augmented_solution, displacement_.Size(), contact_.numPressureDofs())); contact_.update(cycle_, time_, dt); forces_.SetVector(contact_.forces(), 0); + +#ifdef SMITH_USE_TRIBOL + for (const auto& interaction : contact_.getContactInteractions()) { + auto it = contact_interaction_forces_.find(interaction.getInteractionId()); + if (it != contact_interaction_forces_.end()) { + it->second->SetVector(interaction.forces(), 0); + } + } +#endif } /** @@ -430,8 +614,47 @@ class SolidMechanicsContact, auto jacobian = std::unique_ptr(static_cast(&block_J->GetBlock(0, 0))); auto J_T = std::unique_ptr(jacobian->Transpose()); - for (const auto& bc : bcs_.essentials()) { - bc.apply(*J_T, displacement_adjoint_load_, reactions_adjoint_bcs_); + // If a QoI depends on per-interaction contact force duals, the dual-adjoint seeds define an additional contribution + // to the adjoint load: + // dJ/du += (df_i/du)^T * dJ/df_i + // Following SolidMechanics::setAdjointLoad() sign convention, displacement_adjoint_load_ stores the negative of the + // provided dJ/du, so we subtract these contributions here. +#ifdef SMITH_USE_TRIBOL + if (!contact_interaction_force_adjoint_bcs_.empty()) { + FiniteElementDual contact_force_load(displacement_.space(), "contact_force_dual_adjoint_load"); + contact_force_load = 0.0; + + for (const auto& [interaction_id, force_seed] : contact_interaction_force_adjoint_bcs_) { + if (!force_seed) { + continue; + } + + // Only apply if the seed is nonzero. + if (force_seed->Norml2() == 0.0) { + continue; + } + + const auto interaction_J = contactInteraction(interaction_id).jacobianContribution(); + auto* J00 = dynamic_cast(&interaction_J->GetBlock(0, 0)); + SLIC_ERROR_ROOT_IF(!J00, "Expected HypreParMatrix (0,0) block for contact interaction Jacobian."); + + FiniteElementDual tmp(displacement_.space(), "contact_force_dual_adjoint_load_tmp"); + tmp = 0.0; + J00->MultTranspose(*force_seed, tmp); + contact_force_load.Add(1.0, tmp); + } + + displacement_adjoint_load_.Add(-1.0, contact_force_load); + } +#endif + + auto J_e_T = bcs_.eliminateAllEssentialDofsFromMatrix(*J_T); + auto& constrained_dofs = bcs_.allEssentialTrueDofs(); + + mfem::EliminateBC(*J_T, *J_e_T, constrained_dofs, reactions_adjoint_bcs_, displacement_adjoint_load_); + for (int i = 0; i < constrained_dofs.Size(); i++) { + int j = constrained_dofs[i]; + displacement_adjoint_load_[j] = reactions_adjoint_bcs_[j]; } auto& lin_solver = nonlin_solver_->linearSolver(); @@ -514,6 +737,17 @@ class SolidMechanicsContact, /// forces for output FiniteElementDual forces_; + /// per-interaction contact forces for output + std::unordered_map> contact_interaction_forces_; + + /// sorted list of all contact interaction ids + std::vector contact_interaction_ids_sorted_; + + /// per-interaction dual-adjoint (BC) fields for contact force duals + std::unordered_map> contact_interaction_force_adjoint_bcs_; + + /// contact-only shape sensitivities (not stored in StateManager) + std::unordered_map> contact_interaction_shape_sensitivities_; /// contactDOFProlongationOperator std::unique_ptr contact_dof_prolongation_; diff --git a/src/smith/physics/tests/contact_solid_adjoint.cpp b/src/smith/physics/tests/contact_solid_adjoint.cpp index 39bcf3b98..21011c8f9 100644 --- a/src/smith/physics/tests/contact_solid_adjoint.cpp +++ b/src/smith/physics/tests/contact_solid_adjoint.cpp @@ -298,6 +298,126 @@ TEST_F(ContactSensitivityFixture, ReactionShapeSensitivities) EXPECT_NEAR(directional_deriv, directional_deriv_fd, eps); } +TEST_F(ContactSensitivityFixture, SingleContactInteractionForceMagnitudeQoiShapeSensitivities) +{ + // For this test, define a QoI that penalizes the (single-interaction) contact force magnitude: + // J = 0.5 * + // and verify the adjoint-based directional derivative matches a finite-difference check. + constexpr int contact_interaction_id = 0; + + auto solid_solver = createContactSolver(mesh, nonlinear_opts, dyn_opts, mat); + + auto compute_contact_force_qoi = [&](SolidMechT& solver) -> double { + solver.resetStates(); + solver.advanceTimestep(1.0); + const auto f = solver.dual(axom::fmt::format("contact_force_{}", contact_interaction_id)); + return 0.5 * innerProduct(f, f); + }; + + (void)compute_contact_force_qoi(*solid_solver); + + // Compute adjoint/shape sensitivity for this QoI. + const auto f_base = solid_solver->dual(axom::fmt::format("contact_force_{}", contact_interaction_id)); + const auto interaction_jacobian = solid_solver->contactInteraction(contact_interaction_id).jacobianContribution(); + auto* J00 = dynamic_cast(&interaction_jacobian->GetBlock(0, 0)); + SLIC_ERROR_ROOT_IF(!J00, "Expected HypreParMatrix (0,0) block for contact interaction Jacobian."); + + FiniteElementDual displacement_adjoint_load(solid_solver->state("displacement").space(), + "contact_force_qoi_adjoint_load"); + displacement_adjoint_load = 0.0; + J00->MultTranspose(f_base, displacement_adjoint_load); + + // Our shape sensitivity has an explicit term (dJ/d(shape) with u=const) which is only dependent on the contact + // interaction force ( = f_contact,i^T * d(f_contact,i)/d(shape) ) + FiniteElementDual explicit_shape_term(solid_solver->shapeDisplacement().space(), "contact_force_qoi_explicit_shape"); + explicit_shape_term = 0.0; + J00->MultTranspose(f_base, explicit_shape_term); + + solid_solver->setAdjointLoad({{"displacement", displacement_adjoint_load}}); + solid_solver->reverseAdjointTimestep(); + + const FiniteElementDual implicit_shape_term = + static_cast(*solid_solver).computeTimestepShapeSensitivity(); + + FiniteElementDual total_shape_sensitivity(implicit_shape_term.space(), "contact_force_qoi_total_shape"); + total_shape_sensitivity = implicit_shape_term; + total_shape_sensitivity.Add(1.0, explicit_shape_term); + + FiniteElementState derivative_direction(total_shape_sensitivity.space(), "derivative_direction"); + fillDirection(*solid_solver, derivative_direction); + + auto compute_qoi_adjusting_shape = [&](SolidMechT& solver, double scale) -> double { + FiniteElementState shape_disp(solver.shapeDisplacement().space(), "input_shape_displacement"); + shape_disp = 0.0; + shape_disp.Add(scale, derivative_direction); + solver.setShapeDisplacement(shape_disp); + return compute_contact_force_qoi(solver); + }; + + const double qoi_plus = compute_qoi_adjusting_shape(*solid_solver, eps); + const double qoi_minus = compute_qoi_adjusting_shape(*solid_solver, -eps); + + const double directional_deriv = innerProduct(derivative_direction, total_shape_sensitivity); + const double directional_deriv_fd = (qoi_plus - qoi_minus) / (2.0 * eps); + + EXPECT_NEAR(directional_deriv, directional_deriv_fd, 10.0 * eps); +} + +TEST_F(ContactSensitivityFixture, ContactForceDualAdjointBcsMatchesEquivalentAdjointLoad) +{ + // Verify that providing a seed for a contact force dual via setDualAdjointBcs({"contact_force_", seed}) + // produces the same displacement adjoint as providing the equivalent adjoint load + // dJ/du = (df_i/du)^T * dJ/df_i + // directly via setAdjointLoad. + constexpr int contact_interaction_id = 0; + + auto solver = createContactSolver(mesh, nonlinear_opts, dyn_opts, mat); + + solver->resetStates(); + solver->advanceTimestep(1.0); + EXPECT_EQ(1, solver->cycle()); + + FiniteElementState dJ_df(solver->state("displacement").space(), "dJ_df"); + fillDirection(*solver, dJ_df); + + // 1) Solve adjoint with the contact-force seed provided via setDualAdjointBcs. + FiniteElementDual zero_load(solver->state("displacement").space(), "zero_load"); + zero_load = 0.0; + solver->setAdjointLoad({{"displacement", zero_load}}); + solver->setDualAdjointBcs({{axom::fmt::format("contact_force_{}", contact_interaction_id), dJ_df}}); + solver->reverseAdjointTimestep(); + + FiniteElementState lambda_from_seed(solver->adjoint("displacement")); + + // 2) Re-run the forward step and solve the adjoint using the equivalent adjoint load. + solver->resetStates(); + solver->advanceTimestep(1.0); + EXPECT_EQ(1, solver->cycle()); + + const auto interaction_jacobian = solver->contactInteraction(contact_interaction_id).jacobianContribution(); + auto* J00 = dynamic_cast(&interaction_jacobian->GetBlock(0, 0)); + SLIC_ERROR_ROOT_IF(!J00, "Expected HypreParMatrix (0,0) block for contact interaction Jacobian."); + + FiniteElementDual equivalent_load(solver->state("displacement").space(), "equivalent_load"); + equivalent_load = 0.0; + J00->MultTranspose(dJ_df, equivalent_load); + + // Clear the contact-force dual-adjoint seed so the quasiStaticAdjointSolve() path doesn't add it again. + FiniteElementState zero_seed(solver->state("displacement").space(), "zero_seed"); + zero_seed = 0.0; + solver->setDualAdjointBcs({{axom::fmt::format("contact_force_{}", contact_interaction_id), zero_seed}}); + + solver->setAdjointLoad({{"displacement", equivalent_load}}); + solver->reverseAdjointTimestep(); + + const auto& lambda_from_load = solver->adjoint("displacement"); + + FiniteElementState diff(lambda_from_seed); + diff.Add(-1.0, lambda_from_load); + + EXPECT_NEAR(diff.Norml2(), 0.0, 1.0e-10); +} + } // namespace smith int main(int argc, char* argv[])