Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "ComputeWeightedGapLMMechanicalContact.h"
      11             : #include "DisplacedProblem.h"
      12             : #include "Assembly.h"
      13             : #include "MortarContactUtils.h"
      14             : #include "WeightedGapUserObject.h"
      15             : #include "metaphysicl/metaphysicl_version.h"
      16             : #include "metaphysicl/dualsemidynamicsparsenumberarray.h"
      17             : #include "metaphysicl/parallel_dualnumber.h"
      18             : #if METAPHYSICL_MAJOR_VERSION < 2
      19             : #include "metaphysicl/parallel_dynamic_std_array_wrapper.h"
      20             : #else
      21             : #include "metaphysicl/parallel_dynamic_array_wrapper.h"
      22             : #endif
      23             : #include "metaphysicl/parallel_semidynamicsparsenumberarray.h"
      24             : #include "timpi/parallel_sync.h"
      25             : 
      26             : registerMooseObject("ContactApp", ComputeWeightedGapLMMechanicalContact);
      27             : 
      28             : namespace
      29             : {
      30             : const InputParameters &
      31         538 : assignVarsInParams(const InputParameters & params_in)
      32             : {
      33             :   InputParameters & ret = const_cast<InputParameters &>(params_in);
      34         538 :   const auto & disp_x_name = ret.get<std::vector<VariableName>>("disp_x");
      35         538 :   if (disp_x_name.size() != 1)
      36           0 :     mooseError("We require that the disp_x parameter have exactly one coupled name");
      37             : 
      38             :   // We do this so we don't get any variable errors during MortarConstraint(Base) construction
      39        1076 :   ret.set<VariableName>("secondary_variable") = disp_x_name[0];
      40         538 :   ret.set<VariableName>("primary_variable") = disp_x_name[0];
      41             : 
      42         538 :   return ret;
      43             : }
      44             : }
      45             : 
      46             : InputParameters
      47        1073 : ComputeWeightedGapLMMechanicalContact::validParams()
      48             : {
      49        1073 :   InputParameters params = ADMortarConstraint::validParams();
      50        1073 :   params.addClassDescription("Computes the weighted gap that will later be used to enforce the "
      51             :                              "zero-penetration mechanical contact conditions");
      52        1073 :   params.suppressParameter<VariableName>("secondary_variable");
      53        1073 :   params.suppressParameter<VariableName>("primary_variable");
      54        2146 :   params.addRequiredCoupledVar("disp_x", "The x displacement variable");
      55        2146 :   params.addRequiredCoupledVar("disp_y", "The y displacement variable");
      56        2146 :   params.addCoupledVar("disp_z", "The z displacement variable");
      57        2146 :   params.addParam<Real>(
      58        2146 :       "c", 1e6, "Parameter for balancing the size of the gap and contact pressure");
      59        2146 :   params.addParam<bool>(
      60             :       "normalize_c",
      61        2146 :       false,
      62             :       "Whether to normalize c by weighting function norm. When unnormalized "
      63             :       "the value of c effectively depends on element size since in the constraint we compare nodal "
      64             :       "Lagrange Multiplier values to integrated gap values (LM nodal value is independent of "
      65             :       "element size, where integrated values are dependent on element size).");
      66        1073 :   params.set<bool>("ghost_point_neighbors") = true;
      67        1073 :   params.suppressParameter<bool>("ghost_point_neighbors");
      68        1073 :   params.set<bool>("use_displaced_mesh") = true;
      69        1073 :   params.set<bool>("interpolate_normals") = false;
      70        2146 :   params.addRequiredParam<UserObjectName>("weighted_gap_uo", "The weighted gap user object");
      71        1073 :   return params;
      72           0 : }
      73             : 
      74         538 : ComputeWeightedGapLMMechanicalContact::ComputeWeightedGapLMMechanicalContact(
      75         538 :     const InputParameters & parameters)
      76             :   : ADMortarConstraint(assignVarsInParams(parameters)),
      77         538 :     _secondary_disp_x(adCoupledValue("disp_x")),
      78         538 :     _primary_disp_x(adCoupledNeighborValue("disp_x")),
      79         538 :     _secondary_disp_y(adCoupledValue("disp_y")),
      80         538 :     _primary_disp_y(adCoupledNeighborValue("disp_y")),
      81         538 :     _has_disp_z(isCoupled("disp_z")),
      82         538 :     _secondary_disp_z(_has_disp_z ? &adCoupledValue("disp_z") : nullptr),
      83         538 :     _primary_disp_z(_has_disp_z ? &adCoupledNeighborValue("disp_z") : nullptr),
      84        1076 :     _c(getParam<Real>("c")),
      85        1076 :     _normalize_c(getParam<bool>("normalize_c")),
      86         538 :     _nodal(getVar("disp_x", 0)->feType().family == LAGRANGE),
      87         538 :     _disp_x_var(getVar("disp_x", 0)),
      88         538 :     _disp_y_var(getVar("disp_y", 0)),
      89         744 :     _disp_z_var(_has_disp_z ? getVar("disp_z", 0) : nullptr),
      90        1076 :     _weighted_gap_uo(getUserObject<WeightedGapUserObject>("weighted_gap_uo"))
      91             : {
      92         538 :   _weighted_gap_uo.includeNodalNormalDerivatives();
      93             : 
      94        1076 :   if (!getParam<bool>("use_displaced_mesh"))
      95           0 :     paramError(
      96             :         "use_displaced_mesh",
      97             :         "'use_displaced_mesh' must be true for the ComputeWeightedGapLMMechanicalContact object");
      98             : 
      99        1076 :   if (getParam<bool>("interpolate_normals"))
     100           3 :     paramError("interpolate_normals",
     101             :                "Mechanical mortar contact uses normalized secondary nodal normals and cannot be "
     102             :                "combined with quadrature-point normal interpolation.");
     103             : 
     104         535 :   if (!_var->isNodal())
     105           0 :     if (_var->feType().order != static_cast<Order>(0))
     106           0 :       mooseError("Normal contact constraints only support elemental variables of CONSTANT order");
     107         535 : }
     108             : 
     109             : ADReal
     110           0 : ComputeWeightedGapLMMechanicalContact::computeQpResidual(Moose::MortarType)
     111             : {
     112           0 :   mooseError("We should never call computeQpResidual for ComputeWeightedGapLMMechanicalContact");
     113             : }
     114             : 
     115             : void
     116           0 : ComputeWeightedGapLMMechanicalContact::computeQpProperties()
     117             : {
     118           0 : }
     119             : 
     120             : void
     121           0 : ComputeWeightedGapLMMechanicalContact::computeQpIProperties()
     122             : {
     123           0 : }
     124             : 
     125             : void
     126       31178 : ComputeWeightedGapLMMechanicalContact::residualSetup()
     127             : {
     128       31178 : }
     129             : 
     130             : void
     131       15722 : ComputeWeightedGapLMMechanicalContact::jacobianSetup()
     132             : {
     133       15722 : }
     134             : 
     135             : void
     136     5812655 : ComputeWeightedGapLMMechanicalContact::computeResidual(const Moose::MortarType /*mortar_type*/)
     137             : {
     138     5812655 : }
     139             : 
     140             : void
     141     1043372 : ComputeWeightedGapLMMechanicalContact::computeJacobian(const Moose::MortarType mortar_type)
     142             : {
     143             :   // During "computeResidual" and "computeJacobian" we are actually just computing properties on the
     144             :   // mortar segment element mesh. We are *not* actually assembling into the residual/Jacobian. For
     145             :   // the zero-penetration constraint, the property of interest is the map from node to weighted gap.
     146             :   // Computation of the properties proceeds identically for residual and Jacobian evaluation hence
     147             :   // why we simply call computeResidual here. We will assemble into the residual/Jacobian later from
     148             :   // the post() method
     149     1043372 :   computeResidual(mortar_type);
     150     1043372 : }
     151             : 
     152             : void
     153       17754 : ComputeWeightedGapLMMechanicalContact::post()
     154             : {
     155             :   parallel_object_only();
     156             : 
     157       17754 :   const auto & dof_to_weighted_gap = _weighted_gap_uo.dofToWeightedGap();
     158             : 
     159       92902 :   for (const auto & [dof_object, weighted_gap_pr] : dof_to_weighted_gap)
     160             :   {
     161       75148 :     if (dof_object->processor_id() != this->processor_id())
     162         421 :       continue;
     163             : 
     164             :     const auto & [weighted_gap, normalization] = weighted_gap_pr;
     165       74727 :     _weighted_gap_ptr = &weighted_gap;
     166       74727 :     _normalization_ptr = &normalization;
     167             : 
     168       74727 :     enforceConstraintOnDof(dof_object);
     169             :   }
     170       17754 : }
     171             : 
     172             : void
     173       20226 : ComputeWeightedGapLMMechanicalContact::incorrectEdgeDroppingPost(
     174             :     const std::unordered_set<const Node *> & inactive_lm_nodes)
     175             : {
     176       20226 :   const auto & dof_to_weighted_gap = _weighted_gap_uo.dofToWeightedGap();
     177             : 
     178      159015 :   for (const auto & [dof_object, weighted_gap_pr] : dof_to_weighted_gap)
     179             :   {
     180      274862 :     if ((inactive_lm_nodes.find(static_cast<const Node *>(dof_object)) !=
     181      138789 :          inactive_lm_nodes.end()) ||
     182      138789 :         (dof_object->processor_id() != this->processor_id()))
     183        2716 :       continue;
     184             : 
     185             :     const auto & [weighted_gap, normalization] = weighted_gap_pr;
     186      136073 :     _weighted_gap_ptr = &weighted_gap;
     187      136073 :     _normalization_ptr = &normalization;
     188             : 
     189      136073 :     enforceConstraintOnDof(dof_object);
     190             :   }
     191       20226 : }
     192             : 
     193             : void
     194      386402 : ComputeWeightedGapLMMechanicalContact::enforceConstraintOnDof(const DofObject * const dof)
     195             : {
     196      386402 :   const auto & weighted_gap = *_weighted_gap_ptr;
     197      386402 :   const Real c = _normalize_c ? _c / *_normalization_ptr : _c;
     198             : 
     199      386402 :   const auto dof_index = dof->dof_number(_sys.number(), _var->number(), 0);
     200      386402 :   ADReal lm_value = (*_sys.currentSolution())(dof_index);
     201             :   Moose::derivInsert(lm_value.derivatives(), dof_index, 1.);
     202             : 
     203      386402 :   const ADReal dof_residual = std::min(lm_value, weighted_gap * c);
     204             : 
     205      386402 :   addResidualsAndJacobian(_assembly,
     206      772804 :                           std::array<ADReal, 1>{{dof_residual}},
     207      386402 :                           std::array<dof_id_type, 1>{{dof_index}},
     208      386402 :                           _var->scalingFactor());
     209      386402 : }