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_PierceSurface_hh_INCLUDED
49 : #define EnergyPlus_PierceSurface_hh_INCLUDED
50 :
51 : // Purpose: Functions for checking if a ray hits a surface
52 : //
53 : // Author: Stuart Mentzer (Stuart_Mentzer@objexx.com)
54 : //
55 : // History:
56 : // Jun 2015: Last update of legacy version based on DOE-2 DPIERC
57 : // Jan 2016: Initial release
58 : //
59 : // Notes:
60 : // This is filling the role of the former PierceSurface function authored by Fred Winkelmann and based on
61 : // DOE-2.1E subroutine DPIERC and some aspects of this version are analogous
62 : // To match the former behavior rays with origin exactly on the surface are treated as not hitting
63 : // These functions are VERY performance critical for daylighting and solar reflection
64 : // This high-performance implementation was built to complement the octree system for scalability of those systems
65 : // This has been carefully designed for speed but is probably not be optimal yet
66 : // For EnergyPlus most surfaces are rectangular so that is the most important for performance
67 : // Inlining, storing preprocessed values in Surface, 2D projection, & short circuiting are used here for speed
68 : // Agressive inlining options may be needed to get peak performance
69 : // Don't make changes here without validating the performance impact
70 :
71 : // EnergyPlus Headers
72 : #include <EnergyPlus/DataSurfaces.hh>
73 : #include <EnergyPlus/EnergyPlus.hh>
74 : #include <EnergyPlus/Platform.hh>
75 :
76 : // ObjexxFCL Headers
77 : #include <ObjexxFCL/Vector2.hh>
78 : #include <ObjexxFCL/Vector3.hh>
79 : #include <ObjexxFCL/Vector4.hh>
80 :
81 : // C++ Headers
82 : #include <algorithm>
83 : #include <cassert>
84 : #include <limits>
85 :
86 : namespace EnergyPlus {
87 :
88 2 : inline bool PierceSurface_Triangular(DataSurfaces::Surface2D const &s2d, // 2D surface
89 : Vector2<Real64> const &h2d // 2D hit point
90 : )
91 : {
92 : // Purpose: Check if a 2D hit point is in a triangular 2D surface
93 : //
94 : // Author: Stuart Mentzer (Stuart_Mentzer@objexx.com)
95 : //
96 : // History:
97 : // Jan 2016: Initial release
98 : //
99 : // Notes:
100 : // Pulled this case out into separate function to facilitate inlining
101 :
102 : using DataSurfaces::Surface2D;
103 2 : Surface2D::Vertices const &vs(s2d.vertices); // 2D surface vertices
104 2 : Surface2D::Vectors const &es(s2d.edges); // 2D surface edge vectors
105 2 : if (es[0].cross(h2d - vs[0]) < 0.0) {
106 0 : return false;
107 : }
108 2 : if (es[1].cross(h2d - vs[1]) < 0.0) {
109 0 : return false;
110 : }
111 2 : if (es[2].cross(h2d - vs[2]) < 0.0) {
112 0 : return false;
113 : }
114 2 : return true;
115 : } // PierceSurface_Triangular()
116 :
117 5 : inline bool PierceSurface_Convex(DataSurfaces::Surface2D const &s2d, // 2D surface
118 : Vector2<Real64> const &h2d // 2D hit point
119 : )
120 : {
121 : // Purpose: Check if a 2D hit point is in a convex 2D surface
122 : //
123 : // Author: Stuart Mentzer (Stuart_Mentzer@objexx.com)
124 : //
125 : // History:
126 : // Jan 2016: Initial release
127 : //
128 : // Notes:
129 : // Pulled this rare case out into separate function to facilitate inlining
130 : // This is O( n ) complexity so it is isn't used for many-vertex surfaces
131 :
132 : using DataSurfaces::Surface2D;
133 5 : Surface2D::Vertices const &vs(s2d.vertices); // 2D surface vertices
134 5 : Surface2D::Vectors const &es(s2d.edges); // 2D surface edge vectors
135 5 : Surface2D::Vertices::size_type const n(vs.size());
136 5 : assert(n >= 3u);
137 5 : switch (n) {
138 5 : case 8:
139 5 : if (es[7].cross(h2d - vs[7]) < 0.0) {
140 0 : return false;
141 : }
142 : // fallthrough
143 : case 7:
144 5 : if (es[6].cross(h2d - vs[6]) < 0.0) {
145 0 : return false;
146 : }
147 : // fallthrough
148 : case 6:
149 5 : if (es[5].cross(h2d - vs[5]) < 0.0) {
150 0 : return false;
151 : }
152 : // fallthrough
153 : case 5:
154 5 : if (es[4].cross(h2d - vs[4]) < 0.0) {
155 0 : return false;
156 : }
157 : // fallthrough
158 : case 4:
159 5 : if (es[3].cross(h2d - vs[3]) < 0.0) {
160 0 : return false;
161 : }
162 : // fallthrough
163 : case 3:
164 5 : if (es[2].cross(h2d - vs[2]) < 0.0) {
165 0 : return false;
166 : }
167 5 : if (es[1].cross(h2d - vs[1]) < 0.0) {
168 0 : return false;
169 : }
170 5 : if (es[0].cross(h2d - vs[0]) < 0.0) {
171 0 : return false;
172 : }
173 5 : return true;
174 0 : default:
175 0 : for (Surface2D::Vertices::size_type i = 0; i < n; ++i) {
176 0 : if (es[i].cross(h2d - vs[i]) < 0.0) {
177 0 : return false;
178 : }
179 : }
180 0 : return true;
181 : }
182 : } // PierceSurface_Convex()
183 :
184 17 : inline bool PierceSurface_Nonconvex(DataSurfaces::Surface2D const &s2d, // 2D surface
185 : Vector2<Real64> const &h2d // 2D hit point
186 : )
187 : {
188 : // Purpose: Check if a 2D hit point is in a 2D possibly nonconvex surface
189 : //
190 : // Author: Stuart Mentzer (Stuart_Mentzer@objexx.com)
191 : //
192 : // History:
193 : // Jan 2016: Initial release
194 : //
195 : // Notes:
196 : // Pulled this rare case out into separate function to facilitate inlining
197 : // This works for nonconvex "simple" (no edge crossings) polygons
198 : // This is also a fast O( log n ) algorithm for many-vertex convex surfaces
199 :
200 : using DataSurfaces::Surface2D;
201 : using size_type = Surface2D::Vertices::size_type;
202 : using Slab = DataSurfaces::Surface2DSlab;
203 : using Vertex2D = Vector2<Real64>;
204 17 : assert(s2d.vertices.size() >= 3u);
205 17 : Surface2D::Slabs const &slabs(s2d.slabs); // 2D surface y slice slabs
206 17 : Surface2D::SlabYs const &slabYs(s2d.slabYs); // 2D surface slab y coordinates
207 17 : assert(slabYs.size() > 0u);
208 17 : Real64 const yHit(h2d.y); // Hit point y coordinate
209 :
210 : // Find slab with y range containing hit point
211 17 : auto const iHit(std::lower_bound(slabYs.begin(), slabYs.end(), yHit));
212 17 : assert((yHit >= slabYs.front()) && (yHit <= slabYs.back())); // Passed bounding box check so hit point in slabs y range
213 17 : assert(iHit != slabYs.end()); // Hit point can't be above all slabs: passed bounding box check
214 17 : size_type const iSlab(std::min(static_cast<size_type>(iHit - 1 - slabYs.begin()), slabs.size())); // Hit slab index
215 17 : Slab const &slab(slabs[iSlab]);
216 :
217 : // Check hit point within slab bounding box x range
218 17 : Real64 const xHit(h2d.x); // Hit point x coordinate
219 17 : if ((xHit < slab.xl) || (xHit > slab.xu)) {
220 1 : return false; // Hit point outside slab bounding box
221 : }
222 :
223 : // Find edge pair surrounding hit point
224 16 : Slab::Edges const &slabEdges(slab.edges);
225 16 : Slab::EdgesXY const &slabEdgesXY(slab.edgesXY);
226 16 : size_type const nEdges(slabEdges.size());
227 16 : assert(nEdges >= 2u);
228 16 : if (nEdges == 2) { // 2 edges
229 7 : Slab::Edge const se0(slabEdges[0]);
230 7 : Slab::EdgeXY const eXY0(slabEdgesXY[0]);
231 7 : Vertex2D v0(s2d.vertices[se0]);
232 7 : Surface2D::Edge e0(s2d.edges[se0]);
233 7 : Real64 const x0(v0.x + (yHit - v0.y) * eXY0);
234 7 : if (xHit < x0) {
235 1 : return false; // Hit point x is left of left edge
236 : }
237 6 : Slab::Edge const se1(slabEdges[1]);
238 6 : Slab::EdgeXY const eXY1(slabEdgesXY[1]);
239 6 : Vertex2D v1(s2d.vertices[se1]);
240 6 : Surface2D::Edge e1(s2d.edges[se1]);
241 6 : Real64 const x1(v1.x + (yHit - v1.y) * eXY1);
242 6 : if (x1 < xHit) {
243 1 : return false; // Hit point is right of right edge
244 : }
245 11 : } else { // 4+ edges: Binary search for edges surrounding hit point
246 9 : assert(nEdges >= 4u);
247 9 : assert(nEdges % 2 == 0u);
248 9 : size_type l(0u), u(nEdges - 1);
249 9 : Slab::Edge const il(slabEdges[l]);
250 9 : Slab::EdgeXY const eXYl(slabEdgesXY[l]);
251 9 : Vertex2D const &vl(s2d.vertices[il]);
252 9 : Surface2D::Edge const el(s2d.edges[il]);
253 9 : Real64 const xl(vl.x + (yHit - vl.y) * eXYl);
254 9 : if (xHit < xl) {
255 0 : return false; // Hit point x is left of leftmost edge
256 : }
257 9 : Slab::Edge const iu(slabEdges[u]);
258 9 : Slab::EdgeXY const eXYu(slabEdgesXY[u]);
259 9 : Vertex2D const &vu(s2d.vertices[iu]);
260 9 : Surface2D::Edge const eu(s2d.edges[iu]);
261 9 : Real64 const xu(vu.x + (yHit - vu.y) * eXYu);
262 9 : if (xu < xHit) {
263 0 : return false; // Hit point is right of rightmost edge
264 : }
265 24 : while (u - l > 1u) {
266 15 : size_type const m((l + u) / 2);
267 15 : Slab::Edge const im(slabEdges[m]);
268 15 : Slab::EdgeXY const eXYm(slabEdgesXY[m]);
269 15 : Vertex2D const &vm(s2d.vertices[im]);
270 15 : Surface2D::Edge const em(s2d.edges[im]);
271 15 : Real64 xm(vm.x + (yHit - vm.y) * eXYm);
272 15 : if (xHit <= xm) {
273 6 : u = m;
274 : } else {
275 9 : l = m;
276 : }
277 15 : }
278 9 : assert(u - l == 1u);
279 9 : if (u % 2 == 0u) {
280 3 : return false; // Outside of nonconvex surface polygon
281 : }
282 12 : }
283 11 : return true;
284 17 : } // PierceSurface_nonconvex()
285 :
286 : ALWAYS_INLINE
287 : bool PierceSurface_polygon(DataSurfaces::SurfaceData const &surface, // Surface
288 : Vector3<Real64> const &hitPt // Ray-plane intersection point
289 : )
290 : {
291 : // Purpose: Check if hit point on surface plane is in surface polygon
292 : //
293 : // Author: Stuart Mentzer (Stuart_Mentzer@objexx.com)
294 : //
295 : // History:
296 : // Jan 2016: Initial release
297 :
298 : using DataSurfaces::nVerticesBig;
299 : using DataSurfaces::Surface2D;
300 : using Vertex2D = Vector2<Real64>;
301 99381 : Surface2D const &s2d(surface.surface2d);
302 99381 : int const axis(s2d.axis);
303 99381 : Vertex2D const h2d(axis == 0 ? hitPt.y : hitPt.x, axis == 2 ? hitPt.y : hitPt.z); // Hit point in 2D surface's plane
304 99381 : if ((h2d.x < s2d.vl.x) || (s2d.vu.x < h2d.x) || (h2d.y < s2d.vl.y) || (s2d.vu.y < h2d.y)) {
305 99181 : return false; // Misses 2D surface bounding box
306 : }
307 200 : ShapeCat const shapeCat(surface.shapeCat);
308 200 : if (shapeCat == ShapeCat::Rectangular) { // Rectangular is most common: Special case algorithm is faster but assumes these are really rectangular
309 176 : Vertex2D const v0h(h2d - s2d.vertices[0]);
310 176 : Real64 const he1(v0h.dot(s2d.edges[0]));
311 176 : if ((he1 < 0.0) || (he1 > s2d.s1)) {
312 0 : return false;
313 : }
314 176 : Real64 const he3(-v0h.dot(s2d.edges[3]));
315 176 : if ((he3 < 0.0) || (he3 > s2d.s3)) {
316 0 : return false;
317 : }
318 176 : return true;
319 200 : } else if (shapeCat == ShapeCat::Triangular) { // Cross products all nonnegative <=> Hit point in triangle
320 2 : return PierceSurface_Triangular(s2d, h2d);
321 29 : } else if ((shapeCat == ShapeCat::Nonconvex) ||
322 7 : (s2d.vertices.size() >= nVerticesBig)) { // O( log n ) algorithm for nonconvex and many-vertex convex surfaces
323 17 : return PierceSurface_Nonconvex(s2d, h2d);
324 5 : } else if (shapeCat == ShapeCat::Convex) { // O( n ) algorithm for convex surface without too many vertices
325 5 : return PierceSurface_Convex(s2d, h2d);
326 : } else {
327 0 : return false; // Should we assert here also?
328 : }
329 99381 : } // PierceSurface_Polygon()
330 :
331 : ALWAYS_INLINE
332 : bool PierceSurface(DataSurfaces::SurfaceData const &surface, // Surface
333 : Vector3<Real64> const &rayOri, // Ray origin point
334 : Vector3<Real64> const &rayDir, // Ray direction vector
335 : Vector3<Real64> &hitPt // Ray-plane intersection point
336 : )
337 : {
338 : // Purpose: Check if a ray hits a surface and return the point of intersection
339 : // with the surface's plane if they intersect.
340 : // Convex and concave surfaces with 3 or more vertices are supported.
341 : //
342 : // Author: Stuart Mentzer (Stuart_Mentzer@objexx.com)
343 : //
344 : // History:
345 : // Jan 2016: Initial release
346 :
347 : // Find ray intersection with surface plane
348 541095 : DataSurfaces::SurfaceData::Plane const &plane(surface.plane);
349 541095 : Real64 const den((plane.x * rayDir.x) + (plane.y * rayDir.y) + (plane.z * rayDir.z));
350 37 : if (den == 0.0) { // Ray is parallel to plane: This not treated as piercing even if ray lies in plane
351 2893 : return false;
352 : } else { // Ray's line intersects plane
353 538202 : Real64 const num(-((plane.x * rayOri.x) + (plane.y * rayOri.y) + (plane.z * rayOri.z) + plane.w));
354 538202 : if (num * den <=
355 : 0.0) { // Ray points away from surface or ray origin is on surface: This looks odd but is fast way to check for different signs
356 438822 : return false;
357 : } else { // Ray points toward surface: Compute hit point
358 99380 : Real64 const t(num / den); // Ray parameter at plane intersection: hitPt = rayOri + t * rayDir
359 99380 : hitPt.x = rayOri.x + (t * rayDir.x); // Compute by coordinate to avoid Vertex temporaries
360 99380 : hitPt.y = rayOri.y + (t * rayDir.y);
361 99380 : hitPt.z = rayOri.z + (t * rayDir.z);
362 : }
363 : }
364 :
365 : // Check if hit point is in surface polygon
366 99380 : return PierceSurface_polygon(surface, hitPt);
367 : } // PierceSurface()
368 :
369 : ALWAYS_INLINE
370 : bool PierceSurface(EnergyPlusData &state,
371 : int const iSurf, // Surface index
372 : Vector3<Real64> const &rayOri, // Ray origin point
373 : Vector3<Real64> const &rayDir, // Ray direction vector
374 : Vector3<Real64> &hitPt // Ray-plane intersection point
375 : )
376 : {
377 : // Purpose: Overload taking surface index instead of surface
378 : //
379 : // Author: Stuart Mentzer (Stuart_Mentzer@objexx.com)
380 : //
381 : // History:
382 : // Jan 2016: Initial release
383 :
384 1082116 : return PierceSurface(state.dataSurface->Surface(iSurf), rayOri, rayDir, hitPt);
385 : } // PierceSurface()
386 :
387 : ALWAYS_INLINE
388 : bool PierceSurface(DataSurfaces::SurfaceData const &surface, // Surface
389 : Vector3<Real64> const &rayOri, // Ray origin point
390 : Vector3<Real64> const &rayDir, // Ray direction unit vector
391 : Real64 const dMax, // Max distance from rayOri to hit point
392 : Vector3<Real64> &hitPt // Ray-plane intersection point
393 : )
394 : {
395 : // Purpose: Check if a ray hits a surface and return the point of intersection
396 : // with the surface's plane if they intersect.
397 : // Convex and concave surfaces with 3 or more vertices are supported.
398 : // This overload limits the ray-surface distance for a hit.
399 : //
400 : // Author: Stuart Mentzer (Stuart_Mentzer@objexx.com)
401 : //
402 : // History:
403 : // Jan 2016: Initial release
404 :
405 : // Input checks
406 454028 : assert(std::abs(rayDir.mag_squared() - 1.0) <
407 : 6 * std::numeric_limits<Real64>::epsilon()); // Check unit vector (6x is rough estimate. Increase slightly as needed.)
408 454026 : assert(dMax >= 0.0); // Distance must be nonnegative
409 :
410 : // Find ray intersection with surface plane
411 454028 : DataSurfaces::SurfaceData::Plane const &plane(surface.plane);
412 454028 : Real64 const den((plane.x * rayDir.x) + (plane.y * rayDir.y) + (plane.z * rayDir.z));
413 454028 : if (den == 0.0) { // Ray is parallel to plane: This not treated as piercing even if ray lies in plane
414 1315 : return false;
415 : } else { // Ray's line intersects plane
416 452713 : Real64 const num(-((plane.x * rayOri.x) + (plane.y * rayOri.y) + (plane.z * rayOri.z) + plane.w));
417 452713 : if (num * den <=
418 : 0.0) { // Ray points away from surface or ray origin is on surface: This looks odd but is fast way to check for different signs
419 388657 : return false;
420 : } else { // Ray points toward surface: Compute hit point
421 64056 : Real64 const t(num / den); // Ray parameter at plane intersection: hitPt = rayOri + t * rayDir
422 64056 : if (t > dMax) {
423 64055 : return false; // Hit point exceeds distance from rayOri limit
424 : }
425 1 : hitPt.x = rayOri.x + (t * rayDir.x); // Compute by coordinate to avoid Vertex temporaries
426 1 : hitPt.y = rayOri.y + (t * rayDir.y);
427 1 : hitPt.z = rayOri.z + (t * rayDir.z);
428 : }
429 : }
430 :
431 : // Check if hit point is in surface polygon
432 1 : return PierceSurface_polygon(surface, hitPt);
433 : } // PierceSurface()
434 :
435 : ALWAYS_INLINE
436 : bool PierceSurface(EnergyPlusData &state,
437 : int const iSurf, // Surface index
438 : Vector3<Real64> const &rayOri, // Ray origin point
439 : Vector3<Real64> const &rayDir, // Ray direction unit vector
440 : Real64 const dMax, // Max distance from rayOri to hit point
441 : Vector3<Real64> &hitPt // Ray-plane intersection point
442 : )
443 : {
444 : // Purpose: Overload taking surface index instead of surface
445 : //
446 : // Author: Stuart Mentzer (Stuart_Mentzer@objexx.com)
447 : //
448 : // History:
449 : // Jan 2016: Initial release
450 :
451 908052 : return PierceSurface(state.dataSurface->Surface(iSurf), rayOri, rayDir, dMax, hitPt);
452 : }
453 :
454 : } // namespace EnergyPlus
455 :
456 : #endif
|
