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