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