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/dualsemidynamicsparsenumberarray.h"
      16             : #include "metaphysicl/parallel_dualnumber.h"
      17             : #include "metaphysicl/parallel_dynamic_std_array_wrapper.h"
      18             : #include "metaphysicl/parallel_semidynamicsparsenumberarray.h"
      19             : #include "timpi/parallel_sync.h"
      20             : 
      21             : registerMooseObject("ContactApp", ComputeWeightedGapLMMechanicalContact);
      22             : 
      23             : namespace
      24             : {
      25             : const InputParameters &
      26         625 : assignVarsInParams(const InputParameters & params_in)
      27             : {
      28             :   InputParameters & ret = const_cast<InputParameters &>(params_in);
      29         625 :   const auto & disp_x_name = ret.get<std::vector<VariableName>>("disp_x");
      30         625 :   if (disp_x_name.size() != 1)
      31           0 :     mooseError("We require that the disp_x parameter have exactly one coupled name");
      32             : 
      33             :   // We do this so we don't get any variable errors during MortarConstraint(Base) construction
      34        1250 :   ret.set<VariableName>("secondary_variable") = disp_x_name[0];
      35         625 :   ret.set<VariableName>("primary_variable") = disp_x_name[0];
      36             : 
      37         625 :   return ret;
      38             : }
      39             : }
      40             : 
      41             : InputParameters
      42        1250 : ComputeWeightedGapLMMechanicalContact::validParams()
      43             : {
      44        1250 :   InputParameters params = ADMortarConstraint::validParams();
      45        1250 :   params.addClassDescription("Computes the weighted gap that will later be used to enforce the "
      46             :                              "zero-penetration mechanical contact conditions");
      47        1250 :   params.suppressParameter<VariableName>("secondary_variable");
      48        1250 :   params.suppressParameter<VariableName>("primary_variable");
      49        2500 :   params.addRequiredCoupledVar("disp_x", "The x displacement variable");
      50        2500 :   params.addRequiredCoupledVar("disp_y", "The y displacement variable");
      51        2500 :   params.addCoupledVar("disp_z", "The z displacement variable");
      52        2500 :   params.addParam<Real>(
      53        2500 :       "c", 1e6, "Parameter for balancing the size of the gap and contact pressure");
      54        2500 :   params.addParam<bool>(
      55             :       "normalize_c",
      56        2500 :       false,
      57             :       "Whether to normalize c by weighting function norm. When unnormalized "
      58             :       "the value of c effectively depends on element size since in the constraint we compare nodal "
      59             :       "Lagrange Multiplier values to integrated gap values (LM nodal value is independent of "
      60             :       "element size, where integrated values are dependent on element size).");
      61        1250 :   params.set<bool>("use_displaced_mesh") = true;
      62        1250 :   params.set<bool>("interpolate_normals") = false;
      63        2500 :   params.addRequiredParam<UserObjectName>("weighted_gap_uo", "The weighted gap user object");
      64        1250 :   return params;
      65           0 : }
      66             : 
      67         625 : ComputeWeightedGapLMMechanicalContact::ComputeWeightedGapLMMechanicalContact(
      68         625 :     const InputParameters & parameters)
      69             :   : ADMortarConstraint(assignVarsInParams(parameters)),
      70         625 :     _secondary_disp_x(adCoupledValue("disp_x")),
      71         625 :     _primary_disp_x(adCoupledNeighborValue("disp_x")),
      72         625 :     _secondary_disp_y(adCoupledValue("disp_y")),
      73         625 :     _primary_disp_y(adCoupledNeighborValue("disp_y")),
      74         625 :     _has_disp_z(isCoupled("disp_z")),
      75         625 :     _secondary_disp_z(_has_disp_z ? &adCoupledValue("disp_z") : nullptr),
      76         625 :     _primary_disp_z(_has_disp_z ? &adCoupledNeighborValue("disp_z") : nullptr),
      77        1250 :     _c(getParam<Real>("c")),
      78        1250 :     _normalize_c(getParam<bool>("normalize_c")),
      79         625 :     _nodal(getVar("disp_x", 0)->feType().family == LAGRANGE),
      80         625 :     _disp_x_var(getVar("disp_x", 0)),
      81         625 :     _disp_y_var(getVar("disp_y", 0)),
      82         767 :     _disp_z_var(_has_disp_z ? getVar("disp_z", 0) : nullptr),
      83        1250 :     _weighted_gap_uo(getUserObject<WeightedGapUserObject>("weighted_gap_uo"))
      84             : {
      85        1250 :   if (!getParam<bool>("use_displaced_mesh"))
      86           0 :     paramError(
      87             :         "use_displaced_mesh",
      88             :         "'use_displaced_mesh' must be true for the ComputeWeightedGapLMMechanicalContact object");
      89             : 
      90         625 :   if (!_var->isNodal())
      91           0 :     if (_var->feType().order != static_cast<Order>(0))
      92           0 :       mooseError("Normal contact constraints only support elemental variables of CONSTANT order");
      93         625 : }
      94             : 
      95           0 : ADReal ComputeWeightedGapLMMechanicalContact::computeQpResidual(Moose::MortarType)
      96             : {
      97           0 :   mooseError("We should never call computeQpResidual for ComputeWeightedGapLMMechanicalContact");
      98             : }
      99             : 
     100             : void
     101           0 : ComputeWeightedGapLMMechanicalContact::computeQpProperties()
     102             : {
     103           0 : }
     104             : 
     105             : void
     106           0 : ComputeWeightedGapLMMechanicalContact::computeQpIProperties()
     107             : {
     108           0 : }
     109             : 
     110             : void
     111       38015 : ComputeWeightedGapLMMechanicalContact::residualSetup()
     112             : {
     113       38015 : }
     114             : 
     115             : void
     116       22517 : ComputeWeightedGapLMMechanicalContact::jacobianSetup()
     117             : {
     118       22517 : }
     119             : 
     120             : void
     121     4358560 : ComputeWeightedGapLMMechanicalContact::computeResidual(const Moose::MortarType /*mortar_type*/)
     122             : {
     123     4358560 : }
     124             : 
     125             : void
     126      764598 : ComputeWeightedGapLMMechanicalContact::computeJacobian(const Moose::MortarType mortar_type)
     127             : {
     128             :   // During "computeResidual" and "computeJacobian" we are actually just computing properties on the
     129             :   // mortar segment element mesh. We are *not* actually assembling into the residual/Jacobian. For
     130             :   // the zero-penetration constraint, the property of interest is the map from node to weighted gap.
     131             :   // Computation of the properties proceeds identically for residual and Jacobian evaluation hence
     132             :   // why we simply call computeResidual here. We will assemble into the residual/Jacobian later from
     133             :   // the post() method
     134      764598 :   computeResidual(mortar_type);
     135      764598 : }
     136             : 
     137             : void
     138       19714 : ComputeWeightedGapLMMechanicalContact::post()
     139             : {
     140             :   parallel_object_only();
     141             : 
     142       19714 :   const auto & dof_to_weighted_gap = _weighted_gap_uo.dofToWeightedGap();
     143             : 
     144       88114 :   for (const auto & [dof_object, weighted_gap_pr] : dof_to_weighted_gap)
     145             :   {
     146       68400 :     if (dof_object->processor_id() != this->processor_id())
     147         741 :       continue;
     148             : 
     149             :     const auto & [weighted_gap, normalization] = weighted_gap_pr;
     150       67659 :     _weighted_gap_ptr = &weighted_gap;
     151       67659 :     _normalization_ptr = &normalization;
     152             : 
     153       67659 :     enforceConstraintOnDof(dof_object);
     154             :   }
     155       19714 : }
     156             : 
     157             : void
     158       26113 : ComputeWeightedGapLMMechanicalContact::incorrectEdgeDroppingPost(
     159             :     const std::unordered_set<const Node *> & inactive_lm_nodes)
     160             : {
     161       26113 :   const auto & dof_to_weighted_gap = _weighted_gap_uo.dofToWeightedGap();
     162             : 
     163      137328 :   for (const auto & [dof_object, weighted_gap_pr] : dof_to_weighted_gap)
     164             :   {
     165      218074 :     if ((inactive_lm_nodes.find(static_cast<const Node *>(dof_object)) !=
     166      111215 :          inactive_lm_nodes.end()) ||
     167             :         (dof_object->processor_id() != this->processor_id()))
     168        4356 :       continue;
     169             : 
     170             :     const auto & [weighted_gap, normalization] = weighted_gap_pr;
     171      106859 :     _weighted_gap_ptr = &weighted_gap;
     172      106859 :     _normalization_ptr = &normalization;
     173             : 
     174      106859 :     enforceConstraintOnDof(dof_object);
     175             :   }
     176       26113 : }
     177             : 
     178             : void
     179      404897 : ComputeWeightedGapLMMechanicalContact::enforceConstraintOnDof(const DofObject * const dof)
     180             : {
     181      404897 :   const auto & weighted_gap = *_weighted_gap_ptr;
     182      404897 :   const Real c = _normalize_c ? _c / *_normalization_ptr : _c;
     183             : 
     184      404897 :   const auto dof_index = dof->dof_number(_sys.number(), _var->number(), 0);
     185      404897 :   ADReal lm_value = (*_sys.currentSolution())(dof_index);
     186             :   Moose::derivInsert(lm_value.derivatives(), dof_index, 1.);
     187             : 
     188      404897 :   const ADReal dof_residual = std::min(lm_value, weighted_gap * c);
     189             : 
     190      404897 :   addResidualsAndJacobian(_assembly,
     191      404897 :                           std::array<ADReal, 1>{{dof_residual}},
     192      809794 :                           std::array<dof_id_type, 1>{{dof_index}},
     193      404897 :                           _var->scalingFactor());
     194      404897 : }