ComputeWeightedGapLMMechanicalContact.C

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
 
#include "ComputeWeightedGapLMMechanicalContact.h"
#include "DisplacedProblem.h"
#include "Assembly.h"
#include "MortarContactUtils.h"
#include "WeightedGapUserObject.h"
#include "metaphysicl/metaphysicl_version.h"
#include "metaphysicl/dualsemidynamicsparsenumberarray.h"
#include "metaphysicl/parallel_dualnumber.h"
#if METAPHYSICL_MAJOR_VERSION < 2
#include "metaphysicl/parallel_dynamic_std_array_wrapper.h"
#else
#include "metaphysicl/parallel_dynamic_array_wrapper.h"
#endif
#include "metaphysicl/parallel_semidynamicsparsenumberarray.h"
#include "timpi/parallel_sync.h"
 
registerMooseObject("ContactApp", ComputeWeightedGapLMMechanicalContact);
 
namespace
{
const InputParameters &
assignVarsInParams(const InputParameters & params_in)
{
  InputParameters & ret = const_cast<InputParameters &>(params_in);
  const auto & disp_x_name = ret.get<std::vector<VariableName>>("disp_x");
  if (disp_x_name.size() != 1)
    mooseError("We require that the disp_x parameter have exactly one coupled name");
 
  // We do this so we don't get any variable errors during MortarConstraint(Base) construction
  ret.set<VariableName>("secondary_variable") = disp_x_name[0];
  ret.set<VariableName>("primary_variable") = disp_x_name[0];
 
  return ret;
}
}
 
InputParameters
ComputeWeightedGapLMMechanicalContact::validParams()
{
  InputParameters params = ADMortarConstraint::validParams();
  // Mechanical contact residuals use only primal values and normals
  params.set<bool>("_compute_primal_gradients") = false;
  params.addClassDescription("Computes the weighted gap that will later be used to enforce the "
                             "zero-penetration mechanical contact conditions");
  params.suppressParameter<VariableName>("secondary_variable");
  params.suppressParameter<VariableName>("primary_variable");
  params.addRequiredCoupledVar("disp_x", "The x displacement variable");
  params.addRequiredCoupledVar("disp_y", "The y displacement variable");
  params.addCoupledVar("disp_z", "The z displacement variable");
  params.addParam<Real>(
      "c", 1e6, "Parameter for balancing the size of the gap and contact pressure");
  params.addParam<bool>(
      "normalize_c",
      false,
      "Whether to normalize c by weighting function norm. When unnormalized "
      "the value of c effectively depends on element size since in the constraint we compare nodal "
      "Lagrange Multiplier values to integrated gap values (LM nodal value is independent of "
      "element size, where integrated values are dependent on element size).");
  // The nodal normal at a secondary node depends on every face incident to it, so the Jacobian
  // needs the secondary face point-neighbor ghosted and coupled.
  params.set<bool>("ghost_point_neighbors") = true;
  params.set<bool>("use_displaced_mesh") = true;
  params.set<bool>("interpolate_normals") = false;
  params.addRequiredParam<UserObjectName>("weighted_gap_uo", "The weighted gap user object");
  return params;
}
 
ComputeWeightedGapLMMechanicalContact::ComputeWeightedGapLMMechanicalContact(
    const InputParameters & parameters)
  : ADMortarConstraint(assignVarsInParams(parameters)),
    _secondary_disp_x(adCoupledValue("disp_x")),
    _primary_disp_x(adCoupledNeighborValue("disp_x")),
    _secondary_disp_y(adCoupledValue("disp_y")),
    _primary_disp_y(adCoupledNeighborValue("disp_y")),
    _has_disp_z(isCoupled("disp_z")),
    _secondary_disp_z(_has_disp_z ? &adCoupledValue("disp_z") : nullptr),
    _primary_disp_z(_has_disp_z ? &adCoupledNeighborValue("disp_z") : nullptr),
    _c(getParam<Real>("c")),
    _normalize_c(getParam<bool>("normalize_c")),
    _nodal(getVar("disp_x", 0)->feType().family == LAGRANGE),
    _disp_x_var(getVar("disp_x", 0)),
    _disp_y_var(getVar("disp_y", 0)),
    _disp_z_var(_has_disp_z ? getVar("disp_z", 0) : nullptr),
    _weighted_gap_uo(getUserObject<WeightedGapUserObject>("weighted_gap_uo"))
{
  if (!_weighted_gap_uo.usesNodalNormalDerivatives())
    paramError("weighted_gap_uo",
               "Nodal-normal derivatives are not supported by user object '",
               _weighted_gap_uo.name(),
               "'.");
 
  // The user object and constraint must assemble the same interface and displacement variables.
  if (secondarySubdomain() != _weighted_gap_uo.secondarySubdomain() ||
      primarySubdomain() != _weighted_gap_uo.primarySubdomain())
    paramError("weighted_gap_uo",
               "'weighted_gap_uo' must be defined on the same secondary/primary subdomain pair as "
               "this constraint when nodal-normal derivatives are enabled.");
 
  const std::array<std::pair<const MooseVariable *, unsigned int>, 3> displacement_variables{
      {{_disp_x_var, 0}, {_disp_y_var, 1}, {_disp_z_var, 2}}};
  for (const auto & [variable, component] : displacement_variables)
    if (variable != _weighted_gap_uo.dispVar(component))
      paramError("weighted_gap_uo",
                 "'weighted_gap_uo' must use the same displacement variables as this constraint "
                 "when nodal-normal derivatives are enabled.");
 
  if (!getParam<bool>("use_displaced_mesh"))
    paramError(
        "use_displaced_mesh",
        "'use_displaced_mesh' must be true for the ComputeWeightedGapLMMechanicalContact object");
 
  if (getParam<bool>("interpolate_normals"))
    paramError("interpolate_normals",
               "Mechanical mortar contact uses normalized secondary nodal normals and cannot be "
               "combined with quadrature-point normal interpolation.");
 
  if (!_var->isNodal())
    if (_var->feType().order != static_cast<Order>(0))
      mooseError("Normal contact constraints only support elemental variables of CONSTANT order");
}
 
