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