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 319 : SelfShadowSideUserObject::validParams()
19 : {
20 319 : InputParameters params = SideUserObject::validParams();
21 319 : params.addClassDescription("Compute the illumination status for a self shadowing sideset");
22 638 : 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 319 : return params;
29 0 : }
30 :
31 168 : SelfShadowSideUserObject::SelfShadowSideUserObject(const InputParameters & parameters)
32 : : SideUserObject(parameters),
33 336 : _dim(_mesh.dimension()),
34 336 : _raw_direction(coupledPostprocessors("illumination_flux")),
35 168 : _illumination_status(
36 504 : 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 168 : if (_dim != 2 && _dim != 3)
42 0 : mooseError("SelfShadowSideUserObject works only for 2 and 3 dimensional problems.");
43 :
44 : // fetch raw direction
45 672 : for (const auto i : index_range(_raw_direction))
46 504 : _raw_direction[i] = &getPostprocessorValue("illumination_flux", i);
47 168 : }
48 :
49 : void
50 152 : SelfShadowSideUserObject::initialize()
51 : {
52 : // delete list of collected lines or triangles
53 152 : _lines.clear();
54 152 : _triangles.clear();
55 :
56 : // delete local qps
57 152 : _local_qps.clear();
58 :
59 : // update direction and rotation matrix
60 : RealVectorValue direction;
61 608 : for (const auto i : index_range(_raw_direction))
62 456 : direction(i) = *_raw_direction[i];
63 152 : _rotation =
64 152 : _dim == 2 ? RotationMatrix::rotVec2DToX(direction) : RotationMatrix::rotVecToZ(direction);
65 152 : }
66 :
67 : void
68 64152 : SelfShadowSideUserObject::execute()
69 : {
70 64152 : const SideIDType id(_current_elem->id(), _current_side);
71 :
72 : // add triangulated sides to list
73 64152 : if (_dim == 2)
74 2664 : addLines(id);
75 : else
76 61488 : addTriangles(id);
77 :
78 : // save off rotated QP coordinates got local sides
79 64152 : _local_qps.emplace_back(id, _q_point);
80 64152 : rotate(_local_qps.back().second);
81 64152 : }
82 :
83 : void
84 17 : SelfShadowSideUserObject::threadJoin(const UserObject & y)
85 : {
86 : const auto & uo = static_cast<const SelfShadowSideUserObject &>(y);
87 :
88 : // merge lists
89 17 : _lines.insert(_lines.end(), uo._lines.begin(), uo._lines.end());
90 17 : _triangles.insert(_triangles.end(), uo._triangles.begin(), uo._triangles.end());
91 17 : _local_qps.insert(_local_qps.end(), uo._local_qps.begin(), uo._local_qps.end());
92 17 : }
93 :
94 : unsigned int
95 1674112 : SelfShadowSideUserObject::illumination(const SideIDType & id) const
96 : {
97 1674112 : const auto it = _illumination_status.find(id);
98 1674112 : if (it == _illumination_status.end())
99 0 : mooseError("Illumination status was not calculated for current side.");
100 1674112 : return it->second;
101 : }
102 :
103 : void
104 135 : SelfShadowSideUserObject::finalize()
105 : {
106 : // rotate triangulations
107 135 : if (_dim == 2)
108 3928 : for (auto & line : _lines)
109 3864 : rotate(line);
110 : else
111 209303 : for (auto & triangle : _triangles)
112 209232 : rotate(triangle);
113 :
114 : // [compute local projected bounding box (in x or xy)]
115 :
116 : // communicate
117 135 : if (_dim == 2)
118 64 : _communicator.allgather(_lines, /*identical_buffer_sizes=*/false);
119 : else
120 71 : _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 64287 : for (const auto & [id, qps] : _local_qps)
124 : {
125 64152 : auto & illumination = _illumination_status[id];
126 :
127 : // start off assuming no illumination
128 64152 : illumination = 0;
129 :
130 : // current bit in the illumination mask
131 : unsigned int bit = 1;
132 :
133 : // iterate over QPs
134 426684 : 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 362532 : if (_dim == 2 && check2DIllumination(qps[i], id))
138 1986 : illumination |= bit;
139 362532 : if (_dim == 3 && check3DIllumination(qps[i], id))
140 124692 : illumination |= bit;
141 :
142 : // shift to next bit
143 362532 : bit <<= 1;
144 : }
145 : }
146 135 : }
147 :
148 : void
149 2664 : SelfShadowSideUserObject::addLines(const SideIDType & id)
150 : {
151 2664 : const auto & cse = *_current_side_elem;
152 2664 : switch (cse.type())
153 : {
154 1464 : case libMesh::EDGE2:
155 1464 : _lines.emplace_back(cse.node_ref(0), cse.node_ref(1), id);
156 1464 : break;
157 :
158 1200 : case libMesh::EDGE3:
159 1200 : _lines.emplace_back(cse.node_ref(0), cse.node_ref(2), id);
160 1200 : _lines.emplace_back(cse.node_ref(2), cse.node_ref(1), id);
161 1200 : 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 2664 : }
173 :
174 : void
175 61488 : SelfShadowSideUserObject::addTriangles(const SideIDType & id)
176 : {
177 61488 : const auto & cse = *_current_side_elem;
178 61488 : switch (cse.type())
179 : {
180 22392 : case libMesh::TRI3:
181 22392 : _triangles.emplace_back(cse.node_ref(0), cse.node_ref(1), cse.node_ref(2), id);
182 22392 : break;
183 :
184 1104 : case libMesh::TRI6:
185 : case libMesh::TRI7:
186 1104 : _triangles.emplace_back(cse.node_ref(0), cse.node_ref(3), cse.node_ref(5), id);
187 1104 : _triangles.emplace_back(cse.node_ref(3), cse.node_ref(1), cse.node_ref(4), id);
188 1104 : _triangles.emplace_back(cse.node_ref(4), cse.node_ref(2), cse.node_ref(5), id);
189 1104 : _triangles.emplace_back(cse.node_ref(3), cse.node_ref(4), cse.node_ref(5), id);
190 1104 : break;
191 :
192 16704 : case libMesh::QUAD4:
193 16704 : _triangles.emplace_back(cse.node_ref(0), cse.node_ref(1), cse.node_ref(2), id);
194 16704 : _triangles.emplace_back(cse.node_ref(2), cse.node_ref(3), cse.node_ref(0), id);
195 16704 : break;
196 :
197 10644 : case libMesh::QUAD8:
198 10644 : _triangles.emplace_back(cse.node_ref(0), cse.node_ref(4), cse.node_ref(7), id);
199 10644 : _triangles.emplace_back(cse.node_ref(4), cse.node_ref(1), cse.node_ref(5), id);
200 10644 : _triangles.emplace_back(cse.node_ref(5), cse.node_ref(2), cse.node_ref(6), id);
201 10644 : _triangles.emplace_back(cse.node_ref(6), cse.node_ref(3), cse.node_ref(7), id);
202 10644 : _triangles.emplace_back(cse.node_ref(6), cse.node_ref(7), cse.node_ref(4), id);
203 10644 : _triangles.emplace_back(cse.node_ref(4), cse.node_ref(5), cse.node_ref(6), id);
204 10644 : break;
205 :
206 10644 : case libMesh::QUAD9:
207 10644 : _triangles.emplace_back(cse.node_ref(0), cse.node_ref(4), cse.node_ref(7), id);
208 10644 : _triangles.emplace_back(cse.node_ref(4), cse.node_ref(1), cse.node_ref(5), id);
209 10644 : _triangles.emplace_back(cse.node_ref(5), cse.node_ref(2), cse.node_ref(6), id);
210 10644 : _triangles.emplace_back(cse.node_ref(6), cse.node_ref(3), cse.node_ref(7), id);
211 10644 : _triangles.emplace_back(cse.node_ref(8), cse.node_ref(6), cse.node_ref(7), id);
212 10644 : _triangles.emplace_back(cse.node_ref(8), cse.node_ref(5), cse.node_ref(6), id);
213 10644 : _triangles.emplace_back(cse.node_ref(8), cse.node_ref(4), cse.node_ref(5), id);
214 10644 : _triangles.emplace_back(cse.node_ref(8), cse.node_ref(7), cse.node_ref(4), id);
215 10644 : break;
216 :
217 0 : default:
218 0 : mooseError("Unsupported FACE type");
219 : }
220 61488 : }
221 :
222 : bool
223 6828 : SelfShadowSideUserObject::check2DIllumination(const Point & qp, const SideIDType & id)
224 : {
225 6828 : const auto x = qp(0);
226 6828 : const auto y = qp(1);
227 :
228 : // loop over all line segments until one is found that provides shade
229 318709 : for (const auto & line : _lines)
230 : {
231 : const auto & [p1, p2, line_id] = line;
232 :
233 : // make sure a side never shades itself
234 7154 : if (line_id == id)
235 7154 : continue;
236 :
237 309569 : const bool ordered = p1(1) <= p2(1);
238 :
239 309569 : const auto y1 = ordered ? p1(1) : p2(1);
240 309569 : const auto y2 = ordered ? p2(1) : p1(1);
241 :
242 309569 : if (y >= y1 && y <= y2)
243 : {
244 : // line segment is oriented parallel to the irradiation direction
245 9622 : 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 9520 : const auto x1 = ordered ? p1(0) : p2(0);
250 9520 : const auto x2 = ordered ? p2(0) : p1(0);
251 :
252 : // compute intersection location
253 9520 : const auto xs = (x2 - x1) * (y - y1) / (y2 - y1) + x1;
254 9520 : if (x > xs)
255 : return false;
256 : }
257 : }
258 :
259 : return true;
260 : }
261 :
262 : bool
263 355704 : SelfShadowSideUserObject::check3DIllumination(const Point & qp, const SideIDType & id)
264 : {
265 355704 : const auto x = qp(0);
266 355704 : const auto y = qp(1);
267 355704 : const auto z = qp(2);
268 :
269 : // loop over all triangles until one is found that provides shade
270 1018563920 : 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 939460 : if (triangle_id == id)
276 939460 : continue;
277 :
278 1017499768 : const auto x1 = p1(0);
279 1017499768 : const auto x2 = p2(0);
280 1017499768 : const auto x3 = p3(0);
281 :
282 1017499768 : const auto y1 = p1(1);
283 1017499768 : const auto y2 = p2(1);
284 1017499768 : const auto y3 = p3(1);
285 :
286 : // compute barycentric coordinates
287 1017499768 : const auto a = ((y2 - y3) * (x - x3) + (x3 - x2) * (y - y3)) /
288 1017499768 : ((y2 - y3) * (x1 - x3) + (x3 - x2) * (y1 - y3));
289 1017499768 : const auto b = ((y3 - y1) * (x - x3) + (x1 - x3) * (y - y3)) /
290 : ((y2 - y3) * (x1 - x3) + (x3 - x2) * (y1 - y3));
291 1017499768 : const auto c = 1.0 - a - b;
292 :
293 : // are we inside of the projected triangle?
294 1017499768 : 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 420497 : const auto zs = a * p1(2) + b * p2(2) + c * p3(2);
299 420497 : if (z > zs)
300 : return false;
301 : }
302 : }
303 :
304 : return true;
305 : }
306 :
307 : void
308 64152 : SelfShadowSideUserObject::rotate(MooseArray<Point> & points)
309 : {
310 426684 : for (const auto i : index_range(points))
311 362532 : points[i] = _rotation * points[i];
312 64152 : }
313 :
314 : void
315 209232 : SelfShadowSideUserObject::rotate(Triangle & triangle)
316 : {
317 : auto & [p1, p2, p3, triangle_id] = triangle;
318 :
319 209232 : p1 = _rotation * p1;
320 209232 : p2 = _rotation * p2;
321 209232 : p3 = _rotation * p3;
322 : libmesh_ignore(triangle_id);
323 209232 : }
324 :
325 : void
326 3864 : SelfShadowSideUserObject::rotate(LineSegment & line)
327 : {
328 : auto & [p1, p2, line_id] = line;
329 3864 : p1 = _rotation * p1;
330 3864 : p2 = _rotation * p2;
331 : libmesh_ignore(line_id);
332 3864 : }
|
