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 : // C++ Headers
49 : #include <cmath>
50 :
51 : // EnergyPlus Headers
52 : #include <EnergyPlus/Data/EnergyPlusData.hh>
53 : #include <EnergyPlus/DataGlobals.hh>
54 : #include <EnergyPlus/Vectors.hh>
55 :
56 : namespace EnergyPlus::Vectors {
57 : // Module containing the routines dealing with Vector operations
58 :
59 : // MODULE INFORMATION:
60 : // AUTHOR Linda Lawrie
61 : // DATE WRITTEN April 2000
62 :
63 : // PURPOSE OF THIS MODULE:
64 : // This module uses a global "vector" data structure and defines
65 : // operations (using the F90 "operator" features) and other
66 : // manipulations used with Vectors.
67 :
68 : // Vectors should allow for more modern interpretations of the
69 : // calculations used in the shadowing and other surface manipulations.
70 :
71 : // METHODOLOGY EMPLOYED:
72 : // Uses the "Vector" derived type variables to allow for more readable
73 : // calculations.
74 :
75 : // REFERENCES: Original idea from F90 for Scientists and
76 : // Engineers, S. J. Chapman, 1996.
77 : // The module defines 8 operations which can be performed on vectors:
78 : // Operation Operator
79 : // ========= ========
80 : // 1. Creation from a real array =
81 : // 2. Conversion to real array =
82 : // 3. Vector addition +
83 : // 4. Vector subtraction -
84 : // 5. Vector-scalar multiplication (4 cases) *
85 : // 6. Vector-scalar division (2 cases) /
86 : // 7. Dot product .dot.
87 : // 8. Cross product *
88 : // 9. 2d dot product .twoddot.
89 : // 10. 2d Cross product .twodcross.
90 : // It contains a total of 12 procedures to implement those
91 : // operations: array_to_vector, vector_to_array, vector_add,
92 : // vector_subtract, vector_times_real, real_times_vector,
93 : // vector_times_int, int_times_vector, vector_div_real,
94 : // vector_div_int, dot_product, and cross_product.
95 :
96 : // Using/Aliasing
97 : using namespace DataVectorTypes;
98 :
99 : // MODULE PARAMETER DEFINITIONS
100 :
101 : // Object Data
102 : Vector const XUnit(1.0, 0.0, 0.0);
103 : Vector const YUnit(0.0, 1.0, 0.0);
104 : Vector const ZUnit(0.0, 0.0, 1.0);
105 :
106 88754 : Real64 AreaPolygon(int const n, Array1D<Vector> &p)
107 : {
108 :
109 : // PURPOSE OF THIS SUBROUTINE:
110 : // This subroutine calculates the area of a polygon defined by the input vectors.
111 :
112 : // REFERENCE:
113 : // Graphic Gems.
114 :
115 : // Argument array dimensioning
116 88754 : EP_SIZE_CHECK(p, n);
117 :
118 88754 : Vector edge0 = p[1] - p[0];
119 88754 : Vector edge1 = p[2] - p[0];
120 :
121 88754 : Vector edgex = cross(edge0, edge1);
122 88754 : Vector nor = VecNormalize(edgex);
123 :
124 : // Initialize csum
125 88754 : Vector csum;
126 88754 : csum = 0.0;
127 :
128 266262 : for (int i = 0; i <= n - 2; ++i) {
129 177508 : csum += cross(p[i], p[i + 1]);
130 : }
131 88754 : csum += cross(p[n - 1], p[0]);
132 :
133 88754 : Real64 areap = 0.5 * std::abs(dot(nor, csum));
134 :
135 88754 : return areap;
136 88754 : }
137 :
138 674036 : Real64 VecSquaredLength(Vector const &vec)
139 : {
140 :
141 : // PURPOSE OF THIS SUBROUTINE:
142 : // This subroutine calculates the squared length of the input vector.
143 :
144 : // REFERENCE:
145 : // Graphic Gems.
146 :
147 674036 : Real64 vecsqlen = (vec.x * vec.x + vec.y * vec.y + vec.z * vec.z);
148 :
149 674036 : return vecsqlen;
150 : }
151 :
152 674036 : Real64 VecLength(Vector const &vec)
153 : {
154 :
155 : // PURPOSE OF THIS SUBROUTINE:
156 : // This subroutine calculates the length of the input vector.
157 :
158 : // REFERENCE:
159 : // Graphic Gems.
160 :
161 674036 : Real64 veclen = std::sqrt(VecSquaredLength(vec));
162 :
163 674036 : return veclen;
164 : }
165 :
166 0 : Vector VecNegate(Vector const &vec)
167 : {
168 :
169 : // PURPOSE OF THIS SUBROUTINE:
170 : // This subroutine negates the input vector and returns that vector.
171 :
172 : // REFERENCE:
173 : // Graphic Gems.
174 :
175 0 : Vector VecNegate;
176 :
177 0 : VecNegate.x = -vec.x;
178 0 : VecNegate.y = -vec.y;
179 0 : VecNegate.z = -vec.z;
180 :
181 0 : return VecNegate;
182 : }
183 :
184 277498 : Vector VecNormalize(Vector const &vec)
185 : {
186 :
187 : // PURPOSE OF THIS SUBROUTINE:
188 : // This subroutine normalizes the input vector and returns the normalized vector
189 :
190 : // REFERENCE:
191 : // Graphic Gems.
192 :
193 : // Return value
194 277498 : Vector VecNormalize;
195 :
196 277498 : Real64 veclen = VecLength(vec);
197 277498 : if (veclen != 0.0) {
198 277498 : VecNormalize.x = vec.x / veclen;
199 277498 : VecNormalize.y = vec.y / veclen;
200 277498 : VecNormalize.z = vec.z / veclen;
201 : } else {
202 0 : VecNormalize.x = 0.0;
203 0 : VecNormalize.y = 0.0;
204 0 : VecNormalize.z = 0.0;
205 : }
206 :
207 277498 : return VecNormalize;
208 : }
209 :
210 0 : void VecRound(Vector &vec, Real64 const roundto)
211 : {
212 :
213 : // PURPOSE OF THIS SUBROUTINE:
214 : // This subroutine rounds the input vector to a specified "rounding" value.
215 :
216 : // REFERENCE:
217 : // Graphic Gems.
218 :
219 0 : vec.x = nint64(vec.x * roundto) / roundto;
220 0 : vec.y = nint64(vec.y * roundto) / roundto;
221 0 : vec.z = nint64(vec.z * roundto) / roundto;
222 0 : }
223 :
224 51473 : void DetermineAzimuthAndTilt(Array1D<Vector> const &Surf, // Surface Definition
225 : Real64 &Azimuth, // Outward Normal Azimuth Angle
226 : Real64 &Tilt, // Tilt angle of surface
227 : Vector &lcsx,
228 : Vector &lcsy,
229 : Vector &lcsz,
230 : Vector const &NewellSurfaceNormalVector)
231 : {
232 :
233 : // PURPOSE OF THIS SUBROUTINE:
234 : // This subroutine determines the Azimuth (outward normal) angle,
235 : // Tilt angle of a given surface defined by the set of input vectors.
236 :
237 : // REFERENCE:
238 : // Discussions and examples from Bill Carroll, LBNL.
239 :
240 51473 : lcsx = VecNormalize(Surf(3) - Surf(2));
241 51473 : lcsz = NewellSurfaceNormalVector;
242 51473 : lcsy = cross(lcsz, lcsx);
243 :
244 51473 : Real64 costheta = dot(lcsz, ZUnit);
245 :
246 : // if ( fabs(costheta) < 1.0d0) { // normal cases
247 51473 : Real64 constexpr epsilon = 1.12e-16;
248 51473 : Real64 rotang_0 = 0.0;
249 51473 : if (std::abs(costheta) < 1.0 - epsilon) { // Autodesk Added - 1.12e-16 to treat 1 bit from 1.0 as 1.0 to correct different behavior seen in
250 : // release vs debug build due to slight precision differences: May want larger epsilon here
251 : // azimuth
252 35515 : Vector x2 = cross(ZUnit, lcsz);
253 35515 : rotang_0 = std::atan2(dot(x2, YUnit), dot(x2, XUnit));
254 35515 : } else {
255 : // azimuth
256 15958 : rotang_0 = std::atan2(dot(lcsx, YUnit), dot(lcsx, XUnit));
257 : }
258 :
259 51473 : Real64 tlt = std::acos(NewellSurfaceNormalVector.z);
260 51473 : tlt /= Constant::DegToRad;
261 :
262 51473 : Real64 az = rotang_0;
263 :
264 51473 : az /= Constant::DegToRad;
265 51473 : az = mod(450.0 - az, 360.0);
266 51473 : az += 90.0;
267 51473 : if (az < 0.0) {
268 0 : az += 360.0;
269 : }
270 51473 : az = mod(az, 360.0);
271 :
272 : // Clean up angle precision
273 51473 : if (std::abs(az - 360.0) < Constant::OneThousandth) { // Bring small angles to zero
274 1 : az = 0.0;
275 51472 : } else if (std::abs(az - 180.0) <
276 : Constant::OneMillionth) { // Bring angles near 180 to 180 //Autodesk Added to clean up debug--release discrepancies
277 13124 : az = 180.0;
278 : }
279 51473 : if (std::abs(tlt - 180.0) < Constant::OneMillionth) { // Bring angles near 180 to 180 //Autodesk Added to clean up debug--release discrepancies
280 9114 : tlt = 180.0;
281 : }
282 :
283 51473 : Azimuth = az;
284 51473 : Tilt = tlt;
285 51473 : }
286 :
287 50764 : void PlaneEquation(Array1D<Vector> &verts, // Structure of the surface
288 : int const nverts, // Number of vertices in the surface
289 : PlaneEq &plane, // Equation of plane from inputs
290 : bool &error // returns true for degenerate surface
291 : )
292 : {
293 :
294 : // PURPOSE OF THIS SUBROUTINE:
295 : // This subroutine calculates the plane equation for a given surface (which should be planar).
296 :
297 : // REFERENCE:
298 : // Graphic Gems
299 :
300 : // Argument array dimensioning
301 50764 : EP_SIZE_CHECK(verts, nverts);
302 :
303 50764 : Vector normal = Vector(0.0, 0.0, 0.0);
304 50764 : Vector refpt = Vector(0.0, 0.0, 0.0);
305 254747 : for (int i = 0; i <= nverts - 1; ++i) {
306 203983 : Vector const &u(verts[i]);
307 203983 : Vector const &v(i < nverts - 1 ? verts[i + 1] : verts[0]);
308 203983 : normal.x += (u.y - v.y) * (u.z + v.z);
309 203983 : normal.y += (u.z - v.z) * (u.x + v.x);
310 203983 : normal.z += (u.x - v.x) * (u.y + v.y);
311 203983 : refpt += u;
312 : }
313 : // normalize the polygon normal to obtain the first
314 : // three coefficients of the plane equation
315 50764 : Real64 lenvec = VecLength(normal);
316 50764 : error = false;
317 50764 : if (lenvec != 0.0) { // should this be >0
318 50764 : plane.x = normal.x / lenvec;
319 50764 : plane.y = normal.y / lenvec;
320 50764 : plane.z = normal.z / lenvec;
321 : // compute the last coefficient of the plane equation
322 50764 : lenvec *= nverts;
323 50764 : plane.w = -dot(refpt, normal) / lenvec;
324 : } else {
325 0 : error = true;
326 : }
327 50764 : }
328 :
329 183399 : Real64 Pt2Plane(Vector const &pt, // Point for determining the distance
330 : PlaneEq const &pleq // Equation of the plane
331 : )
332 : {
333 :
334 : // PURPOSE OF THIS SUBROUTINE:
335 : // This subroutine calculates the distance from a point to the plane (of a surface). Used to determine the reveal of a heat transfer subsurface.
336 :
337 : // REFERENCE:
338 : // Graphic Gems
339 :
340 : Real64 PtDist; // Distance of the point to the plane
341 :
342 183399 : PtDist = (pleq.x * pt.x) + (pleq.y * pt.y) + (pleq.z * pt.z) + pleq.w;
343 :
344 183399 : return PtDist;
345 : }
346 :
347 80974 : void CreateNewellAreaVector(Array1D<Vector> const &VList, int const NSides, Vector &OutNewellAreaVector)
348 : {
349 :
350 : // SUBROUTINE INFORMATION:
351 : // AUTHOR Linda Lawrie
352 : // DATE WRITTEN May 2004
353 :
354 : // PURPOSE OF THIS SUBROUTINE:
355 : // This subroutine creates a "Newell" vector from the vector list for a surface face. Also the Newell Area vector.
356 :
357 : // REFERENCES:
358 : // Collaboration with Bill Carroll, LBNL.
359 :
360 80974 : OutNewellAreaVector = 0.0;
361 :
362 80974 : Vector V1 = VList(2) - VList(1);
363 244650 : for (int Vert = 3; Vert <= NSides; ++Vert) {
364 163676 : Vector V2 = VList(Vert) - VList(1);
365 163676 : OutNewellAreaVector += cross(V1, V2);
366 163676 : V1 = V2;
367 163676 : }
368 :
369 80974 : OutNewellAreaVector /= 2.0;
370 80974 : }
371 :
372 51473 : void CreateNewellSurfaceNormalVector(Array1D<Vector> const &VList, int const NSides, Vector &OutNewellSurfaceNormalVector)
373 : {
374 :
375 : // SUBROUTINE INFORMATION:
376 : // AUTHOR Linda Lawrie
377 : // DATE WRITTEN Jan 2011
378 :
379 : // PURPOSE OF THIS SUBROUTINE:
380 : // This subroutine creates a "Newell" surface normal vector from the vector list for a surface face.
381 :
382 : // REFERENCES:
383 : // September 2010: from OpenGL.org
384 : // Begin Function CalculateSurfaceNormal (Input Polygon) Returns Vector
385 : // Set Vertex Normal to (0, 0, 0)
386 : // Begin Cycle for Index in [0, Polygon.vertexNumber)
387 : // Set Vertex Current to Polygon.verts[Index]
388 : // Set Vertex Next to Polygon.verts[(Index plus 1) mod Polygon.vertexNumber]
389 : // Set Normal.x to Sum of Normal.x and (multiply (Current.y minus Next.y) by (Current.z plus Next.z)
390 : // Set Normal.y to Sum of Normal.y and (multiply (Current.z minus Next.z) by (Current.x plus Next.x)
391 : // Set Normal.z to Sum of Normal.z and (multiply (Current.x minus Next.x) by (Current.y plus Next.y)
392 : // End Cycle
393 : // Returning Normalize(Normal)
394 : // End Function
395 :
396 51473 : OutNewellSurfaceNormalVector = 0.0;
397 51473 : Real64 xvalue = 0.0;
398 51473 : Real64 yvalue = 0.0;
399 51473 : Real64 zvalue = 0.0;
400 :
401 : // IF (NSides > 3) THEN
402 258206 : for (int Side = 1; Side <= NSides; ++Side) {
403 206733 : int curVert = Side;
404 206733 : int nextVert = Side + 1;
405 206733 : if (nextVert > NSides) {
406 51473 : nextVert = 1;
407 : }
408 206733 : xvalue += (VList(curVert).y - VList(nextVert).y) * (VList(curVert).z + VList(nextVert).z);
409 206733 : yvalue += (VList(curVert).z - VList(nextVert).z) * (VList(curVert).x + VList(nextVert).x);
410 206733 : zvalue += (VList(curVert).x - VList(nextVert).x) * (VList(curVert).y + VList(nextVert).y);
411 : }
412 :
413 51473 : OutNewellSurfaceNormalVector.x = xvalue;
414 51473 : OutNewellSurfaceNormalVector.y = yvalue;
415 51473 : OutNewellSurfaceNormalVector.z = zvalue;
416 51473 : OutNewellSurfaceNormalVector = VecNormalize(OutNewellSurfaceNormalVector);
417 51473 : }
418 :
419 6842 : void CompareTwoVectors(Vector const &vector1, // standard vector
420 : Vector const &vector2, // standard vector
421 : bool &areSame, // true if the two vectors are the same within specified tolerance
422 : Real64 const tolerance // specified tolerance
423 : )
424 : {
425 :
426 : // SUBROUTINE INFORMATION:
427 : // AUTHOR Linda Lawrie
428 : // DATE WRITTEN February 2012
429 :
430 : // PURPOSE OF THIS SUBROUTINE:
431 : // This routine will provide the ability to compare two vectors (e.g. surface normals) to be the same within a specified tolerance.
432 :
433 : // METHODOLOGY EMPLOYED:
434 : // compare each element (x,y,z)
435 :
436 6842 : areSame = true;
437 6842 : if (std::abs(vector1.x - vector2.x) > tolerance) {
438 0 : areSame = false;
439 : }
440 6842 : if (std::abs(vector1.y - vector2.y) > tolerance) {
441 1 : areSame = false;
442 : }
443 6842 : if (std::abs(vector1.z - vector2.z) > tolerance) {
444 0 : areSame = false;
445 : }
446 6842 : }
447 :
448 43908 : void CalcCoPlanarNess(Array1D<Vector> &Surf, int const NSides, bool &IsCoPlanar, Real64 &MaxDist, int &ErrorVertex)
449 : {
450 :
451 : // SUBROUTINE INFORMATION:
452 : // AUTHOR Linda Lawrie
453 : // DATE WRITTEN June 2004
454 :
455 : // PURPOSE OF THIS SUBROUTINE:
456 : // This subroutine provides the calculation to determine if the surface is planar or not.
457 :
458 : // REFERENCES:
459 : // Eric W. Weisstein. "Coplanar." From MathWorld--A Wolfram Web Resource.
460 : // http://mathworld.wolfram.com/Coplanar.html
461 :
462 : // Argument array dimensioning
463 43908 : EP_SIZE_CHECK(Surf, NSides);
464 :
465 : bool plerror;
466 43908 : PlaneEq NewellPlane;
467 :
468 43908 : IsCoPlanar = true;
469 43908 : MaxDist = 0.0;
470 43908 : ErrorVertex = 0;
471 :
472 : // Use first three to determine plane
473 43908 : PlaneEquation(Surf, NSides, NewellPlane, plerror);
474 :
475 220403 : for (int vert = 1; vert <= NSides; ++vert) {
476 176495 : Real64 dist = Pt2Plane(Surf(vert), NewellPlane);
477 176495 : if (std::abs(dist) > MaxDist) {
478 10845 : MaxDist = std::abs(dist);
479 10845 : ErrorVertex = vert;
480 : }
481 : }
482 :
483 43908 : if (std::abs(MaxDist) > Constant::SmallDistance) {
484 0 : IsCoPlanar = false;
485 : }
486 43908 : }
487 :
488 : std::vector<int>
489 16 : PointsInPlane(Array1D<Vector> &BaseSurf, int const BaseSides, Array1D<Vector> const &QuerySurf, int const QuerySides, bool &ErrorFound)
490 : {
491 16 : std::vector<int> pointIndices;
492 :
493 16 : PlaneEq NewellPlane;
494 16 : PlaneEquation(BaseSurf, BaseSides, NewellPlane, ErrorFound);
495 :
496 80 : for (int vert = 1; vert <= QuerySides; ++vert) {
497 64 : Real64 dist = Pt2Plane(QuerySurf(vert), NewellPlane);
498 64 : if (std::abs(dist) < Constant::SmallDistance) { // point on query surface is co-planar with base surface
499 32 : pointIndices.push_back(vert);
500 : }
501 : }
502 32 : return pointIndices;
503 0 : }
504 :
505 5174 : Real64 CalcPolyhedronVolume(EnergyPlusData const &state, Polyhedron const &Poly)
506 : {
507 :
508 : // SUBROUTINE INFORMATION:
509 : // AUTHOR Linda Lawrie
510 : // DATE WRITTEN June 2004
511 :
512 : // PURPOSE OF THIS SUBROUTINE:
513 : // This subroutine provides the volume calculation for a polyhedron (i.e. Zone).
514 :
515 : // REFERENCES:
516 : // Conversations with Bill Carroll, LBNL.
517 :
518 : // Object Data
519 5174 : Vector p3FaceOrigin;
520 :
521 5174 : Real64 Volume = 0.0;
522 :
523 41129 : for (int NFace = 1; NFace <= Poly.NumSurfaceFaces; ++NFace) {
524 35955 : p3FaceOrigin = Poly.SurfaceFace(NFace).FacePoints(2);
525 35955 : Real64 PyramidVolume = dot(Poly.SurfaceFace(NFace).NewellAreaVector, (p3FaceOrigin - state.dataVectors->p0));
526 35955 : Volume += PyramidVolume / 3.0;
527 : }
528 5174 : return Volume;
529 5174 : }
530 :
531 : } // namespace EnergyPlus::Vectors
|
