Line data    Source code

       1             : /*************************************************/
       2             : /*           DO NOT MODIFY THIS HEADER           */
       3             : /*                                               */
       4             : /*                     MASTODON                  */
       5             : /*                                               */
       6             : /*    (c) 2015 Battelle Energy Alliance, LLC     */
       7             : /*            ALL RIGHTS RESERVED                */
       8             : /*                                               */
       9             : /*   Prepared by Battelle Energy Alliance, LLC   */
      10             : /*     With the U. S. Department of Energy       */
      11             : /*                                               */
      12             : /*     See COPYRIGHT for full restrictions       */
      13             : /*************************************************/
      14             : 
      15             : // MASTODON includes
      16             : 
      17             : #include "StressDivergenceSpring.h"
      18             : 
      19             : // MOOSE includes
      20             : #include "Assembly.h"
      21             : #include "Material.h"
      22             : #include "MooseVariable.h"
      23             : #include "SystemBase.h"
      24             : #include "RankTwoTensor.h"
      25             : #include "NonlinearSystem.h"
      26             : #include "MooseMesh.h"
      27             : 
      28             : // libmesh includes
      29             : #include "libmesh/quadrature.h"
      30             : 
      31             : registerMooseObject("MastodonApp", StressDivergenceSpring);
      32             : 
      33             : InputParameters
      34          72 : StressDivergenceSpring::validParams()
      35             : {
      36          72 :   InputParameters params = Kernel::validParams();
      37          72 :   params.addClassDescription("Kernel for spring element");
      38         144 :   params.addRequiredParam<unsigned int>(
      39             :       "component",
      40             :       "An integer corresponding to the direction "
      41             :       "the variable this kernel acts in. (0 for x, "
      42             :       "1 for y, 2 for z, 3 for rot_x, 4 for rot_y and 5 for rot_z).");
      43         144 :   params.addRequiredCoupledVar("displacements", "The displacement variables for spring.");
      44         144 :   params.addRequiredCoupledVar("rotations", "The rotation variables for the spring.");
      45          72 :   params.set<bool>("use_displaced_mesh") = true;
      46          72 :   return params;
      47           0 : }
      48             : 
      49          36 : StressDivergenceSpring::StressDivergenceSpring(const InputParameters & parameters)
      50             :   : Kernel(parameters),
      51          36 :     _component(getParam<unsigned int>("component")),
      52          36 :     _ndisp(coupledComponents("displacements")),
      53          36 :     _disp_var(_ndisp),
      54          36 :     _nrot(coupledComponents("rotations")),
      55          36 :     _rot_var(_nrot),
      56          72 :     _spring_forces_global(getMaterialPropertyByName<ColumnMajorMatrix>("global_forces")),
      57          36 :     _spring_moments_global(getMaterialPropertyByName<ColumnMajorMatrix>("global_moments")),
      58          36 :     _kdd(getMaterialPropertyByName<ColumnMajorMatrix>("displacement_stiffness_matrix")),
      59          36 :     _krr(getMaterialPropertyByName<ColumnMajorMatrix>("rotation_stiffness_matrix")),
      60          36 :     _total_global_to_local_rotation(
      61          36 :         getMaterialPropertyByName<ColumnMajorMatrix>("total_global_to_local_rotation"))
      62             : {
      63          36 :   if (_component > 5)
      64           0 :     mooseError("Error in StressDivergenceSpring block ",
      65             :                name(),
      66             :                ". Please enter an integer value between 0 and 5 for the 'component' parameter.");
      67             : 
      68          36 :   if (_ndisp != _nrot)
      69           0 :     mooseError("Error in StressDivergenceSpring block ",
      70             :                name(),
      71             :                ". The number of displacement and rotation variables should be the same.");
      72             : 
      73         144 :   for (unsigned int i = 0; i < _ndisp; ++i)
      74         108 :     _disp_var[i] = coupled("displacements", i);
      75             : 
      76         144 :   for (unsigned int i = 0; i < _nrot; ++i)
      77         108 :     _rot_var[i] = coupled("rotations", i);
      78          36 : }
      79             : 
      80             : void
      81       57600 : StressDivergenceSpring::computeResidual()
      82             : {
      83             :   // Accessing residual vector, re, from MOOSE assembly
      84       57600 :   prepareVectorTag(_assembly, _var.number());
      85             :   mooseAssert(_local_re.size() == 2, "Spring element has and only has two nodes.");
      86             : 
      87             :   // Calculating residual for node 0 (external forces on node 0)
      88       57600 :   if (_component < 3)
      89       28800 :     _local_re(0) = -_spring_forces_global[0](_component);
      90             :   else
      91       28800 :     _local_re(0) = -_spring_moments_global[0](_component - 3);
      92             : 
      93             :   // External force on node 1 = -1 * external force on node 0
      94       57600 :   _local_re(1) = -_local_re(0);
      95             : 
      96       57600 :   accumulateTaggedLocalResidual();
      97             : 
      98       57600 :   if (_has_save_in)
      99           0 :     for (unsigned int i = 0; i < _save_in.size(); ++i)
     100           0 :       _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
     101       57600 : }
     102             : 
     103             : void
     104       19200 : StressDivergenceSpring::computeJacobian()
     105             : {
     106             :   // Access Jacobian; size is n x n (n is number of nodes)
     107       19200 :   prepareMatrixTag(_assembly, _var.number(), _var.number());
     108             : 
     109       57600 :   for (unsigned int i = 0; i < _test.size(); ++i)
     110      115200 :     for (unsigned int j = 0; j < _phi.size(); ++j)
     111             :     {
     112       76800 :       if (_component < 3)
     113       57600 :         _local_ke(i, j) += (i == j ? 1 : -1) * _kdd[0](_component, _component);
     114             :       else
     115       57600 :         _local_ke(i, j) += (i == j ? 1 : -1) * _krr[0](_component - 3, _component - 3);
     116             :     }
     117             : 
     118       19200 :   accumulateTaggedLocalMatrix();
     119             : 
     120       19200 :   if (_has_diag_save_in)
     121             :   {
     122             :     unsigned int rows = _local_ke.m();
     123           0 :     DenseVector<Number> diag(rows);
     124           0 :     for (unsigned int i = 0; i < rows; ++i)
     125           0 :       diag(i) = _local_ke(i, i);
     126             : 
     127           0 :     for (unsigned int i = 0; i < _diag_save_in.size(); ++i)
     128           0 :       _diag_save_in[i]->sys().solution().add_vector(diag, _diag_save_in[i]->dofIndices());
     129             :   }
     130       19200 : }
     131             : 
     132             : void
     133      115200 : StressDivergenceSpring::computeOffDiagJacobian(const unsigned int jvar_num)
     134             : // coupling one variable to another (disp x to disp y, etc)
     135             : {
     136      115200 :   if (jvar_num == _var.number())
     137             :     // jacobian calculation if jvar is the same as the current variable i.e.,
     138             :     // diagonal elements
     139       19200 :     computeJacobian();
     140             : 
     141             :   // Off-diagonal elements are zero for linear elastic spring
     142      115200 : }