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 "DiscreteLineSegmentInterface.h"
11 : #include "MooseObject.h"
12 : #include "MooseEnum.h"
13 : #include "MooseUtils.h"
14 : #include "Numerics.h"
15 :
16 : InputParameters
17 11849 : DiscreteLineSegmentInterface::validParams()
18 : {
19 11849 : InputParameters params = emptyInputParameters();
20 :
21 23698 : params.addRequiredParam<Point>("position", "Start position of axis in 3-D space [m]");
22 23698 : params.addRequiredParam<RealVectorValue>(
23 : "orientation",
24 : "Direction of axis from start position to end position (no need to normalize)");
25 23698 : params.addParam<Real>("rotation", 0.0, "Angle of rotation about the x-axis [degrees]");
26 23698 : params.addRequiredParam<std::vector<Real>>("length", "Length of each axial section [m]");
27 23698 : params.addRequiredParam<std::vector<unsigned int>>("n_elems",
28 : "Number of elements in each axial section");
29 :
30 11849 : return params;
31 0 : }
32 :
33 5930 : DiscreteLineSegmentInterface::DiscreteLineSegmentInterface(const MooseObject * moose_object)
34 5930 : : _position(moose_object->parameters().get<Point>("position")),
35 5930 : _dir_unnormalized(moose_object->parameters().get<RealVectorValue>("orientation")),
36 5930 : _dir(initializeDirectionVector(_dir_unnormalized)),
37 5930 : _rotation(moose_object->parameters().get<Real>("rotation")),
38 5930 : _lengths(moose_object->parameters().get<std::vector<Real>>("length")),
39 5930 : _length(std::accumulate(_lengths.begin(), _lengths.end(), 0.0)),
40 5930 : _n_elems(moose_object->parameters().get<std::vector<unsigned int>>("n_elems")),
41 5930 : _n_elem(std::accumulate(_n_elems.begin(), _n_elems.end(), 0)),
42 5930 : _n_sections(_lengths.size()),
43 5930 : _section_end(_n_sections),
44 5930 : _x_centers(_n_elem),
45 5930 : _R(computeDirectionTransformationTensor(_dir)),
46 5930 : _Rx(computeXRotationTransformationTensor(_rotation)),
47 5930 : _R_inv(_R.inverse()),
48 5930 : _Rx_inv(_Rx.inverse()),
49 11860 : _moose_object_name_dlsi(moose_object->name())
50 : {
51 5930 : std::partial_sum(_lengths.begin(), _lengths.end(), _section_end.begin());
52 :
53 5930 : if (_lengths.size() != _n_elems.size())
54 2 : mooseError("The parameters 'length' and 'n_elems' must have the same number of entries.");
55 :
56 : // Compute the axial coordinates of the centers of each element
57 : unsigned int k_section_begin = 0;
58 : Real x_begin = 0.0;
59 12382 : for (unsigned int j = 0; j < _n_sections; j++)
60 : {
61 6454 : const Real dx = _lengths[j] / _n_elems[j];
62 148391 : for (unsigned int i = 0; i < _n_elems[j]; i++)
63 : {
64 141937 : const unsigned int k = k_section_begin + i;
65 141937 : _x_centers[k] = x_begin + 0.5 * dx;
66 141937 : x_begin += dx;
67 : }
68 6454 : k_section_begin += _n_elems[j];
69 : }
70 5928 : }
71 :
72 : RealVectorValue
73 5930 : DiscreteLineSegmentInterface::initializeDirectionVector(const RealVectorValue & dir_unnormalized)
74 : {
75 5930 : if (!MooseUtils::absoluteFuzzyEqual(dir_unnormalized.norm(), 0.0))
76 5930 : return dir_unnormalized / dir_unnormalized.norm();
77 : else
78 0 : mooseError("The parameter 'orientation' must not be the zero vector.");
79 : }
80 :
81 : RealTensorValue
82 5984 : DiscreteLineSegmentInterface::computeDirectionTransformationTensor(const RealVectorValue & dir)
83 : {
84 5984 : if (MooseUtils::absoluteFuzzyEqual(dir.norm(), 0.0))
85 0 : mooseError("The direction vector must not be the zero vector.");
86 :
87 5984 : const auto dir_normalized = dir / dir.norm();
88 5984 : const Real theta = acos(dir_normalized(2));
89 5984 : const Real aphi = atan2(dir_normalized(1), dir_normalized(0));
90 : return RealTensorValue(
91 5984 : RealVectorValue(cos(aphi) * sin(theta), -sin(aphi), -cos(aphi) * cos(theta)),
92 5984 : RealVectorValue(sin(aphi) * sin(theta), cos(aphi), -sin(aphi) * cos(theta)),
93 5984 : RealVectorValue(cos(theta), 0.0, sin(theta)));
94 : }
95 :
96 : RealTensorValue
97 5984 : DiscreteLineSegmentInterface::computeXRotationTransformationTensor(const Real & rotation)
98 : {
99 5984 : const Real rotation_rad = M_PI * rotation / 180.;
100 : const RealVectorValue Rx_x(1., 0., 0.);
101 5984 : const RealVectorValue Rx_y(0., cos(rotation_rad), -sin(rotation_rad));
102 : const RealVectorValue Rx_z(0., sin(rotation_rad), cos(rotation_rad));
103 5984 : return RealTensorValue(Rx_x, Rx_y, Rx_z);
104 : }
105 :
106 : Real
107 6086 : DiscreteLineSegmentInterface::computeAxialCoordinate(const Point & p) const
108 : {
109 6086 : const Real ax_coord = _dir * (p - _position);
110 6086 : if (MooseUtils::absoluteFuzzyLessThan(ax_coord, 0.0) ||
111 : MooseUtils::absoluteFuzzyGreaterThan(ax_coord, _length))
112 2 : mooseError(_moose_object_name_dlsi,
113 : ": The point ",
114 : p,
115 : " has an invalid axial coordinate (",
116 : ax_coord,
117 : "). Valid axial coordinates are in the range (0,",
118 2 : _length,
119 : ").");
120 : else
121 6084 : return ax_coord;
122 : }
123 :
124 : Real
125 1995 : DiscreteLineSegmentInterface::computeRadialCoordinate(const Point & p) const
126 : {
127 1995 : const RealVectorValue v = (p - _position);
128 1995 : return v.cross(_dir).norm();
129 : }
130 :
131 : unsigned int
132 1995 : DiscreteLineSegmentInterface::getAxialSectionIndex(const Point & p) const
133 : {
134 1995 : const Real axial_coordinate = computeAxialCoordinate(p);
135 4465 : for (unsigned int i = 0; i < _n_sections; ++i)
136 4465 : if (MooseUtils::absoluteFuzzyLessEqual(axial_coordinate, _section_end[i]))
137 1995 : return i;
138 :
139 0 : mooseError("No axial section index was found.");
140 : }
141 :
142 : unsigned int
143 1600 : DiscreteLineSegmentInterface::getAxialElementIndex(const Point & p_center) const
144 : {
145 1600 : const Real axial_coordinate = computeAxialCoordinate(p_center);
146 16800 : for (unsigned int i = 0; i < _n_elem; ++i)
147 16800 : if (MooseUtils::absoluteFuzzyEqual(axial_coordinate, _x_centers[i]))
148 1600 : return i;
149 :
150 0 : mooseError("No axial element index was found.");
151 : }
152 :
153 : Point
154 374114 : DiscreteLineSegmentInterface::computeRealPointFromReferencePoint(const Point & p,
155 : const RealVectorValue & position,
156 : const RealTensorValue & R,
157 : const RealTensorValue & Rx)
158 : {
159 374114 : return R * (Rx * p) + position;
160 : }
161 :
162 : Point
163 373520 : DiscreteLineSegmentInterface::computeRealPointFromReferencePoint(const Point & p) const
164 : {
165 373520 : return computeRealPointFromReferencePoint(p, _position, _R, _Rx);
166 : }
167 :
168 : Point
169 0 : DiscreteLineSegmentInterface::computeReferencePointFromRealPoint(const Point & p,
170 : const RealVectorValue & position,
171 : const RealTensorValue & R_inv,
172 : const RealTensorValue & Rx_inv)
173 : {
174 0 : return Rx_inv * R_inv * (p - position);
175 : }
176 :
177 : Point
178 0 : DiscreteLineSegmentInterface::computeReferencePointFromRealPoint(const Point & p) const
179 : {
180 0 : return computeReferencePointFromRealPoint(p, _position, _R_inv, _Rx_inv);
181 : }
182 :
183 : MooseEnum
184 272 : DiscreteLineSegmentInterface::getAlignmentAxis(const RealVectorValue & dir)
185 : {
186 544 : MooseEnum axis("x y z");
187 272 : if (THM::areParallelVectors(dir, RealVectorValue(1, 0, 0)))
188 0 : axis = "x";
189 272 : else if (THM::areParallelVectors(dir, RealVectorValue(0, 1, 0)))
190 540 : axis = "y";
191 2 : else if (THM::areParallelVectors(dir, RealVectorValue(0, 0, 1)))
192 0 : axis = "z";
193 :
194 272 : return axis;
195 0 : }
196 :
197 : MooseEnum
198 0 : DiscreteLineSegmentInterface::getAlignmentAxis() const
199 : {
200 0 : return getAlignmentAxis(_dir);
201 : }
202 :
203 : std::vector<Real>
204 54 : DiscreteLineSegmentInterface::getElementBoundaryCoordinates(
205 : const RealVectorValue & position,
206 : const RealVectorValue & orientation,
207 : const Real & rotation,
208 : const std::vector<Real> & lengths,
209 : const std::vector<unsigned int> & n_elems)
210 : {
211 54 : const auto axis = getAlignmentAxis(orientation);
212 :
213 : unsigned int d;
214 54 : if (axis == "x")
215 : d = 0;
216 54 : else if (axis == "y")
217 : d = 1;
218 0 : else if (axis == "z")
219 : d = 2;
220 : else
221 0 : mooseError("Invalid axis.");
222 :
223 54 : const auto R = computeDirectionTransformationTensor(orientation);
224 54 : const auto Rx = computeXRotationTransformationTensor(rotation);
225 :
226 54 : const unsigned int n_elems_total = std::accumulate(n_elems.begin(), n_elems.end(), 0);
227 54 : const unsigned int n_sections = lengths.size();
228 :
229 : unsigned int i_section_start = 0;
230 : Real section_start_ref_position = 0.0;
231 54 : std::vector<Real> element_boundary_coordinates(n_elems_total + 1);
232 54 : element_boundary_coordinates[0] =
233 54 : computeRealPointFromReferencePoint(Point(0, 0, 0), position, R, Rx)(d);
234 108 : for (unsigned int j = 0; j < n_sections; ++j)
235 : {
236 54 : const Real section_dx = lengths[j] / n_elems[j];
237 594 : for (unsigned int k = 0; k < n_elems[j]; ++k)
238 : {
239 540 : const unsigned int i = i_section_start + k + 1;
240 540 : const Real ref_position = section_start_ref_position + section_dx * k;
241 : const Point ref_point = Point(ref_position, 0, 0);
242 540 : element_boundary_coordinates[i] =
243 540 : computeRealPointFromReferencePoint(ref_point, position, R, Rx)(d);
244 : }
245 :
246 54 : section_start_ref_position += lengths[j];
247 54 : i_section_start += n_elems[j];
248 : }
249 :
250 54 : return element_boundary_coordinates;
251 54 : }
252 :
253 : std::vector<Real>
254 0 : DiscreteLineSegmentInterface::getElementBoundaryCoordinates() const
255 : {
256 0 : return getElementBoundaryCoordinates(_position, _dir, _rotation, _lengths, _n_elems);
257 : }
|
