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 "WeightedGapUserObject.h"
      11             : #include "MooseVariableField.h"
      12             : #include "SubProblem.h"
      13             : #include "MortarUtils.h"
      14             : #include "MooseUtils.h"
      15             : #include "MortarContactUtils.h"
      16             : #include "AutomaticMortarGeneration.h"
      17             : 
      18             : #include "libmesh/quadrature.h"
      19             : 
      20             : #include <limits>
      21             : 
      22             : InputParameters
      23        2730 : WeightedGapUserObject::validParams()
      24             : {
      25        2730 :   InputParameters params = MortarUserObject::validParams();
      26        2730 :   params += MortarConsumerInterface::validParams();
      27        2730 :   params += TwoMaterialPropertyInterface::validParams();
      28        5460 :   params.addRequiredCoupledVar("disp_x", "The x displacement variable");
      29        5460 :   params.addRequiredCoupledVar("disp_y", "The y displacement variable");
      30        5460 :   params.addCoupledVar("disp_z", "The z displacement variable");
      31        2730 :   params.set<bool>("use_displaced_mesh") = true;
      32        2730 :   params.set<bool>("interpolate_normals") = false;
      33       10920 :   params.set<ExecFlagEnum>("execute_on") = {EXEC_LINEAR, EXEC_NONLINEAR};
      34        2730 :   params.suppressParameter<ExecFlagEnum>("execute_on");
      35        2730 :   return params;
      36        2730 : }
      37             : 
      38        1207 : WeightedGapUserObject::WeightedGapUserObject(const InputParameters & parameters)
      39             :   : MortarUserObject(parameters),
      40        2414 :     _fe_problem(*getCheckedPointerParam<FEProblemBase *>("_fe_problem_base")),
      41        1207 :     _nodal(getVar("disp_x", 0)->feType().family == LAGRANGE),
      42        1207 :     _disp_x_var(getVar("disp_x", 0)),
      43        1207 :     _disp_y_var(getVar("disp_y", 0)),
      44        1207 :     _has_disp_z(isCoupled("disp_z")),
      45        1207 :     _disp_z_var(_has_disp_z ? getVar("disp_z", 0) : nullptr),
      46        1207 :     _secondary_disp_x(_disp_x_var->adSln()),
      47        1207 :     _primary_disp_x(_disp_x_var->adSlnNeighbor()),
      48        1207 :     _secondary_disp_y(_disp_y_var->adSln()),
      49        1207 :     _primary_disp_y(_disp_y_var->adSlnNeighbor()),
      50        1207 :     _secondary_disp_z(_has_disp_z ? &_disp_z_var->adSln() : nullptr),
      51        1207 :     _primary_disp_z(_has_disp_z ? &_disp_z_var->adSlnNeighbor() : nullptr),
      52        2414 :     _coord(_assembly.mortarCoordTransformation())
      53             : {
      54        2414 :   if (!getParam<bool>("use_displaced_mesh"))
      55           0 :     paramError("use_displaced_mesh",
      56             :                "'use_displaced_mesh' must be true for the WeightedGapUserObject object");
      57        1207 : }
      58             : 
      59             : void
      60        1191 : WeightedGapUserObject::initialSetup()
      61             : {
      62        1191 :   MortarUserObject::initialSetup();
      63        1191 :   _test = &test();
      64        1191 : }
      65             : 
      66             : void
      67     8733668 : WeightedGapUserObject::computeQpProperties()
      68             : {
      69             :   // Trim interior node variable derivatives
      70             :   const auto & primary_ip_lowerd_map = amg().getPrimaryIpToLowerElementMap(
      71     8733668 :       *_lower_primary_elem, *_lower_primary_elem->interior_parent(), *_lower_secondary_elem);
      72             :   const auto & secondary_ip_lowerd_map =
      73     8733668 :       amg().getSecondaryIpToLowerElementMap(*_lower_secondary_elem);
      74             : 
      75     8733668 :   std::array<const MooseVariable *, 3> var_array{{_disp_x_var, _disp_y_var, _disp_z_var}};
      76             :   std::array<ADReal, 3> primary_disp{
      77     8733668 :       {_primary_disp_x[_qp], _primary_disp_y[_qp], _has_disp_z ? (*_primary_disp_z)[_qp] : 0}};
      78     8733668 :   std::array<ADReal, 3> secondary_disp{{_secondary_disp_x[_qp],
      79     8733668 :                                         _secondary_disp_y[_qp],
      80     8733668 :                                         _has_disp_z ? (*_secondary_disp_z)[_qp] : 0}};
      81             : 
      82     8733668 :   trimInteriorNodeDerivatives(primary_ip_lowerd_map, var_array, primary_disp, false);
      83     8733668 :   trimInteriorNodeDerivatives(secondary_ip_lowerd_map, var_array, secondary_disp, true);
      84             : 
      85             :   const ADReal & prim_x = primary_disp[0];
      86             :   const ADReal & prim_y = primary_disp[1];
      87             :   const ADReal * prim_z = nullptr;
      88     8733668 :   if (_has_disp_z)
      89             :     prim_z = &primary_disp[2];
      90             : 
      91             :   const ADReal & sec_x = secondary_disp[0];
      92             :   const ADReal & sec_y = secondary_disp[1];
      93             :   const ADReal * sec_z = nullptr;
      94     8733668 :   if (_has_disp_z)
      95             :     sec_z = &secondary_disp[2];
      96             : 
      97             :   // Compute gap vector
      98     8733668 :   ADRealVectorValue gap_vec = _phys_points_primary[_qp] - _phys_points_secondary[_qp];
      99             : 
     100             :   // Generic displacement for interface problems
     101     8733668 :   _qp_displacement_nodal(0) = prim_x - sec_x;
     102     8733668 :   _qp_displacement_nodal(1) = prim_y - sec_y;
     103     8733668 :   if (_has_disp_z)
     104     7287760 :     _qp_displacement_nodal(2) = *prim_z - *sec_z;
     105             : 
     106     8733668 :   _qp_displacement_nodal *= _JxW_msm[_qp] * _coord[_qp];
     107             : 
     108     8733668 :   gap_vec(0).derivatives() = prim_x.derivatives() - sec_x.derivatives();
     109     8733668 :   gap_vec(1).derivatives() = prim_y.derivatives() - sec_y.derivatives();
     110     8733668 :   if (_has_disp_z)
     111     7287760 :     gap_vec(2).derivatives() = prim_z->derivatives() - sec_z->derivatives();
     112             : 
     113             :   // Compute integration point quantities: Normals (geometry) is averaged at the node, but not
     114             :   // interpolated within the weak integration.
     115     8733668 :   _qp_gap_nodal = gap_vec * (_JxW_msm[_qp] * _coord[_qp]);
     116             : 
     117             :   // To do normalization of constraint coefficient (c_n)
     118     8733668 :   _qp_factor = _JxW_msm[_qp] * _coord[_qp];
     119     8733668 : }
     120             : 
     121             : void
     122    32042856 : WeightedGapUserObject::computeQpIProperties()
     123             : {
     124             :   mooseAssert(_normals.size() == _lower_secondary_elem->n_nodes(),
     125             :               "Making sure that _normals is the expected size");
     126             : 
     127             :   // Get the _dof_to_weighted_gap map
     128    32042856 :   const auto * const dof = static_cast<const DofObject *>(_lower_secondary_elem->node_ptr(_i));
     129             : 
     130             :   auto & [weighted_gap, normalization] = _dof_to_weighted_gap[dof];
     131             : 
     132    32042856 :   weighted_gap += (*_test)[_i][_qp] * _qp_gap_nodal * _normals[_i];
     133    32042856 :   normalization += (*_test)[_i][_qp] * _qp_factor;
     134             : 
     135    64085712 :   _dof_to_weighted_displacements[dof] += (*_test)[_i][_qp] * _qp_displacement_nodal;
     136    32042856 : }
     137             : 
     138             : void
     139      176060 : WeightedGapUserObject::initialize()
     140             : {
     141             :   _dof_to_weighted_gap.clear();
     142             :   _dof_to_weighted_displacements.clear();
     143      176060 : }
     144             : 
     145             : void
     146      176060 : WeightedGapUserObject::finalize()
     147             : {
     148             :   // If the constraint is performed by the owner, then we don't need any data sent back; the owner
     149             :   // will take care of it. But if the constraint is not performed by the owner and we might have to
     150             :   // do some of the constraining ourselves, then we need data sent back to us
     151      176060 :   const bool send_data_back = !constrainedByOwner();
     152      176060 :   Moose::Mortar::Contact::communicateGaps(_dof_to_weighted_gap,
     153      176060 :                                           _subproblem.mesh(),
     154      176060 :                                           _nodal,
     155             :                                           /*normalize_c*/ true,
     156             :                                           _communicator,
     157             :                                           send_data_back);
     158      176060 : }
     159             : 
     160             : void
     161     2535982 : WeightedGapUserObject::execute()
     162             : {
     163    11269650 :   for (_qp = 0; _qp < _qrule_msm->n_points(); _qp++)
     164             :   {
     165     8733668 :     computeQpProperties();
     166    40776524 :     for (_i = 0; _i < _test->size(); ++_i)
     167    32042856 :       computeQpIProperties();
     168             :   }
     169     2535982 : }
     170             : 
     171             : Real
     172      171851 : WeightedGapUserObject::getNormalGap(const Node * const node) const
     173             : {
     174      171851 :   const auto it = _dof_to_weighted_gap.find(_subproblem.mesh().nodePtr(node->id()));
     175             : 
     176             :   // We are returning the physical weighted gap for analysis purposes
     177      171851 :   if (it != _dof_to_weighted_gap.end())
     178       21797 :     return physicalGap(it->second);
     179             :   else
     180             :     return 0.0;
     181             : }