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 0 : 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 0 : Surface2D::Vertices const &vs(s2d.vertices); // 2D surface vertices
104 0 : Surface2D::Vectors const &es(s2d.edges); // 2D surface edge vectors
105 0 : if (es[0].cross(h2d - vs[0]) < 0.0) {
106 0 : return false;
107 : }
108 0 : if (es[1].cross(h2d - vs[1]) < 0.0) {
109 0 : return false;
110 : }
111 0 : if (es[2].cross(h2d - vs[2]) < 0.0) {
112 0 : return false;
113 : }
114 0 : return true;
115 : } // PierceSurface_Triangular()
116 :
117 23976 : 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 23976 : Surface2D::Vertices const &vs(s2d.vertices); // 2D surface vertices
134 23976 : Surface2D::Vectors const &es(s2d.edges); // 2D surface edge vectors
135 23976 : Surface2D::Vertices::size_type const n(vs.size());
136 23976 : assert(n >= 3u);
137 23976 : switch (n) {
138 0 : case 8:
139 0 : if (es[7].cross(h2d - vs[7]) < 0.0) {
140 0 : return false;
141 : }
142 : // fallthrough
143 : case 7:
144 0 : if (es[6].cross(h2d - vs[6]) < 0.0) {
145 0 : return false;
146 : }
147 : // fallthrough
148 : case 6:
149 0 : if (es[5].cross(h2d - vs[5]) < 0.0) {
150 0 : return false;
151 : }
152 : // fallthrough
153 : case 5:
154 0 : if (es[4].cross(h2d - vs[4]) < 0.0) {
155 0 : return false;
156 : }
157 : // fallthrough
158 : case 4:
159 23976 : if (es[3].cross(h2d - vs[3]) < 0.0) {
160 5978 : return false;
161 : }
162 : // fallthrough
163 : case 3:
164 17998 : if (es[2].cross(h2d - vs[2]) < 0.0) {
165 18 : return false;
166 : }
167 17980 : if (es[1].cross(h2d - vs[1]) < 0.0) {
168 18 : return false;
169 : }
170 17962 : if (es[0].cross(h2d - vs[0]) < 0.0) {
171 5902 : return false;
172 : }
173 12060 : 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 0 : 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 0 : assert(s2d.vertices.size() >= 3u);
205 0 : Surface2D::Slabs const &slabs(s2d.slabs); // 2D surface y slice slabs
206 0 : Surface2D::SlabYs const &slabYs(s2d.slabYs); // 2D surface slab y coordinates
207 0 : assert(slabYs.size() > 0u);
208 0 : Real64 const yHit(h2d.y); // Hit point y coordinate
209 :
210 : // Find slab with y range containing hit point
211 0 : auto const iHit(std::lower_bound(slabYs.begin(), slabYs.end(), yHit));
212 0 : assert((yHit >= slabYs.front()) && (yHit <= slabYs.back())); // Passed bounding box check so hit point in slabs y range
213 0 : assert(iHit != slabYs.end()); // Hit point can't be above all slabs: passed bounding box check
214 0 : size_type const iSlab(std::min(static_cast<size_type>(iHit - 1 - slabYs.begin()), slabs.size())); // Hit slab index
215 0 : Slab const &slab(slabs[iSlab]);
216 :
217 : // Check hit point within slab bounding box x range
218 0 : Real64 const xHit(h2d.x); // Hit point x coordinate
219 0 : if ((xHit < slab.xl) || (xHit > slab.xu)) {
220 0 : return false; // Hit point outside slab bounding box
221 : }
222 :
223 : // Find edge pair surrounding hit point
224 0 : Slab::Edges const &slabEdges(slab.edges);
225 0 : Slab::EdgesXY const &slabEdgesXY(slab.edgesXY);
226 0 : size_type const nEdges(slabEdges.size());
227 0 : assert(nEdges >= 2u);
228 0 : if (nEdges == 2) { // 2 edges
229 0 : Slab::Edge const se0(slabEdges[0]);
230 0 : Slab::EdgeXY const eXY0(slabEdgesXY[0]);
231 0 : Vertex2D v0(s2d.vertices[se0]);
232 0 : Surface2D::Edge e0(s2d.edges[se0]);
233 0 : Real64 const x0(v0.x + (yHit - v0.y) * eXY0);
234 0 : if (xHit < x0) {
235 0 : return false; // Hit point x is left of left edge
236 : }
237 0 : Slab::Edge const se1(slabEdges[1]);
238 0 : Slab::EdgeXY const eXY1(slabEdgesXY[1]);
239 0 : Vertex2D v1(s2d.vertices[se1]);
240 0 : Surface2D::Edge e1(s2d.edges[se1]);
241 0 : Real64 const x1(v1.x + (yHit - v1.y) * eXY1);
242 0 : if (x1 < xHit) {
243 0 : return false; // Hit point is right of right edge
244 : }
245 0 : } else { // 4+ edges: Binary search for edges surrounding hit point
246 0 : assert(nEdges >= 4u);
247 0 : assert(nEdges % 2 == 0u);
248 0 : size_type l(0u), u(nEdges - 1);
249 0 : Slab::Edge const il(slabEdges[l]);
250 0 : Slab::EdgeXY const eXYl(slabEdgesXY[l]);
251 0 : Vertex2D const &vl(s2d.vertices[il]);
252 0 : Surface2D::Edge const el(s2d.edges[il]);
253 0 : Real64 const xl(vl.x + (yHit - vl.y) * eXYl);
254 0 : if (xHit < xl) {
255 0 : return false; // Hit point x is left of leftmost edge
256 : }
257 0 : Slab::Edge const iu(slabEdges[u]);
258 0 : Slab::EdgeXY const eXYu(slabEdgesXY[u]);
259 0 : Vertex2D const &vu(s2d.vertices[iu]);
260 0 : Surface2D::Edge const eu(s2d.edges[iu]);
261 0 : Real64 const xu(vu.x + (yHit - vu.y) * eXYu);
262 0 : if (xu < xHit) {
263 0 : return false; // Hit point is right of rightmost edge
264 : }
265 0 : while (u - l > 1u) {
266 0 : size_type const m((l + u) / 2);
267 0 : Slab::Edge const im(slabEdges[m]);
268 0 : Slab::EdgeXY const eXYm(slabEdgesXY[m]);
269 0 : Vertex2D const &vm(s2d.vertices[im]);
270 0 : Surface2D::Edge const em(s2d.edges[im]);
271 0 : Real64 xm(vm.x + (yHit - vm.y) * eXYm);
272 0 : if (xHit <= xm) {
273 0 : u = m;
274 : } else {
275 0 : l = m;
276 : }
277 0 : }
278 0 : assert(u - l == 1u);
279 0 : if (u % 2 == 0u) {
280 0 : return false; // Outside of nonconvex surface polygon
281 : }
282 0 : }
283 0 : return true;
284 0 : } // 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 52430718 : Surface2D const &s2d(surface.surface2d);
302 52430718 : int const axis(s2d.axis);
303 52430718 : Vertex2D const h2d(axis == 0 ? hitPt.y : hitPt.x, axis == 2 ? hitPt.y : hitPt.z); // Hit point in 2D surface's plane
304 52430718 : if ((h2d.x < s2d.vl.x) || (s2d.vu.x < h2d.x) || (h2d.y < s2d.vl.y) || (s2d.vu.y < h2d.y)) {
305 51354136 : return false; // Misses 2D surface bounding box
306 : }
307 1076582 : ShapeCat const shapeCat(surface.shapeCat);
308 1076582 : if (shapeCat == ShapeCat::Rectangular) { // Rectangular is most common: Special case algorithm is faster but assumes these are really rectangular
309 1052606 : Vertex2D const v0h(h2d - s2d.vertices[0]);
310 1052606 : Real64 const he1(v0h.dot(s2d.edges[0]));
311 1052606 : if ((he1 < 0.0) || (he1 > s2d.s1)) {
312 34677 : return false;
313 : }
314 1017929 : Real64 const he3(-v0h.dot(s2d.edges[3]));
315 1017929 : if ((he3 < 0.0) || (he3 > s2d.s3)) {
316 7512 : return false;
317 : }
318 1010417 : return true;
319 1076582 : } else if (shapeCat == ShapeCat::Triangular) { // Cross products all nonnegative <=> Hit point in triangle
320 0 : return PierceSurface_Triangular(s2d, h2d);
321 47952 : } else if ((shapeCat == ShapeCat::Nonconvex) ||
322 23976 : (s2d.vertices.size() >= nVerticesBig)) { // O( log n ) algorithm for nonconvex and many-vertex convex surfaces
323 0 : return PierceSurface_Nonconvex(s2d, h2d);
324 23976 : } else if (shapeCat == ShapeCat::Convex) { // O( n ) algorithm for convex surface without too many vertices
325 23976 : return PierceSurface_Convex(s2d, h2d);
326 : } else {
327 0 : return false; // Should we assert here also?
328 : }
329 52430718 : } // 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 128188365 : DataSurfaces::SurfaceData::Plane const &plane(surface.plane);
349 128188365 : Real64 const den((plane.x * rayDir.x) + (plane.y * rayDir.y) + (plane.z * rayDir.z));
350 76035379 : if (den == 0.0) { // Ray is parallel to plane: This not treated as piercing even if ray lies in plane
351 8539 : return false;
352 : } else { // Ray's line intersects plane
353 128179826 : Real64 const num(-((plane.x * rayOri.x) + (plane.y * rayOri.y) + (plane.z * rayOri.z) + plane.w));
354 128179826 : 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 75916562 : return false;
357 : } else { // Ray points toward surface: Compute hit point
358 52263264 : Real64 const t(num / den); // Ray parameter at plane intersection: hitPt = rayOri + t * rayDir
359 52263264 : hitPt.x = rayOri.x + (t * rayDir.x); // Compute by coordinate to avoid Vertex temporaries
360 52263264 : hitPt.y = rayOri.y + (t * rayDir.y);
361 52263264 : hitPt.z = rayOri.z + (t * rayDir.z);
362 : }
363 : }
364 :
365 : // Check if hit point is in surface polygon
366 52263264 : 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 112246618 : 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 15473387 : 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 15473387 : assert(dMax >= 0.0); // Distance must be nonnegative
409 :
410 : // Find ray intersection with surface plane
411 15473387 : DataSurfaces::SurfaceData::Plane const &plane(surface.plane);
412 15473387 : Real64 const den((plane.x * rayDir.x) + (plane.y * rayDir.y) + (plane.z * rayDir.z));
413 15473387 : if (den == 0.0) { // Ray is parallel to plane: This not treated as piercing even if ray lies in plane
414 1426 : return false;
415 : } else { // Ray's line intersects plane
416 15471961 : Real64 const num(-((plane.x * rayOri.x) + (plane.y * rayOri.y) + (plane.z * rayOri.z) + plane.w));
417 15471961 : 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 8768288 : return false;
420 : } else { // Ray points toward surface: Compute hit point
421 6703673 : Real64 const t(num / den); // Ray parameter at plane intersection: hitPt = rayOri + t * rayDir
422 6703673 : if (t > dMax) {
423 6536219 : return false; // Hit point exceeds distance from rayOri limit
424 : }
425 167454 : hitPt.x = rayOri.x + (t * rayDir.x); // Compute by coordinate to avoid Vertex temporaries
426 167454 : hitPt.y = rayOri.y + (t * rayDir.y);
427 167454 : hitPt.z = rayOri.z + (t * rayDir.z);
428 : }
429 : }
430 :
431 : // Check if hit point is in surface polygon
432 167454 : 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 24867122 : return PierceSurface(state.dataSurface->Surface(iSurf), rayOri, rayDir, dMax, hitPt);
452 : }
453 :
454 : } // namespace EnergyPlus
455 :
456 : #endif
|
