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