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