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 "HomogenizedTotalLagrangianStressDivergence.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "Function.h"
      14             : #include "MooseVariableScalar.h"
      15             : 
      16             : registerMooseObject("SolidMechanicsApp", HomogenizedTotalLagrangianStressDivergence);
      17             : 
      18             : InputParameters
      19         624 : HomogenizedTotalLagrangianStressDivergence::validParams()
      20             : {
      21         624 :   InputParameters params = HomogenizationInterface<TotalLagrangianStressDivergence>::validParams();
      22         624 :   params.addClassDescription("Total Lagrangian stress equilibrium kernel with "
      23             :                              "homogenization constraint Jacobian terms");
      24        1248 :   params.renameCoupledVar(
      25             :       "scalar_variable", "macro_var", "Optional scalar field with the macro gradient");
      26             : 
      27        1248 :   params.addParam<bool>(
      28        1248 :       "off_diagonal_jacobian", true, "Include the off diagonal parts of the constraint Jacobian");
      29             : 
      30         624 :   return params;
      31           0 : }
      32             : 
      33         312 : HomogenizedTotalLagrangianStressDivergence::HomogenizedTotalLagrangianStressDivergence(
      34         312 :     const InputParameters & parameters)
      35             :   : HomogenizationInterface<TotalLagrangianStressDivergence>(parameters),
      36         624 :     _use_off_diagonal(getParam<bool>("off_diagonal_jacobian"))
      37             : {
      38         312 : }
      39             : 
      40             : std::set<std::string>
      41         312 : HomogenizedTotalLagrangianStressDivergence::additionalROVariables()
      42             : {
      43             :   // Add the scalar variable to the list of variables that this kernel contributes to
      44             :   std::set<std::string> vars = TotalLagrangianStressDivergence::additionalROVariables();
      45         312 :   vars.insert(_kappa_var_ptr->name());
      46         312 :   return vars;
      47             : }
      48             : 
      49             : void
      50     2376640 : HomogenizedTotalLagrangianStressDivergence::computeScalarResidual()
      51             : {
      52     2376640 :   if (_alpha != 0)
      53     1574157 :     return;
      54             : 
      55      802483 :   std::vector<Real> scalar_residuals(_k_order);
      56             : 
      57             :   // only assemble scalar residual once; i.e. when handling the first displacement component
      58     7104231 :   for (_qp = 0; _qp < _qrule->n_points(); _qp++)
      59             :   {
      60     6301748 :     initScalarQpResidual();
      61     6301748 :     const auto dV = _JxW[_qp] * _coord[_qp];
      62             : 
      63             :     // index for residual vector
      64             :     unsigned int h = 0;
      65             : 
      66    52122856 :     for (const auto & [indices, constraint] : cmap())
      67             :     {
      68    45821108 :       const auto [i, j] = indices;
      69    45821108 :       const auto [ctype, ctarget] = constraint;
      70    45821108 :       const auto cval = ctarget->value(_t, _q_point[_qp]);
      71             : 
      72             :       // value to be constrained
      73             :       Real val;
      74    45821108 :       if (_large_kinematics)
      75             :       {
      76    25179316 :         if (ctype == Homogenization::ConstraintType::Stress)
      77     8266256 :           val = _pk1[_qp](i, j);
      78    16913060 :         else if (ctype == Homogenization::ConstraintType::Strain)
      79    29513748 :           val = _F[_qp](i, j) - (Real(i == j));
      80             :         else
      81           0 :           mooseError("Unknown constraint type in the integral!");
      82             :       }
      83             :       else
      84             :       {
      85    20641792 :         if (ctype == Homogenization::ConstraintType::Stress)
      86    10287360 :           val = _pk1[_qp](i, j);
      87    10354432 :         else if (ctype == Homogenization::ConstraintType::Strain)
      88    15506688 :           val = 0.5 * (_F[_qp](i, j) + _F[_qp](j, i)) - (Real(i == j));
      89             :         else
      90           0 :           mooseError("Unknown constraint type in the integral!");
      91             :       }
      92             : 
      93    45821108 :       scalar_residuals[h++] += (val - cval) * dV;
      94             :     }
      95             :   }
      96             : 
      97      802483 :   addResiduals(
      98      802483 :       _assembly, scalar_residuals, _kappa_var_ptr->dofIndices(), _kappa_var_ptr->scalingFactor());
      99      802483 : }
     100             : 
     101             : void
     102      166892 : HomogenizedTotalLagrangianStressDivergence::computeScalarJacobian()
     103             : {
     104      166892 :   if (_alpha != 0)
     105             :     return;
     106             : 
     107       62270 :   _local_ke.resize(_k_order, _k_order);
     108             : 
     109             :   // only assemble scalar residual once; i.e. when handling the first displacement component
     110      483798 :   for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     111             :   {
     112      421528 :     initScalarQpJacobian(_kappa_var);
     113      421528 :     const auto dV = _JxW[_qp] * _coord[_qp];
     114             : 
     115             :     // index for Jacobian row
     116             :     unsigned int h = 0;
     117             : 
     118     3532592 :     for (const auto & [indices1, constraint1] : cmap())
     119             :     {
     120     3111064 :       const auto [i, j] = indices1;
     121     3111064 :       const auto ctype = constraint1.first;
     122             : 
     123             :       // index for Jacobian col
     124             :       unsigned int m = 0;
     125             : 
     126    28262192 :       for (const auto & [indices2, constraint2] : cmap())
     127             :       {
     128    25151128 :         const auto [k, l] = indices2;
     129    25151128 :         if (ctype == Homogenization::ConstraintType::Stress)
     130     8114672 :           _local_ke(h, m++) += dV * (_dpk1[_qp](i, j, k, l));
     131    17036456 :         else if (ctype == Homogenization::ConstraintType::Strain)
     132             :         {
     133    17036456 :           if (_large_kinematics)
     134    28704688 :             _local_ke(h, m++) += dV * (Real(i == k && j == l));
     135             :           else
     136     4914168 :             _local_ke(h, m++) += dV * (0.5 * Real(i == k && j == l) + 0.5 * Real(i == l && j == k));
     137             :         }
     138             :         else
     139           0 :           mooseError("Unknown constraint type in Jacobian calculator!");
     140             :       }
     141     3111064 :       h++;
     142             :     }
     143             :   }
     144             : 
     145       62270 :   addJacobian(_assembly,
     146             :               _local_ke,
     147       62270 :               _kappa_var_ptr->dofIndices(),
     148       62270 :               _kappa_var_ptr->dofIndices(),
     149       62270 :               _kappa_var_ptr->scalingFactor());
     150             : }
     151             : 
     152             : void
     153      463784 : HomogenizedTotalLagrangianStressDivergence::computeScalarOffDiagJacobian(
     154             :     const unsigned int jvar_num)
     155             : {
     156      463784 :   if (!_use_off_diagonal)
     157             :     return;
     158             : 
     159             :   // ONLY assemble the contribution from _alpha component, which is connected with _var
     160             :   // The other components are handled by other kernel instances with other _alpha
     161      380840 :   if (jvar_num != _var.number())
     162             :     return;
     163             : 
     164             :   const auto & jvar = getVariable(jvar_num);
     165      138950 :   const auto jvar_size = jvar.phiSize();
     166      138950 :   _local_ke.resize(_k_order, jvar_size);
     167             : 
     168     1108062 :   for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     169             :   {
     170      969112 :     const auto dV = _JxW[_qp] * _coord[_qp];
     171             : 
     172             :     // index for Jacobian row
     173             :     unsigned int h = 0;
     174             : 
     175     8304944 :     for (const auto & [indices, constraint] : cmap())
     176             :     {
     177     7335832 :       std::tie(_m, _n) = indices;
     178     7335832 :       _ctype = constraint.first;
     179     7335832 :       initScalarQpOffDiagJacobian(jvar);
     180    64239096 :       for (_j = 0; _j < jvar_size; _j++)
     181    56903264 :         _local_ke(h, _j) += dV * computeScalarQpOffDiagJacobian(jvar_num);
     182     7335832 :       h++;
     183             :     }
     184             :   }
     185             : 
     186      138950 :   addJacobian(_assembly,
     187             :               _local_ke,
     188      138950 :               _kappa_var_ptr->dofIndices(),
     189      138950 :               jvar.dofIndices(),
     190      138950 :               _kappa_var_ptr->scalingFactor());
     191             : }
     192             : 
     193             : void
     194      166598 : HomogenizedTotalLagrangianStressDivergence::computeOffDiagJacobianScalarLocal(
     195             :     const unsigned int svar_num)
     196             : {
     197      166598 :   if (!_use_off_diagonal)
     198             :     return;
     199             : 
     200             :   // Just in case, skip any other scalar variables
     201      138950 :   if (svar_num != _kappa_var)
     202             :     return;
     203             : 
     204      138950 :   _local_ke.resize(_test.size(), _k_order);
     205             : 
     206     1108062 :   for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     207             :   {
     208             :     unsigned int l = 0;
     209      969112 :     const auto dV = _JxW[_qp] * _coord[_qp];
     210     8304944 :     for (const auto & [indices, constraint] : cmap())
     211             :     {
     212             :       // copy constraint indices to protected variables to pass to Qp routine
     213     7335832 :       std::tie(_m, _n) = indices;
     214     7335832 :       _ctype = constraint.first;
     215     7335832 :       initScalarQpJacobian(svar_num);
     216    64239096 :       for (_i = 0; _i < _test.size(); _i++)
     217    56903264 :         _local_ke(_i, l) += dV * computeQpOffDiagJacobianScalar(svar_num);
     218     7335832 :       l++;
     219             :     }
     220             :   }
     221             : 
     222      138950 :   addJacobian(
     223      138950 :       _assembly, _local_ke, _var.dofIndices(), _kappa_var_ptr->dofIndices(), _var.scalingFactor());
     224             : }
     225             : 
     226             : Real
     227    56903264 : HomogenizedTotalLagrangianStressDivergence::computeQpOffDiagJacobianScalar(
     228             :     unsigned int /*svar_num*/)
     229             : {
     230    56903264 :   return _dpk1[_qp].contractionKl(_m, _n, gradTest(_alpha));
     231             : }
     232             : 
     233             : Real
     234    56903264 : HomogenizedTotalLagrangianStressDivergence::computeScalarQpOffDiagJacobian(
     235             :     unsigned int /*jvar_num*/)
     236             : {
     237    56903264 :   if (_ctype == Homogenization::ConstraintType::Stress)
     238    23587328 :     return _dpk1[_qp].contractionIj(_m, _n, gradTrial(_alpha));
     239    33315936 :   else if (_ctype == Homogenization::ConstraintType::Strain)
     240    33315936 :     if (_large_kinematics)
     241    49461728 :       return Real(_m == _alpha) * gradTrial(_alpha)(_m, _n);
     242             :     else
     243     9460224 :       return 0.5 * (Real(_m == _alpha) * gradTrial(_alpha)(_m, _n) +
     244     5899520 :                     Real(_n == _alpha) * gradTrial(_alpha)(_n, _m));
     245             :   else
     246           0 :     mooseError("Unknown constraint type in kernel calculation!");
     247             : }