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,
};
 
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([[maybe_unused]] double time, [[maybe_unused]] double dt,
                        [[maybe_unused]] const std::vector<ConstFieldPtr>& fields) const
  {
    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);
    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
    mfem::Vector gaps = gaps_hpv;
    return gaps;
  };
 
  std::unique_ptr<mfem::HypreParMatrix> jacobian([[maybe_unused]] double time, [[maybe_unused]] double dt,
                                                 [[maybe_unused]] const std::vector<ConstFieldPtr>& fields,
                                                 [[maybe_unused]] int direction) const
  {
    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;
 
    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);
        }
      }
    }
 
    SLIC_ERROR_IF(J_contact->IsZeroBlock(iblock, jblock), "attempting to extract a null block");
    auto dgdu = dynamic_cast<mfem::HypreParMatrix*>(&J_contact->GetBlock(iblock, jblock));
    std::unique_ptr<mfem::HypreParMatrix> dgdu_unique(dgdu);
    return dgdu_unique;
  };
 
 protected:
  mutable ContactData contact_;
 
  ContactOptions contact_opts_;
 
  int interaction_id_;
};
 
}  // namespace smith
 
#endif