Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://www.mooseframework.org
       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 "StressDivergenceTensorsTruss.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "Assembly.h"
      14             : #include "Material.h"
      15             : #include "MooseVariable.h"
      16             : #include "SystemBase.h"
      17             : 
      18             : registerMooseObject("TensorMechanicsApp", StressDivergenceTensorsTruss);
      19             : 
      20             : InputParameters
      21         112 : StressDivergenceTensorsTruss::validParams()
      22             : {
      23         112 :   InputParameters params = Kernel::validParams();
      24         112 :   params.addClassDescription("Kernel for truss element");
      25         224 :   params.addRequiredParam<unsigned int>("component",
      26             :                                         "An integer corresponding to the direction "
      27             :                                         "the variable this kernel acts in. (0 for x, "
      28             :                                         "1 for y, 2 for z)");
      29         224 :   params.addCoupledVar("displacements",
      30             :                        "The string of displacements suitable for the problem statement");
      31         224 :   params.addCoupledVar("temperature", "The temperature");
      32         224 :   params.addCoupledVar("area", "Cross-sectional area of truss element");
      33         224 :   params.addParam<std::string>("base_name", "Material property base name");
      34         112 :   params.set<bool>("use_displaced_mesh") = true;
      35         112 :   return params;
      36           0 : }
      37             : 
      38          62 : StressDivergenceTensorsTruss::StressDivergenceTensorsTruss(const InputParameters & parameters)
      39             :   : Kernel(parameters),
      40          62 :     _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
      41         124 :     _axial_stress(getMaterialPropertyByName<Real>(_base_name + "axial_stress")),
      42         124 :     _e_over_l(getMaterialPropertyByName<Real>(_base_name + "e_over_l")),
      43         124 :     _component(getParam<unsigned int>("component")),
      44          62 :     _ndisp(coupledComponents("displacements")),
      45          62 :     _temp_coupled(isCoupled("temperature")),
      46          62 :     _temp_var(_temp_coupled ? coupled("temperature") : 0),
      47          62 :     _area(coupledValue("area")),
      48         124 :     _orientation(NULL)
      49             : {
      50         204 :   for (unsigned int i = 0; i < _ndisp; ++i)
      51         284 :     _disp_var.push_back(coupled("displacements", i));
      52          62 : }
      53             : 
      54             : void
      55          62 : StressDivergenceTensorsTruss::initialSetup()
      56             : {
      57          62 :   _orientation = &_subproblem.assembly(_tid, _sys.number()).getFE(FEType(), 1)->get_dxyzdxi();
      58          62 : }
      59             : 
      60             : void
      61        4144 : StressDivergenceTensorsTruss::computeResidual()
      62             : {
      63        4144 :   prepareVectorTag(_assembly, _var.number());
      64             : 
      65             :   mooseAssert(_local_re.size() == 2, "Truss element has and only has two nodes.");
      66             : 
      67        4144 :   RealGradient orientation((*_orientation)[0]);
      68        4144 :   orientation /= orientation.norm();
      69             : 
      70        4144 :   VectorValue<Real> force_local = _axial_stress[0] * _area[0] * orientation;
      71             : 
      72        4144 :   _local_re(0) = -force_local(_component);
      73        4144 :   _local_re(1) = -_local_re(0);
      74             : 
      75        4144 :   accumulateTaggedLocalResidual();
      76             : 
      77        4144 :   if (_has_save_in)
      78             :   {
      79             :     Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
      80        8288 :     for (unsigned int i = 0; i < _save_in.size(); ++i)
      81        4144 :       _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
      82             :   }
      83        4144 : }
      84             : 
      85             : Real
      86        3424 : StressDivergenceTensorsTruss::computeStiffness(unsigned int i, unsigned int j)
      87             : {
      88        3424 :   RealGradient orientation((*_orientation)[0]);
      89        3424 :   orientation /= orientation.norm();
      90             : 
      91        3424 :   return orientation(i) * orientation(j) * _e_over_l[0] * _area[0];
      92             : }
      93             : 
      94             : void
      95         392 : StressDivergenceTensorsTruss::computeJacobian()
      96             : {
      97         392 :   prepareMatrixTag(_assembly, _var.number(), _var.number());
      98             : 
      99        1176 :   for (unsigned int i = 0; i < _test.size(); ++i)
     100        2352 :     for (unsigned int j = 0; j < _phi.size(); ++j)
     101        2352 :       _local_ke(i, j) += (i == j ? 1 : -1) * computeStiffness(_component, _component);
     102             : 
     103         392 :   accumulateTaggedLocalMatrix();
     104             : 
     105         392 :   if (_has_diag_save_in)
     106             :   {
     107             :     unsigned int rows = _local_ke.m();
     108           0 :     DenseVector<Number> diag(rows);
     109           0 :     for (unsigned int i = 0; i < rows; ++i)
     110           0 :       diag(i) = _local_ke(i, i);
     111             : 
     112             :     Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
     113           0 :     for (unsigned int i = 0; i < _diag_save_in.size(); ++i)
     114           0 :       _diag_save_in[i]->sys().solution().add_vector(diag, _diag_save_in[i]->dofIndices());
     115             :   }
     116         392 : }
     117             : 
     118             : void
     119         784 : StressDivergenceTensorsTruss::computeOffDiagJacobian(const unsigned int jvar_num)
     120             : {
     121         784 :   if (jvar_num == _var.number())
     122         320 :     computeJacobian();
     123             :   else
     124             :   {
     125             :     const auto & jvar = getVariable(jvar_num);
     126             : 
     127             :     // This (undisplaced) jvar could potentially yield the wrong phi size if this object is acting
     128             :     // on the displaced mesh
     129         464 :     auto phi_size = jvar.dofIndices().size();
     130             : 
     131             :     unsigned int coupled_component = 0;
     132             :     bool disp_coupled = false;
     133             : 
     134         840 :     for (unsigned int i = 0; i < _ndisp; ++i)
     135         840 :       if (jvar_num == _disp_var[i])
     136             :       {
     137             :         coupled_component = i;
     138             :         disp_coupled = true;
     139             :         break;
     140             :       }
     141             : 
     142             :     if (disp_coupled)
     143             :     {
     144         464 :       prepareMatrixTag(_assembly, _var.number(), jvar_num);
     145             : 
     146        1392 :       for (unsigned int i = 0; i < _test.size(); ++i)
     147        2784 :         for (unsigned int j = 0; j < phi_size; ++j)
     148        2784 :           _local_ke(i, j) += (i == j ? 1 : -1) * computeStiffness(_component, coupled_component);
     149             : 
     150         464 :       accumulateTaggedLocalMatrix();
     151             :     }
     152             :     else if (false) // Need some code here for coupling with temperature
     153             :     {
     154             :     }
     155             :   }
     156         784 : }