diff --git a/ContinuationSolvers b/ContinuationSolvers index 14048e7c7..b7c26c473 160000 --- a/ContinuationSolvers +++ b/ContinuationSolvers @@ -1 +1 @@ -Subproject commit 14048e7c714236127666712a07d32d5eb8b0f435 +Subproject commit b7c26c473ffb1aaaefc82832b6e1f625b239992c diff --git a/README.md b/README.md index 62afbdab2..e43b3ff23 100644 --- a/README.md +++ b/README.md @@ -70,6 +70,10 @@ PackageName: BLT PackageHomePage: https://github.com/LLNL/blt PackageLicenseDeclared: BSD-3-Clause +PackageName: ContinuationSolvers +PackageHomePage: https://github.com/LLNL/ContinuationSolvers +PackageLicenseDeclared: BSD-3-Clause + PackageName: MFEM PackageHomePage: https://github.com/mfem/mfem PackageLicenseDeclared: BSD-3-Clause diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index b98f72a5c..f78aa7014 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -26,6 +26,9 @@ add_subdirectory(conduction) # Add the contact examples add_subdirectory(contact) +# Add the contact-homotopy examples +add_subdirectory(contact/homotopy) + # Add the buckling examples add_subdirectory(buckling) diff --git a/examples/contact/CMakeLists.txt b/examples/contact/CMakeLists.txt index a74a24ba5..8263292a3 100644 --- a/examples/contact/CMakeLists.txt +++ b/examples/contact/CMakeLists.txt @@ -10,7 +10,6 @@ if(TRIBOL_FOUND AND STRUMPACK_DIR) ironing.cpp sphere.cpp twist.cpp - constraint_twist.cpp ) foreach(filename ${CONTACT_EXAMPLES_SOURCES}) diff --git a/examples/contact/constraint_twist.cpp b/examples/contact/constraint_twist.cpp deleted file mode 100644 index efb1bd155..000000000 --- a/examples/contact/constraint_twist.cpp +++ /dev/null @@ -1,78 +0,0 @@ -// Copyright (c) Lawrence Livermore National Security, LLC and -// other Smith Project Developers. See the top-level LICENSE file for -// details. -// -// SPDX-License-Identifier: (BSD-3-Clause) - -#include - -#include -#include - -#include "axom/slic.hpp" - -#include "mfem.hpp" - -#include "smith/smith.hpp" - -// this example is intended to eventually replace twist.cpp -int main(int argc, char* argv[]) -{ - // Initialize and automatically finalize MPI and other libraries - smith::ApplicationManager applicationManager(argc, argv); - - // NOTE: p must be equal to 1 to work with Tribol's mortar method - constexpr int p = 1; - // NOTE: dim must be equal to 3 - constexpr int dim = 3; - - using VectorSpace = smith::H1; - - // Create DataStore - std::string name = "contact_twist_example"; - axom::sidre::DataStore datastore; - smith::StateManager::initialize(datastore, name + "_data"); - - // Construct the appropriate dimension mesh and give it to the data store - std::string filename = SMITH_REPO_DIR "/data/meshes/twohex_for_contact.mesh"; - auto mesh = std::make_shared(smith::buildMeshFromFile(filename), "twist_mesh", 3, 0); - - mesh->addDomainOfBoundaryElements("fixed_surface", smith::by_attr(3)); - mesh->addDomainOfBoundaryElements("driven_surface", smith::by_attr(6)); - - smith::ContactOptions contact_options{.method = smith::ContactMethod::SingleMortar, - .enforcement = smith::ContactEnforcement::LagrangeMultiplier, - .type = smith::ContactType::Frictionless, - .jacobian = smith::ContactJacobian::Exact}; - - std::string contact_constraint_name = "default_contact"; - - // Specify the contact interaction - auto contact_interaction_id = 0; - std::set surface_1_boundary_attributes({4}); - std::set surface_2_boundary_attributes({5}); - smith::ContactConstraint contact_constraint(contact_interaction_id, mesh->mfemParMesh(), - surface_1_boundary_attributes, surface_2_boundary_attributes, - contact_options, contact_constraint_name); - - smith::FiniteElementState shape = smith::StateManager::newState(VectorSpace{}, "shape", mesh->tag()); - smith::FiniteElementState disp = smith::StateManager::newState(VectorSpace{}, "displacement", mesh->tag()); - - std::vector contact_states; - contact_states = {shape, disp}; - // initialize displacement - contact_states[smith::ContactFields::DISP].setFromFieldFunction([](smith::tensor x) { - auto u = 0.1 * x; - return u; - }); - - contact_states[smith::ContactFields::SHAPE] = 0.0; - - double time = 0.0, dt = 1.0; - int direction = 0; - auto input_states = getConstFieldPointers(contact_states); - auto gap = contact_constraint.evaluate(time, dt, input_states); - auto gap_Jacobian = contact_constraint.jacobian(time, dt, input_states, direction); - - return 0; -} diff --git a/examples/contact/homotopy/CMakeLists.txt b/examples/contact/homotopy/CMakeLists.txt new file mode 100644 index 000000000..1209453ee --- /dev/null +++ b/examples/contact/homotopy/CMakeLists.txt @@ -0,0 +1,22 @@ +# Copyright (c) Lawrence Livermore National Security, LLC and +# other Serac Project Developers. See the top-level LICENSE file for +# details. +# +# SPDX-License-Identifier: (BSD-3-Clause) + +if (TRIBOL_FOUND AND STRUMPACK_DIR AND SMITH_ENABLE_CONTINUATION) + blt_add_executable( NAME two_blocks_example + SOURCES two_blocks.cpp + OUTPUT_DIR ${EXAMPLE_OUTPUT_DIRECTORY} + DEPENDS_ON smith_physics smith_mesh_utils continuation_solver + ) + smith_add_example_test(NAME two_blocks + COMMAND two_blocks_example + NUM_MPI_TASKS 4) + install( + FILES + two_blocks.cpp + DESTINATION + examples/smith/two_blocks + ) +endif() diff --git a/examples/contact/homotopy/two_blocks.cpp b/examples/contact/homotopy/two_blocks.cpp new file mode 100644 index 000000000..34492cd54 --- /dev/null +++ b/examples/contact/homotopy/two_blocks.cpp @@ -0,0 +1,394 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and +// other Serac Project Developers. See the top-level LICENSE file for +// details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#include + +#include +#include + +#include "axom/slic.hpp" + +#include "mfem.hpp" +#include "smith/physics/contact/contact_config.hpp" +#include "shared/mesh/MeshBuilder.hpp" + +// ContinuationSolver headers +#include "problems/Problems.hpp" +#include "solvers/HomotopySolver.hpp" +#include "utilities.hpp" + +#include "smith/smith.hpp" + +static constexpr int dim = 3; +static constexpr int disp_order = 1; + +using VectorSpace = smith::H1; +using DensitySpace = smith::L2; + +using SolidMaterial = smith::solid_mechanics::NeoHookeanWithFieldDensity; +using SolidWeakFormT = smith::SolidWeakForm>; + +enum FIELD +{ + DISP = SolidWeakFormT::DISPLACEMENT, + VELO = SolidWeakFormT::VELOCITY, + ACCEL = SolidWeakFormT::ACCELERATION, + DENSITY = SolidWeakFormT::NUM_STATES +}; + +/* Nonlinear problem of the form + * F(X) = [ r(u) + (dc/du)^T l ] = [ 0 ] + * [ -c(u) ] [ 0 ] + * X = [ u ] + * [ l ] + * + * wherein r(u) is the elasticity nonlinear residual + * c(u) are the tied gap contacts + * u are the displacement dofs + * l are the Lagrange multipliers + * + * This problem inherits from EqualityConstrainedHomotopyProblem + * for compatibility with the HomotopySolver. + */ +template +class TiedContactProblem : public EqualityConstrainedHomotopyProblem { + protected: + std::unique_ptr drdu_; + std::unique_ptr dcdu_; + std::vector contact_states_; + std::vector residual_states_; + std::shared_ptr weak_form_; + std::unique_ptr shape_disp_; + std::shared_ptr mesh_; + std::shared_ptr constraints_; + double time_ = 0.0; + double dt_ = 0.0; + std::vector jacobian_weights_ = {1.0, 0.0, 0.0, 0.0}; + + public: + TiedContactProblem(std::vector contact_states, std::vector residual_states, + std::shared_ptr mesh, std::shared_ptr weak_form, + std::shared_ptr constraints, mfem::Array fixed_tdof_list, + mfem::Array disp_tdof_list, mfem::Vector uDC); + mfem::Vector residual(const mfem::Vector& u, bool new_point) const; + mfem::HypreParMatrix* residualJacobian(const mfem::Vector& u, bool new_point); + mfem::Vector constraint(const mfem::Vector& u, bool new_point) const; + mfem::HypreParMatrix* constraintJacobian(const mfem::Vector& u, bool new_point); + mfem::Vector constraintJacobianTvp(const mfem::Vector& u, const mfem::Vector& l, bool new_point) const; + void fullDisplacement(const mfem::Vector& x, mfem::Vector& u); + virtual ~TiedContactProblem(); +}; + +class ParaviewWriter { + public: + using StateVecs = std::vector>; + using DualVecs = std::vector>; + + ParaviewWriter(std::unique_ptr pv_, const StateVecs& states_) + : pv(std::move(pv_)), states(states_) + { + } + + ParaviewWriter(std::unique_ptr pv_, const StateVecs& states_, const StateVecs& duals_) + : pv(std::move(pv_)), states(states_), dual_states(duals_) + { + } + + void write(int step, double time, const std::vector& current_states) + { + SMITH_MARK_FUNCTION; + SLIC_ERROR_ROOT_IF(current_states.size() != states.size(), "wrong number of output states to write"); + + for (size_t n = 0; n < states.size(); ++n) { + auto& state = states[n]; + *state = *current_states[n]; + state->gridFunction(); + } + + pv->SetCycle(step); + pv->SetTime(time); + pv->Save(); + } + + private: + std::unique_ptr pv; + StateVecs states; + StateVecs dual_states; +}; + +auto createParaviewOutput(const mfem::ParMesh& mesh, const std::vector& states, + std::string output_name) +{ + if (output_name == "") { + output_name = "default"; + } + + ParaviewWriter::StateVecs output_states; + for (const auto& s : states) { + output_states.push_back(std::make_shared(s->space(), s->name())); + } + + auto non_const_mesh = const_cast(&mesh); + auto paraview_dc = std::make_unique(output_name, non_const_mesh); + int max_order_in_fields = 0; + + // Find the maximum polynomial order in the physics module's states + for (const auto& state : output_states) { + paraview_dc->RegisterField(state->name(), &state->gridFunction()); + max_order_in_fields = std::max(max_order_in_fields, state->space().GetOrder(0)); + } + + // Set the options for the paraview output files + paraview_dc->SetLevelsOfDetail(max_order_in_fields); + paraview_dc->SetHighOrderOutput(true); + paraview_dc->SetDataFormat(mfem::VTKFormat::BINARY); + paraview_dc->SetCompression(true); + + return ParaviewWriter(std::move(paraview_dc), output_states, {}); +} + +// this example is intended to eventually replace twist.cpp +int main(int argc, char* argv[]) +{ + // Initialize and automatically finalize MPI and other libraries + smith::ApplicationManager applicationManager(argc, argv); + + int visualize = 1; + int visualize_all_iterates = 1; + // command line arguments + axom::CLI::App app{"Two block contact."}; + app.add_option("--visualize", visualize, "solution visualization") + ->default_val("1") // Matches value set above + ->check(axom::CLI::Range(0, 1)); + app.set_help_flag("--help"); + + CLI11_PARSE(app, argc, argv); + + // Create DataStore + std::string name = "two_block_example"; + axom::sidre::DataStore datastore; + smith::StateManager::initialize(datastore, name + "_data"); + + // Construct the appropriate dimension mesh and give it to the data store + std::string filename = SMITH_REPO_DIR "/data/meshes/twohex_for_contact.mesh"; + auto mesh = std::make_shared(smith::buildMeshFromFile(filename), "two_block_mesh", 3, 0); + + mesh->addDomainOfBoundaryElements("fixed_surface", smith::by_attr(3)); + mesh->addDomainOfBoundaryElements("driven_surface", smith::by_attr(6)); + + smith::ContactOptions contact_options{.method = smith::ContactMethod::SingleMortar, + .enforcement = smith::ContactEnforcement::LagrangeMultiplier, + .type = smith::ContactType::TiedNormal, + .jacobian = smith::ContactJacobian::Exact}; + + std::string contact_constraint_name = "default_contact"; + + // Specify the contact interaction + auto contact_interaction_id = 0; + std::set surface_1_boundary_attributes({4}); + std::set surface_2_boundary_attributes({5}); + auto contact_constraint = std::make_shared( + contact_interaction_id, mesh->mfemParMesh(), surface_1_boundary_attributes, surface_2_boundary_attributes, + contact_options, contact_constraint_name); + + smith::FiniteElementState shape = smith::StateManager::newState(VectorSpace{}, "shape", mesh->tag()); + smith::FiniteElementState disp = smith::StateManager::newState(VectorSpace{}, "displacement", mesh->tag()); + smith::FiniteElementState velo = smith::StateManager::newState(VectorSpace{}, "velocity", mesh->tag()); + smith::FiniteElementState accel = smith::StateManager::newState(VectorSpace{}, "acceleration", mesh->tag()); + smith::FiniteElementState density = smith::StateManager::newState(DensitySpace{}, "density", mesh->tag()); + + std::vector contact_states; + std::vector states; + std::vector params; + contact_states = {shape, disp}; + states = {disp, velo, accel}; + params = {density}; + + // initialize displacement + contact_states[smith::ContactFields::DISP].setFromFieldFunction([](smith::tensor x) { + auto u = 0.0 * x; + return u; + }); + + contact_states[smith::ContactFields::SHAPE] = 0.0; + states[FIELD::VELO] = 0.0; + states[FIELD::ACCEL] = 0.0; + params[0] = 1.0; + + std::string physics_name = "solid"; + + // construct residual + auto solid_mechanics_weak_form = + std::make_shared(physics_name, mesh, states[FIELD::DISP].space(), getSpaces(params)); + + SolidMaterial mat; + mat.K = 1.0; + mat.G = 0.5; + solid_mechanics_weak_form->setMaterial(smith::DependsOn<0>{}, mesh->entireBodyName(), mat); + + smith::tensor constant_force{}; + for (int i = 0; i < dim; i++) { + constant_force[i] = 0.0; + } + constant_force[dim - 1] = 0.0; + + solid_mechanics_weak_form->addBodyIntegral(mesh->entireBodyName(), [constant_force](double /* t */, auto x) { + return smith::tuple{constant_force, 0.0 * smith::get(x)}; + }); + + auto residual_state_ptrs = smith::getFieldPointers(states, params); + auto contact_state_ptrs = smith::getFieldPointers(contact_states); + + // homogeneous Dirichlet boundary conditions + mfem::Array ess_fixed_tdof_list; + mfem::Array ess_disp_tdof_list; + mfem::Array ess_bdr_marker(mesh->mfemParMesh().bdr_attributes.Max()); + ess_bdr_marker = 0; + ess_bdr_marker[2] = 1; + ess_bdr_marker[5] = 0; + states[FIELD::DISP].space().GetEssentialTrueDofs(ess_bdr_marker, ess_fixed_tdof_list); + ess_bdr_marker = 0; + ess_bdr_marker[5] = 1; + states[FIELD::DISP].space().GetEssentialTrueDofs(ess_bdr_marker, ess_disp_tdof_list); + mfem::Vector uDC(states[FIELD::DISP].space().GetTrueVSize()); + uDC = 0.0; + auto applied_displacement = [](smith::tensor /*x*/) { + smith::tensor u{}; + u[0] = 0.0; + u[1] = 0.0; + u[2] = -0.06; + return u; + }; + states[FIELD::DISP].setFromFieldFunction(applied_displacement); + uDC.Set(1.0, states[FIELD::DISP]); + + TiedContactProblem problem(contact_state_ptrs, residual_state_ptrs, mesh, solid_mechanics_weak_form, + contact_constraint, ess_fixed_tdof_list, ess_disp_tdof_list, uDC); + double nonlinear_absolute_tol = 1.e-6; + int nonlinear_max_iterations = 30; + int nonlinear_print_level = 1; + // optimization variables + auto X0 = problem.GetOptimizationVariable(); + auto Xf = problem.GetOptimizationVariable(); + + HomotopySolver solver(&problem); + solver.SetTol(nonlinear_absolute_tol); + solver.SetMaxIter(nonlinear_max_iterations); + solver.SetPrintLevel(nonlinear_print_level); + solver.EnableRegularizedNewtonMode(); + solver.EnableSaveIterates(); + solver.Mult(X0, Xf); + bool converged = solver.GetConverged(); + SLIC_WARNING_ROOT_IF(!converged, "Homotopy solver did not converge"); + + // visualize + auto writer = createParaviewOutput(mesh->mfemParMesh(), smith::getConstFieldPointers(states), "contact"); + if (visualize) { + mfem::Vector u(states[FIELD::DISP].space().GetTrueVSize()); + u = problem.GetDisplacement(X0); + states[FIELD::DISP].Set(1.0, u); + writer.write(0, 0.0, smith::getConstFieldPointers(states)); + if (!visualize_all_iterates) { + u = problem.GetDisplacement(Xf); + states[FIELD::DISP].Set(1.0, u); + writer.write(1, 1.0, smith::getConstFieldPointers(states)); + } else { + auto iterates = solver.GetIterates(); + for (int i = 0; i < iterates.Size(); i++) { + u = problem.GetDisplacement(*iterates[i]); + states[FIELD::DISP].Set(1.0, u); + writer.write((i + 1), static_cast(i + 1), smith::getConstFieldPointers(states)); + } + } + } + return 0; +} + +template +TiedContactProblem::TiedContactProblem(std::vector contact_states, + std::vector residual_states, + std::shared_ptr mesh, + std::shared_ptr weak_form, + std::shared_ptr constraints, + mfem::Array fixed_tdof_list, + mfem::Array disp_tdof_list, mfem::Vector uDC) + : EqualityConstrainedHomotopyProblem(fixed_tdof_list, disp_tdof_list, uDC), + weak_form_(weak_form), + mesh_(mesh), + constraints_(constraints) +{ + // copy residual states + residual_states_.resize(residual_states.size()); + std::copy(residual_states.begin(), residual_states.end(), residual_states_.begin()); + + // copy contact states + contact_states_.resize(contact_states.size()); + std::copy(contact_states.begin(), contact_states.end(), contact_states_.begin()); + + // number of "internal" non-essential dofs + const int dimu = + residual_states[FIELD::DISP]->space().GetTrueVSize() - fixed_tdof_list.Size() - disp_tdof_list.Size(); + // number of contact constraints + const int dimc = constraints->numPressureDofs(); + // EqualityConstrainedHomotopyProblem utility function + SetSizes(dimu, dimc); + + // shape_disp field + shape_disp_ = std::make_unique(mesh->newShapeDisplacement()); +} + +template +mfem::Vector TiedContactProblem::residual(const mfem::Vector& u, bool /*new_point*/) const +{ + residual_states_[FIELD::DISP]->Set(1.0, u); + auto res = weak_form_->residual(time_, dt_, shape_disp_.get(), smith::getConstFieldPointers(residual_states_)); + return res; +}; + +template +mfem::HypreParMatrix* TiedContactProblem::residualJacobian(const mfem::Vector& u, bool /*new_point*/) +{ + residual_states_[FIELD::DISP]->Set(1.0, u); + drdu_ = weak_form_->jacobian(time_, dt_, shape_disp_.get(), smith::getConstFieldPointers(residual_states_), + jacobian_weights_); + return drdu_.get(); +} + +template +mfem::Vector TiedContactProblem::constraint(const mfem::Vector& u, bool /*new_point*/) const +{ + bool new_point = true; + contact_states_[smith::ContactFields::DISP]->Set(1.0, u); + auto gap = constraints_->evaluate(time_, dt_, smith::getConstFieldPointers(contact_states_), new_point); + return gap; +} + +template +mfem::HypreParMatrix* TiedContactProblem::constraintJacobian(const mfem::Vector& u, + bool /*new_point*/) +{ + bool new_point = true; + contact_states_[smith::ContactFields::DISP]->Set(1.0, u); + dcdu_ = constraints_->jacobian(time_, dt_, smith::getConstFieldPointers(contact_states_), smith::ContactFields::DISP, + new_point); + return dcdu_.get(); +} + +template +mfem::Vector TiedContactProblem::constraintJacobianTvp(const mfem::Vector& u, const mfem::Vector& l, + bool /*new_point*/) const +{ + bool new_point = true; + contact_states_[smith::ContactFields::DISP]->Set(1.0, u); + auto res_contribution = constraints_->residual_contribution(time_, dt_, smith::getConstFieldPointers(contact_states_), + l, smith::ContactFields::DISP, new_point); + return res_contribution; +} + +template +TiedContactProblem::~TiedContactProblem() +{ +} diff --git a/examples/inertia_relief/inertia_relief_example.cpp b/examples/inertia_relief/inertia_relief_example.cpp index e99e796f9..fb02919dd 100644 --- a/examples/inertia_relief/inertia_relief_example.cpp +++ b/examples/inertia_relief/inertia_relief_example.cpp @@ -111,49 +111,46 @@ auto createParaviewOutput(const mfem::ParMesh& mesh, const std::vector= 0 - * Q(x, y) = 0 - * Here, F and x are both 0-dimensional - * and Q(x, y) = [ r(u) + (dc/du)^T l] - * [-c(u)] - * y = [ u ] - * [ l ] - * we use the approximate Jacobian - * dQ/dy \approx [ dr/du (dc/du)^T] - * [-dc/du 0 ] - * we note that the sign-convention with regard to "c" is important - * as the approximate Jacobian is positive semi-definite when dr/du is - * and thus the NLMC problem is guaranteed to be semi-monotone. +/* Nonlinear problem of the form + * F(X) = [ r(u) + (dc/du)^T l ] = [ 0 ] + * [ -c(u) ] [ 0 ] + * X = [ u ] + * [ l ] + * + * wherein r(u) is the elasticity nonlinear residual + * c(u) are the tied gap contacts + * u are the displacement dofs + * l are the Lagrange multipliers + * + * This problem inherits from EqualityConstrainedHomotopyProblem + * for compatibility with the HomotopySolver. */ class InertialReliefProblem : public EqualityConstrainedHomotopyProblem { protected: - mfem::HypreParMatrix* drdu_ = nullptr; - mfem::HypreParMatrix* dcdu_ = nullptr; - int dimu_; - int dimc_; - int dimuglb_; - int dimcglb_; - mfem::Array y_partition_; - std::vector obj_states_; - std::vector all_states_; - std::shared_ptr weak_form_; - std::unique_ptr shape_disp_; + std::unique_ptr drdu_; // Jacobian of residual + std::unique_ptr dcdu_; // Jacobian of constraint + std::vector obj_states_; // states for objective evaluation + std::vector all_states_; // states for weak_form evaluation + std::shared_ptr weak_form_; // weak_form + std::unique_ptr shape_disp_; // shape displacement std::shared_ptr mesh_; - std::vector> constraints_; - double time_ = 0.0; - double dt_ = 0.0; - std::vector jacobian_weights_ = {1.0, 0.0, 0.0, 0.0}; + std::vector> constraints_; // vector of constraints + double time_ = 0.0; // parameter for constraint and weak_form member function calls + double dt_ = 0.0; // parameter for constraint and weak_form member function calls + std::vector jacobian_weights_ = {1.0, 0.0, 0.0, 0.0}; // weights for weak_form_->jacobian calls + mutable mfem::Vector constraint_cached_; + mutable mfem::Vector residual_cached_; + mutable mfem::Vector JTvp_cached_; public: InertialReliefProblem(std::vector obj_states, std::vector all_states, std::shared_ptr mesh, std::shared_ptr weak_form, std::vector> constraints); - mfem::Vector residual(const mfem::Vector& u) const; - mfem::Vector constraintJacobianTvp(const mfem::Vector& u, const mfem::Vector& l) const; - mfem::Vector constraint(const mfem::Vector& u) const; - mfem::HypreParMatrix* constraintJacobian(const mfem::Vector& u); - mfem::HypreParMatrix* residualJacobian(const mfem::Vector& u); + mfem::Vector residual(const mfem::Vector& u, bool new_point) const; + mfem::Vector constraintJacobianTvp(const mfem::Vector& u, const mfem::Vector& l, bool new_point) const; + mfem::Vector constraint(const mfem::Vector& u, bool new_point) const; + mfem::HypreParMatrix* constraintJacobian(const mfem::Vector& u, bool new_point); + mfem::HypreParMatrix* residualJacobian(const mfem::Vector& u, bool new_point); virtual ~InertialReliefProblem(); }; @@ -227,7 +224,7 @@ int main(int argc, char* argv[]) // construct residual auto solid_mechanics_weak_form = - std::make_shared(physics_name, mesh, states[DISP].space(), getSpaces(params)); + std::make_shared(physics_name, mesh, states[FIELD::DISP].space(), getSpaces(params)); SolidMaterial mat; mat.K = 1.0; @@ -258,9 +255,9 @@ int main(int argc, char* argv[]) double time = 0.0; double dt = 1.0; auto all_states = getConstFieldPointers(states, params); - auto objective_states = {all_states[DISP], all_states[DENSITY]}; + auto objective_states = {all_states[FIELD::DISP], all_states[FIELD::DENSITY]}; - ObjectiveT::SpacesT param_space_ptrs{&all_states[DISP]->space(), &all_states[DENSITY]->space()}; + ObjectiveT::SpacesT param_space_ptrs{&all_states[FIELD::DISP]->space(), &all_states[FIELD::DENSITY]->space()}; ObjectiveT mass_objective("mass constraining", mesh, param_space_ptrs); @@ -268,7 +265,7 @@ int main(int argc, char* argv[]) [](double /*t*/, auto /*X*/, auto RHO) { return get(RHO); }); double mass = mass_objective.evaluate(time, dt, shape_disp.get(), objective_states); - smith::tensor initial_cg; + smith::tensor initial_cg; // center of gravity for (int i = 0; i < dim; ++i) { auto cg_objective = std::make_shared("translation " + std::to_string(i), mesh, param_space_ptrs); @@ -309,8 +306,8 @@ int main(int argc, char* argv[]) } auto non_const_states = getFieldPointers(states, params); // create an inertial relief problem - InertialReliefProblem problem({non_const_states[DISP], non_const_states[DENSITY]}, non_const_states, mesh, - solid_mechanics_weak_form, constraints); + InertialReliefProblem problem({non_const_states[FIELD::DISP], non_const_states[FIELD::DENSITY]}, non_const_states, + mesh, solid_mechanics_weak_form, constraints); // optimization variables auto X0 = problem.GetOptimizationVariable(); @@ -321,8 +318,9 @@ int main(int argc, char* argv[]) // set solver options solver.SetTol(nonlinear_absolute_tol); solver.SetMaxIter(nonlinear_max_iterations); - + solver.EnableRegularizedNewtonMode(); // solve the inertia relief problem + solver.SetPrintLevel(2); solver.Mult(X0, Xf); // extract displacement and Lagrange multipliers mfem::Vector displacement_sol = problem.GetDisplacement(Xf); @@ -333,6 +331,12 @@ int main(int argc, char* argv[]) double multiplier_norm = mfem::GlobalLpNorm(2, multiplier_sol.Norml2(), MPI_COMM_WORLD); SLIC_INFO_ROOT(axom::fmt::format("||displacement|| = {}", displacement_norm)); SLIC_INFO_ROOT(axom::fmt::format("||multiplier|| = {}", multiplier_norm)); + auto adjoint = problem.GetOptimizationVariable(); + auto adjoint_load = problem.GetOptimizationVariable(); + adjoint = 0.0; + adjoint_load = 1.0; + problem.AdjointSolve(displacement_sol, adjoint_load, adjoint); + if (visualize) { writer.write(1, 1.0, objective_states); } @@ -358,136 +362,132 @@ InertialReliefProblem::InertialReliefProblem(std::vectorspace().GetTrueDofOffsets(); - dimc_ = static_cast(constraints_.size()); + int dim_displacement = all_states_[FIELD::DISP]->space().GetTrueVSize(); + int dim_constraints = static_cast(constraints_.size()); int myid = mfem::Mpi::WorldRank(); - cOffsets[0] = 0; - cOffsets[1] = dimc_; if (myid > 0) { - dimc_ = 0; - cOffsets[0] = cOffsets[1]; + dim_constraints = 0; } - - dimu_ = all_states_[FIELD::DISP]->space().GetTrueVSize(); - MPI_Allreduce(&dimc_, &dimcglb_, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(&dimu_, &dimuglb_, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); - - SetSizes(uOffsets, cOffsets); + SetSizes(dim_displacement, dim_constraints); + + constraint_cached_.SetSize(dim_constraints); + constraint_cached_ = 0.0; + residual_cached_.SetSize(dim_displacement); + residual_cached_ = 0.0; + JTvp_cached_.SetSize(dim_displacement); + JTvp_cached_ = 0.0; } // residual callback -mfem::Vector InertialReliefProblem::residual(const mfem::Vector& u) const +mfem::Vector InertialReliefProblem::residual(const mfem::Vector& u, bool new_point) const { - obj_states_[DISP]->Set(1.0, u); - auto res_vector = weak_form_->residual(time_, dt_, shape_disp_.get(), smith::getConstFieldPointers(all_states_)); - return res_vector; + if (new_point) { + obj_states_[FIELD::DISP]->Set(1.0, u); + residual_cached_.Set( + 1.0, weak_form_->residual(time_, dt_, shape_disp_.get(), smith::getConstFieldPointers(all_states_))); + } + return residual_cached_; } // constraint Jacobian transpose vector product -mfem::Vector InertialReliefProblem::constraintJacobianTvp(const mfem::Vector& u, const mfem::Vector& l) const +mfem::Vector InertialReliefProblem::constraintJacobianTvp(const mfem::Vector& u, const mfem::Vector& l, + bool new_point) const { - obj_states_[DISP]->Set(1.0, u); - std::vector multipliers(constraints_.size()); - for (int i = 0; i < dimc_; i++) { - multipliers[static_cast(i)] = l(i); - } - const int nconstraints = static_cast(constraints_.size()); - MPI_Bcast(multipliers.data(), nconstraints, MPI_DOUBLE, 0, MPI_COMM_WORLD); - - mfem::Vector constraint_gradient(dimu_); - constraint_gradient = 0.0; - mfem::Vector output_vec(dimu_); - output_vec = 0.0; - - for (size_t i = 0; i < constraints_.size(); i++) { - mfem::Vector grad_temp = - constraints_[i]->gradient(time_, dt_, shape_disp_.get(), smith::getConstFieldPointers(obj_states_), DISP); - constraint_gradient.Set(1.0, grad_temp); - output_vec.Add(multipliers[i], constraint_gradient); + if (new_point) { + int dim_constraints = GetMultiplierDim(); + int dim_displacement = GetDisplacementDim(); + obj_states_[FIELD::DISP]->Set(1.0, u); + std::vector multipliers(constraints_.size()); + for (int i = 0; i < dim_constraints; i++) { + multipliers[static_cast(i)] = l(i); + } + const int nconstraints = static_cast(constraints_.size()); + MPI_Bcast(multipliers.data(), nconstraints, MPI_DOUBLE, 0, MPI_COMM_WORLD); + + mfem::Vector constraint_gradient(dim_displacement); + constraint_gradient = 0.0; + JTvp_cached_ = 0.0; + for (size_t i = 0; i < constraints_.size(); i++) { + mfem::Vector grad_temp = constraints_[i]->gradient(time_, dt_, shape_disp_.get(), + smith::getConstFieldPointers(obj_states_), FIELD::DISP); + constraint_gradient.Set(1.0, grad_temp); + JTvp_cached_.Add(multipliers[i], constraint_gradient); + } } - return output_vec; + return JTvp_cached_; } // Jacobian of the residual -mfem::HypreParMatrix* InertialReliefProblem::residualJacobian(const mfem::Vector& u) +mfem::HypreParMatrix* InertialReliefProblem::residualJacobian(const mfem::Vector& u, bool new_point) { - obj_states_[DISP]->Set(1.0, u); - auto drdu_unique = - weak_form_->jacobian(time_, dt_, shape_disp_.get(), getConstFieldPointers(all_states_), jacobian_weights_); - - if (drdu_) { - delete drdu_; + if (new_point) { + obj_states_[FIELD::DISP]->Set(1.0, u); + drdu_.reset(); + drdu_ = weak_form_->jacobian(time_, dt_, shape_disp_.get(), getConstFieldPointers(all_states_), jacobian_weights_); } - drdu_ = drdu_unique.release(); - SLIC_ERROR_ROOT_IF(drdu_->Height() != dimu_ || drdu_->Width() != dimu_, "residual Jacobian of an unexpected shape"); - return drdu_; + int dim_displacement = GetDisplacementDim(); + SLIC_ERROR_ROOT_IF(drdu_->Height() != dim_displacement || drdu_->Width() != dim_displacement, + "residual Jacobian of an unexpected shape"); + return drdu_.get(); } // constraint callback -mfem::Vector InertialReliefProblem::constraint(const mfem::Vector& u) const +mfem::Vector InertialReliefProblem::constraint(const mfem::Vector& u, bool new_point) const { - obj_states_[DISP]->Set(1.0, u); - mfem::Vector output_vec(dimc_); - output_vec = 0.0; - - for (size_t i = 0; i < constraints_.size(); i++) { - const int idx = static_cast(i); - const size_t i2 = static_cast(idx); - SLIC_ERROR_ROOT_IF(i2 != i, "Constraint index is out of range, bad cast from size_t to int"); - - double constraint_i = - constraints_[i]->evaluate(time_, dt_, shape_disp_.get(), smith::getConstFieldPointers(obj_states_)); - if (dimc_ > 0) { - output_vec(idx) = constraint_i; + if (new_point) { + int dim_constraints = GetMultiplierDim(); + obj_states_[FIELD::DISP]->Set(1.0, u); + + for (size_t i = 0; i < constraints_.size(); i++) { + const int idx = static_cast(i); + const size_t i2 = static_cast(idx); + SLIC_ERROR_ROOT_IF(i2 != i, "Constraint index is out of range, bad cast from size_t to int"); + + double constraint_i = + constraints_[i]->evaluate(time_, dt_, shape_disp_.get(), smith::getConstFieldPointers(obj_states_)); + if (dim_constraints > 0) { + constraint_cached_(idx) = constraint_i; + } } } - return output_vec; + return constraint_cached_; } // Jacobian of the constraint -mfem::HypreParMatrix* InertialReliefProblem::constraintJacobian(const mfem::Vector& u) +mfem::HypreParMatrix* InertialReliefProblem::constraintJacobian(const mfem::Vector& u, bool new_point) { - obj_states_[DISP]->Set(1.0, u); - int myid = mfem::Mpi::WorldRank(); - int nentries = dimuglb_; - if (myid > 0) { - nentries = 0; - } - mfem::SparseMatrix dcdumat(dimc_, dimuglb_, nentries); - mfem::Array cols; - cols.SetSize(dimuglb_); - for (int i = 0; i < dimuglb_; i++) { - cols[i] = i; - } - for (size_t i = 0; i < constraints_.size(); i++) { - const int idx = static_cast(i); - const size_t i2 = static_cast(idx); - SLIC_ERROR_ROOT_IF(i2 != i, "Constraint index is out of range, bad cast from size_t to int"); - mfem::HypreParVector gradVector(MPI_COMM_WORLD, dimuglb_, uOffsets_); - gradVector.Set( - 1.0, constraints_[i]->gradient(time_, dt_, shape_disp_.get(), smith::getConstFieldPointers(obj_states_), DISP)); - mfem::Vector* globalGradVector = gradVector.GlobalVector(); - if (myid == 0) { - dcdumat.SetRow(idx, cols, *globalGradVector); + int dim_constraints = GetMultiplierDim(); + int glbdim_displacement = GetGlobalDisplacementDim(); + if (new_point) { + obj_states_[FIELD::DISP]->Set(1.0, u); + // dense rows + int nentries = glbdim_displacement; + if (dimc_ == 0) { + nentries = 0; } - delete globalGradVector; - } - dcdumat.Threshold(1.e-20); - dcdumat.Finalize(); - if (dcdu_) { - delete dcdu_; + mfem::SparseMatrix dcdumat(dim_constraints, glbdim_displacement, nentries); + mfem::Array cols; + cols.SetSize(glbdim_displacement); + for (int i = 0; i < glbdim_displacement; i++) { + cols[i] = i; + } + std::unique_ptr globalGradVector; + for (size_t i = 0; i < constraints_.size(); i++) { + const int idx = static_cast(i); + const size_t i2 = static_cast(idx); + SLIC_ERROR_ROOT_IF(i2 != i, "Constraint index is out of range, bad cast from size_t to int"); + mfem::HypreParVector gradVector(MPI_COMM_WORLD, glbdim_displacement, uOffsets_); + gradVector.Set(1.0, constraints_[i]->gradient(time_, dt_, shape_disp_.get(), + smith::getConstFieldPointers(obj_states_), FIELD::DISP)); + globalGradVector.reset(gradVector.GlobalVector()); + if (dim_constraints > 0) { + dcdumat.SetRow(idx, cols, *globalGradVector.get()); + } + } + dcdumat.Finalize(); + dcdu_.reset(GenerateHypreParMatrixFromSparseMatrix(cOffsets_, uOffsets_, &dcdumat)); } - dcdu_ = GenerateHypreParMatrixFromSparseMatrix(uOffsets_, cOffsets_, &dcdumat); - return dcdu_; + return dcdu_.get(); } -InertialReliefProblem::~InertialReliefProblem() -{ - if (drdu_) { - delete drdu_; - } - if (dcdu_) { - delete dcdu_; - } -} +InertialReliefProblem::~InertialReliefProblem() {} diff --git a/src/smith/physics/constraint.hpp b/src/smith/physics/constraint.hpp index f0954e735..1c081e91e 100644 --- a/src/smith/physics/constraint.hpp +++ b/src/smith/physics/constraint.hpp @@ -41,9 +41,11 @@ class Constraint { * @param time time * @param dt time step * @param fields vector of smith::FiniteElementState* + * @param new_point boolean indicating if this is a new evaluation point or not * @return mfem::Vector which is the constraint evaluation */ - virtual mfem::Vector evaluate(double time, double dt, const std::vector& fields) const = 0; + virtual mfem::Vector evaluate(double time, double dt, const std::vector& fields, + bool new_point = true) const = 0; /** @brief Virtual interface for computing constraint Jacobian from a vector of smith::FiniteElementState* * @@ -51,11 +53,12 @@ class Constraint { * @param dt time step * @param fields vector of smith::FiniteElementState* * @param direction index for which field to take the gradient with respect to + * @param new_point boolean indicating if this is a new evaluation point or not * @return std::unique_ptr */ virtual std::unique_ptr jacobian(double time, double dt, - const std::vector& fields, - int direction) const = 0; + const std::vector& fields, int direction, + bool new_point = true) const = 0; /** @brief Virtual interface for computing constraint Jacobian_tilde from a vector of smith::FiniteElementState* * Jacobian_tilde is an optional approximation of the true Jacobian @@ -64,13 +67,14 @@ class Constraint { * @param dt time step * @param fields vector of smith::FiniteElementState* * @param direction index for which field to take the gradient with respect to + * @param new_point boolean indicating if this is a new evaluation point or not * @return std::unique_ptr */ virtual std::unique_ptr jacobian_tilde(double time, double dt, - const std::vector& fields, - int direction) const + const std::vector& fields, int direction, + bool new_point = true) const { - return jacobian(time, dt, fields, direction); + return jacobian(time, dt, fields, direction, new_point); }; /** @brief Virtual interface for computing residual contribution Jacobian_tilde^T multiplier from a vector of @@ -81,12 +85,14 @@ class Constraint { * @param fields vector of smith::FiniteElementState* * @param multipliers mfem::Vector of Lagrange multipliers * @param direction index for which field to take the gradient with respect to + * @param new_point boolean indicating if this is a new evaluation point or not * @return std::Vector */ virtual mfem::Vector residual_contribution(double time, double dt, const std::vector& fields, - const mfem::Vector& multipliers, int direction) const + const mfem::Vector& multipliers, int direction, + bool new_point = true) const { - std::unique_ptr jac = jacobian_tilde(time, dt, fields, direction); + std::unique_ptr jac = jacobian_tilde(time, dt, fields, direction, new_point); mfem::Vector y(jac->Width()); y = 0.0; SLIC_ERROR_ROOT_IF(jac->Height() != multipliers.Size(), "Incompatible matrix and vector sizes."); @@ -102,12 +108,13 @@ class Constraint { * @param fields vector of smith::FiniteElementState* * @param multipliers mfem::Vector of Lagrange multipliers * @param direction index for which field to take the gradient with respect to + * @param new_point boolean indicating if this is a new evaluation point or not * @return std::unique_ptr */ virtual std::unique_ptr residual_contribution_jacobian( [[maybe_unused]] double time, [[maybe_unused]] double dt, [[maybe_unused]] const std::vector& fields, [[maybe_unused]] const mfem::Vector& multipliers, - [[maybe_unused]] int direction) const + [[maybe_unused]] int direction, [[maybe_unused]] bool new_point = true) const { SLIC_ERROR_ROOT(axom::fmt::format("Base class must override residual_contribution_jacobian before usage")); std::unique_ptr res_contr_jacobian = nullptr; diff --git a/src/smith/physics/contact_constraint.hpp b/src/smith/physics/contact_constraint.hpp index fed0b7540..bb36b1408 100644 --- a/src/smith/physics/contact_constraint.hpp +++ b/src/smith/physics/contact_constraint.hpp @@ -42,12 +42,51 @@ enum ContactFields DISP, }; +/** @brief Interface for extracting the iblock, jblock block from a std::unique_ptr + * said block is returned as a std::unique_ptr if possible + * All other blocks are deleted + * + * @param block_operator the block operator + * @param iblock row block index + * @param jblock column block index + * @return std::unique_ptr The requested block of block_operator + */ +static std::unique_ptr safelyObtainBlock(mfem::BlockOperator* block_operator, int iblock, + int jblock, bool own_blocks = false) +{ + SLIC_ERROR_IF(iblock < 0 || jblock < 0, "block indicies must be non-negative"); + SLIC_ERROR_IF(iblock > block_operator->NumRowBlocks() || jblock > block_operator->NumColBlocks(), + "one or more block indicies are too large and the requested block does not exist"); + SLIC_ERROR_IF(block_operator->IsZeroBlock(iblock, jblock), "attempting to extract a null block"); + block_operator->owns_blocks = false; + for (int i = 0; i < block_operator->NumRowBlocks(); i++) { + for (int j = 0; j < block_operator->NumColBlocks(); j++) { + if (i == iblock && j == jblock) { + continue; + } + if (!block_operator->IsZeroBlock(i, j) && !own_blocks) { + delete &block_operator->GetBlock(i, j); + } + } + } + auto Ablk = dynamic_cast(&block_operator->GetBlock(iblock, jblock)); + SLIC_ERROR_IF(!Ablk, "failed cast block to mfem::HypreParMatrix"); + if (own_blocks) { + // deep copy --> unique_ptr + auto Ablk_unique = std::make_unique(*Ablk); + return Ablk_unique; + } else { + std::unique_ptr Ablk_unique(Ablk); + return Ablk_unique; + } +}; + class FiniteElementState; /** * @brief A ContactConstraint defines a gap constraint associated to contact problem * - * This class stores the details of a single contact interaction between two surfaces. It also interfaces provides a + * This class stores the details of a single contact interaction between two surfaces. It also provides a * description of a contact constraint given by a single contact interaction. A ContactConstraint can have a single * ContactInteraction and will default to LagrangeMultiplier as it will be up to the solver that takes this * ContactConstraint to determine how it will enforce the constraint. @@ -83,17 +122,21 @@ class ContactConstraint : public Constraint { * @param time time * @param dt time step * @param fields vector of smith::FiniteElementState* + * @param new_point boolean indicating if this is a new point or not * @return mfem::Vector which is the constraint evaluation */ - mfem::Vector evaluate([[maybe_unused]] double time, [[maybe_unused]] double dt, - [[maybe_unused]] const std::vector& fields) const + mfem::Vector evaluate(double time, double dt, const std::vector& fields, + bool new_point = true) const override { contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]); tribol::setLagrangeMultiplierOptions(interaction_id_, tribol::ImplicitEvalMode::MORTAR_GAP); - // note: Tribol does not use cycle. - int cycle = 0; - contact_.update(cycle, time, dt); + if (new_point) { + // note: Tribol does not use cycle. + int cycle = 0; + contact_.update(cycle, time, dt); + pressures_set_ = false; + } 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 @@ -107,39 +150,148 @@ class ContactConstraint : public Constraint { * @param dt time step * @param fields vector of smith::FiniteElementState* * @param direction index for which field to take the gradient with respect to - * @return std::unique_ptr + * @param new_point boolean indicating if this is a new point or not + * @return std::unique_ptr The true Jacobian */ - std::unique_ptr jacobian([[maybe_unused]] double time, [[maybe_unused]] double dt, - [[maybe_unused]] const std::vector& fields, - [[maybe_unused]] int direction) const + std::unique_ptr jacobian(double time, double dt, const std::vector& fields, + int direction, bool new_point = true) const override { + SLIC_ERROR_IF(direction != ContactFields::DISP, "requesting a non displacement-field derivative"); contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]); tribol::setLagrangeMultiplierOptions(interaction_id_, tribol::ImplicitEvalMode::MORTAR_JACOBIAN); - int cycle = 0; - contact_.update(cycle, time, dt); - auto J_contact = contact_.mergedJacobian(); - J_contact->owns_blocks = false; - + if (new_point) { + int cycle = 0; + contact_.update(cycle, time, dt); + J_contact_ = contact_.mergedJacobian(); + pressures_set_ = false; + } int iblock = 1; int jblock = 0; - for (int i = 0; i < 2; i++) { - for (int j = 0; j < 2; j++) { - if (i == iblock && j == jblock) { - continue; - } - if (!J_contact->IsZeroBlock(i, j)) { - delete &J_contact->GetBlock(i, j); - } + auto dgdu = safelyObtainBlock(J_contact_.get(), iblock, jblock, own_blocks_); + return dgdu; + }; + + /** @brief Interface for computing residual contribution Jacobian_tilde^T multiplier from a vector of + * smith::FiniteElementState* + * + * @param time time + * @param dt time step + * @param fields vector of smith::FiniteElementState* + * @param multipliers mfem::Vector of Lagrange multipliers + * @param direction index for which field to take the gradient with respect to + * @param new_point boolean indicating if this is a new point or not + * @return std::Vector + */ + mfem::Vector residual_contribution(double time, double dt, const std::vector& fields, + const mfem::Vector& multipliers, int direction, + bool new_point = true) const override + { + SLIC_ERROR_IF(direction != ContactFields::DISP, "requesting a non displacement-field derivative"); + contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]); + tribol::setLagrangeMultiplierOptions(interaction_id_, tribol::ImplicitEvalMode::MORTAR_GAP); + + int cycle = 0; + if (new_point) { + // we need to call update first to update gaps + for (auto& interaction : contact_.getContactInteractions()) { + interaction.evalJacobian(false); } + contact_.update(cycle, time, dt); + pressures_set_ = false; + } + if (!pressures_set_) { + // 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); + } + contact_.update(cycle, time, dt); + pressures_set_ = true; } - SLIC_ERROR_IF(J_contact->IsZeroBlock(iblock, jblock), "attempting to extract a null block"); - auto dgdu = dynamic_cast(&J_contact->GetBlock(iblock, jblock)); - std::unique_ptr dgdu_unique(dgdu); - return dgdu_unique; + return contact_.forces(); }; + /** @brief Interface for computing Jacobians of the residual contribution from a vector of + * smith::FiniteElementState* + * + * @param time time + * @param dt time step + * @param fields vector of smith::FiniteElementState* + * @param multipliers mfem::Vector of Lagrange multipliers + * @param direction index for which field to take the gradient with respect to + * @param new_point boolean indicating if this is a new point or not + * @return std::unique_ptr + */ + std::unique_ptr residual_contribution_jacobian(double time, double dt, + const std::vector& fields, + const mfem::Vector& multipliers, int direction, + bool new_point = true) const override + { + SLIC_ERROR_IF(direction != ContactFields::DISP, "requesting a non displacement-field derivative"); + contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]); + tribol::setLagrangeMultiplierOptions(interaction_id_, tribol::ImplicitEvalMode::MORTAR_JACOBIAN); + + int cycle = 0; + if (new_point) { + // we need to call update first to update gaps + for (auto& interaction : contact_.getContactInteractions()) { + interaction.evalJacobian(false); + } + contact_.update(cycle, time, dt); + pressures_set_ = false; + } + if (!pressures_set_) { + // 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); + } + contact_.update(cycle, time, dt); + J_contact_ = contact_.mergedJacobian(); + pressures_set_ = true; + } + int iblock = 0; + int jblock = 0; + auto Hessian = safelyObtainBlock(J_contact_.get(), iblock, jblock, own_blocks_); + return Hessian; + }; + + /** @brief Interface for computing contact constraint Jacobian_tilde from a vector of smith::FiniteElementState* + * + * @param time time + * @param dt time step + * @param fields vector of smith::FiniteElementState* + * @param direction index for which field to take the gradient with respect to + * @param new_point boolean indicating if this is a new point or not + * @return std::unique_ptr + */ + std::unique_ptr jacobian_tilde(double time, double dt, const std::vector& fields, + int direction, bool new_point = true) const override + { + SLIC_ERROR_IF(direction != ContactFields::DISP, "requesting a non displacement-field derivative"); + contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]); + tribol::setLagrangeMultiplierOptions(interaction_id_, tribol::ImplicitEvalMode::MORTAR_JACOBIAN); + + int cycle = 0; + if (new_point) { + contact_.update(cycle, time, dt); + J_contact_.reset(); + J_contact_ = contact_.mergedJacobian(); + pressures_set_ = false; + } + int iblock = 0; + int jblock = 1; + auto dgduT = safelyObtainBlock(J_contact_.get(), iblock, jblock, own_blocks_); + std::unique_ptr dgdu(dgduT->Transpose()); + return dgdu; + }; + + int numPressureDofs() const { return contact_.numPressureDofs(); } + protected: /** * @brief ContactData which has various contact calls @@ -155,6 +307,21 @@ class ContactConstraint : public Constraint { * @brief interaction_id Unique identifier for the ContactInteraction (used in Tribol) */ int interaction_id_; + + /** + * @brief J_contact_ to hold contact derivatives + */ + mutable std::unique_ptr J_contact_; + + /** + * @brief own_blocks_ temporary boolean + */ + const bool own_blocks_ = true; + + /** + * @brief pressures_set_ are the Lagrange multipliers set + */ + mutable bool pressures_set_ = false; }; } // namespace smith