Line data    Source code

       1             : // MASTODON includes
       2             : #include "StressDivergenceDamper.h"
       3             : 
       4             : // MOOSE includes
       5             : #include "Assembly.h"
       6             : #include "Material.h"
       7             : #include "MooseVariable.h"
       8             : #include "SystemBase.h"
       9             : #include "NonlinearSystem.h"
      10             : #include "MooseMesh.h"
      11             : 
      12             : // libmesh includes
      13             : #include "libmesh/quadrature.h"
      14             : 
      15             : registerMooseObject("MastodonApp", StressDivergenceDamper);
      16             : 
      17             : InputParameters
      18          90 : StressDivergenceDamper::validParams()
      19             : {
      20          90 :   InputParameters params = Kernel::validParams();
      21          90 :   params.addClassDescription("Kernel for two-noded nonlinear fluid viscous damper element");
      22         180 :   params.addRequiredParam<unsigned int>("component",
      23             :                                         "An integer corresponding to the direction "
      24             :                                         "the variable this kernel acts in. (0 for x, "
      25             :                                         "1 for y, 2 for z).");
      26         180 :   params.addRequiredCoupledVar("displacements", "The displacement variables for Damper.");
      27          90 :   params.set<bool>("use_displaced_mesh") = true;
      28          90 :   return params;
      29           0 : }
      30             : 
      31          45 : StressDivergenceDamper::StressDivergenceDamper(const InputParameters & parameters)
      32             :   : Kernel(parameters),
      33          45 :     _component(getParam<unsigned int>("component")),
      34          45 :     _ndisp(coupledComponents("displacements")),
      35          45 :     _disp_var(_ndisp),
      36          90 :     _Fg(getMaterialPropertyByName<ColumnMajorMatrix>("global_forces")),
      37          90 :     _Kg(getMaterialPropertyByName<ColumnMajorMatrix>("global_stiffness_matrix"))
      38             : 
      39             : {
      40          45 :   if (_component > 2)
      41           0 :     mooseError("Error in StressDivergenceDamper block ",
      42             :                name(),
      43             :                ". Please enter an integer value between 0 and 2 for the 'component' parameter.");
      44             : 
      45         180 :   for (unsigned int i = 0; i < _ndisp; ++i)
      46         135 :     _disp_var[i] = coupled("displacements", i);
      47          45 : }
      48             : 
      49             : void
      50       48420 : StressDivergenceDamper::computeResidual()
      51             : {
      52             :   // Accessing residual vector, re, from MOOSE assembly
      53       48420 :   prepareVectorTag(_assembly, _var.number());
      54             :   mooseAssert(_local_re.size() == 2, "Damper element only has two nodes.");
      55             : 
      56             :   // Calculating residual for node 0 (external forces on node 0)
      57       48420 :   _local_re(0) = _Fg[0](_component);
      58             : 
      59             :   // External force on node 1 (external forces on node 1)
      60       48420 :   _local_re(1) = _Fg[0](_component + 3);
      61             : 
      62       48420 :   accumulateTaggedLocalResidual();
      63             : 
      64       48420 :   if (_has_save_in)
      65       96840 :     for (unsigned int i = 0; i < _save_in.size(); ++i)
      66       48420 :       _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
      67       48420 : }
      68             : 
      69             : void
      70        4836 : StressDivergenceDamper::computeJacobian()
      71             : {
      72             :   // Access Jacobian; size is n x n (n is number of nodes)
      73        4836 :   prepareMatrixTag(_assembly, _var.number(), _var.number());
      74             : 
      75             :   // i and j are looping over nodes
      76       14508 :   for (unsigned int i = 0; i < _test.size(); ++i)
      77       29016 :     for (unsigned int j = 0; j < _phi.size(); ++j)
      78             :     {
      79       19344 :       _local_ke(i, j) += _Kg[0](i * 3 + _component, j * 3 + _component);
      80             :     }
      81             : 
      82        4836 :   accumulateTaggedLocalMatrix();
      83             : 
      84        4836 :   if (_has_diag_save_in)
      85             :   {
      86             :     unsigned int rows = _local_ke.m();
      87           0 :     DenseVector<Number> diag(rows);
      88           0 :     for (unsigned int i = 0; i < rows; ++i)
      89           0 :       diag(i) = _local_ke(i, i);
      90             : 
      91           0 :     for (unsigned int i = 0; i < _diag_save_in.size(); ++i)
      92           0 :       _diag_save_in[i]->sys().solution().add_vector(diag, _diag_save_in[i]->dofIndices());
      93             :   }
      94        4836 : }
      95             : 
      96             : void
      97       14508 : StressDivergenceDamper::computeOffDiagJacobian(const unsigned int jvar_num)
      98             : // coupling one variable to another (disp x to disp y, etc)
      99             : {
     100       14508 :   if (jvar_num == _var.number())
     101             :     // jacobian calculation if jvar is the same as the current variable i.e.,
     102             :     // diagonal elements
     103        4836 :     computeJacobian();
     104             : 
     105             :   else
     106             :   // jacobian calculation for off-diagonal elements
     107             :   {
     108             :     unsigned int coupled_component = 0;
     109             :     bool coupled = false;
     110             :     // finding which variable jvar is
     111       19344 :     for (unsigned int i = 0; i < _ndisp; ++i)
     112             :     {
     113       19344 :       if (jvar_num == _disp_var[i])
     114             :       {
     115             :         coupled_component = i;
     116             :         coupled = true;
     117             :         break;
     118             :       }
     119             :     }
     120             :     // getting the jacobian from assembly
     121        9672 :     prepareMatrixTag(_assembly, _var.number(), jvar_num);
     122        9672 :     if (coupled)
     123             :     {
     124       29016 :       for (unsigned int i = 0; i < _test.size(); ++i)
     125       58032 :         for (unsigned int j = 0; j < _phi.size(); ++j)
     126       38688 :           _local_ke(i, j) += _Kg[0](i * 3 + _component, j * 3 + coupled_component);
     127        9672 :       accumulateTaggedLocalMatrix();
     128             :     }
     129             :   }
     130       14508 : }