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             : #ifdef NEML2_ENABLED
      11             : 
      12             : // NEML2 includes
      13             : #include "neml2/tensors/functions/stack.h"
      14             : #include "neml2/tensors/functions/sum.h"
      15             : 
      16             : // MOOSE includes
      17             : #include "NEML2StressDivergence.h"
      18             : 
      19             : registerMooseObject("SolidMechanicsApp", NEML2StressDivergence);
      20             : 
      21             : InputParameters
      22           0 : NEML2StressDivergence::validParams()
      23             : {
      24           0 :   InputParameters params = NEML2PostKernel::validParams();
      25           0 :   params.addClassDescription(
      26             :       "This user object calculates the stress divergence for a given set of displacement "
      27             :       "variables and assemble the contribution into the residual vector.");
      28           0 :   params.addRequiredParam<std::vector<NonlinearVariableName>>(
      29             :       "displacements",
      30             :       "The displacements variables whose function space will be used to define the test "
      31             :       "functions.");
      32           0 :   params.addRequiredParam<std::string>("stress",
      33             :                                        "The name of the NEML2 variable to use as the stress.");
      34           0 :   params.addRequiredParam<std::string>("residual", "The tag for the residual");
      35           0 :   return params;
      36           0 : }
      37             : 
      38           0 : NEML2StressDivergence::NEML2StressDivergence(const InputParameters & parameters)
      39             :   : NEML2PostKernel(parameters),
      40           0 :     _residual(dynamic_cast<PetscVector<Real> *>(&_sys.system().add_vector(
      41           0 :         getParam<std::string>("residual"), false, libMesh::ParallelType::GHOSTED))),
      42           0 :     _stress(_constitutive.getOutput(getParam<std::string>("stress"))),
      43           0 :     _disp_vars(getParam<std::vector<NonlinearVariableName>>("displacements")),
      44           0 :     _ndisp(_disp_vars.size())
      45             : {
      46             :   mooseAssert(_residual, "Failed to cast residual to PetscVector<Real>");
      47             : 
      48           0 :   if (_disp_vars.size() < 1 || _disp_vars.size() > 3)
      49           0 :     mooseError("NEML2StressDivergence requires 1 to 3 displacement variables, got ",
      50           0 :                _disp_vars.size());
      51             : 
      52           0 :   _grad_test_x = &_fe.getPhiGradient(_disp_vars[0]);
      53           0 :   _grad_test_y = _ndisp >= 2 ? &_fe.getPhiGradient(_disp_vars[1]) : nullptr;
      54           0 :   _grad_test_z = _ndisp >= 3 ? &_fe.getPhiGradient(_disp_vars[2]) : nullptr;
      55           0 : }
      56             : 
      57             : void
      58           0 : NEML2StressDivergence::forward()
      59             : {
      60           0 :   if (!_constitutive.outputReady())
      61           0 :     return;
      62             : 
      63           0 :   torch::InferenceMode guard;
      64             : 
      65           0 :   auto dphix = neml2::Tensor(*_grad_test_x, 3);
      66           0 :   auto dphiy = _ndisp >= 2 ? neml2::Tensor(*_grad_test_y, 3) : neml2::Tensor::zeros_like(dphix);
      67           0 :   auto dphiz = _ndisp >= 3 ? neml2::Tensor(*_grad_test_z, 3) : neml2::Tensor::zeros_like(dphix);
      68           0 :   auto dphi = neml2::base_stack({dphix, dphiy, dphiz}, 0);            // (nelem, ndofe, nqp; 3, 3)
      69           0 :   auto stress = neml2::R2(neml2::SR2(_stress)).dynamic_unsqueeze(-2); // (nelem, 1,     nqp; 3, 3)
      70             : 
      71             :   // weak form
      72           0 :   auto re_qp = neml2::base_sum(dphi * neml2::Tensor(stress), -1); // (nelem, ndofe, nqp; 3)
      73             : 
      74             :   // element integration
      75             :   auto JxWxT =
      76           0 :       _neml2_assembly.JxWxT().dynamic_unsqueeze(-2).base_unsqueeze(0); // (nelem, 1, nqp; 1)
      77           0 :   auto re = neml2::dynamic_sum(re_qp * JxWxT, -1)
      78           0 :                 .base_index({neml2::indexing::Slice(0, _ndisp)}); // (nelem, ndofe; ndisp)
      79             : 
      80             :   // assemble residual
      81           0 :   for (auto i : make_range(_ndisp))
      82             :   {
      83           0 :     auto re_i = re.base_index({i}).contiguous().cpu();
      84           0 :     const auto & dofmap_i = _fe.getGlobalDofMap(_disp_vars[i]);
      85           0 :     _residual->add_vector(re_i.data_ptr<Real>(), dofmap_i);
      86             :   }
      87           0 :   _residual->close();
      88           0 : }
      89             : 
      90             : #endif