contact_constraint.hpp

Go to the documentation of this file.

// Copyright Lawrence Livermore National Security, LLC and
// other Smith Project Developers. See the top-level LICENSE file for
// details.
//
// SPDX-License-Identifier: (BSD-3-Clause)
 
#pragma once
 
#include <vector>
 
#include "smith/smith_config.hpp"
 
#ifdef SMITH_USE_TRIBOL
 
#include "tribol/interface/tribol.hpp"
 
#include "smith/physics/constraint.hpp"
#include "smith/physics/field_types.hpp"
#include "smith/physics/contact/contact_config.hpp"
#include "smith/physics/contact/contact_data.hpp"
 
namespace mfem {
class Vector;
class HypreParMatrix;
}  // namespace mfem
 
namespace smith {
 
enum ContactFields
{
  SHAPE,
  DISP,
};
 
static std::unique_ptr<mfem::HypreParMatrix> obtainBlock(mfem::BlockOperator* block_operator, int iblock, int jblock)
{
  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");
  auto Ablk = dynamic_cast<mfem::HypreParMatrix*>(&block_operator->GetBlock(iblock, jblock));
  SLIC_ERROR_IF(!Ablk, "failed cast block to mfem::HypreParMatrix");
  // deep copy --> unique_ptr
  auto Ablk_unique = std::make_unique<mfem::HypreParMatrix>(*Ablk);
  return Ablk_unique;
};
 
class FiniteElementState;
 
class ContactConstraint : public Constraint {
 public:
  ContactConstraint(int interaction_id, const mfem::ParMesh& mesh, const std::set<int>& bdry_attr_surf1,
                    const std::set<int>& bdry_attr_surf2, ContactOptions contact_opts,
                    const std::string& name = "contact_constraint")
      : Constraint(name), contact_(mesh), contact_opts_{contact_opts}
  {
    contact_opts_.enforcement = ContactEnforcement::LagrangeMultiplier;
    contact_.addContactInteraction(interaction_id, bdry_attr_surf1, bdry_attr_surf2, contact_opts_);
    interaction_id_ = interaction_id;
  }
 
  virtual ~ContactConstraint() {}
 
  mfem::Vector evaluate(double time, double dt, const std::vector<ConstFieldPtr>& fields,
                        bool update_fields = true) const override
  {
    if (update_fields) {
      // note: Tribol does not use cycle.
      int cycle = 0;
      contact_.updateGeometry(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
      cached_gap_.SetSize(gaps_hpv.Size());
      cached_gap_.Set(1.0, gaps_hpv);
    }
    return cached_gap_;
  };
 
  std::unique_ptr<mfem::HypreParMatrix> jacobian(double time, double dt, const std::vector<ConstFieldPtr>& fields,
                                                 int direction, bool update_fields = true,
                                                 [[maybe_unused]] bool fresh_derivative = true) const override
  {
    SLIC_ERROR_IF(direction != ContactFields::DISP, "requesting a non displacement-field derivative");
    // 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) {
      int cycle = 0;
      if (update_fields) {
        contact_.updateGeometry(cycle, time, dt, *fields[ContactFields::SHAPE], *fields[ContactFields::DISP], true);
      } else {
        contact_.updateGeometry(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
    auto dgdu = obtainBlock(J_contact_.get(), 1, 0);
    return dgdu;
  };
 
  mfem::Vector residual_contribution(double time, double dt, const std::vector<ConstFieldPtr>& fields,
                                     const mfem::Vector& multipliers, int direction, bool update_fields = true,
                                     bool fresh_derivative = true) const override
  {
    SLIC_ERROR_IF(direction != ContactFields::DISP, "requesting a non displacement-field derivative");
    int cycle = 0;
    if (update_fields || fresh_derivative) {
      if (update_fields) {
        contact_.updateForcesAndJacobian(cycle, time, dt, *fields[ContactFields::SHAPE], *fields[ContactFields::DISP],
                                         multipliers);
      } else {
        contact_.updateForcesAndJacobian(cycle, time, dt, std::nullopt, std::nullopt, multipliers);
      }
    }
    return contact_.forces();
  };
 
  std::unique_ptr<mfem::HypreParMatrix> residual_contribution_jacobian(double time, double dt,
                                                                       const std::vector<ConstFieldPtr>& fields,
                                                                       const mfem::Vector& multipliers, int direction,
                                                                       bool update_fields = true,
                                                                       bool fresh_derivative = true) const override
  {
    SLIC_ERROR_IF(direction != ContactFields::DISP, "requesting a non displacement-field derivative");
 
    int cycle = 0;
    if (update_fields || fresh_derivative) {
      if (update_fields) {
        contact_.updateForcesAndJacobian(cycle, time, dt, *fields[ContactFields::SHAPE], *fields[ContactFields::DISP],
                                         multipliers);
      } else {
        contact_.updateForcesAndJacobian(cycle, time, dt, std::nullopt, std::nullopt, multipliers);
      }
      J_contact_ = contact_.mergedJacobian();
    }
    // obtain (0, 0) block entry from the 2 x 2 block contact linear system
    auto Hessian = obtainBlock(J_contact_.get(), 0, 0);
    return Hessian;
  };
 
  std::unique_ptr<mfem::HypreParMatrix> jacobian_tilde(double time, double dt, const std::vector<ConstFieldPtr>& fields,
                                                       int direction, bool update_fields = true,
                                                       [[maybe_unused]] bool fresh_derivative = true) const override
  {
    SLIC_ERROR_IF(direction != ContactFields::DISP, "requesting a non displacement-field derivative");
 
    int cycle = 0;
    if (update_fields || fresh_derivative) {
      mfem::Vector p = contact_.mergedPressures();
      if (update_fields) {
        contact_.updateForcesAndJacobian(cycle, time, dt, *fields[ContactFields::SHAPE], *fields[ContactFields::DISP],
                                         p);
      } else {
        contact_.updateForcesAndJacobian(cycle, time, dt, std::nullopt, std::nullopt, p);
      }
      J_contact_ = contact_.mergedJacobian();
    }
    // obtain (0, 1) block entry from the 2 x 2 block contact linear system
    auto dgduT = obtainBlock(J_contact_.get(), 0, 1);
    std::unique_ptr<mfem::HypreParMatrix> dgdu(dgduT->Transpose());
    return dgdu;
  };
 
  int numPressureDofs() const { return contact_.numPressureDofs(); }
 
 protected:
  mutable ContactData contact_;
 
  ContactOptions contact_opts_;
 
  int interaction_id_;
 
  mutable std::unique_ptr<mfem::BlockOperator> J_contact_;
 
  mutable mfem::Vector cached_gap_;
};
 
}  // namespace smith
 
#endif