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