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 4415 : 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 4415 : EP_SIZE_CHECK(p, n);
117 :
118 4415 : Vector edge0 = p[1] - p[0];
119 4415 : Vector edge1 = p[2] - p[0];
120 :
121 4415 : Vector edgex = cross(edge0, edge1);
122 4415 : Vector nor = VecNormalize(edgex);
123 :
124 : // Initialize csum
125 4415 : Vector csum;
126 4415 : csum = 0.0;
127 :
128 13246 : for (int i = 0; i <= n - 2; ++i) {
129 8831 : csum += cross(p[i], p[i + 1]);
130 : }
131 4415 : csum += cross(p[n - 1], p[0]);
132 :
133 4415 : Real64 areap = 0.5 * std::abs(dot(nor, csum));
134 :
135 4415 : return areap;
136 4415 : }
137 :
138 29025 : 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 29025 : Real64 vecsqlen = (vec.x * vec.x + vec.y * vec.y + vec.z * vec.z);
148 :
149 29025 : return vecsqlen;
150 : }
151 :
152 29025 : 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 29025 : Real64 veclen = std::sqrt(VecSquaredLength(vec));
162 :
163 29025 : 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 12355 : 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 12355 : Vector VecNormalize;
195 :
196 12355 : Real64 veclen = VecLength(vec);
197 12355 : if (veclen != 0.0) {
198 12352 : VecNormalize.x = vec.x / veclen;
199 12352 : VecNormalize.y = vec.y / veclen;
200 12352 : VecNormalize.z = vec.z / veclen;
201 : } else {
202 3 : VecNormalize.x = 0.0;
203 3 : VecNormalize.y = 0.0;
204 3 : VecNormalize.z = 0.0;
205 : }
206 :
207 12355 : return VecNormalize;
208 : }
209 :
210 1 : 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 1 : vec.x = nint64(vec.x * roundto) / roundto;
220 1 : vec.y = nint64(vec.y * roundto) / roundto;
221 1 : vec.z = nint64(vec.z * roundto) / roundto;
222 1 : }
223 :
224 2353 : 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 2353 : lcsx = VecNormalize(Surf(3) - Surf(2));
241 2353 : lcsz = NewellSurfaceNormalVector;
242 2353 : lcsy = cross(lcsz, lcsx);
243 :
244 2353 : Real64 costheta = dot(lcsz, ZUnit);
245 :
246 : // if ( fabs(costheta) < 1.0d0) { // normal cases
247 2353 : Real64 constexpr epsilon = 1.12e-16;
248 2353 : Real64 rotang_0 = 0.0;
249 2353 : 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 1544 : Vector x2 = cross(ZUnit, lcsz);
253 1544 : rotang_0 = std::atan2(dot(x2, YUnit), dot(x2, XUnit));
254 1544 : } else {
255 : // azimuth
256 809 : rotang_0 = std::atan2(dot(lcsx, YUnit), dot(lcsx, XUnit));
257 : }
258 :
259 2353 : Real64 tlt = std::acos(NewellSurfaceNormalVector.z);
260 2353 : tlt /= Constant::DegToRad;
261 :
262 2353 : Real64 az = rotang_0;
263 :
264 2353 : az /= Constant::DegToRad;
265 2353 : az = mod(450.0 - az, 360.0);
266 2353 : az += 90.0;
267 2353 : if (az < 0.0) az += 360.0;
268 2353 : az = mod(az, 360.0);
269 :
270 : // Clean up angle precision
271 2353 : if (std::abs(az - 360.0) < Constant::OneThousandth) { // Bring small angles to zero
272 1 : az = 0.0;
273 2352 : } else if (std::abs(az - 180.0) <
274 : Constant::OneMillionth) { // Bring angles near 180 to 180 //Autodesk Added to clean up debug--release discrepancies
275 870 : az = 180.0;
276 : }
277 2353 : if (std::abs(tlt - 180.0) < Constant::OneMillionth) { // Bring angles near 180 to 180 //Autodesk Added to clean up debug--release discrepancies
278 456 : tlt = 180.0;
279 : }
280 :
281 2353 : Azimuth = az;
282 2353 : Tilt = tlt;
283 2353 : }
284 :
285 1789 : void PlaneEquation(Array1D<Vector> &verts, // Structure of the surface
286 : int const nverts, // Number of vertices in the surface
287 : PlaneEq &plane, // Equation of plane from inputs
288 : bool &error // returns true for degenerate surface
289 : )
290 : {
291 :
292 : // PURPOSE OF THIS SUBROUTINE:
293 : // This subroutine calculates the plane equation for a given surface (which should be planar).
294 :
295 : // REFERENCE:
296 : // Graphic Gems
297 :
298 : // Argument array dimensioning
299 1789 : EP_SIZE_CHECK(verts, nverts);
300 :
301 1789 : Vector normal = Vector(0.0, 0.0, 0.0);
302 1789 : Vector refpt = Vector(0.0, 0.0, 0.0);
303 8949 : for (int i = 0; i <= nverts - 1; ++i) {
304 7160 : Vector const &u(verts[i]);
305 7160 : Vector const &v(i < nverts - 1 ? verts[i + 1] : verts[0]);
306 7160 : normal.x += (u.y - v.y) * (u.z + v.z);
307 7160 : normal.y += (u.z - v.z) * (u.x + v.x);
308 7160 : normal.z += (u.x - v.x) * (u.y + v.y);
309 7160 : refpt += u;
310 : }
311 : // normalize the polygon normal to obtain the first
312 : // three coefficients of the plane equation
313 1789 : Real64 lenvec = VecLength(normal);
314 1789 : error = false;
315 1789 : if (lenvec != 0.0) { // should this be >0
316 1789 : plane.x = normal.x / lenvec;
317 1789 : plane.y = normal.y / lenvec;
318 1789 : plane.z = normal.z / lenvec;
319 : // compute the last coefficient of the plane equation
320 1789 : lenvec *= nverts;
321 1789 : plane.w = -dot(refpt, normal) / lenvec;
322 : } else {
323 0 : error = true;
324 : }
325 1789 : }
326 :
327 6767 : Real64 Pt2Plane(Vector const &pt, // Point for determining the distance
328 : PlaneEq const &pleq // Equation of the plane
329 : )
330 : {
331 :
332 : // PURPOSE OF THIS SUBROUTINE:
333 : // This subroutine calculates the distance from a point to the plane (of a surface). Used to determine the reveal of a heat transfer subsurface.
334 :
335 : // REFERENCE:
336 : // Graphic Gems
337 :
338 : Real64 PtDist; // Distance of the point to the plane
339 :
340 6767 : PtDist = (pleq.x * pt.x) + (pleq.y * pt.y) + (pleq.z * pt.z) + pleq.w;
341 :
342 6767 : return PtDist;
343 : }
344 :
345 3777 : void CreateNewellAreaVector(Array1D<Vector> const &VList, int const NSides, Vector &OutNewellAreaVector)
346 : {
347 :
348 : // SUBROUTINE INFORMATION:
349 : // AUTHOR Linda Lawrie
350 : // DATE WRITTEN May 2004
351 :
352 : // PURPOSE OF THIS SUBROUTINE:
353 : // This subroutine creates a "Newell" vector from the vector list for a surface face. Also the Newell Area vector.
354 :
355 : // REFERENCES:
356 : // Collaboration with Bill Carroll, LBNL.
357 :
358 3777 : OutNewellAreaVector = 0.0;
359 :
360 3777 : Vector V1 = VList(2) - VList(1);
361 11330 : for (int Vert = 3; Vert <= NSides; ++Vert) {
362 7553 : Vector V2 = VList(Vert) - VList(1);
363 7553 : OutNewellAreaVector += cross(V1, V2);
364 7553 : V1 = V2;
365 7553 : }
366 :
367 3777 : OutNewellAreaVector /= 2.0;
368 3777 : }
369 :
370 2353 : void CreateNewellSurfaceNormalVector(Array1D<Vector> const &VList, int const NSides, Vector &OutNewellSurfaceNormalVector)
371 : {
372 :
373 : // SUBROUTINE INFORMATION:
374 : // AUTHOR Linda Lawrie
375 : // DATE WRITTEN Jan 2011
376 :
377 : // PURPOSE OF THIS SUBROUTINE:
378 : // This subroutine creates a "Newell" surface normal vector from the vector list for a surface face.
379 :
380 : // REFERENCES:
381 : // September 2010: from OpenGL.org
382 : // Begin Function CalculateSurfaceNormal (Input Polygon) Returns Vector
383 : // Set Vertex Normal to (0, 0, 0)
384 : // Begin Cycle for Index in [0, Polygon.vertexNumber)
385 : // Set Vertex Current to Polygon.verts[Index]
386 : // Set Vertex Next to Polygon.verts[(Index plus 1) mod Polygon.vertexNumber]
387 : // Set Normal.x to Sum of Normal.x and (multiply (Current.y minus Next.y) by (Current.z plus Next.z)
388 : // Set Normal.y to Sum of Normal.y and (multiply (Current.z minus Next.z) by (Current.x plus Next.x)
389 : // Set Normal.z to Sum of Normal.z and (multiply (Current.x minus Next.x) by (Current.y plus Next.y)
390 : // End Cycle
391 : // Returning Normalize(Normal)
392 : // End Function
393 :
394 2353 : OutNewellSurfaceNormalVector = 0.0;
395 2353 : Real64 xvalue = 0.0;
396 2353 : Real64 yvalue = 0.0;
397 2353 : Real64 zvalue = 0.0;
398 :
399 : // IF (NSides > 3) THEN
400 11762 : for (int Side = 1; Side <= NSides; ++Side) {
401 9409 : int curVert = Side;
402 9409 : int nextVert = Side + 1;
403 9409 : if (nextVert > NSides) nextVert = 1;
404 9409 : xvalue += (VList(curVert).y - VList(nextVert).y) * (VList(curVert).z + VList(nextVert).z);
405 9409 : yvalue += (VList(curVert).z - VList(nextVert).z) * (VList(curVert).x + VList(nextVert).x);
406 9409 : zvalue += (VList(curVert).x - VList(nextVert).x) * (VList(curVert).y + VList(nextVert).y);
407 : }
408 :
409 2353 : OutNewellSurfaceNormalVector.x = xvalue;
410 2353 : OutNewellSurfaceNormalVector.y = yvalue;
411 2353 : OutNewellSurfaceNormalVector.z = zvalue;
412 2353 : OutNewellSurfaceNormalVector = VecNormalize(OutNewellSurfaceNormalVector);
413 2353 : }
414 :
415 238 : void CompareTwoVectors(Vector const &vector1, // standard vector
416 : Vector const &vector2, // standard vector
417 : bool &areSame, // true if the two vectors are the same within specified tolerance
418 : Real64 const tolerance // specified tolerance
419 : )
420 : {
421 :
422 : // SUBROUTINE INFORMATION:
423 : // AUTHOR Linda Lawrie
424 : // DATE WRITTEN February 2012
425 :
426 : // PURPOSE OF THIS SUBROUTINE:
427 : // This routine will provide the ability to compare two vectors (e.g. surface normals) to be the same within a specified tolerance.
428 :
429 : // METHODOLOGY EMPLOYED:
430 : // compare each element (x,y,z)
431 :
432 238 : areSame = true;
433 238 : if (std::abs(vector1.x - vector2.x) > tolerance) areSame = false;
434 238 : if (std::abs(vector1.y - vector2.y) > tolerance) areSame = false;
435 238 : if (std::abs(vector1.z - vector2.z) > tolerance) areSame = false;
436 238 : }
437 :
438 1657 : void CalcCoPlanarNess(Array1D<Vector> &Surf, int const NSides, bool &IsCoPlanar, Real64 &MaxDist, int &ErrorVertex)
439 : {
440 :
441 : // SUBROUTINE INFORMATION:
442 : // AUTHOR Linda Lawrie
443 : // DATE WRITTEN June 2004
444 :
445 : // PURPOSE OF THIS SUBROUTINE:
446 : // This subroutine provides the calculation to determine if the surface is planar or not.
447 :
448 : // REFERENCES:
449 : // Eric W. Weisstein. "Coplanar." From MathWorld--A Wolfram Web Resource.
450 : // http://mathworld.wolfram.com/Coplanar.html
451 :
452 : // Argument array dimensioning
453 1657 : EP_SIZE_CHECK(Surf, NSides);
454 :
455 : bool plerror;
456 1657 : PlaneEq NewellPlane;
457 :
458 1657 : IsCoPlanar = true;
459 1657 : MaxDist = 0.0;
460 1657 : ErrorVertex = 0;
461 :
462 : // Use first three to determine plane
463 1657 : PlaneEquation(Surf, NSides, NewellPlane, plerror);
464 :
465 8289 : for (int vert = 1; vert <= NSides; ++vert) {
466 6632 : Real64 dist = Pt2Plane(Surf(vert), NewellPlane);
467 6632 : if (std::abs(dist) > MaxDist) {
468 308 : MaxDist = std::abs(dist);
469 308 : ErrorVertex = vert;
470 : }
471 : }
472 :
473 1657 : if (std::abs(MaxDist) > Constant::SmallDistance) IsCoPlanar = false;
474 1657 : }
475 :
476 : std::vector<int>
477 1 : PointsInPlane(Array1D<Vector> &BaseSurf, int const BaseSides, Array1D<Vector> const &QuerySurf, int const QuerySides, bool &ErrorFound)
478 : {
479 1 : std::vector<int> pointIndices;
480 :
481 1 : PlaneEq NewellPlane;
482 1 : PlaneEquation(BaseSurf, BaseSides, NewellPlane, ErrorFound);
483 :
484 5 : for (int vert = 1; vert <= QuerySides; ++vert) {
485 4 : Real64 dist = Pt2Plane(QuerySurf(vert), NewellPlane);
486 4 : if (std::abs(dist) < Constant::SmallDistance) { // point on query surface is co-planar with base surface
487 2 : pointIndices.push_back(vert);
488 : }
489 : }
490 2 : return pointIndices;
491 0 : }
492 :
493 185 : Real64 CalcPolyhedronVolume(EnergyPlusData const &state, Polyhedron const &Poly)
494 : {
495 :
496 : // SUBROUTINE INFORMATION:
497 : // AUTHOR Linda Lawrie
498 : // DATE WRITTEN June 2004
499 :
500 : // PURPOSE OF THIS SUBROUTINE:
501 : // This subroutine provides the volume calculation for a polyhedron (i.e. Zone).
502 :
503 : // REFERENCES:
504 : // Conversations with Bill Carroll, LBNL.
505 :
506 : // Object Data
507 185 : Vector p3FaceOrigin;
508 :
509 185 : Real64 Volume = 0.0;
510 :
511 1303 : for (int NFace = 1; NFace <= Poly.NumSurfaceFaces; ++NFace) {
512 1118 : p3FaceOrigin = Poly.SurfaceFace(NFace).FacePoints(2);
513 1118 : Real64 PyramidVolume = dot(Poly.SurfaceFace(NFace).NewellAreaVector, (p3FaceOrigin - state.dataVectors->p0));
514 1118 : Volume += PyramidVolume / 3.0;
515 : }
516 185 : return Volume;
517 185 : }
518 :
519 : } // namespace EnergyPlus::Vectors
|
