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 "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("SolidMechanicsApp", StressDivergenceTensorsTruss);
      19             : 
      20             : InputParameters
      21         224 : StressDivergenceTensorsTruss::validParams()
      22             : {
      23         224 :   InputParameters params = Kernel::validParams();
      24         224 :   params.addClassDescription("Kernel for truss element");
      25         448 :   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         448 :   params.addCoupledVar("displacements",
      30             :                        "The string of displacements suitable for the problem statement");
      31         448 :   params.addCoupledVar("temperature", "The temperature");
      32         448 :   params.addCoupledVar("area", "Cross-sectional area of truss element");
      33         448 :   params.addParam<std::string>("base_name", "Material property base name");
      34         224 :   params.set<bool>("use_displaced_mesh") = true;
      35         224 :   return params;
      36           0 : }
      37             : 
      38         124 : StressDivergenceTensorsTruss::StressDivergenceTensorsTruss(const InputParameters & parameters)
      39             :   : Kernel(parameters),
      40         124 :     _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
      41         248 :     _axial_stress(getMaterialPropertyByName<Real>(_base_name + "axial_stress")),
      42         248 :     _e_over_l(getMaterialPropertyByName<Real>(_base_name + "e_over_l")),
      43         248 :     _component(getParam<unsigned int>("component")),
      44         124 :     _ndisp(coupledComponents("displacements")),
      45         124 :     _temp_coupled(isCoupled("temperature")),
      46         124 :     _temp_var(_temp_coupled ? coupled("temperature") : 0),
      47         124 :     _area(coupledValue("area")),
      48         248 :     _orientation(NULL)
      49             : {
      50         408 :   for (unsigned int i = 0; i < _ndisp; ++i)
      51         568 :     _disp_var.push_back(coupled("displacements", i));
      52         124 : }
      53             : 
      54             : void
      55         124 : StressDivergenceTensorsTruss::initialSetup()
      56             : {
      57         124 :   _orientation = &_subproblem.assembly(_tid, _sys.number()).getFE(FEType(), 1)->get_dxyzdxi();
      58         124 : }
      59             : 
      60             : void
      61        7632 : StressDivergenceTensorsTruss::computeResidual()
      62             : {
      63        7632 :   prepareVectorTag(_assembly, _var.number());
      64             : 
      65             :   mooseAssert(_local_re.size() == 2, "Truss element has and only has two nodes.");
      66             : 
      67        7632 :   RealGradient orientation((*_orientation)[0]);
      68        7632 :   orientation /= orientation.norm();
      69             : 
      70        7632 :   VectorValue<Real> force_local = _axial_stress[0] * _area[0] * orientation;
      71             : 
      72        7632 :   _local_re(0) = -force_local(_component);
      73        7632 :   _local_re(1) = -_local_re(0);
      74             : 
      75        7632 :   accumulateTaggedLocalResidual();
      76             : 
      77        7632 :   if (_has_save_in)
      78             :   {
      79             :     Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
      80       15264 :     for (unsigned int i = 0; i < _save_in.size(); ++i)
      81        7632 :       _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
      82             :   }
      83        7632 : }
      84             : 
      85             : Real
      86        6848 : StressDivergenceTensorsTruss::computeStiffness(unsigned int i, unsigned int j)
      87             : {
      88        6848 :   RealGradient orientation((*_orientation)[0]);
      89        6848 :   orientation /= orientation.norm();
      90             : 
      91        6848 :   return orientation(i) * orientation(j) * _e_over_l[0] * _area[0];
      92             : }
      93             : 
      94             : void
      95         784 : StressDivergenceTensorsTruss::computeJacobian()
      96             : {
      97         784 :   prepareMatrixTag(_assembly, _var.number(), _var.number());
      98             : 
      99        2352 :   for (unsigned int i = 0; i < _test.size(); ++i)
     100        4704 :     for (unsigned int j = 0; j < _phi.size(); ++j)
     101        4704 :       _local_ke(i, j) += (i == j ? 1 : -1) * computeStiffness(_component, _component);
     102             : 
     103         784 :   accumulateTaggedLocalMatrix();
     104             : 
     105         784 :   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         784 : }
     117             : 
     118             : void
     119        1568 : StressDivergenceTensorsTruss::computeOffDiagJacobian(const unsigned int jvar_num)
     120             : {
     121        1568 :   if (jvar_num == _var.number())
     122         640 :     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         928 :     auto phi_size = jvar.dofIndices().size();
     130             : 
     131             :     unsigned int coupled_component = 0;
     132             :     bool disp_coupled = false;
     133             : 
     134        1680 :     for (unsigned int i = 0; i < _ndisp; ++i)
     135        1680 :       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         928 :       prepareMatrixTag(_assembly, _var.number(), jvar_num);
     145             : 
     146        2784 :       for (unsigned int i = 0; i < _test.size(); ++i)
     147        5568 :         for (unsigned int j = 0; j < phi_size; ++j)
     148        5568 :           _local_ke(i, j) += (i == j ? 1 : -1) * computeStiffness(_component, coupled_component);
     149             : 
     150         928 :       accumulateTaggedLocalMatrix();
     151             :     }
     152             :     else if (false) // Need some code here for coupling with temperature
     153             :     {
     154             :     }
     155             :   }
     156        1568 : }