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_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) return false;
106 0 : if (es[1].cross(h2d - vs[1]) < 0.0) return false;
107 0 : if (es[2].cross(h2d - vs[2]) < 0.0) return false;
108 0 : return true;
109 : } // PierceSurface_Triangular()
110 :
111 23940 : inline bool PierceSurface_Convex(DataSurfaces::Surface2D const &s2d, // 2D surface
112 : Vector2<Real64> const &h2d // 2D hit point
113 : )
114 : {
115 : // Purpose: Check if a 2D hit point is in a convex 2D surface
116 : //
117 : // Author: Stuart Mentzer (Stuart_Mentzer@objexx.com)
118 : //
119 : // History:
120 : // Jan 2016: Initial release
121 : //
122 : // Notes:
123 : // Pulled this rare case out into separate function to facilitate inlining
124 : // This is O( n ) complexity so it is isn't used for many-vertex surfaces
125 :
126 : using DataSurfaces::Surface2D;
127 23940 : Surface2D::Vertices const &vs(s2d.vertices); // 2D surface vertices
128 23940 : Surface2D::Vectors const &es(s2d.edges); // 2D surface edge vectors
129 23940 : Surface2D::Vertices::size_type const n(vs.size());
130 23940 : assert(n >= 3u);
131 23940 : switch (n) {
132 0 : case 8:
133 0 : if (es[7].cross(h2d - vs[7]) < 0.0) {
134 0 : return false;
135 : }
136 : // fallthrough
137 : case 7:
138 0 : if (es[6].cross(h2d - vs[6]) < 0.0) {
139 0 : return false;
140 : }
141 : // fallthrough
142 : case 6:
143 0 : if (es[5].cross(h2d - vs[5]) < 0.0) {
144 0 : return false;
145 : }
146 : // fallthrough
147 : case 5:
148 0 : if (es[4].cross(h2d - vs[4]) < 0.0) {
149 0 : return false;
150 : }
151 : // fallthrough
152 : case 4:
153 23940 : if (es[3].cross(h2d - vs[3]) < 0.0) {
154 5978 : return false;
155 : }
156 : // fallthrough
157 : case 3:
158 17962 : if (es[2].cross(h2d - vs[2]) < 0.0) {
159 18 : return false;
160 : }
161 17944 : if (es[1].cross(h2d - vs[1]) < 0.0) {
162 18 : return false;
163 : }
164 17926 : if (es[0].cross(h2d - vs[0]) < 0.0) {
165 5902 : return false;
166 : }
167 12024 : return true;
168 0 : default:
169 0 : for (Surface2D::Vertices::size_type i = 0; i < n; ++i) {
170 0 : if (es[i].cross(h2d - vs[i]) < 0.0) return false;
171 : }
172 0 : return true;
173 : }
174 : } // PierceSurface_Convex()
175 :
176 0 : inline bool PierceSurface_Nonconvex(DataSurfaces::Surface2D const &s2d, // 2D surface
177 : Vector2<Real64> const &h2d // 2D hit point
178 : )
179 : {
180 : // Purpose: Check if a 2D hit point is in a 2D possibly nonconvex surface
181 : //
182 : // Author: Stuart Mentzer (Stuart_Mentzer@objexx.com)
183 : //
184 : // History:
185 : // Jan 2016: Initial release
186 : //
187 : // Notes:
188 : // Pulled this rare case out into separate function to facilitate inlining
189 : // This works for nonconvex "simple" (no edge crossings) polygons
190 : // This is also a fast O( log n ) algorithm for many-vertex convex surfaces
191 :
192 : using DataSurfaces::Surface2D;
193 : using size_type = Surface2D::Vertices::size_type;
194 : using Slab = DataSurfaces::Surface2DSlab;
195 : using Vertex2D = Vector2<Real64>;
196 0 : assert(s2d.vertices.size() >= 3u);
197 0 : Surface2D::Slabs const &slabs(s2d.slabs); // 2D surface y slice slabs
198 0 : Surface2D::SlabYs const &slabYs(s2d.slabYs); // 2D surface slab y coordinates
199 0 : assert(slabYs.size() > 0u);
200 0 : Real64 const yHit(h2d.y); // Hit point y coordinate
201 :
202 : // Find slab with y range containing hit point
203 0 : auto const iHit(std::lower_bound(slabYs.begin(), slabYs.end(), yHit));
204 0 : assert((yHit >= slabYs.front()) && (yHit <= slabYs.back())); // Passed bounding box check so hit point in slabs y range
205 0 : assert(iHit != slabYs.end()); // Hit point can't be above all slabs: passed bounding box check
206 0 : size_type const iSlab(std::min(static_cast<size_type>(iHit - 1 - slabYs.begin()), slabs.size())); // Hit slab index
207 0 : Slab const &slab(slabs[iSlab]);
208 :
209 : // Check hit point within slab bounding box x range
210 0 : Real64 const xHit(h2d.x); // Hit point x coordinate
211 0 : if ((xHit < slab.xl) || (xHit > slab.xu)) return false; // Hit point outside slab bounding box
212 :
213 : // Find edge pair surrounding hit point
214 0 : Slab::Edges const &slabEdges(slab.edges);
215 0 : Slab::EdgesXY const &slabEdgesXY(slab.edgesXY);
216 0 : size_type const nEdges(slabEdges.size());
217 0 : assert(nEdges >= 2u);
218 0 : if (nEdges == 2) { // 2 edges
219 0 : Slab::Edge const se0(slabEdges[0]);
220 0 : Slab::EdgeXY const eXY0(slabEdgesXY[0]);
221 0 : Vertex2D v0(s2d.vertices[se0]);
222 0 : Surface2D::Edge e0(s2d.edges[se0]);
223 0 : Real64 const x0(v0.x + (yHit - v0.y) * eXY0);
224 0 : if (xHit < x0) return false; // Hit point x is left of left edge
225 0 : Slab::Edge const se1(slabEdges[1]);
226 0 : Slab::EdgeXY const eXY1(slabEdgesXY[1]);
227 0 : Vertex2D v1(s2d.vertices[se1]);
228 0 : Surface2D::Edge e1(s2d.edges[se1]);
229 0 : Real64 const x1(v1.x + (yHit - v1.y) * eXY1);
230 0 : if (x1 < xHit) return false; // Hit point is right of right edge
231 0 : } else { // 4+ edges: Binary search for edges surrounding hit point
232 0 : assert(nEdges >= 4u);
233 0 : assert(nEdges % 2 == 0u);
234 0 : size_type l(0u), u(nEdges - 1);
235 0 : Slab::Edge const il(slabEdges[l]);
236 0 : Slab::EdgeXY const eXYl(slabEdgesXY[l]);
237 0 : Vertex2D const &vl(s2d.vertices[il]);
238 0 : Surface2D::Edge const el(s2d.edges[il]);
239 0 : Real64 const xl(vl.x + (yHit - vl.y) * eXYl);
240 0 : if (xHit < xl) return false; // Hit point x is left of leftmost edge
241 0 : Slab::Edge const iu(slabEdges[u]);
242 0 : Slab::EdgeXY const eXYu(slabEdgesXY[u]);
243 0 : Vertex2D const &vu(s2d.vertices[iu]);
244 0 : Surface2D::Edge const eu(s2d.edges[iu]);
245 0 : Real64 const xu(vu.x + (yHit - vu.y) * eXYu);
246 0 : if (xu < xHit) return false; // Hit point is right of rightmost edge
247 0 : while (u - l > 1u) {
248 0 : size_type const m((l + u) / 2);
249 0 : Slab::Edge const im(slabEdges[m]);
250 0 : Slab::EdgeXY const eXYm(slabEdgesXY[m]);
251 0 : Vertex2D const &vm(s2d.vertices[im]);
252 0 : Surface2D::Edge const em(s2d.edges[im]);
253 0 : Real64 xm(vm.x + (yHit - vm.y) * eXYm);
254 0 : if (xHit <= xm) {
255 0 : u = m;
256 : } else {
257 0 : l = m;
258 : }
259 0 : }
260 0 : assert(u - l == 1u);
261 0 : if (u % 2 == 0u) return false; // Outside of nonconvex surface polygon
262 0 : }
263 0 : return true;
264 0 : } // PierceSurface_nonconvex()
265 :
266 : ALWAYS_INLINE
267 : bool PierceSurface_polygon(DataSurfaces::SurfaceData const &surface, // Surface
268 : Vector3<Real64> const &hitPt // Ray-plane intersection point
269 : )
270 : {
271 : // Purpose: Check if hit point on surface plane is in surface polygon
272 : //
273 : // Author: Stuart Mentzer (Stuart_Mentzer@objexx.com)
274 : //
275 : // History:
276 : // Jan 2016: Initial release
277 :
278 : using DataSurfaces::nVerticesBig;
279 : using DataSurfaces::Surface2D;
280 : using Vertex2D = Vector2<Real64>;
281 48937276 : Surface2D const &s2d(surface.surface2d);
282 48937276 : int const axis(s2d.axis);
283 48937276 : Vertex2D const h2d(axis == 0 ? hitPt.y : hitPt.x, axis == 2 ? hitPt.y : hitPt.z); // Hit point in 2D surface's plane
284 48937276 : if ((h2d.x < s2d.vl.x) || (s2d.vu.x < h2d.x) || (h2d.y < s2d.vl.y) || (s2d.vu.y < h2d.y)) return false; // Misses 2D surface bounding box
285 1042540 : ShapeCat const shapeCat(surface.shapeCat);
286 1042540 : if (shapeCat == ShapeCat::Rectangular) { // Rectangular is most common: Special case algorithm is faster but assumes these are really rectangular
287 1018600 : Vertex2D const v0h(h2d - s2d.vertices[0]);
288 1018600 : Real64 const he1(v0h.dot(s2d.edges[0]));
289 1018600 : if ((he1 < 0.0) || (he1 > s2d.s1)) return false;
290 983923 : Real64 const he3(-v0h.dot(s2d.edges[3]));
291 983923 : if ((he3 < 0.0) || (he3 > s2d.s3)) return false;
292 976411 : return true;
293 1042540 : } else if (shapeCat == ShapeCat::Triangular) { // Cross products all nonnegative <=> Hit point in triangle
294 0 : return PierceSurface_Triangular(s2d, h2d);
295 47880 : } else if ((shapeCat == ShapeCat::Nonconvex) ||
296 23940 : (s2d.vertices.size() >= nVerticesBig)) { // O( log n ) algorithm for nonconvex and many-vertex convex surfaces
297 0 : return PierceSurface_Nonconvex(s2d, h2d);
298 23940 : } else if (shapeCat == ShapeCat::Convex) { // O( n ) algorithm for convex surface without too many vertices
299 23940 : return PierceSurface_Convex(s2d, h2d);
300 : } else {
301 0 : return false; // Should we assert here also?
302 : }
303 48937276 : } // PierceSurface_Polygon()
304 :
305 : ALWAYS_INLINE
306 : bool PierceSurface(DataSurfaces::SurfaceData const &surface, // Surface
307 : Vector3<Real64> const &rayOri, // Ray origin point
308 : Vector3<Real64> const &rayDir, // Ray direction vector
309 : Vector3<Real64> &hitPt // Ray-plane intersection point
310 : )
311 : {
312 : // Purpose: Check if a ray hits a surface and return the point of intersection
313 : // with the surface's plane if they intersect.
314 : // Convex and concave surfaces with 3 or more vertices are supported.
315 : //
316 : // Author: Stuart Mentzer (Stuart_Mentzer@objexx.com)
317 : //
318 : // History:
319 : // Jan 2016: Initial release
320 :
321 : // Find ray intersection with surface plane
322 120863126 : DataSurfaces::SurfaceData::Plane const &plane(surface.plane);
323 120863126 : Real64 const den((plane.x * rayDir.x) + (plane.y * rayDir.y) + (plane.z * rayDir.z));
324 68714700 : if (den == 0.0) { // Ray is parallel to plane: This not treated as piercing even if ray lies in plane
325 7330 : return false;
326 : } else { // Ray's line intersects plane
327 120855796 : Real64 const num(-((plane.x * rayOri.x) + (plane.y * rayOri.y) + (plane.z * rayOri.z) + plane.w));
328 120855796 : if (num * den <=
329 : 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
330 72085961 : return false;
331 : } else { // Ray points toward surface: Compute hit point
332 48769835 : Real64 const t(num / den); // Ray parameter at plane intersection: hitPt = rayOri + t * rayDir
333 48769835 : hitPt.x = rayOri.x + (t * rayDir.x); // Compute by coordinate to avoid Vertex temporaries
334 48769835 : hitPt.y = rayOri.y + (t * rayDir.y);
335 48769835 : hitPt.z = rayOri.z + (t * rayDir.z);
336 : }
337 : }
338 :
339 : // Check if hit point is in surface polygon
340 48769835 : return PierceSurface_polygon(surface, hitPt);
341 : } // PierceSurface()
342 :
343 : ALWAYS_INLINE
344 : bool PierceSurface(EnergyPlusData &state,
345 : int const iSurf, // Surface index
346 : Vector3<Real64> const &rayOri, // Ray origin point
347 : Vector3<Real64> const &rayDir, // Ray direction vector
348 : Vector3<Real64> &hitPt // Ray-plane intersection point
349 : )
350 : {
351 : // Purpose: Overload taking surface index instead of surface
352 : //
353 : // Author: Stuart Mentzer (Stuart_Mentzer@objexx.com)
354 : //
355 : // History:
356 : // Jan 2016: Initial release
357 :
358 112237498 : return PierceSurface(state.dataSurface->Surface(iSurf), rayOri, rayDir, hitPt);
359 : } // PierceSurface()
360 :
361 : ALWAYS_INLINE
362 : bool PierceSurface(DataSurfaces::SurfaceData const &surface, // Surface
363 : Vector3<Real64> const &rayOri, // Ray origin point
364 : Vector3<Real64> const &rayDir, // Ray direction unit vector
365 : Real64 const dMax, // Max distance from rayOri to hit point
366 : Vector3<Real64> &hitPt // Ray-plane intersection point
367 : )
368 : {
369 : // Purpose: Check if a ray hits a surface and return the point of intersection
370 : // with the surface's plane if they intersect.
371 : // Convex and concave surfaces with 3 or more vertices are supported.
372 : // This overload limits the ray-surface distance for a hit.
373 : //
374 : // Author: Stuart Mentzer (Stuart_Mentzer@objexx.com)
375 : //
376 : // History:
377 : // Jan 2016: Initial release
378 :
379 : // Input checks
380 14850212 : assert(std::abs(rayDir.mag_squared() - 1.0) <
381 : 6 * std::numeric_limits<Real64>::epsilon()); // Check unit vector (6x is rough estimate. Increase slightly as needed.)
382 14850212 : assert(dMax >= 0.0); // Distance must be nonnegative
383 :
384 : // Find ray intersection with surface plane
385 14850212 : DataSurfaces::SurfaceData::Plane const &plane(surface.plane);
386 14850212 : Real64 const den((plane.x * rayDir.x) + (plane.y * rayDir.y) + (plane.z * rayDir.z));
387 14850212 : if (den == 0.0) { // Ray is parallel to plane: This not treated as piercing even if ray lies in plane
388 1335 : return false;
389 : } else { // Ray's line intersects plane
390 14848877 : Real64 const num(-((plane.x * rayOri.x) + (plane.y * rayOri.y) + (plane.z * rayOri.z) + plane.w));
391 14848877 : if (num * den <=
392 : 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
393 8564668 : return false;
394 : } else { // Ray points toward surface: Compute hit point
395 6284209 : Real64 const t(num / den); // Ray parameter at plane intersection: hitPt = rayOri + t * rayDir
396 6284209 : if (t > dMax) return false; // Hit point exceeds distance from rayOri limit
397 167441 : hitPt.x = rayOri.x + (t * rayDir.x); // Compute by coordinate to avoid Vertex temporaries
398 167441 : hitPt.y = rayOri.y + (t * rayDir.y);
399 167441 : hitPt.z = rayOri.z + (t * rayDir.z);
400 : }
401 : }
402 :
403 : // Check if hit point is in surface polygon
404 167441 : return PierceSurface_polygon(surface, hitPt);
405 : } // PierceSurface()
406 :
407 : ALWAYS_INLINE
408 : bool PierceSurface(EnergyPlusData &state,
409 : int const iSurf, // Surface index
410 : Vector3<Real64> const &rayOri, // Ray origin point
411 : Vector3<Real64> const &rayDir, // Ray direction unit vector
412 : Real64 const dMax, // Max distance from rayOri to hit point
413 : Vector3<Real64> &hitPt // Ray-plane intersection point
414 : )
415 : {
416 : // Purpose: Overload taking surface index instead of surface
417 : //
418 : // Author: Stuart Mentzer (Stuart_Mentzer@objexx.com)
419 : //
420 : // History:
421 : // Jan 2016: Initial release
422 :
423 24867122 : return PierceSurface(state.dataSurface->Surface(iSurf), rayOri, rayDir, dMax, hitPt);
424 : }
425 :
426 : } // namespace EnergyPlus
427 :
428 : #endif
|
