Line data Source code
1 : // EnergyPlus, Copyright (c) 1996-2024, The Board of Trustees of the University of Illinois,
2 : // The Regents of the University of California, through Lawrence Berkeley National Laboratory
3 : // (subject to receipt of any required approvals from the U.S. Dept. of Energy), Oak Ridge
4 : // National Laboratory, managed by UT-Battelle, Alliance for Sustainable Energy, LLC, and other
5 : // contributors. All rights reserved.
6 : //
7 : // NOTICE: This Software was developed under funding from the U.S. Department of Energy and the
8 : // U.S. Government consequently retains certain rights. As such, the U.S. Government has been
9 : // granted for itself and others acting on its behalf a paid-up, nonexclusive, irrevocable,
10 : // worldwide license in the Software to reproduce, distribute copies to the public, prepare
11 : // derivative works, and perform publicly and display publicly, and to permit others to do so.
12 : //
13 : // Redistribution and use in source and binary forms, with or without modification, are permitted
14 : // provided that the following conditions are met:
15 : //
16 : // (1) Redistributions of source code must retain the above copyright notice, this list of
17 : // conditions and the following disclaimer.
18 : //
19 : // (2) Redistributions in binary form must reproduce the above copyright notice, this list of
20 : // conditions and the following disclaimer in the documentation and/or other materials
21 : // provided with the distribution.
22 : //
23 : // (3) Neither the name of the University of California, Lawrence Berkeley National Laboratory,
24 : // the University of Illinois, U.S. Dept. of Energy nor the names of its contributors may be
25 : // used to endorse or promote products derived from this software without specific prior
26 : // written permission.
27 : //
28 : // (4) Use of EnergyPlus(TM) Name. If Licensee (i) distributes the software in stand-alone form
29 : // without changes from the version obtained under this License, or (ii) Licensee makes a
30 : // reference solely to the software portion of its product, Licensee must refer to the
31 : // software as "EnergyPlus version X" software, where "X" is the version number Licensee
32 : // obtained under this License and may not use a different name for the software. Except as
33 : // specifically required in this Section (4), Licensee shall not use in a company name, a
34 : // product name, in advertising, publicity, or other promotional activities any name, trade
35 : // name, trademark, logo, or other designation of "EnergyPlus", "E+", "e+" or confusingly
36 : // similar designation, without the U.S. Department of Energy's prior written consent.
37 : //
38 : // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
39 : // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY
40 : // AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
41 : // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
42 : // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
43 : // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
44 : // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
45 : // OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
46 : // POSSIBILITY OF SUCH DAMAGE.
47 :
48 : #ifndef EnergyPlus_SurfaceOctree_hh_INCLUDED
49 : #define EnergyPlus_SurfaceOctree_hh_INCLUDED
50 :
51 : // EnergyPlus Headers
52 : #include <EnergyPlus/EPVector.hh>
53 : #include <EnergyPlus/EnergyPlus.hh>
54 :
55 : // ObjexxFCL Headers
56 : #include <ObjexxFCL/Array1.fwd.hh>
57 : #include <ObjexxFCL/Vector3.hh>
58 :
59 : // C++ Headers
60 : #include <cassert>
61 : #include <cstdint>
62 : #include <vector>
63 :
64 : namespace EnergyPlus {
65 :
66 : // Forward
67 : struct EnergyPlusData;
68 :
69 : namespace DataSurfaces {
70 : struct SurfaceData;
71 : }
72 :
73 : // Package: Surface Octree System
74 : //
75 : // Purpose: Spatial sort of surfaces for fast, scalable identification of active surfaces for algorithms
76 : // making spatial queries such as solar shading, solar reflection, and daylighting obstruction
77 : //
78 : // Author: Stuart Mentzer (Stuart_Mentzer@objexx.com)
79 : //
80 : // History:
81 : // Sep 2015: Experimental code
82 : // Jan 2016: Initial release
83 : //
84 : // Notes: See the .cc file
85 :
86 : class SurfaceOctreeCube
87 : {
88 :
89 : public: // Types
90 : using Real = Real64;
91 : using Surface = DataSurfaces::SurfaceData;
92 : using Vertex = ObjexxFCL::Vector3<Real>;
93 : using Surfaces = std::vector<Surface *>;
94 : using size_type = Surfaces::size_type;
95 :
96 : public: // Creation
97 : // Default Constructor
98 796 : SurfaceOctreeCube() : d_(0u), n_(0u), l_(Vertex(0.0)), u_(Vertex(0.0)), c_(Vertex(0.0)), w_(0.0), r_(0.0)
99 : {
100 7164 : for (std::uint8_t i = 0; i < 8; ++i)
101 6368 : cubes_[i] = nullptr; // VC++ 2013 compatible initialization
102 796 : }
103 :
104 : // Surfaces Outer Cube Constructor
105 : SurfaceOctreeCube(EPVector<Surface> &surfaces) : d_(0u), n_(0u), l_(Vertex(0.0)), u_(Vertex(0.0)), c_(Vertex(0.0)), w_(0.0), r_(0.0)
106 : {
107 : for (std::uint8_t i = 0; i < 8; ++i)
108 : cubes_[i] = nullptr; // VC++ 2013 compatible initialization
109 : init(surfaces);
110 : }
111 :
112 : // Box Constructor
113 480 : SurfaceOctreeCube(std::uint8_t const d, Vertex const &l, Vertex const &u, Real const w)
114 480 : : d_(d), n_(0u), l_(l), u_(u), c_(cen(l, u)), w_(w), r_(0.75 * (w * w))
115 : {
116 4320 : for (std::uint8_t i = 0; i < 8; ++i)
117 3840 : cubes_[i] = nullptr; // VC++ 2013 compatible initialization
118 480 : assert(valid());
119 480 : }
120 :
121 : // Destructor
122 1276 : ~SurfaceOctreeCube()
123 : {
124 1756 : for (std::uint8_t i = 0; i < n_; ++i)
125 480 : delete cubes_[i];
126 1276 : }
127 :
128 : public: // Properties
129 : // Depth
130 : std::uint8_t d() const
131 : {
132 : return d_;
133 : }
134 :
135 : // Depth
136 : std::uint8_t depth() const
137 : {
138 : return d_;
139 : }
140 :
141 : // Number of Children
142 : std::uint8_t nChildren() const
143 : {
144 : return n_;
145 : }
146 :
147 : // Number of Sub-Cubes
148 : std::uint8_t nSubCube() const
149 : {
150 : return n_;
151 : }
152 :
153 : // Lower Corner
154 : Vertex const &l() const
155 : {
156 : return l_;
157 : }
158 :
159 : // Upper Corner
160 : Vertex const &u() const
161 : {
162 : return u_;
163 : }
164 :
165 : // Center Point
166 : Vertex const &c() const
167 : {
168 : return c_;
169 : }
170 :
171 : // Center Point
172 : Vertex const ¢er() const
173 : {
174 : return c_;
175 : }
176 :
177 : // Width
178 : Real w() const
179 : {
180 : return w_;
181 : }
182 :
183 : // Width
184 : Real width() const
185 : {
186 : return w_;
187 : }
188 :
189 : // Surfaces
190 : Surfaces const &surfaces() const
191 : {
192 : return surfaces_;
193 : }
194 :
195 : // Surfaces
196 : Surfaces::size_type surfaces_size() const
197 : {
198 : return surfaces_.size();
199 : }
200 :
201 : // Surfaces Begin Iterator
202 : Surfaces::const_iterator surfaces_begin() const
203 : {
204 : return surfaces_.begin();
205 : }
206 :
207 : // Surfaces Begin Iterator
208 : Surfaces::iterator surfaces_begin()
209 : {
210 : return surfaces_.begin();
211 : }
212 :
213 : // Surfaces End Iterator
214 : Surfaces::const_iterator surfaces_end() const
215 : {
216 : return surfaces_.end();
217 : }
218 :
219 : // Surfaces End Iterator
220 : Surfaces::iterator surfaces_end()
221 : {
222 : return surfaces_.end();
223 : }
224 :
225 : public: // Methods
226 : // Vertex in Cube?
227 40428024 : bool contains(Vertex const &v) const
228 : {
229 40428024 : return (l_.x <= v.x) && (v.x <= u_.x) && (l_.y <= v.y) && (v.y <= u_.y) && (l_.z <= v.z) && (v.z <= u_.z);
230 : }
231 :
232 : // Surface in Cube?
233 : bool contains(Surface const &surface) const;
234 :
235 : // Line Segment Intersects Enclosing Sphere?
236 : bool segmentIntersectsSphere(Vertex const &a, Vertex const &b) const
237 : {
238 : Vertex const ab(b - a);
239 : Real const ab_mag_squared(ab.mag_squared());
240 : if (ab_mag_squared == 0.0) { // Segment is a point
241 : return ObjexxFCL::distance_squared(a, c_) <= r_;
242 : } else { // Might pay to check if a or b in sphere first in some applications
243 : Vertex const ac(c_ - a);
244 : Real const projection_fac(((ac.x * ab.x) + (ac.y * ab.y) + (ac.z * ab.z)) / ab_mag_squared);
245 : if ((0.0 <= projection_fac) && (projection_fac <= 1.0)) { // Projected (closest) point is on ab segment
246 : return ObjexxFCL::distance_squared(ac, projection_fac * ab) <= r_;
247 : } else { // Projection (closest) point is outside of ab segment: Intersects iff a or b are in sphere
248 : return (ObjexxFCL::distance_squared(a, c_) <= r_) || (ObjexxFCL::distance_squared(b, c_) <= r_);
249 : }
250 : }
251 : }
252 :
253 : // Ray Intersects Enclosing Sphere?
254 : bool rayIntersectsSphere(Vertex const &a, Vertex const &dir) const
255 : {
256 : assert(std::abs(dir.mag_squared() - 1.0) < 4 * std::numeric_limits<Real>::epsilon()); // Check unit vector
257 : // Might pay to check if a in sphere first in some applications
258 : Vertex const ac(c_ - a);
259 : Real const projection_fac((ac.x * dir.x) + (ac.y * dir.y) + (ac.z * dir.z));
260 : if (0.0 <= projection_fac) { // Projected (closest) point is on ray
261 : return ObjexxFCL::distance_squared(ac, projection_fac * dir) <= r_;
262 : } else { // Projection (closest) point is outside of ray: Intersects iff a is in sphere
263 : return ObjexxFCL::distance_squared(a, c_) <= r_;
264 : }
265 : }
266 :
267 : // Line Intersects Enclosing Sphere?
268 : bool lineIntersectsSphere(Vertex const &a, Vertex const &dir) const
269 : {
270 : assert(std::abs(dir.mag_squared() - 1.0) < 4 * std::numeric_limits<Real>::epsilon()); // Check unit vector
271 : Vertex const ac(c_ - a);
272 : return ac.mag_squared() - ObjexxFCL::square((ac.x * dir.x) + (ac.y * dir.y) + (ac.z * dir.z)) <= r_;
273 : }
274 :
275 : // Line Segment Intersects Cube?
276 5853602 : bool segmentIntersectsCube(Vertex const &a, Vertex const &b) const
277 : {
278 : // Check if either segment endpoint is in cube
279 5853602 : if (contains(a) || contains(b)) return true;
280 :
281 : // Use separating axis theorem (faster variants exist)
282 4193948 : Vertex const m(mid(a, b) - c_); // Mid-point relative to cube center
283 4193948 : Vertex const mb(b - c_ - m); // ab mid-point to b half segment vector
284 4193948 : Vertex const e(std::abs(mb.x), std::abs(mb.y), std::abs(mb.z)); // Extent of half ab segment
285 4193948 : Real const h(0.5 * w_); // Half-width
286 : // Check if x,y,z axes are separating
287 4193948 : if (std::abs(m.x) > h + e.x) return false;
288 2141194 : if (std::abs(m.y) > h + e.y) return false;
289 707352 : if (std::abs(m.z) > h + e.z) return false;
290 : // Check if cross product axes are separating
291 136029 : if (std::abs(m.y * mb.z - m.z * mb.y) > h * (e.z + e.y)) return false;
292 115864 : if (std::abs(m.x * mb.z - m.z * mb.x) > h * (e.z + e.x)) return false;
293 88814 : if (std::abs(m.x * mb.y - m.y * mb.x) > h * (e.y + e.x)) return false;
294 87565 : return true;
295 4193948 : }
296 :
297 : // Ray Intersects Cube?
298 30003993 : bool rayIntersectsCube(Vertex const &a, Vertex const &dir, Vertex const &dir_inv) const
299 : {
300 : // Note: dir_inv coordinates corresponding to a zero dir coordinate are not used and can be set to zero
301 :
302 30003993 : assert(std::abs(dir.mag_squared() - 1.0) < 4 * std::numeric_limits<Real>::epsilon()); // Check unit vector
303 30003993 : assert((dir.x == 0.0) || (std::abs(dir_inv.x - (1.0 / dir.x)) < 2 * std::numeric_limits<Real>::epsilon() * std::abs(dir_inv.x)));
304 30003993 : assert((dir.y == 0.0) || (std::abs(dir_inv.y - (1.0 / dir.y)) < 2 * std::numeric_limits<Real>::epsilon() * std::abs(dir_inv.y)));
305 30003993 : assert((dir.z == 0.0) || (std::abs(dir_inv.z - (1.0 / dir.z)) < 2 * std::numeric_limits<Real>::epsilon() * std::abs(dir_inv.z)));
306 :
307 : // Check if ray origin is in cube
308 30003993 : if (contains(a)) return true;
309 :
310 : // Refined Smits' method: Largest distance to 3 visible cube face planes along ray is the candidate entry point
311 : Real tx, ty, tz; // Ray position parameter for intersections with box face planes
312 23582539 : if (dir.x > 0.0) {
313 11038830 : tx = (l_.x - a.x) * dir_inv.x;
314 12543709 : } else if (dir.x < 0.0) {
315 12542177 : tx = (u_.x - a.x) * dir_inv.x;
316 : } else { // dir.x == 0
317 1532 : tx = std::numeric_limits<Real>::lowest();
318 : }
319 23582539 : if (dir.y > 0.0) {
320 11826131 : ty = (l_.y - a.y) * dir_inv.y;
321 11756408 : } else if (dir.y < 0.0) {
322 11756067 : ty = (u_.y - a.y) * dir_inv.y;
323 : } else { // dir.y == 0
324 341 : ty = std::numeric_limits<Real>::lowest();
325 : }
326 23582539 : if (dir.z > 0.0) {
327 18448594 : tz = (l_.z - a.z) * dir_inv.z;
328 5133945 : } else if (dir.z < 0.0) {
329 5133945 : tz = (u_.z - a.z) * dir_inv.z;
330 : } else { // dir.z == 0
331 0 : tz = std::numeric_limits<Real>::lowest();
332 : }
333 :
334 23582539 : Real const tmax(ObjexxFCL::max(tx, ty, tz));
335 23582539 : if (tmax >= 0.0) { // Entry point is on ray: See if it is within cube extent
336 13730413 : if (tx == tmax) {
337 5650489 : Real const y(a.y + (tmax * dir.y));
338 5650489 : if ((y < l_.y) || (y > u_.y)) return false;
339 2199248 : Real const z(a.z + (tmax * dir.z));
340 2199248 : if ((z < l_.z) || (z > u_.z)) return false;
341 8079924 : } else if (ty == tmax) {
342 4809117 : Real const x(a.x + (tmax * dir.x));
343 4809117 : if ((x < l_.x) || (x > u_.x)) return false;
344 1643242 : Real const z(a.z + (tmax * dir.z));
345 1643242 : if ((z < l_.z) || (z > u_.z)) return false;
346 : } else { // tz == tmax
347 3270807 : Real const x(a.x + (tmax * dir.x));
348 3270807 : if ((x < l_.x) || (x > u_.x)) return false;
349 1197572 : Real const y(a.y + (tmax * dir.y));
350 1197572 : if ((y < l_.y) || (y > u_.y)) return false;
351 : }
352 2663352 : return true;
353 : } else { // Entry point is on backwards projection of ray
354 9852126 : return false;
355 : }
356 : }
357 :
358 : // Ray Intersects Cube?
359 : bool rayIntersectsCube(Vertex const &a, Vertex const &dir) const
360 : {
361 : return rayIntersectsCube(a, dir, safe_inverse(dir)); // Inefficient if called in loop with same dir
362 : }
363 :
364 : // Line Intersects Cube?
365 : bool lineIntersectsCube(Vertex const &a, Vertex const &dir, Vertex const &dir_inv) const
366 : {
367 : // Note: dir_inv coordinates corresponding to a zero dir coordinate are not used and can be set to zero
368 :
369 : assert(std::abs(dir.mag_squared() - 1.0) < 4 * std::numeric_limits<Real>::epsilon()); // Check unit vector
370 : assert((dir.x == 0.0) || (std::abs(dir_inv.x - (1.0 / dir.x)) < 2 * std::numeric_limits<Real>::epsilon() * std::abs(dir_inv.x)));
371 : assert((dir.y == 0.0) || (std::abs(dir_inv.y - (1.0 / dir.y)) < 2 * std::numeric_limits<Real>::epsilon() * std::abs(dir_inv.y)));
372 : assert((dir.z == 0.0) || (std::abs(dir_inv.z - (1.0 / dir.z)) < 2 * std::numeric_limits<Real>::epsilon() * std::abs(dir_inv.z)));
373 :
374 : // Smits' method: Largest distance to 3 visible cube faces along ray is the candidate entry point
375 : Real txl, txu, tyl, tyu, tzl, tzu; // Ray position parameter for intersections with box face planes
376 : if (dir.x > 0.0) {
377 : txl = (l_.x - a.x) * dir_inv.x;
378 : txu = (u_.x - a.x) * dir_inv.x;
379 : } else if (dir.x < 0.0) {
380 : txl = (u_.x - a.x) * dir_inv.x;
381 : txu = (l_.x - a.x) * dir_inv.x;
382 : } else { // dir.x == 0
383 : if ((l_.x - a.x <= 0.0) && (u_.x - a.x >= 0.0)) {
384 : txl = std::numeric_limits<Real>::lowest();
385 : txu = std::numeric_limits<Real>::max();
386 : } else {
387 : return false;
388 : }
389 : }
390 : if (dir.y > 0.0) {
391 : tyl = (l_.y - a.y) * dir_inv.y;
392 : tyu = (u_.y - a.y) * dir_inv.y;
393 : } else if (dir.y < 0.0) {
394 : tyl = (u_.y - a.y) * dir_inv.y;
395 : tyu = (l_.y - a.y) * dir_inv.y;
396 : } else { // dir.y == 0
397 : if ((l_.y - a.y <= 0.0) && (u_.y - a.y >= 0.0)) {
398 : tyl = std::numeric_limits<Real>::lowest();
399 : tyu = std::numeric_limits<Real>::max();
400 : } else {
401 : return false;
402 : }
403 : }
404 : if ((txl > tyu) || (tyl > txu)) return false;
405 : if (dir.z > 0.0) {
406 : tzl = (l_.z - a.z) * dir_inv.z;
407 : tzu = (u_.z - a.z) * dir_inv.z;
408 : } else if (dir.z < 0.0) {
409 : tzl = (u_.z - a.z) * dir_inv.z;
410 : tzu = (l_.z - a.z) * dir_inv.z;
411 : } else { // dir.z == 0
412 : if ((l_.z - a.z <= 0.0) && (u_.z - a.z >= 0.0)) {
413 : tzl = std::numeric_limits<Real>::lowest();
414 : tzu = std::numeric_limits<Real>::max();
415 : } else {
416 : return false;
417 : }
418 : }
419 : return ((txl <= tzu) && (tzl <= txu) && (tyl <= tzu) && (tzl <= tyu));
420 : }
421 :
422 : // Line Intersects Cube?
423 : bool lineIntersectsCube(Vertex const &a, Vertex const &dir) const
424 : {
425 : return lineIntersectsCube(a, dir, safe_inverse(dir)); // Inefficient if called in loop with same dir
426 : }
427 :
428 : // Surfaces Outer Cube Initilization
429 : void init(EPVector<Surface> &surfaces);
430 :
431 : // Surfaces that Line Segment Intersects Cube's Enclosing Sphere
432 : void surfacesSegmentIntersectsSphere(Vertex const &a, Vertex const &b, Surfaces &surfaces) const
433 : {
434 : if (segmentIntersectsSphere(a, b)) {
435 : if (!surfaces_.empty()) surfaces.insert(surfaces.end(), surfaces_.begin(), surfaces_.end()); // Add this cube's surfaces
436 : for (std::uint8_t i = 0; i < n_; ++i) { // Recurse
437 : cubes_[i]->surfacesSegmentIntersectsSphere(a, b, surfaces);
438 : }
439 : }
440 : }
441 :
442 : // Surfaces that Ray Intersects Cube's Enclosing Sphere
443 : void surfacesRayIntersectsSphere(Vertex const &a, Vertex const &dir, Surfaces &surfaces) const
444 : {
445 : if (rayIntersectsSphere(a, dir)) {
446 : if (!surfaces_.empty()) surfaces.insert(surfaces.end(), surfaces_.begin(), surfaces_.end()); // Add this cube's surfaces
447 : for (std::uint8_t i = 0; i < n_; ++i) { // Recurse
448 : cubes_[i]->surfacesRayIntersectsSphere(a, dir, surfaces);
449 : }
450 : }
451 : }
452 :
453 : // Surfaces that Line Intersects Cube's Enclosing Sphere
454 : void surfacesLineIntersectsSphere(Vertex const &a, Vertex const &dir, Surfaces &surfaces) const
455 : {
456 : if (lineIntersectsSphere(a, dir)) {
457 : if (!surfaces_.empty()) surfaces.insert(surfaces.end(), surfaces_.begin(), surfaces_.end()); // Add this cube's surfaces
458 : for (std::uint8_t i = 0; i < n_; ++i) { // Recurse
459 : cubes_[i]->surfacesLineIntersectsSphere(a, dir, surfaces);
460 : }
461 : }
462 : }
463 :
464 : // Surfaces that Line Segment Intersects Cube
465 : void surfacesSegmentIntersectsCube(Vertex const &a, Vertex const &b, Surfaces &surfaces) const
466 : {
467 : if (segmentIntersectsCube(a, b)) {
468 : if (!surfaces_.empty()) surfaces.insert(surfaces.end(), surfaces_.begin(), surfaces_.end()); // Add this cube's surfaces
469 : for (std::uint8_t i = 0; i < n_; ++i) { // Recurse
470 : cubes_[i]->surfacesSegmentIntersectsCube(a, b, surfaces);
471 : }
472 : }
473 : }
474 :
475 : // Surfaces that Ray Intersects Cube
476 : void surfacesRayIntersectsCube(Vertex const &a, Vertex const &dir, Vertex const &dir_inv, Surfaces &surfaces) const
477 : {
478 : if (rayIntersectsCube(a, dir, dir_inv)) {
479 : if (!surfaces_.empty()) surfaces.insert(surfaces.end(), surfaces_.begin(), surfaces_.end()); // Add this cube's surfaces
480 : for (std::uint8_t i = 0; i < n_; ++i) { // Recurse
481 : cubes_[i]->surfacesRayIntersectsCube(a, dir, dir_inv, surfaces);
482 : }
483 : }
484 : }
485 :
486 : // Surfaces that Ray Intersects Cube
487 : void surfacesRayIntersectsCube(Vertex const &a, Vertex const &dir, Surfaces &surfaces) const
488 : {
489 : surfacesRayIntersectsCube(a, dir, safe_inverse(dir), surfaces); // Inefficient if called in loop with same dir
490 : }
491 :
492 : // Surfaces that Line Intersects Cube
493 : void surfacesLineIntersectsCube(Vertex const &a, Vertex const &dir, Vertex const &dir_inv, Surfaces &surfaces) const
494 : {
495 : if (lineIntersectsCube(a, dir, dir_inv)) {
496 : if (!surfaces_.empty()) surfaces.insert(surfaces.end(), surfaces_.begin(), surfaces_.end()); // Add this cube's surfaces
497 : for (std::uint8_t i = 0; i < n_; ++i) { // Recurse
498 : cubes_[i]->surfacesLineIntersectsCube(a, dir, dir_inv, surfaces);
499 : }
500 : }
501 : }
502 :
503 : // Surfaces that Line Intersects Cube
504 : void surfacesLineIntersectsCube(Vertex const &a, Vertex const &dir, Surfaces &surfaces) const
505 : {
506 : surfacesLineIntersectsCube(a, dir, safe_inverse(dir), surfaces); // Inefficient if called in loop with same dir
507 : }
508 :
509 : // Seek a Surface in Cube that Line Segment Intersects and Satisfies Predicate
510 5853602 : template <typename Predicate> bool hasSurfaceSegmentIntersectsCube(Vertex const &a, Vertex const &b, Predicate const &predicate) const
511 : {
512 5853602 : if (segmentIntersectsCube(a, b)) {
513 78482239 : for (auto const *surface_p : surfaces_) { // Try this cube's surfaces
514 76735020 : if (predicate(*surface_p)) return true;
515 : }
516 7090403 : for (std::uint8_t i = 0; i < n_; ++i) { // Recurse
517 5343184 : if (cubes_[i]->hasSurfaceSegmentIntersectsCube(a, b, predicate)) return true;
518 : }
519 : }
520 5853602 : return false;
521 : }
522 :
523 : // Seek a Surface in Cube that Ray Intersects and Satisfies Predicate
524 : template <typename Predicate>
525 0 : bool hasSurfaceRayIntersectsCube(Vertex const &a, Vertex const &dir, Vertex const &dir_inv, Predicate const &predicate) const
526 : {
527 0 : if (rayIntersectsCube(a, dir, dir_inv)) {
528 0 : for (auto const *surface_p : surfaces_) { // Try this cube's surfaces
529 0 : if (predicate(*surface_p)) return true;
530 : }
531 0 : for (std::uint8_t i = 0; i < n_; ++i) { // Recurse
532 0 : if (cubes_[i]->hasSurfaceRayIntersectsCube(a, dir, dir_inv, predicate)) return true;
533 : }
534 : }
535 0 : return false;
536 : }
537 :
538 : // Seek a Surface in Cube that Ray Intersects and Satisfies Predicate
539 : template <typename Predicate> bool hasSurfaceRayIntersectsCube(Vertex const &a, Vertex const &dir, Predicate const &predicate) const
540 : {
541 : return hasSurfaceRayIntersectsCube(a, dir, safe_inverse(dir), predicate); // Inefficient if called in loop with same dir
542 : }
543 :
544 : // Process Surfaces in Cube that Ray Intersects with Function
545 : template <typename Function>
546 0 : void processSurfaceRayIntersectsCube(Vertex const &a, Vertex const &dir, Vertex const &dir_inv, Function const &function) const
547 : {
548 0 : if (rayIntersectsCube(a, dir, dir_inv)) {
549 0 : for (auto const *surface_p : surfaces_) { // Process this cube's surfaces
550 0 : function(*surface_p);
551 : }
552 0 : for (std::uint8_t i = 0; i < n_; ++i) { // Recurse
553 0 : cubes_[i]->processSurfaceRayIntersectsCube(a, dir, dir_inv, function);
554 : }
555 : }
556 0 : }
557 :
558 : // Process Surfaces in Cube that Ray Intersects with Function
559 : template <typename Function> void processSurfaceRayIntersectsCube(Vertex const &a, Vertex const &dir, Function const &function) const
560 : {
561 : processSurfaceRayIntersectsCube(a, dir, safe_inverse(dir), function); // Inefficient if called in loop with same dir
562 : }
563 :
564 : // Process Surfaces in Cube that Ray Intersects Stopping if Predicate Satisfied
565 : template <typename Predicate>
566 30003993 : bool processSomeSurfaceRayIntersectsCube(
567 : EnergyPlusData &state, Vertex const &a, Vertex const &dir, Vertex const &dir_inv, Predicate const &predicate) const
568 : {
569 30003993 : if (rayIntersectsCube(a, dir, dir_inv)) {
570 308933240 : for (auto const *surface_p : surfaces_) { // Process this cube's surfaces
571 299776131 : if (predicate(*surface_p)) return true; // Don't need to process more surfaces
572 : }
573 36854542 : for (std::uint8_t i = 0; i < n_; ++i) { // Recurse
574 27884295 : if (cubes_[i]->processSomeSurfaceRayIntersectsCube(state, a, dir, dir_inv, predicate))
575 42256 : return true; // Don't need to process more surfaces
576 : }
577 : }
578 29889434 : return false;
579 : }
580 :
581 : // Process Surfaces in Cube that Ray Intersects Stopping if Predicate Satisfied
582 : template <typename Predicate>
583 : bool processSomeSurfaceRayIntersectsCube(EnergyPlusData &state, Vertex const &a, Vertex const &dir, Predicate const &predicate) const
584 : {
585 : return processSomeSurfaceRayIntersectsCube(state, a, dir, safe_inverse(dir), predicate); // Inefficient if called in loop with same dir
586 : }
587 :
588 : public: // Static Methods
589 : // Octree-Safe Vector Inverse
590 2119698 : static Vertex safe_inverse(Vertex const &v)
591 : {
592 2119698 : return Vertex((v.x != 0.0 ? 1.0 / v.x : 0.0), (v.y != 0.0 ? 1.0 / v.y : 0.0), (v.z != 0.0 ? 1.0 / v.z : 0.0));
593 : }
594 :
595 : private: // Methods
596 : // Valid?
597 : bool valid() const;
598 :
599 : // Add a Surface
600 6197 : void add(Surface &surface)
601 : {
602 6197 : surfaces_.push_back(&surface);
603 6197 : }
604 :
605 : // Branch to Sub-Tree
606 : void branch();
607 :
608 : // Surface Branch Processing
609 : void surfaceBranch(Surface &surface);
610 :
611 : private: // Static Methods
612 : // Vertex in Cube Defined by Lower and Upper Corners?
613 0 : static bool contains(Vertex const &l, Vertex const &u, Vertex const &v)
614 : {
615 0 : return (l.x <= v.x) && (v.x <= u.x) && (l.y <= v.y) && (v.y <= u.y) && (l.z <= v.z) && (v.z <= u.z);
616 : }
617 :
618 : // Surface in Cube Defined by Lower and Upper Corners?
619 : static bool contains(Vertex const &l, Vertex const &u, Surface const &surface);
620 :
621 : private: // Static Data
622 : static std::uint8_t const maxDepth_; // Max tree depth
623 : static size_type const maxSurfaces_; // Max surfaces in a cube before subdividing
624 :
625 : private: // Data
626 : std::uint8_t d_; // Depth in tree
627 : std::uint8_t n_; // Number of active sub-cubes ([0,8])
628 : Vertex l_; // Lower corner
629 : Vertex u_; // Upper corner
630 : Vertex c_; // Center point
631 : Real w_; // Side width
632 : Real r_; // Enclosing sphere radius
633 : SurfaceOctreeCube *cubes_[8]; // Sub-cubes
634 : Surfaces surfaces_; // Surfaces in this cube
635 :
636 : }; // SurfaceOctreeCube
637 :
638 : } // namespace EnergyPlus
639 :
640 : #endif
|
