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 "SelfShadowSideUserObject.h"
11 : #include "RotationMatrix.h"
12 :
13 : #include "libmesh/parallel_algebra.h"
14 :
15 : registerMooseObject("HeatTransferApp", SelfShadowSideUserObject);
16 :
17 : InputParameters
18 677 : SelfShadowSideUserObject::validParams()
19 : {
20 677 : InputParameters params = SideUserObject::validParams();
21 677 : params.addClassDescription("Compute the illumination status for a self shadowing sideset");
22 1354 : params.addRequiredRangeCheckedParam<std::vector<PostprocessorName>>(
23 : "illumination_flux",
24 : "illumination_flux_size>0 && illumination_flux_size<=3",
25 : "Radiation direction vector. Each component of the vector can be a constant number or a "
26 : "postprocessor name (the latter enables time varying radiation directions - spatial "
27 : "variation is not supported)");
28 677 : return params;
29 0 : }
30 :
31 364 : SelfShadowSideUserObject::SelfShadowSideUserObject(const InputParameters & parameters)
32 : : SideUserObject(parameters),
33 728 : _dim(_mesh.dimension()),
34 728 : _raw_direction(coupledPostprocessors("illumination_flux")),
35 364 : _illumination_status(
36 1092 : declareRestartableData<std::map<SideIDType, unsigned int>>("illumination_status"))
37 : {
38 : // we should check the coordinate system (i.e. permit only 0,0,+-1 for RZ)
39 :
40 : // check problem dimension
41 364 : if (_dim != 2 && _dim != 3)
42 0 : mooseError("SelfShadowSideUserObject works only for 2 and 3 dimensional problems.");
43 :
44 : // fetch raw direction
45 1456 : for (const auto i : index_range(_raw_direction))
46 1092 : _raw_direction[i] = &getPostprocessorValue("illumination_flux", i);
47 364 : }
48 :
49 : void
50 284 : SelfShadowSideUserObject::initialize()
51 : {
52 : // delete list of collected lines or triangles
53 284 : _lines.clear();
54 284 : _triangles.clear();
55 :
56 : // delete local qps
57 284 : _local_qps.clear();
58 :
59 : // update direction and rotation matrix
60 : RealVectorValue direction;
61 1136 : for (const auto i : index_range(_raw_direction))
62 852 : direction(i) = *_raw_direction[i];
63 284 : _rotation =
64 284 : _dim == 2 ? RotationMatrix::rotVec2DToX(direction) : RotationMatrix::rotVecToZ(direction);
65 284 : }
66 :
67 : void
68 96048 : SelfShadowSideUserObject::execute()
69 : {
70 96048 : const SideIDType id(_current_elem->id(), _current_side);
71 :
72 : // add triangulated sides to list
73 96048 : if (_dim == 2)
74 3996 : addLines(id);
75 : else
76 92052 : addTriangles(id);
77 :
78 : // save off rotated QP coordinates got local sides
79 96048 : _local_qps.emplace_back(id, _q_point);
80 96048 : rotate(_local_qps.back().second);
81 96048 : }
82 :
83 : void
84 51 : SelfShadowSideUserObject::threadJoin(const UserObject & y)
85 : {
86 : const auto & uo = static_cast<const SelfShadowSideUserObject &>(y);
87 :
88 : // merge lists
89 51 : _lines.insert(_lines.end(), uo._lines.begin(), uo._lines.end());
90 51 : _triangles.insert(_triangles.end(), uo._triangles.begin(), uo._triangles.end());
91 51 : _local_qps.insert(_local_qps.end(), uo._local_qps.begin(), uo._local_qps.end());
92 51 : }
93 :
94 : unsigned int
95 2458476 : SelfShadowSideUserObject::illumination(const SideIDType & id) const
96 : {
97 2458476 : const auto it = _illumination_status.find(id);
98 2458476 : if (it == _illumination_status.end())
99 0 : mooseError("Illumination status was not calculated for current side.");
100 2458476 : return it->second;
101 : }
102 :
103 : void
104 233 : SelfShadowSideUserObject::finalize()
105 : {
106 : // rotate triangulations
107 233 : if (_dim == 2)
108 5908 : for (auto & line : _lines)
109 5796 : rotate(line);
110 : else
111 313609 : for (auto & triangle : _triangles)
112 313488 : rotate(triangle);
113 :
114 : // [compute local projected bounding box (in x or xy)]
115 :
116 : // communicate
117 233 : if (_dim == 2)
118 112 : _communicator.allgather(_lines, /*identical_buffer_sizes=*/false);
119 : else
120 121 : _communicator.allgather(_triangles, /*identical_buffer_sizes=*/false);
121 :
122 : // otherwise we iterate over QPs and check if any other side is in the way of the radiation
123 96281 : for (const auto & [id, qps] : _local_qps)
124 : {
125 96048 : auto & illumination = _illumination_status[id];
126 :
127 : // start off assuming no illumination
128 96048 : illumination = 0;
129 :
130 : // current bit in the illumination mask
131 : unsigned int bit = 1;
132 :
133 : // iterate over QPs
134 639126 : for (const auto i : index_range(qps))
135 : {
136 : // we set the bit at position _qp in the bitmask if the QP is illuminated
137 543078 : if (_dim == 2 && check2DIllumination(qps[i], id))
138 2979 : illumination |= bit;
139 543078 : if (_dim == 3 && check3DIllumination(qps[i], id))
140 186696 : illumination |= bit;
141 :
142 : // shift to next bit
143 543078 : bit <<= 1;
144 : }
145 : }
146 233 : }
147 :
148 : void
149 3996 : SelfShadowSideUserObject::addLines(const SideIDType & id)
150 : {
151 3996 : const auto & cse = *_current_side_elem;
152 3996 : switch (cse.type())
153 : {
154 2196 : case libMesh::EDGE2:
155 2196 : _lines.emplace_back(cse.node_ref(0), cse.node_ref(1), id);
156 2196 : break;
157 :
158 1800 : case libMesh::EDGE3:
159 1800 : _lines.emplace_back(cse.node_ref(0), cse.node_ref(2), id);
160 1800 : _lines.emplace_back(cse.node_ref(2), cse.node_ref(1), id);
161 1800 : break;
162 :
163 0 : case libMesh::EDGE4:
164 0 : _lines.emplace_back(cse.node_ref(0), cse.node_ref(2), id);
165 0 : _lines.emplace_back(cse.node_ref(2), cse.node_ref(3), id);
166 0 : _lines.emplace_back(cse.node_ref(3), cse.node_ref(1), id);
167 0 : break;
168 :
169 0 : default:
170 0 : mooseError("Unsupported EDGE type");
171 : }
172 3996 : }
173 :
174 : void
175 92052 : SelfShadowSideUserObject::addTriangles(const SideIDType & id)
176 : {
177 92052 : const auto & cse = *_current_side_elem;
178 92052 : switch (cse.type())
179 : {
180 33588 : case libMesh::TRI3:
181 33588 : _triangles.emplace_back(cse.node_ref(0), cse.node_ref(1), cse.node_ref(2), id);
182 33588 : break;
183 :
184 1656 : case libMesh::TRI6:
185 : case libMesh::TRI7:
186 1656 : _triangles.emplace_back(cse.node_ref(0), cse.node_ref(3), cse.node_ref(5), id);
187 1656 : _triangles.emplace_back(cse.node_ref(3), cse.node_ref(1), cse.node_ref(4), id);
188 1656 : _triangles.emplace_back(cse.node_ref(4), cse.node_ref(2), cse.node_ref(5), id);
189 1656 : _triangles.emplace_back(cse.node_ref(3), cse.node_ref(4), cse.node_ref(5), id);
190 1656 : break;
191 :
192 24876 : case libMesh::QUAD4:
193 24876 : _triangles.emplace_back(cse.node_ref(0), cse.node_ref(1), cse.node_ref(2), id);
194 24876 : _triangles.emplace_back(cse.node_ref(2), cse.node_ref(3), cse.node_ref(0), id);
195 24876 : break;
196 :
197 15966 : case libMesh::QUAD8:
198 15966 : _triangles.emplace_back(cse.node_ref(0), cse.node_ref(4), cse.node_ref(7), id);
199 15966 : _triangles.emplace_back(cse.node_ref(4), cse.node_ref(1), cse.node_ref(5), id);
200 15966 : _triangles.emplace_back(cse.node_ref(5), cse.node_ref(2), cse.node_ref(6), id);
201 15966 : _triangles.emplace_back(cse.node_ref(6), cse.node_ref(3), cse.node_ref(7), id);
202 15966 : _triangles.emplace_back(cse.node_ref(6), cse.node_ref(7), cse.node_ref(4), id);
203 15966 : _triangles.emplace_back(cse.node_ref(4), cse.node_ref(5), cse.node_ref(6), id);
204 15966 : break;
205 :
206 15966 : case libMesh::QUAD9:
207 15966 : _triangles.emplace_back(cse.node_ref(0), cse.node_ref(4), cse.node_ref(7), id);
208 15966 : _triangles.emplace_back(cse.node_ref(4), cse.node_ref(1), cse.node_ref(5), id);
209 15966 : _triangles.emplace_back(cse.node_ref(5), cse.node_ref(2), cse.node_ref(6), id);
210 15966 : _triangles.emplace_back(cse.node_ref(6), cse.node_ref(3), cse.node_ref(7), id);
211 15966 : _triangles.emplace_back(cse.node_ref(8), cse.node_ref(6), cse.node_ref(7), id);
212 15966 : _triangles.emplace_back(cse.node_ref(8), cse.node_ref(5), cse.node_ref(6), id);
213 15966 : _triangles.emplace_back(cse.node_ref(8), cse.node_ref(4), cse.node_ref(5), id);
214 15966 : _triangles.emplace_back(cse.node_ref(8), cse.node_ref(7), cse.node_ref(4), id);
215 15966 : break;
216 :
217 0 : default:
218 0 : mooseError("Unsupported FACE type");
219 : }
220 92052 : }
221 :
222 : bool
223 10242 : SelfShadowSideUserObject::check2DIllumination(const Point & qp, const SideIDType & id)
224 : {
225 10242 : const auto x = qp(0);
226 10242 : const auto y = qp(1);
227 :
228 : // loop over all line segments until one is found that provides shade
229 483898 : for (const auto & line : _lines)
230 : {
231 : const auto & [p1, p2, line_id] = line;
232 :
233 : // make sure a side never shades itself
234 11078 : if (line_id == id)
235 11078 : continue;
236 :
237 469841 : const bool ordered = p1(1) <= p2(1);
238 :
239 469841 : const auto y1 = ordered ? p1(1) : p2(1);
240 469841 : const auto y2 = ordered ? p2(1) : p1(1);
241 :
242 469841 : if (y >= y1 && y <= y2)
243 : {
244 : // line segment is oriented parallel to the irradiation direction
245 14790 : if (std::abs(y2 - y1) < libMesh::TOLERANCE)
246 : return false;
247 :
248 : // segment is in line with the QP in radiation direction. Is it in front or behind?
249 14637 : const auto x1 = ordered ? p1(0) : p2(0);
250 14637 : const auto x2 = ordered ? p2(0) : p1(0);
251 :
252 : // compute intersection location
253 14637 : const auto xs = (x2 - x1) * (y - y1) / (y2 - y1) + x1;
254 14637 : if (x > xs)
255 : return false;
256 : }
257 : }
258 :
259 : return true;
260 : }
261 :
262 : bool
263 532836 : SelfShadowSideUserObject::check3DIllumination(const Point & qp, const SideIDType & id)
264 : {
265 532836 : const auto x = qp(0);
266 532836 : const auto y = qp(1);
267 532836 : const auto z = qp(2);
268 :
269 : // loop over all triangles until one is found that provides shade
270 1535064876 : for (const auto & triangle : _triangles)
271 : {
272 : const auto & [p1, p2, p3, triangle_id] = triangle;
273 :
274 : // make sure a side never shades itself
275 1455509 : if (triangle_id == id)
276 1455509 : continue;
277 :
278 1533422671 : const auto x1 = p1(0);
279 1533422671 : const auto x2 = p2(0);
280 1533422671 : const auto x3 = p3(0);
281 :
282 1533422671 : const auto y1 = p1(1);
283 1533422671 : const auto y2 = p2(1);
284 1533422671 : const auto y3 = p3(1);
285 :
286 : // compute barycentric coordinates
287 1533422671 : const auto a = ((y2 - y3) * (x - x3) + (x3 - x2) * (y - y3)) /
288 1533422671 : ((y2 - y3) * (x1 - x3) + (x3 - x2) * (y1 - y3));
289 1533422671 : const auto b = ((y3 - y1) * (x - x3) + (x1 - x3) * (y - y3)) /
290 : ((y2 - y3) * (x1 - x3) + (x3 - x2) * (y1 - y3));
291 1533422671 : const auto c = 1.0 - a - b;
292 :
293 : // are we inside of the projected triangle?
294 1533422671 : if (0 <= a && a <= 1 && 0 <= b && b <= 1 && 0 <= c && c <= 1)
295 : {
296 : // Is the intersection it in front or behind the QP? (interpolate z using the barycentric
297 : // coordinates)
298 638039 : const auto zs = a * p1(2) + b * p2(2) + c * p3(2);
299 638039 : if (z > zs)
300 : return false;
301 : }
302 : }
303 :
304 : return true;
305 : }
306 :
307 : void
308 96048 : SelfShadowSideUserObject::rotate(MooseArray<Point> & points)
309 : {
310 639126 : for (const auto i : index_range(points))
311 543078 : points[i] = _rotation * points[i];
312 96048 : }
313 :
314 : void
315 313488 : SelfShadowSideUserObject::rotate(Triangle & triangle)
316 : {
317 : auto & [p1, p2, p3, triangle_id] = triangle;
318 :
319 313488 : p1 = _rotation * p1;
320 313488 : p2 = _rotation * p2;
321 313488 : p3 = _rotation * p3;
322 : libmesh_ignore(triangle_id);
323 313488 : }
324 :
325 : void
326 5796 : SelfShadowSideUserObject::rotate(LineSegment & line)
327 : {
328 : auto & [p1, p2, line_id] = line;
329 5796 : p1 = _rotation * p1;
330 5796 : p2 = _rotation * p2;
331 : libmesh_ignore(line_id);
332 5796 : }
|
