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