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