Line data Source code
1 : // EnergyPlus, Copyright (c) 1996-2025, 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 801 : 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 7209 : for (std::uint8_t i = 0; i < 8; ++i) {
101 6408 : cubes_[i] = nullptr; // VC++ 2013 compatible initialization
102 : }
103 801 : }
104 :
105 : // Surfaces Outer Cube Constructor
106 0 : 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)
107 : {
108 0 : for (std::uint8_t i = 0; i < 8; ++i) {
109 0 : cubes_[i] = nullptr; // VC++ 2013 compatible initialization
110 : }
111 0 : init(surfaces);
112 0 : }
113 :
114 : // Box Constructor
115 557 : SurfaceOctreeCube(std::uint8_t const d, Vertex const &l, Vertex const &u, Real const w)
116 557 : : d_(d), n_(0u), l_(l), u_(u), c_(cen(l, u)), w_(w), r_(0.75 * (w * w))
117 : {
118 5013 : for (std::uint8_t i = 0; i < 8; ++i) {
119 4456 : cubes_[i] = nullptr; // VC++ 2013 compatible initialization
120 : }
121 557 : assert(valid());
122 557 : }
123 :
124 : // Destructor
125 1358 : ~SurfaceOctreeCube()
126 : {
127 1915 : for (std::uint8_t i = 0; i < n_; ++i) {
128 557 : delete cubes_[i];
129 : }
130 1358 : }
131 :
132 : public: // Properties
133 : // Depth
134 : std::uint8_t d() const
135 : {
136 : return d_;
137 : }
138 :
139 : // Depth
140 : std::uint8_t depth() const
141 : {
142 : return d_;
143 : }
144 :
145 : // Number of Children
146 : std::uint8_t nChildren() const
147 : {
148 : return n_;
149 : }
150 :
151 : // Number of Sub-Cubes
152 : std::uint8_t nSubCube() const
153 : {
154 : return n_;
155 : }
156 :
157 : // Lower Corner
158 0 : Vertex const &l() const
159 : {
160 0 : return l_;
161 : }
162 :
163 : // Upper Corner
164 0 : Vertex const &u() const
165 : {
166 0 : return u_;
167 : }
168 :
169 : // Center Point
170 0 : Vertex const &c() const
171 : {
172 0 : return c_;
173 : }
174 :
175 : // Center Point
176 : Vertex const ¢er() const
177 : {
178 : return c_;
179 : }
180 :
181 : // Width
182 0 : Real w() const
183 : {
184 0 : return w_;
185 : }
186 :
187 : // Width
188 : Real width() const
189 : {
190 : return w_;
191 : }
192 :
193 : // Surfaces
194 : Surfaces const &surfaces() const
195 : {
196 : return surfaces_;
197 : }
198 :
199 : // Surfaces
200 : Surfaces::size_type surfaces_size() const
201 : {
202 : return surfaces_.size();
203 : }
204 :
205 : // Surfaces Begin Iterator
206 : Surfaces::const_iterator surfaces_begin() const
207 : {
208 : return surfaces_.begin();
209 : }
210 :
211 : // Surfaces Begin Iterator
212 : Surfaces::iterator surfaces_begin()
213 : {
214 : return surfaces_.begin();
215 : }
216 :
217 : // Surfaces End Iterator
218 : Surfaces::const_iterator surfaces_end() const
219 : {
220 : return surfaces_.end();
221 : }
222 :
223 : // Surfaces End Iterator
224 : Surfaces::iterator surfaces_end()
225 : {
226 : return surfaces_.end();
227 : }
228 :
229 : public: // Methods
230 : // Vertex in Cube?
231 42856471 : bool contains(Vertex const &v) const
232 : {
233 42856471 : 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);
234 : }
235 :
236 : // Surface in Cube?
237 : bool contains(Surface const &surface) const;
238 :
239 : // Line Segment Intersects Enclosing Sphere?
240 0 : bool segmentIntersectsSphere(Vertex const &a, Vertex const &b) const
241 : {
242 0 : Vertex const ab(b - a);
243 0 : Real const ab_mag_squared(ab.mag_squared());
244 0 : if (ab_mag_squared == 0.0) { // Segment is a point
245 0 : return ObjexxFCL::distance_squared(a, c_) <= r_;
246 : } else { // Might pay to check if a or b in sphere first in some applications
247 0 : Vertex const ac(c_ - a);
248 0 : Real const projection_fac(((ac.x * ab.x) + (ac.y * ab.y) + (ac.z * ab.z)) / ab_mag_squared);
249 0 : if ((0.0 <= projection_fac) && (projection_fac <= 1.0)) { // Projected (closest) point is on ab segment
250 0 : return ObjexxFCL::distance_squared(ac, projection_fac * ab) <= r_;
251 : } else { // Projection (closest) point is outside of ab segment: Intersects iff a or b are in sphere
252 0 : return (ObjexxFCL::distance_squared(a, c_) <= r_) || (ObjexxFCL::distance_squared(b, c_) <= r_);
253 : }
254 0 : }
255 0 : }
256 :
257 : // Ray Intersects Enclosing Sphere?
258 0 : bool rayIntersectsSphere(Vertex const &a, Vertex const &dir) const
259 : {
260 0 : assert(std::abs(dir.mag_squared() - 1.0) < 4 * std::numeric_limits<Real>::epsilon()); // Check unit vector
261 : // Might pay to check if a in sphere first in some applications
262 0 : Vertex const ac(c_ - a);
263 0 : Real const projection_fac((ac.x * dir.x) + (ac.y * dir.y) + (ac.z * dir.z));
264 0 : if (0.0 <= projection_fac) { // Projected (closest) point is on ray
265 0 : return ObjexxFCL::distance_squared(ac, projection_fac * dir) <= r_;
266 : } else { // Projection (closest) point is outside of ray: Intersects iff a is in sphere
267 0 : return ObjexxFCL::distance_squared(a, c_) <= r_;
268 : }
269 0 : }
270 :
271 : // Line Intersects Enclosing Sphere?
272 0 : bool lineIntersectsSphere(Vertex const &a, Vertex const &dir) const
273 : {
274 0 : assert(std::abs(dir.mag_squared() - 1.0) < 4 * std::numeric_limits<Real>::epsilon()); // Check unit vector
275 0 : Vertex const ac(c_ - a);
276 0 : return ac.mag_squared() - ObjexxFCL::square((ac.x * dir.x) + (ac.y * dir.y) + (ac.z * dir.z)) <= r_;
277 0 : }
278 :
279 : // Line Segment Intersects Cube?
280 6281823 : bool segmentIntersectsCube(Vertex const &a, Vertex const &b) const
281 : {
282 : // Check if either segment endpoint is in cube
283 6281823 : if (contains(a) || contains(b)) {
284 1797866 : return true;
285 : }
286 :
287 : // Use separating axis theorem (faster variants exist)
288 4483957 : Vertex const m(mid(a, b) - c_); // Mid-point relative to cube center
289 4483957 : Vertex const mb(b - c_ - m); // ab mid-point to b half segment vector
290 4483957 : Vertex const e(std::abs(mb.x), std::abs(mb.y), std::abs(mb.z)); // Extent of half ab segment
291 4483957 : Real const h(0.5 * w_); // Half-width
292 : // Check if x,y,z axes are separating
293 4483957 : if (std::abs(m.x) > h + e.x) {
294 2227916 : return false;
295 : }
296 2256041 : if (std::abs(m.y) > h + e.y) {
297 1517279 : return false;
298 : }
299 738762 : if (std::abs(m.z) > h + e.z) {
300 601122 : return false;
301 : }
302 : // Check if cross product axes are separating
303 137640 : if (std::abs(m.y * mb.z - m.z * mb.y) > h * (e.z + e.y)) {
304 20165 : return false;
305 : }
306 117475 : if (std::abs(m.x * mb.z - m.z * mb.x) > h * (e.z + e.x)) {
307 27050 : return false;
308 : }
309 90425 : if (std::abs(m.x * mb.y - m.y * mb.x) > h * (e.y + e.x)) {
310 1790 : return false;
311 : }
312 88635 : return true;
313 4483957 : }
314 :
315 : // Ray Intersects Cube?
316 31691428 : bool rayIntersectsCube(Vertex const &a, Vertex const &dir, Vertex const &dir_inv) const
317 : {
318 : // Note: dir_inv coordinates corresponding to a zero dir coordinate are not used and can be set to zero
319 :
320 31691428 : assert(std::abs(dir.mag_squared() - 1.0) < 4 * std::numeric_limits<Real>::epsilon()); // Check unit vector
321 31691428 : 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)));
322 31691428 : 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)));
323 31691428 : 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)));
324 :
325 : // Check if ray origin is in cube
326 31691428 : if (contains(a)) {
327 6910290 : return true;
328 : }
329 :
330 : // Refined Smits' method: Largest distance to 3 visible cube face planes along ray is the candidate entry point
331 : Real tx, ty, tz; // Ray position parameter for intersections with box face planes
332 24781138 : if (dir.x > 0.0) {
333 11609518 : tx = (l_.x - a.x) * dir_inv.x;
334 13171620 : } else if (dir.x < 0.0) {
335 13169919 : tx = (u_.x - a.x) * dir_inv.x;
336 : } else { // dir.x == 0
337 1701 : tx = std::numeric_limits<Real>::lowest();
338 : }
339 24781138 : if (dir.y > 0.0) {
340 12406152 : ty = (l_.y - a.y) * dir_inv.y;
341 12374986 : } else if (dir.y < 0.0) {
342 12373774 : ty = (u_.y - a.y) * dir_inv.y;
343 : } else { // dir.y == 0
344 1212 : ty = std::numeric_limits<Real>::lowest();
345 : }
346 24781138 : if (dir.z > 0.0) {
347 19203853 : tz = (l_.z - a.z) * dir_inv.z;
348 5577285 : } else if (dir.z < 0.0) {
349 5577285 : tz = (u_.z - a.z) * dir_inv.z;
350 : } else { // dir.z == 0
351 0 : tz = std::numeric_limits<Real>::lowest();
352 : }
353 :
354 24781138 : Real const tmax(ObjexxFCL::max(tx, ty, tz));
355 24781138 : if (tmax >= 0.0) { // Entry point is on ray: See if it is within cube extent
356 14439455 : if (tx == tmax) {
357 5942153 : Real const y(a.y + (tmax * dir.y));
358 5942153 : if ((y < l_.y) || (y > u_.y)) {
359 3637679 : return false;
360 : }
361 2304474 : Real const z(a.z + (tmax * dir.z));
362 2304474 : if ((z < l_.z) || (z > u_.z)) {
363 1094616 : return false;
364 : }
365 8497302 : } else if (ty == tmax) {
366 5050752 : Real const x(a.x + (tmax * dir.x));
367 5050752 : if ((x < l_.x) || (x > u_.x)) {
368 3335165 : return false;
369 : }
370 1715587 : Real const z(a.z + (tmax * dir.z));
371 1715587 : if ((z < l_.z) || (z > u_.z)) {
372 823436 : return false;
373 : }
374 : } else { // tz == tmax
375 3446550 : Real const x(a.x + (tmax * dir.x));
376 3446550 : if ((x < l_.x) || (x > u_.x)) {
377 2166424 : return false;
378 : }
379 1280126 : Real const y(a.y + (tmax * dir.y));
380 1280126 : if ((y < l_.y) || (y > u_.y)) {
381 574442 : return false;
382 : }
383 : }
384 2807693 : return true;
385 : } else { // Entry point is on backwards projection of ray
386 10341683 : return false;
387 : }
388 : }
389 :
390 : // Ray Intersects Cube?
391 0 : bool rayIntersectsCube(Vertex const &a, Vertex const &dir) const
392 : {
393 0 : return rayIntersectsCube(a, dir, safe_inverse(dir)); // Inefficient if called in loop with same dir
394 : }
395 :
396 : // Line Intersects Cube?
397 0 : bool lineIntersectsCube(Vertex const &a, Vertex const &dir, Vertex const &dir_inv) const
398 : {
399 : // Note: dir_inv coordinates corresponding to a zero dir coordinate are not used and can be set to zero
400 :
401 0 : assert(std::abs(dir.mag_squared() - 1.0) < 4 * std::numeric_limits<Real>::epsilon()); // Check unit vector
402 0 : 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)));
403 0 : 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)));
404 0 : 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)));
405 :
406 : // Smits' method: Largest distance to 3 visible cube faces along ray is the candidate entry point
407 : Real txl, txu, tyl, tyu, tzl, tzu; // Ray position parameter for intersections with box face planes
408 0 : if (dir.x > 0.0) {
409 0 : txl = (l_.x - a.x) * dir_inv.x;
410 0 : txu = (u_.x - a.x) * dir_inv.x;
411 0 : } else if (dir.x < 0.0) {
412 0 : txl = (u_.x - a.x) * dir_inv.x;
413 0 : txu = (l_.x - a.x) * dir_inv.x;
414 : } else { // dir.x == 0
415 0 : if ((l_.x - a.x <= 0.0) && (u_.x - a.x >= 0.0)) {
416 0 : txl = std::numeric_limits<Real>::lowest();
417 0 : txu = std::numeric_limits<Real>::max();
418 : } else {
419 0 : return false;
420 : }
421 : }
422 0 : if (dir.y > 0.0) {
423 0 : tyl = (l_.y - a.y) * dir_inv.y;
424 0 : tyu = (u_.y - a.y) * dir_inv.y;
425 0 : } else if (dir.y < 0.0) {
426 0 : tyl = (u_.y - a.y) * dir_inv.y;
427 0 : tyu = (l_.y - a.y) * dir_inv.y;
428 : } else { // dir.y == 0
429 0 : if ((l_.y - a.y <= 0.0) && (u_.y - a.y >= 0.0)) {
430 0 : tyl = std::numeric_limits<Real>::lowest();
431 0 : tyu = std::numeric_limits<Real>::max();
432 : } else {
433 0 : return false;
434 : }
435 : }
436 0 : if ((txl > tyu) || (tyl > txu)) {
437 0 : return false;
438 : }
439 0 : if (dir.z > 0.0) {
440 0 : tzl = (l_.z - a.z) * dir_inv.z;
441 0 : tzu = (u_.z - a.z) * dir_inv.z;
442 0 : } else if (dir.z < 0.0) {
443 0 : tzl = (u_.z - a.z) * dir_inv.z;
444 0 : tzu = (l_.z - a.z) * dir_inv.z;
445 : } else { // dir.z == 0
446 0 : if ((l_.z - a.z <= 0.0) && (u_.z - a.z >= 0.0)) {
447 0 : tzl = std::numeric_limits<Real>::lowest();
448 0 : tzu = std::numeric_limits<Real>::max();
449 : } else {
450 0 : return false;
451 : }
452 : }
453 0 : return ((txl <= tzu) && (tzl <= txu) && (tyl <= tzu) && (tzl <= tyu));
454 : }
455 :
456 : // Line Intersects Cube?
457 0 : bool lineIntersectsCube(Vertex const &a, Vertex const &dir) const
458 : {
459 0 : return lineIntersectsCube(a, dir, safe_inverse(dir)); // Inefficient if called in loop with same dir
460 : }
461 :
462 : // Surfaces Outer Cube Initilization
463 : void init(EPVector<Surface> &surfaces);
464 :
465 : // Surfaces that Line Segment Intersects Cube's Enclosing Sphere
466 0 : void surfacesSegmentIntersectsSphere(Vertex const &a, Vertex const &b, Surfaces &surfaces) const
467 : {
468 0 : if (segmentIntersectsSphere(a, b)) {
469 0 : if (!surfaces_.empty()) {
470 0 : surfaces.insert(surfaces.end(), surfaces_.begin(), surfaces_.end()); // Add this cube's surfaces
471 : }
472 0 : for (std::uint8_t i = 0; i < n_; ++i) { // Recurse
473 0 : cubes_[i]->surfacesSegmentIntersectsSphere(a, b, surfaces);
474 : }
475 : }
476 0 : }
477 :
478 : // Surfaces that Ray Intersects Cube's Enclosing Sphere
479 0 : void surfacesRayIntersectsSphere(Vertex const &a, Vertex const &dir, Surfaces &surfaces) const
480 : {
481 0 : if (rayIntersectsSphere(a, dir)) {
482 0 : if (!surfaces_.empty()) {
483 0 : surfaces.insert(surfaces.end(), surfaces_.begin(), surfaces_.end()); // Add this cube's surfaces
484 : }
485 0 : for (std::uint8_t i = 0; i < n_; ++i) { // Recurse
486 0 : cubes_[i]->surfacesRayIntersectsSphere(a, dir, surfaces);
487 : }
488 : }
489 0 : }
490 :
491 : // Surfaces that Line Intersects Cube's Enclosing Sphere
492 0 : void surfacesLineIntersectsSphere(Vertex const &a, Vertex const &dir, Surfaces &surfaces) const
493 : {
494 0 : if (lineIntersectsSphere(a, dir)) {
495 0 : if (!surfaces_.empty()) {
496 0 : surfaces.insert(surfaces.end(), surfaces_.begin(), surfaces_.end()); // Add this cube's surfaces
497 : }
498 0 : for (std::uint8_t i = 0; i < n_; ++i) { // Recurse
499 0 : cubes_[i]->surfacesLineIntersectsSphere(a, dir, surfaces);
500 : }
501 : }
502 0 : }
503 :
504 : // Surfaces that Line Segment Intersects Cube
505 0 : void surfacesSegmentIntersectsCube(Vertex const &a, Vertex const &b, Surfaces &surfaces) const
506 : {
507 0 : if (segmentIntersectsCube(a, b)) {
508 0 : if (!surfaces_.empty()) {
509 0 : surfaces.insert(surfaces.end(), surfaces_.begin(), surfaces_.end()); // Add this cube's surfaces
510 : }
511 0 : for (std::uint8_t i = 0; i < n_; ++i) { // Recurse
512 0 : cubes_[i]->surfacesSegmentIntersectsCube(a, b, surfaces);
513 : }
514 : }
515 0 : }
516 :
517 : // Surfaces that Ray Intersects Cube
518 0 : void surfacesRayIntersectsCube(Vertex const &a, Vertex const &dir, Vertex const &dir_inv, Surfaces &surfaces) const
519 : {
520 0 : if (rayIntersectsCube(a, dir, dir_inv)) {
521 0 : if (!surfaces_.empty()) {
522 0 : surfaces.insert(surfaces.end(), surfaces_.begin(), surfaces_.end()); // Add this cube's surfaces
523 : }
524 0 : for (std::uint8_t i = 0; i < n_; ++i) { // Recurse
525 0 : cubes_[i]->surfacesRayIntersectsCube(a, dir, dir_inv, surfaces);
526 : }
527 : }
528 0 : }
529 :
530 : // Surfaces that Ray Intersects Cube
531 0 : void surfacesRayIntersectsCube(Vertex const &a, Vertex const &dir, Surfaces &surfaces) const
532 : {
533 0 : surfacesRayIntersectsCube(a, dir, safe_inverse(dir), surfaces); // Inefficient if called in loop with same dir
534 0 : }
535 :
536 : // Surfaces that Line Intersects Cube
537 0 : void surfacesLineIntersectsCube(Vertex const &a, Vertex const &dir, Vertex const &dir_inv, Surfaces &surfaces) const
538 : {
539 0 : if (lineIntersectsCube(a, dir, dir_inv)) {
540 0 : if (!surfaces_.empty()) {
541 0 : surfaces.insert(surfaces.end(), surfaces_.begin(), surfaces_.end()); // Add this cube's surfaces
542 : }
543 0 : for (std::uint8_t i = 0; i < n_; ++i) { // Recurse
544 0 : cubes_[i]->surfacesLineIntersectsCube(a, dir, dir_inv, surfaces);
545 : }
546 : }
547 0 : }
548 :
549 : // Surfaces that Line Intersects Cube
550 0 : void surfacesLineIntersectsCube(Vertex const &a, Vertex const &dir, Surfaces &surfaces) const
551 : {
552 0 : surfacesLineIntersectsCube(a, dir, safe_inverse(dir), surfaces); // Inefficient if called in loop with same dir
553 0 : }
554 :
555 : // Seek a Surface in Cube that Line Segment Intersects and Satisfies Predicate
556 6281823 : template <typename Predicate> bool hasSurfaceSegmentIntersectsCube(Vertex const &a, Vertex const &b, Predicate const &predicate) const
557 : {
558 6281823 : if (segmentIntersectsCube(a, b)) {
559 87286908 : for (auto const *surface_p : surfaces_) { // Try this cube's surfaces
560 85400407 : if (predicate(*surface_p)) {
561 0 : return true;
562 : }
563 : }
564 7626399 : for (std::uint8_t i = 0; i < n_; ++i) { // Recurse
565 5739898 : if (cubes_[i]->hasSurfaceSegmentIntersectsCube(a, b, predicate)) {
566 0 : return true;
567 : }
568 : }
569 : }
570 6281823 : return false;
571 : }
572 :
573 : // Seek a Surface in Cube that Ray Intersects and Satisfies Predicate
574 : template <typename Predicate>
575 0 : bool hasSurfaceRayIntersectsCube(Vertex const &a, Vertex const &dir, Vertex const &dir_inv, Predicate const &predicate) const
576 : {
577 0 : if (rayIntersectsCube(a, dir, dir_inv)) {
578 0 : for (auto const *surface_p : surfaces_) { // Try this cube's surfaces
579 0 : if (predicate(*surface_p)) {
580 0 : return true;
581 : }
582 : }
583 0 : for (std::uint8_t i = 0; i < n_; ++i) { // Recurse
584 0 : if (cubes_[i]->hasSurfaceRayIntersectsCube(a, dir, dir_inv, predicate)) {
585 0 : return true;
586 : }
587 : }
588 : }
589 0 : return false;
590 : }
591 :
592 : // Seek a Surface in Cube that Ray Intersects and Satisfies Predicate
593 0 : template <typename Predicate> bool hasSurfaceRayIntersectsCube(Vertex const &a, Vertex const &dir, Predicate const &predicate) const
594 : {
595 0 : return hasSurfaceRayIntersectsCube(a, dir, safe_inverse(dir), predicate); // Inefficient if called in loop with same dir
596 : }
597 :
598 : // Process Surfaces in Cube that Ray Intersects with Function
599 : template <typename Function>
600 0 : void processSurfaceRayIntersectsCube(Vertex const &a, Vertex const &dir, Vertex const &dir_inv, Function const &function) const
601 : {
602 0 : if (rayIntersectsCube(a, dir, dir_inv)) {
603 0 : for (auto const *surface_p : surfaces_) { // Process this cube's surfaces
604 0 : function(*surface_p);
605 : }
606 0 : for (std::uint8_t i = 0; i < n_; ++i) { // Recurse
607 0 : cubes_[i]->processSurfaceRayIntersectsCube(a, dir, dir_inv, function);
608 : }
609 : }
610 0 : }
611 :
612 : // Process Surfaces in Cube that Ray Intersects with Function
613 0 : template <typename Function> void processSurfaceRayIntersectsCube(Vertex const &a, Vertex const &dir, Function const &function) const
614 : {
615 0 : processSurfaceRayIntersectsCube(a, dir, safe_inverse(dir), function); // Inefficient if called in loop with same dir
616 0 : }
617 :
618 : // Process Surfaces in Cube that Ray Intersects Stopping if Predicate Satisfied
619 : template <typename Predicate>
620 31691428 : bool processSomeSurfaceRayIntersectsCube(
621 : EnergyPlusData &state, Vertex const &a, Vertex const &dir, Vertex const &dir_inv, Predicate const &predicate) const
622 : {
623 31691428 : if (rayIntersectsCube(a, dir, dir_inv)) {
624 346368486 : for (auto const *surface_p : surfaces_) { // Process this cube's surfaces
625 336544556 : if (predicate(*surface_p)) {
626 105947 : return true; // Don't need to process more surfaces
627 : }
628 : }
629 38964035 : for (std::uint8_t i = 0; i < n_; ++i) { // Recurse
630 29436223 : if (cubes_[i]->processSomeSurfaceRayIntersectsCube(state, a, dir, dir_inv, predicate)) {
631 84224 : return true; // Don't need to process more surfaces
632 : }
633 : }
634 : }
635 31501257 : return false;
636 : }
637 :
638 : // Process Surfaces in Cube that Ray Intersects Stopping if Predicate Satisfied
639 : template <typename Predicate>
640 0 : bool processSomeSurfaceRayIntersectsCube(EnergyPlusData &state, Vertex const &a, Vertex const &dir, Predicate const &predicate) const
641 : {
642 0 : return processSomeSurfaceRayIntersectsCube(state, a, dir, safe_inverse(dir), predicate); // Inefficient if called in loop with same dir
643 : }
644 :
645 : public: // Static Methods
646 : // Octree-Safe Vector Inverse
647 2255205 : static Vertex safe_inverse(Vertex const &v)
648 : {
649 2255205 : 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));
650 : }
651 :
652 : private: // Methods
653 : // Valid?
654 : bool valid() const;
655 :
656 : // Add a Surface
657 7908 : void add(Surface &surface)
658 : {
659 7908 : surfaces_.push_back(&surface);
660 7908 : }
661 :
662 : // Branch to Sub-Tree
663 : void branch();
664 :
665 : // Surface Branch Processing
666 : void surfaceBranch(Surface &surface);
667 :
668 : private: // Static Methods
669 : // Vertex in Cube Defined by Lower and Upper Corners?
670 0 : static bool contains(Vertex const &l, Vertex const &u, Vertex const &v)
671 : {
672 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);
673 : }
674 :
675 : // Surface in Cube Defined by Lower and Upper Corners?
676 : static bool contains(Vertex const &l, Vertex const &u, Surface const &surface);
677 :
678 : private: // Static Data
679 : static std::uint8_t const maxDepth_; // Max tree depth
680 : static size_type const maxSurfaces_; // Max surfaces in a cube before subdividing
681 :
682 : private: // Data
683 : std::uint8_t d_; // Depth in tree
684 : std::uint8_t n_; // Number of active sub-cubes ([0,8])
685 : Vertex l_; // Lower corner
686 : Vertex u_; // Upper corner
687 : Vertex c_; // Center point
688 : Real w_; // Side width
689 : Real r_; // Enclosing sphere radius
690 : SurfaceOctreeCube *cubes_[8]; // Sub-cubes
691 : Surfaces surfaces_; // Surfaces in this cube
692 :
693 : }; // SurfaceOctreeCube
694 :
695 : } // namespace EnergyPlus
696 :
697 : #endif
|
