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         730 : assignVarsInParams(const InputParameters & params_in)
      32             : {
      33             :   InputParameters & ret = const_cast<InputParameters &>(params_in);
      34         730 :   const auto & disp_x_name = ret.get<std::vector<VariableName>>("disp_x");
      35         730 :   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        1460 :   ret.set<VariableName>("secondary_variable") = disp_x_name[0];
      40         730 :   ret.set<VariableName>("primary_variable") = disp_x_name[0];
      41             : 
      42         730 :   return ret;
      43             : }
      44             : }
      45             : 
      46             : InputParameters
      47        1460 : ComputeWeightedGapLMMechanicalContact::validParams()
      48             : {
      49        1460 :   InputParameters params = ADMortarConstraint::validParams();
      50        1460 :   params.addClassDescription("Computes the weighted gap that will later be used to enforce the "
      51             :                              "zero-penetration mechanical contact conditions");
      52        1460 :   params.suppressParameter<VariableName>("secondary_variable");
      53        1460 :   params.suppressParameter<VariableName>("primary_variable");
      54        2920 :   params.addRequiredCoupledVar("disp_x", "The x displacement variable");
      55        2920 :   params.addRequiredCoupledVar("disp_y", "The y displacement variable");
      56        2920 :   params.addCoupledVar("disp_z", "The z displacement variable");
      57        2920 :   params.addParam<Real>(
      58        2920 :       "c", 1e6, "Parameter for balancing the size of the gap and contact pressure");
      59        2920 :   params.addParam<bool>(
      60             :       "normalize_c",
      61        2920 :       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        1460 :   params.set<bool>("use_displaced_mesh") = true;
      67        1460 :   params.set<bool>("interpolate_normals") = false;
      68        2920 :   params.addRequiredParam<UserObjectName>("weighted_gap_uo", "The weighted gap user object");
      69        1460 :   return params;
      70           0 : }
      71             : 
      72         730 : ComputeWeightedGapLMMechanicalContact::ComputeWeightedGapLMMechanicalContact(
      73         730 :     const InputParameters & parameters)
      74             :   : ADMortarConstraint(assignVarsInParams(parameters)),
      75         730 :     _secondary_disp_x(adCoupledValue("disp_x")),
      76         730 :     _primary_disp_x(adCoupledNeighborValue("disp_x")),
      77         730 :     _secondary_disp_y(adCoupledValue("disp_y")),
      78         730 :     _primary_disp_y(adCoupledNeighborValue("disp_y")),
      79         730 :     _has_disp_z(isCoupled("disp_z")),
      80         730 :     _secondary_disp_z(_has_disp_z ? &adCoupledValue("disp_z") : nullptr),
      81         730 :     _primary_disp_z(_has_disp_z ? &adCoupledNeighborValue("disp_z") : nullptr),
      82        1460 :     _c(getParam<Real>("c")),
      83        1460 :     _normalize_c(getParam<bool>("normalize_c")),
      84         730 :     _nodal(getVar("disp_x", 0)->feType().family == LAGRANGE),
      85         730 :     _disp_x_var(getVar("disp_x", 0)),
      86         730 :     _disp_y_var(getVar("disp_y", 0)),
      87         890 :     _disp_z_var(_has_disp_z ? getVar("disp_z", 0) : nullptr),
      88        1460 :     _weighted_gap_uo(getUserObject<WeightedGapUserObject>("weighted_gap_uo"))
      89             : {
      90        1460 :   if (!getParam<bool>("use_displaced_mesh"))
      91           0 :     paramError(
      92             :         "use_displaced_mesh",
      93             :         "'use_displaced_mesh' must be true for the ComputeWeightedGapLMMechanicalContact object");
      94             : 
      95         730 :   if (!_var->isNodal())
      96           0 :     if (_var->feType().order != static_cast<Order>(0))
      97           0 :       mooseError("Normal contact constraints only support elemental variables of CONSTANT order");
      98         730 : }
      99             : 
     100             : ADReal
     101           0 : ComputeWeightedGapLMMechanicalContact::computeQpResidual(Moose::MortarType)
     102             : {
     103           0 :   mooseError("We should never call computeQpResidual for ComputeWeightedGapLMMechanicalContact");
     104             : }
     105             : 
     106             : void
     107           0 : ComputeWeightedGapLMMechanicalContact::computeQpProperties()
     108             : {
     109           0 : }
     110             : 
     111             : void
     112           0 : ComputeWeightedGapLMMechanicalContact::computeQpIProperties()
     113             : {
     114           0 : }
     115             : 
     116             : void
     117       43151 : ComputeWeightedGapLMMechanicalContact::residualSetup()
     118             : {
     119       43151 : }
     120             : 
     121             : void
     122       24607 : ComputeWeightedGapLMMechanicalContact::jacobianSetup()
     123             : {
     124       24607 : }
     125             : 
     126             : void
     127     5742757 : ComputeWeightedGapLMMechanicalContact::computeResidual(const Moose::MortarType /*mortar_type*/)
     128             : {
     129     5742757 : }
     130             : 
     131             : void
     132     1011343 : ComputeWeightedGapLMMechanicalContact::computeJacobian(const Moose::MortarType mortar_type)
     133             : {
     134             :   // During "computeResidual" and "computeJacobian" we are actually just computing properties on the
     135             :   // mortar segment element mesh. We are *not* actually assembling into the residual/Jacobian. For
     136             :   // the zero-penetration constraint, the property of interest is the map from node to weighted gap.
     137             :   // Computation of the properties proceeds identically for residual and Jacobian evaluation hence
     138             :   // why we simply call computeResidual here. We will assemble into the residual/Jacobian later from
     139             :   // the post() method
     140     1011343 :   computeResidual(mortar_type);
     141     1011343 : }
     142             : 
     143             : void
     144       23375 : ComputeWeightedGapLMMechanicalContact::post()
     145             : {
     146             :   parallel_object_only();
     147             : 
     148       23375 :   const auto & dof_to_weighted_gap = _weighted_gap_uo.dofToWeightedGap();
     149             : 
     150      113832 :   for (const auto & [dof_object, weighted_gap_pr] : dof_to_weighted_gap)
     151             :   {
     152       90457 :     if (dof_object->processor_id() != this->processor_id())
     153         741 :       continue;
     154             : 
     155             :     const auto & [weighted_gap, normalization] = weighted_gap_pr;
     156       89716 :     _weighted_gap_ptr = &weighted_gap;
     157       89716 :     _normalization_ptr = &normalization;
     158             : 
     159       89716 :     enforceConstraintOnDof(dof_object);
     160             :   }
     161       23375 : }
     162             : 
     163             : void
     164       28786 : ComputeWeightedGapLMMechanicalContact::incorrectEdgeDroppingPost(
     165             :     const std::unordered_set<const Node *> & inactive_lm_nodes)
     166             : {
     167       28786 :   const auto & dof_to_weighted_gap = _weighted_gap_uo.dofToWeightedGap();
     168             : 
     169      164330 :   for (const auto & [dof_object, weighted_gap_pr] : dof_to_weighted_gap)
     170             :   {
     171      266732 :     if ((inactive_lm_nodes.find(static_cast<const Node *>(dof_object)) !=
     172      135544 :          inactive_lm_nodes.end()) ||
     173      135544 :         (dof_object->processor_id() != this->processor_id()))
     174        4356 :       continue;
     175             : 
     176             :     const auto & [weighted_gap, normalization] = weighted_gap_pr;
     177      131188 :     _weighted_gap_ptr = &weighted_gap;
     178      131188 :     _normalization_ptr = &normalization;
     179             : 
     180      131188 :     enforceConstraintOnDof(dof_object);
     181             :   }
     182       28786 : }
     183             : 
     184             : void
     185      480087 : ComputeWeightedGapLMMechanicalContact::enforceConstraintOnDof(const DofObject * const dof)
     186             : {
     187      480087 :   const auto & weighted_gap = *_weighted_gap_ptr;
     188      480087 :   const Real c = _normalize_c ? _c / *_normalization_ptr : _c;
     189             : 
     190      480087 :   const auto dof_index = dof->dof_number(_sys.number(), _var->number(), 0);
     191      480087 :   ADReal lm_value = (*_sys.currentSolution())(dof_index);
     192             :   Moose::derivInsert(lm_value.derivatives(), dof_index, 1.);
     193             : 
     194      480087 :   const ADReal dof_residual = std::min(lm_value, weighted_gap * c);
     195             : 
     196      480087 :   addResidualsAndJacobian(_assembly,
     197      960174 :                           std::array<ADReal, 1>{{dof_residual}},
     198      480087 :                           std::array<dof_id_type, 1>{{dof_index}},
     199      480087 :                           _var->scalingFactor());
     200      480087 : }