ADReal
ComputeWeightedGapLMMechanicalContact::computeQpResidual(Moose::MortarType)
{
  mooseError("We should never call computeQpResidual for ComputeWeightedGapLMMechanicalContact");
}
 
void
ComputeWeightedGapLMMechanicalContact::computeQpProperties()
{
}
 
void
ComputeWeightedGapLMMechanicalContact::computeQpIProperties()
{
}
 
void
ComputeWeightedGapLMMechanicalContact::residualSetup()
{
}
 
void
ComputeWeightedGapLMMechanicalContact::jacobianSetup()
{
}
 
void
ComputeWeightedGapLMMechanicalContact::computeResidual(const Moose::MortarType /*mortar_type*/)
{
}
 
void
ComputeWeightedGapLMMechanicalContact::computeJacobian(const Moose::MortarType mortar_type)
{
  // During "computeResidual" and "computeJacobian" we are actually just computing properties on the
  // mortar segment element mesh. We are *not* actually assembling into the residual/Jacobian. For
  // the zero-penetration constraint, the property of interest is the map from node to weighted gap.
  // Computation of the properties proceeds identically for residual and Jacobian evaluation hence
  // why we simply call computeResidual here. We will assemble into the residual/Jacobian later from
  // the post() method
  computeResidual(mortar_type);
}
 
void
ComputeWeightedGapLMMechanicalContact::post()
{
  parallel_object_only();
 
  const auto & dof_to_weighted_gap = _weighted_gap_uo.dofToWeightedGap();
 
  for (const auto & [dof_object, weighted_gap_pr] : dof_to_weighted_gap)
  {
    if (dof_object->processor_id() != this->processor_id())
      continue;
 
    const auto & [weighted_gap, normalization] = weighted_gap_pr;
    _weighted_gap_ptr = &weighted_gap;
    _normalization_ptr = &normalization;
 
    enforceConstraintOnDof(dof_object);
  }
}
 
void
ComputeWeightedGapLMMechanicalContact::incorrectEdgeDroppingPost(
    const std::unordered_set<const Node *> & inactive_lm_nodes)
{
  const auto & dof_to_weighted_gap = _weighted_gap_uo.dofToWeightedGap();
 
  for (const auto & [dof_object, weighted_gap_pr] : dof_to_weighted_gap)
  {
    if ((inactive_lm_nodes.find(static_cast<const Node *>(dof_object)) !=
         inactive_lm_nodes.end()) ||
        (dof_object->processor_id() != this->processor_id()))
      continue;
 
    const auto & [weighted_gap, normalization] = weighted_gap_pr;
    _weighted_gap_ptr = &weighted_gap;
    _normalization_ptr = &normalization;
 
    enforceConstraintOnDof(dof_object);
  }
}
 
void
ComputeWeightedGapLMMechanicalContact::enforceConstraintOnDof(const DofObject * const dof)
{
  const auto & weighted_gap = *_weighted_gap_ptr;
  // normalizeCDivisor() returns *_normalization_ptr unchanged unless node-based scaling is active,
  // in which case it substitutes a coverage-independent (full-element) integral so the two
  // mechanisms don't double-count coverage (PR review discussion; not from Popp 2013).
  const Real c =
      _normalize_c ? _c / _weighted_gap_uo.normalizeCDivisor(dof, *_normalization_ptr) : _c;
 
  // Scaling factor kappa_j (Popp 2013 eq. 34; 1 when disabled/fully covered): dividing the gap by
  // kappa_j gives g_bar/kappa_j, preserving the complementarity root with zhat_j = kappa_j
  // lambda_j.
  const Real kappa = _weighted_gap_uo.nodalScale(dof);
 
  const auto dof_index = dof->dof_number(_sys.number(), _var->number(), 0);
  ADReal lm_value = (*_sys.currentSolution())(dof_index);
  Moose::derivInsert(lm_value.derivatives(), dof_index, 1.);
 
  const ADReal dof_residual = std::min(lm_value, weighted_gap * c / kappa);
 
  addResidualsAndJacobian(_assembly,
                          std::array<ADReal, 1>{{dof_residual}},
                          std::array<dof_id_type, 1>{{dof_index}},
                          _var->scalingFactor());
}