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 "DisplacementAboutAxis.h"
11 : #include "Function.h"
12 :
13 : registerMooseObject("SolidMechanicsApp", DisplacementAboutAxis);
14 :
15 : InputParameters
16 422 : DisplacementAboutAxis::validParams()
17 : {
18 422 : InputParameters params = DirichletBCBase::validParams();
19 422 : params.addClassDescription("Implements a boundary condition that enforces rotational"
20 : "displacement around an axis on a boundary");
21 422 : addDisplacementAboutAxisParams(params);
22 844 : params.addRequiredParam<int>("component", "The component for the rotational displacement");
23 422 : params.set<bool>("use_displaced_mesh") = false;
24 422 : params.set<bool>("preset") = true;
25 422 : return params;
26 0 : }
27 :
28 : void
29 422 : addDisplacementAboutAxisParams(InputParameters & params)
30 : {
31 844 : MooseEnum units("degrees radians");
32 844 : params.addRequiredParam<FunctionName>(
33 : "function", "The function providing the total angle of rotation or the angular velocity.");
34 844 : params.addRequiredParam<MooseEnum>("angle_units",
35 : units,
36 422 : "The units of the angle of rotation. Choices are:" +
37 422 : units.getRawNames());
38 844 : params.addRequiredParam<RealVectorValue>("axis_origin", "Origin of the axis of rotation");
39 844 : params.addRequiredParam<RealVectorValue>("axis_direction", "Direction of the axis of rotation");
40 844 : params.addParam<bool>(
41 844 : "angular_velocity", false, "If true interprets the function value as an angular velocity");
42 844 : params.addRequiredCoupledVar("displacements",
43 : "The string of displacements suitable for the problem statement");
44 422 : }
45 :
46 212 : DisplacementAboutAxis::DisplacementAboutAxis(const InputParameters & parameters)
47 : : DirichletBCBase(parameters),
48 212 : _component(getParam<int>("component")),
49 212 : _func(getFunction("function")),
50 424 : _angle_units(getParam<MooseEnum>("angle_units")),
51 212 : _axis_origin(getParam<RealVectorValue>("axis_origin")),
52 212 : _axis_direction(getParam<RealVectorValue>("axis_direction")),
53 212 : _ndisp(coupledComponents("displacements")),
54 212 : _disp_old(_ndisp),
55 636 : _angular_velocity(getParam<bool>("angular_velocity"))
56 : {
57 212 : if (_component < 0 || _component > 2)
58 1 : mooseError("Invalid component given for ", name(), ": ", _component, ".");
59 :
60 211 : if (_axis_direction.norm() == 0.)
61 1 : mooseError("Please specify a non-zero direction vector for the axis_direction in ", name());
62 210 : }
63 :
64 : void
65 208 : DisplacementAboutAxis::initialSetup()
66 : {
67 208 : calculateUnitDirectionVector();
68 208 : calculateTransformationMatrices();
69 :
70 208 : if (_angular_velocity)
71 592 : for (unsigned int i = 0; i < _ndisp; ++i)
72 444 : _disp_old[i] = &coupledDofValuesOld("displacements", i);
73 : else
74 240 : for (unsigned int i = 0; i < _ndisp; ++i)
75 180 : _disp_old[i] = nullptr;
76 208 : }
77 :
78 : Real
79 203536 : DisplacementAboutAxis::computeQpValue()
80 : {
81 203536 : Point p(*_current_node);
82 :
83 203536 : Real angle(_func.value(_t, *_current_node));
84 203536 : if (_angle_units == "degrees")
85 203536 : angle = angle * libMesh::pi / 180.0;
86 :
87 203536 : if (_angular_velocity)
88 84256 : angle *= _dt;
89 :
90 203536 : ColumnMajorMatrix p_old(4, 1);
91 203536 : p_old(0, 0) = p(0);
92 203536 : p_old(1, 0) = p(1);
93 203536 : p_old(2, 0) = p(2);
94 203536 : p_old(3, 0) = 1;
95 :
96 203536 : if (_angular_velocity)
97 337024 : for (unsigned int i = 0; i < _ndisp; i++)
98 252768 : p_old(i, 0) += (*_disp_old[i])[_qp];
99 :
100 203536 : ColumnMajorMatrix p_new = rotateAroundAxis(p_old, angle);
101 :
102 407072 : return p_new(_component, 0) - p(_component);
103 : }
104 :
105 : ColumnMajorMatrix
106 203536 : DisplacementAboutAxis::rotateAroundAxis(const ColumnMajorMatrix & p0, const Real angle)
107 : {
108 203536 : ColumnMajorMatrix rotate_about_z(4, 4);
109 203536 : rotate_about_z(0, 0) = cos(angle);
110 203536 : rotate_about_z(0, 1) = -sin(angle);
111 203536 : rotate_about_z(0, 2) = 0;
112 203536 : rotate_about_z(0, 3) = 0;
113 203536 : rotate_about_z(1, 0) = sin(angle);
114 203536 : rotate_about_z(1, 1) = cos(angle);
115 203536 : rotate_about_z(1, 2) = 0;
116 203536 : rotate_about_z(1, 3) = 0;
117 203536 : rotate_about_z(2, 0) = 0;
118 203536 : rotate_about_z(2, 1) = 0;
119 203536 : rotate_about_z(2, 2) = 1;
120 203536 : rotate_about_z(2, 3) = 0;
121 203536 : rotate_about_z(3, 0) = 0;
122 203536 : rotate_about_z(3, 1) = 0;
123 203536 : rotate_about_z(3, 2) = 0;
124 203536 : rotate_about_z(3, 3) = 1;
125 :
126 : ColumnMajorMatrix transform =
127 203536 : _transformation_matrix_inv * rotate_about_z * _transformation_matrix;
128 407072 : return transform * p0;
129 : }
130 :
131 : void
132 208 : DisplacementAboutAxis::calculateUnitDirectionVector()
133 : {
134 208 : Real magnitude = _axis_direction.norm();
135 : _axis_direction /= magnitude;
136 208 : }
137 :
138 : void
139 208 : DisplacementAboutAxis::calculateTransformationMatrices()
140 : {
141 : // These parts of the transformation matrix only depend on the axis of rotation:
142 :
143 : Real length = _axis_direction.norm_sq();
144 208 : Real v = _axis_direction(1) * _axis_direction(1) + _axis_direction(2) * _axis_direction(2);
145 :
146 208 : ColumnMajorMatrix transl(4, 4);
147 208 : transl(0, 0) = 1;
148 208 : transl(0, 1) = 0;
149 208 : transl(0, 2) = 0;
150 208 : transl(0, 3) = -_axis_origin(0);
151 208 : transl(1, 0) = 0;
152 208 : transl(1, 1) = 1;
153 208 : transl(1, 2) = 0;
154 208 : transl(1, 3) = -_axis_origin(1);
155 208 : transl(2, 0) = 0;
156 208 : transl(2, 1) = 0;
157 208 : transl(2, 2) = 1;
158 208 : transl(2, 3) = -_axis_origin(2);
159 208 : transl(3, 0) = 0;
160 208 : transl(3, 1) = 0;
161 208 : transl(3, 2) = 0;
162 208 : transl(3, 3) = 1;
163 :
164 208 : ColumnMajorMatrix rotate_about_x(4, 4);
165 208 : rotate_about_x(0, 0) = 1;
166 208 : rotate_about_x(0, 1) = 0;
167 208 : rotate_about_x(0, 2) = 0;
168 208 : rotate_about_x(0, 3) = 0;
169 208 : rotate_about_x(1, 0) = 0;
170 208 : rotate_about_x(1, 1) = _axis_direction(2) / v;
171 208 : rotate_about_x(1, 2) = -_axis_direction(1) / v;
172 208 : rotate_about_x(1, 3) = 0;
173 208 : rotate_about_x(2, 0) = 0;
174 208 : rotate_about_x(2, 1) = _axis_direction(1) / v;
175 208 : rotate_about_x(2, 2) = _axis_direction(2) / v;
176 208 : rotate_about_x(2, 3) = 0;
177 208 : rotate_about_x(3, 0) = 0;
178 208 : rotate_about_x(3, 1) = 0;
179 208 : rotate_about_x(3, 2) = 0;
180 208 : rotate_about_x(3, 3) = 1;
181 :
182 208 : ColumnMajorMatrix rotate_about_y(4, 4);
183 208 : rotate_about_y(0, 0) = v / length;
184 208 : rotate_about_y(0, 1) = 0;
185 208 : rotate_about_y(0, 2) = -_axis_direction(0) / length;
186 208 : rotate_about_y(0, 3) = 0;
187 208 : rotate_about_y(1, 0) = 0;
188 208 : rotate_about_y(1, 1) = 1;
189 208 : rotate_about_y(1, 2) = 0;
190 208 : rotate_about_y(1, 3) = 0;
191 208 : rotate_about_y(2, 0) = _axis_direction(0) / length;
192 208 : rotate_about_y(2, 1) = 0;
193 208 : rotate_about_y(2, 2) = v / length;
194 208 : rotate_about_y(2, 3) = 0;
195 208 : rotate_about_y(3, 0) = 0;
196 208 : rotate_about_y(3, 1) = 0;
197 208 : rotate_about_y(3, 2) = 0;
198 208 : rotate_about_y(3, 3) = 1;
199 :
200 208 : ColumnMajorMatrix transl_inv(4, 4);
201 208 : transl.inverse(transl_inv);
202 208 : ColumnMajorMatrix rotx_inv(4, 4);
203 208 : rotate_about_x.inverse(rotx_inv);
204 208 : ColumnMajorMatrix roty_inv(4, 4);
205 208 : rotate_about_y.inverse(roty_inv);
206 :
207 208 : _transformation_matrix = rotate_about_y * rotate_about_x * transl;
208 416 : _transformation_matrix_inv = transl_inv * rotx_inv * roty_inv;
209 208 : }
|
