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 : #include "StressDivergenceIsolator.h"
17 :
18 : // MOOSE includes
19 : #include "Assembly.h"
20 : #include "Material.h"
21 : #include "MooseVariable.h"
22 : #include "SystemBase.h"
23 : #include "NonlinearSystem.h"
24 : #include "MooseMesh.h"
25 :
26 : // libmesh includes
27 : #include "libmesh/quadrature.h"
28 :
29 : registerMooseObject("MastodonApp", StressDivergenceIsolator);
30 :
31 : InputParameters
32 492 : StressDivergenceIsolator::validParams()
33 : {
34 492 : InputParameters params = Kernel::validParams();
35 492 : params.addClassDescription("Kernel for isolator element");
36 984 : params.addRequiredParam<unsigned int>(
37 : "component",
38 : "An integer corresponding to the direction "
39 : "the variable this kernel acts in. (0 for x, "
40 : "1 for y, 2 for z, 3 for rot_x, 4 for rot_y and 5 for rot_z).");
41 984 : params.addRequiredCoupledVar("displacements", "The displacement variables for isolator.");
42 984 : params.addRequiredCoupledVar("rotations", "The rotation variables for the isolator.");
43 984 : params.addParam<Real>(
44 : "zeta",
45 984 : 0.0,
46 : "Name of material property or a constant real number defining the zeta parameter for the "
47 : "Rayleigh damping.");
48 1476 : params.addRangeCheckedParam<Real>(
49 984 : "alpha", 0.0, "alpha >= -0.3333 & alpha <= 0.0", "alpha parameter for HHT time integration");
50 984 : params.addParam<bool>("static_initialization",
51 984 : false,
52 : "Set to true to get the system to "
53 : "equilibrium under gravity by running a "
54 : "quasi-static analysis (by solving Ku = F) "
55 : "in the first time step");
56 492 : params.set<bool>("use_displaced_mesh") = true;
57 492 : return params;
58 0 : }
59 :
60 240 : StressDivergenceIsolator::StressDivergenceIsolator(const InputParameters & parameters)
61 : : Kernel(parameters),
62 240 : _component(getParam<unsigned int>("component")),
63 240 : _ndisp(coupledComponents("displacements")),
64 240 : _disp_var(_ndisp),
65 240 : _nrot(coupledComponents("rotations")),
66 240 : _rot_var(_nrot),
67 480 : _Fg(getMaterialPropertyByName<ColumnMajorMatrix>("global_forces")),
68 240 : _Kg(getMaterialPropertyByName<ColumnMajorMatrix>("global_stiffness_matrix")),
69 480 : _zeta(getParam<Real>("zeta")),
70 480 : _alpha(getParam<Real>("alpha")),
71 240 : _isDamped(_zeta != 0.0 || std::abs(_alpha) > 0.0),
72 294 : _Fg_old(_isDamped ? &getMaterialPropertyOld<ColumnMajorMatrix>("global_forces") : nullptr),
73 294 : _Fg_older(_isDamped ? &getMaterialPropertyOlder<ColumnMajorMatrix>("global_forces") : nullptr),
74 720 : _static_initialization(getParam<bool>("static_initialization"))
75 : {
76 240 : if (_component > 5)
77 0 : mooseError("Error in StressDivergenceIsolator block ",
78 : name(),
79 : ". Please enter an integer value between 0 and 5 for the 'component' parameter.");
80 :
81 240 : if (_ndisp != _nrot)
82 0 : mooseError("Error in StressDivergenceIsolator block ",
83 : name(),
84 : ". The number of displacement and rotation variables should be the same.");
85 :
86 960 : for (unsigned int i = 0; i < _ndisp; ++i)
87 720 : _disp_var[i] = coupled("displacements", i);
88 :
89 960 : for (unsigned int i = 0; i < _nrot; ++i)
90 720 : _rot_var[i] = coupled("rotations", i);
91 240 : }
92 :
93 : void
94 207432 : StressDivergenceIsolator::computeResidual()
95 : {
96 207432 : prepareVectorTag(_assembly, _var.number());
97 :
98 : mooseAssert(_local_re.size() == 2, "Isolator element only has two nodes.");
99 :
100 207432 : ColumnMajorMatrix global_force_res = _Fg[0];
101 :
102 : // add contributions from stiffness proportional damping (non-zero _zeta) or HHT time integration
103 : // (non-zero _alpha)
104 207432 : if (_isDamped && _dt > 0.0 && !(_static_initialization && _t <= 2 * _dt))
105 0 : global_force_res = global_force_res * (1.0 + _alpha + (1.0 + _alpha) * _zeta / _dt) -
106 0 : (*_Fg_old)[0] * (_alpha + (1.0 + 2.0 * _alpha) * _zeta / _dt) +
107 0 : (*_Fg_older)[0] * (_alpha * _zeta / _dt);
108 :
109 207432 : _local_re(0) = global_force_res(_component);
110 :
111 207432 : _local_re(1) = global_force_res(_component + 6);
112 :
113 207432 : accumulateTaggedLocalResidual();
114 :
115 207432 : if (_has_save_in)
116 414864 : for (unsigned int i = 0; i < _save_in.size(); ++i)
117 207432 : _save_in[i]->sys().solution().add_vector(_local_re, _save_in[i]->dofIndices());
118 207432 : }
119 :
120 : void
121 59112 : StressDivergenceIsolator::computeJacobian()
122 : {
123 59112 : prepareMatrixTag(_assembly, _var.number(), _var.number());
124 :
125 : // i and j are looping over nodes
126 177336 : for (unsigned int i = 0; i < _test.size(); ++i)
127 354672 : for (unsigned int j = 0; j < _phi.size(); ++j)
128 236448 : _local_ke(i, j) += _Kg[0](i * 6 + _component, j * 6 + _component);
129 :
130 : // scaling factor for Rayliegh damping and HHT time integration
131 59112 : if (_isDamped && _dt > 0.0 && !(_static_initialization && _t == _dt))
132 0 : _local_ke *= (1.0 + _alpha + (1.0 + _alpha) * _zeta / _dt);
133 :
134 59112 : accumulateTaggedLocalMatrix();
135 :
136 59112 : if (_has_diag_save_in)
137 : {
138 : unsigned int rows = _local_ke.m();
139 0 : DenseVector<Number> diag(rows);
140 0 : for (unsigned int i = 0; i < rows; ++i)
141 0 : diag(i) = _local_ke(i, i);
142 :
143 0 : for (unsigned int i = 0; i < _diag_save_in.size(); ++i)
144 0 : _diag_save_in[i]->sys().solution().add_vector(diag, _diag_save_in[i]->dofIndices());
145 : }
146 59112 : }
147 :
148 : void
149 354672 : StressDivergenceIsolator::computeOffDiagJacobian(const unsigned int jvar_num)
150 : // coupling one variable to another (disp x to disp y, etc)
151 : {
152 354672 : if (jvar_num == _var.number())
153 : // jacobian calculation if jvar is the same as the current variable i.e.,
154 : // diagonal elements
155 59112 : computeJacobian();
156 :
157 : else
158 : // jacobian calculation for off-diagonal elements
159 : {
160 : unsigned int coupled_component = 0;
161 : bool coupled = false;
162 : // finding which variable jvar is
163 591120 : for (unsigned int i = 0; i < _ndisp; ++i)
164 : {
165 591120 : if (jvar_num == _disp_var[i])
166 : {
167 : coupled_component = i;
168 : coupled = true;
169 : break;
170 : }
171 443340 : else if (jvar_num == _rot_var[i])
172 : {
173 147780 : coupled_component = i + 3;
174 : coupled = true;
175 147780 : break;
176 : }
177 : }
178 :
179 295560 : prepareMatrixTag(_assembly, _var.number(), jvar_num);
180 :
181 : // getting the jacobian from assembly
182 : // DenseMatrix<Number> & ke = _assembly.jacobianBlock(_var.number(), jvar_num);
183 295560 : if (coupled)
184 : {
185 886680 : for (unsigned int i = 0; i < _test.size(); ++i)
186 1773360 : for (unsigned int j = 0; j < _phi.size(); ++j)
187 1182240 : _local_ke(i, j) += _Kg[0](i * 6 + _component, j * 6 + coupled_component);
188 : }
189 :
190 : // scaling factor for Rayliegh damping and HHT time integration
191 295560 : if (_isDamped && _dt > 0.0 && !(_static_initialization && _t == _dt))
192 0 : _local_ke *= (1.0 + _alpha + (1.0 + _alpha) * _zeta / _dt);
193 :
194 295560 : accumulateTaggedLocalMatrix();
195 : }
196 354672 : }
|
