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 : // 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 : // MODIFIED na
66 : // RE-ENGINEERED na
67 :
68 : // PURPOSE OF THIS MODULE:
69 : // This module uses a global "vector" data structure and defines
70 : // operations (using the F90 "operator" features) and other
71 : // manipulations used with Vectors.
72 :
73 : // Vectors should allow for more modern interpretations of the
74 : // calculations used in the shadowing and other surface manipulations.
75 :
76 : // METHODOLOGY EMPLOYED:
77 : // Uses the "Vector" derived type variables to allow for more readable
78 : // calculations.
79 :
80 : // REFERENCES: Original idea from F90 for Scientists and
81 : // Engineers, S. J. Chapman, 1996.
82 : // The module defines 8 operations which can be performed on vectors:
83 : // Operation Operator
84 : // ========= ========
85 : // 1. Creation from a real array =
86 : // 2. Conversion to real array =
87 : // 3. Vector addition +
88 : // 4. Vector subtraction -
89 : // 5. Vector-scalar multiplication (4 cases) *
90 : // 6. Vector-scalar division (2 cases) /
91 : // 7. Dot product .dot.
92 : // 8. Cross product *
93 : // 9. 2d dot product .twoddot.
94 : // 10. 2d Cross product .twodcross.
95 : // It contains a total of 12 procedures to implement those
96 : // operations: array_to_vector, vector_to_array, vector_add,
97 : // vector_subtract, vector_times_real, real_times_vector,
98 : // vector_times_int, int_times_vector, vector_div_real,
99 : // vector_div_int, dot_product, and cross_product.
100 :
101 : // OTHER NOTES: none
102 :
103 : // Using/Aliasing
104 : using namespace DataVectorTypes;
105 :
106 : // MODULE PARAMETER DEFINITIONS
107 :
108 : // Object Data
109 771 : Vector const XUnit(1.0, 0.0, 0.0);
110 771 : Vector const YUnit(0.0, 1.0, 0.0);
111 771 : Vector const ZUnit(0.0, 0.0, 1.0);
112 :
113 : // DERIVED TYPE DEFINITIONS
114 : // na
115 :
116 : // MODULE VARIABLE DECLARATIONS:
117 : // na
118 :
119 : // SUBROUTINE SPECIFICATIONS FOR MODULE <module_name>
120 :
121 : // Functions
122 :
123 82296 : Real64 AreaPolygon(int const n, Array1D<Vector> &p)
124 : {
125 :
126 : // PURPOSE OF THIS SUBROUTINE:
127 : // This subroutine calculates the area of a polygon defined by the
128 : // input vectors.
129 :
130 : // REFERENCE:
131 : // Graphic Gems.
132 :
133 : // Return value
134 : Real64 areap;
135 :
136 : // Argument array dimensioning
137 82296 : EP_SIZE_CHECK(p, n);
138 :
139 : // Locals
140 : // SUBROUTINE ARGUMENT DEFINITIONS:
141 :
142 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
143 : int i;
144 :
145 : // Object Data
146 164592 : Vector edge0;
147 164592 : Vector edge1;
148 164592 : Vector nor;
149 164592 : Vector edgex;
150 164592 : Vector csum;
151 :
152 82296 : edge0 = p[1] - p[0];
153 82296 : edge1 = p[2] - p[0];
154 :
155 82296 : edgex = cross(edge0, edge1);
156 82296 : nor = VecNormalize(edgex);
157 :
158 : // Initialize csum
159 82296 : csum = 0.0;
160 :
161 246888 : for (i = 0; i <= n - 2; ++i) {
162 164592 : csum += cross(p[i], p[i + 1]);
163 : }
164 82296 : csum += cross(p[n - 1], p[0]);
165 :
166 82296 : areap = 0.5 * std::abs(dot(nor, csum));
167 :
168 164592 : return areap;
169 : }
170 :
171 722126 : Real64 VecSquaredLength(Vector const &vec)
172 : {
173 :
174 : // PURPOSE OF THIS SUBROUTINE:
175 : // This subroutine calculates the squared length of the input vector.
176 :
177 : // REFERENCE:
178 : // Graphic Gems.
179 :
180 : // Return value
181 : Real64 vecsqlen;
182 :
183 : // Locals
184 : // SUBROUTINE ARGUMENT DEFINITIONS:
185 :
186 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
187 : // na
188 :
189 722126 : vecsqlen = (vec.x * vec.x + vec.y * vec.y + vec.z * vec.z);
190 :
191 722126 : return vecsqlen;
192 : }
193 :
194 722126 : Real64 VecLength(Vector const &vec)
195 : {
196 :
197 : // PURPOSE OF THIS SUBROUTINE:
198 : // This subroutine calculates the length of the input vector.
199 :
200 : // REFERENCE:
201 : // Graphic Gems.
202 :
203 : // Return value
204 : Real64 veclen;
205 :
206 : // Locals
207 : // SUBROUTINE ARGUMENT DEFINITIONS:
208 :
209 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
210 : // na
211 :
212 722126 : veclen = std::sqrt(VecSquaredLength(vec));
213 :
214 722126 : return veclen;
215 : }
216 :
217 0 : Vector VecNegate(Vector const &vec)
218 : {
219 :
220 : // PURPOSE OF THIS SUBROUTINE:
221 : // This subroutine negates the input vector and returns that vector.
222 :
223 : // REFERENCE:
224 : // Graphic Gems.
225 :
226 : // Return value
227 0 : Vector VecNegate;
228 :
229 : // Locals
230 : // SUBROUTINE ARGUMENT DEFINITIONS:
231 :
232 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
233 : // na
234 :
235 0 : VecNegate.x = -vec.x;
236 0 : VecNegate.y = -vec.y;
237 0 : VecNegate.z = -vec.z;
238 :
239 0 : return VecNegate;
240 : }
241 :
242 354576 : Vector VecNormalize(Vector const &vec)
243 : {
244 :
245 : // PURPOSE OF THIS SUBROUTINE:
246 : // This subroutine normalizes the input vector and returns the normalized
247 : // vector
248 :
249 : // REFERENCE:
250 : // Graphic Gems.
251 :
252 : // Return value
253 354576 : Vector VecNormalize;
254 :
255 : // Locals
256 : // SUBROUTINE ARGUMENT DEFINITIONS:
257 :
258 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
259 : Real64 veclen;
260 :
261 354576 : veclen = VecLength(vec);
262 354576 : if (veclen != 0.0) {
263 354576 : VecNormalize.x = vec.x / veclen;
264 354576 : VecNormalize.y = vec.y / veclen;
265 354576 : VecNormalize.z = vec.z / veclen;
266 : } else {
267 0 : VecNormalize.x = 0.0;
268 0 : VecNormalize.y = 0.0;
269 0 : VecNormalize.z = 0.0;
270 : }
271 :
272 354576 : return VecNormalize;
273 : }
274 :
275 0 : void VecRound(Vector &vec, Real64 const roundto)
276 : {
277 :
278 : // PURPOSE OF THIS SUBROUTINE:
279 : // This subroutine rounds the input vector to a specified "rounding" value.
280 :
281 : // REFERENCE:
282 : // Graphic Gems.
283 :
284 : // Locals
285 : // SUBROUTINE ARGUMENT DEFINITIONS:
286 :
287 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
288 : // na
289 :
290 0 : vec.x = nint64(vec.x * roundto) / roundto;
291 0 : vec.y = nint64(vec.y * roundto) / roundto;
292 0 : vec.z = nint64(vec.z * roundto) / roundto;
293 0 : }
294 :
295 48198 : void DetermineAzimuthAndTilt(Array1D<Vector> const &Surf, // Surface Definition
296 : [[maybe_unused]] int const NSides, // Number of sides to surface
297 : Real64 &Azimuth, // Outward Normal Azimuth Angle
298 : Real64 &Tilt, // Tilt angle of surface
299 : Vector &lcsx,
300 : Vector &lcsy,
301 : Vector &lcsz,
302 : [[maybe_unused]] Real64 const surfaceArea,
303 : Vector const &NewellSurfaceNormalVector)
304 : {
305 :
306 : // PURPOSE OF THIS SUBROUTINE:
307 : // This subroutine determines the Azimuth (outward normal) angle,
308 : // Tilt angle of a given surface defined by the set of input vectors.
309 :
310 : // REFERENCE:
311 : // Discussions and examples from Bill Carroll, LBNL.
312 :
313 : // Argument array dimensioning
314 :
315 : // Locals
316 : // SUBROUTINE ARGUMENT DEFINITIONS:
317 :
318 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
319 : // TYPE(Vector) :: x3,y3,z3,v12
320 : // TYPE(Vector) :: y2
321 : Real64 costheta;
322 : Real64 rotang_0;
323 : // REAL(r64) rotang_2
324 :
325 : Real64 az;
326 : // REAL(r64) azm
327 : Real64 tlt;
328 : // REAL(r64) newtlt
329 : // REAL(r64) roundval
330 : // REAL(r64) xcomp
331 : // REAL(r64) ycomp
332 : // REAL(r64) zcomp
333 : // REAL(r64) proj
334 : // integer :: scount
335 : // integer :: nvert1
336 : // REAL(r64) :: tltcos
337 :
338 : // Object Data
339 96396 : Vector x2;
340 96396 : Vector x3a;
341 96396 : Vector v12a;
342 96396 : Vector v0;
343 96396 : Vector v1;
344 96396 : Vector v2;
345 96396 : Vector cs3_2;
346 96396 : Vector cs3_0;
347 96396 : Vector cs3_1;
348 96396 : Vector z3;
349 :
350 : //!! x3=VecNormalize(Surf(2)-Surf(1))
351 : //!! v12=Surf(3)-Surf(2)
352 :
353 : //!! z3=VecNormalize(x3*v12)
354 : //!! y3=z3*x3
355 : //!! roundval=10000.d0
356 : //!! CALL VecRound(x3,roundval)
357 : //!! CALL VecRound(y3,roundval)
358 : //!! CALL VecRound(z3,roundval)
359 :
360 : //!!! Direction cosines, local coordinates.
361 : //!!! write(OUTPUT,*) 'lcs:'
362 : //!!! write(OUTPUT,*) 'x=',x3
363 : //!!! write(OUTPUT,*) 'y=',y3
364 : //!!! write(OUTPUT,*) 'z=',z3
365 48198 : x3a = VecNormalize(Surf(3) - Surf(2));
366 48198 : v12a = Surf(1) - Surf(2);
367 :
368 : //!! lcsx=x3a
369 : //!! lcsz=VecNormalize(x3a*v12a)
370 : //!! lcsy=lcsz*x3a
371 :
372 48198 : lcsx = x3a;
373 48198 : lcsz = NewellSurfaceNormalVector;
374 48198 : lcsy = cross(lcsz, x3a);
375 :
376 : //!!
377 :
378 : // Vec3d v0(p1 - p0); ! BGL has different conventions...p0=surf(2), etc
379 48198 : v0 = Surf(3) - Surf(2);
380 : // Vec3d v1(p2 - p0);
381 48198 : v1 = Surf(1) - Surf(2);
382 :
383 : // Vec3d v2 = cross(v0,v1);
384 48198 : v2 = cross(v0, v1);
385 : // cs3[2] = norm(v2); // z
386 48198 : cs3_2 = VecNormalize(v2);
387 : // cs3[0] = norm(v0); // x
388 48198 : cs3_0 = VecNormalize(v0);
389 : // cs3[1] = cross(cs3[2],cs3[0]); // y
390 48198 : cs3_1 = cross(cs3_2, cs3_0);
391 : // Vec3d z3 = cs3[2];
392 48198 : z3 = cs3_2;
393 : // double costheta = dot(z3,Ref_CS[2]);
394 48198 : costheta = dot(z3, ZUnit);
395 :
396 : // if ( fabs(costheta) < 1.0d0) { // normal cases
397 48198 : if (std::abs(costheta) < 1.0 - 1.12e-16) { // Autodesk Added - 1.12e-16 to treat 1 bit from 1.0 as 1.0 to correct different behavior seen in
398 : // release vs debug build due to slight precision differences: May want larger epsilon here
399 : // // azimuth
400 : // Vec3d x2 = cross(Ref_CS[2],z3); // order is important; x2 = x1
401 : // RotAng[0] = ATAN2(dot(x2,Ref_CS[1]),dot(x2,Ref_CS[0]));
402 33358 : x2 = cross(ZUnit, z3);
403 33358 : rotang_0 = std::atan2(dot(x2, YUnit), dot(x2, XUnit));
404 :
405 : } else {
406 :
407 : // }
408 : // else { // special cases: tilt angle theta = 0, PI
409 : // // azimuth
410 : // RotAng[0] = ATAN2(dot(cs3[0],Ref_CS[1]),dot(cs3[0],Ref_CS[0]) );
411 14840 : rotang_0 = std::atan2(dot(cs3_0, YUnit), dot(cs3_0, XUnit));
412 : }
413 :
414 48198 : tlt = std::acos(NewellSurfaceNormalVector.z);
415 48198 : tlt /= DataGlobalConstants::DegToRadians;
416 :
417 48198 : az = rotang_0;
418 :
419 48198 : az /= DataGlobalConstants::DegToRadians;
420 48198 : az = mod(450.0 - az, 360.0);
421 48198 : az += 90.0;
422 48198 : if (az < 0.0) az += 360.0;
423 48198 : az = mod(az, 360.0);
424 :
425 : // Clean up angle precision
426 48198 : if (std::abs(az - 360.0) < 1.0e-3) { // Bring small angles to zero
427 1 : az = 0.0;
428 48197 : } else if (std::abs(az - 180.0) < 1.0e-6) { // Bring angles near 180 to 180 //Autodesk Added to clean up debug--release discrepancies
429 12353 : az = 180.0;
430 : }
431 48198 : if (std::abs(tlt - 180.0) < 1.0e-6) { // Bring angles near 180 to 180 //Autodesk Added to clean up debug--release discrepancies
432 8527 : tlt = 180.0;
433 : }
434 :
435 48198 : Azimuth = az;
436 48198 : Tilt = tlt;
437 48198 : }
438 :
439 47162 : void PlaneEquation(Array1D<Vector> &verts, // Structure of the surface
440 : int const nverts, // Number of vertices in the surface
441 : PlaneEq &plane, // Equation of plane from inputs
442 : bool &error // returns true for degenerate surface
443 : )
444 : {
445 :
446 : // PURPOSE OF THIS SUBROUTINE:
447 : // This subroutine calculates the plane equation for a given
448 : // surface (which should be planar).
449 :
450 : // REFERENCE:
451 : // Graphic Gems
452 :
453 : // Argument array dimensioning
454 47162 : EP_SIZE_CHECK(verts, nverts);
455 :
456 : // Locals
457 : // SUBROUTINE ARGUMENT DEFINITIONS:
458 :
459 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
460 : int i;
461 : Real64 lenvec;
462 :
463 : // Object Data
464 94324 : Vector normal;
465 94324 : Vector refpt;
466 :
467 : // - - - begin - - -
468 47162 : normal = Vector(0.0, 0.0, 0.0);
469 47162 : refpt = Vector(0.0, 0.0, 0.0);
470 236569 : for (i = 0; i <= nverts - 1; ++i) {
471 189407 : Vector const &u(verts[i]);
472 189407 : Vector const &v(i < nverts - 1 ? verts[i + 1] : verts[0]);
473 189407 : normal.x += (u.y - v.y) * (u.z + v.z);
474 189407 : normal.y += (u.z - v.z) * (u.x + v.x);
475 189407 : normal.z += (u.x - v.x) * (u.y + v.y);
476 189407 : refpt += u;
477 : }
478 : // normalize the polygon normal to obtain the first
479 : // three coefficients of the plane equation
480 47162 : lenvec = VecLength(normal);
481 47162 : error = false;
482 47162 : if (lenvec != 0.0) { // should this be >0
483 47162 : plane.x = normal.x / lenvec;
484 47162 : plane.y = normal.y / lenvec;
485 47162 : plane.z = normal.z / lenvec;
486 : // compute the last coefficient of the plane equation
487 47162 : lenvec *= nverts;
488 47162 : plane.w = -dot(refpt, normal) / lenvec;
489 : } else {
490 0 : error = true;
491 : }
492 47162 : }
493 :
494 170023 : Real64 Pt2Plane(Vector const &pt, // Point for determining the distance
495 : PlaneEq const &pleq // Equation of the plane
496 : )
497 : {
498 :
499 : // PURPOSE OF THIS SUBROUTINE:
500 : // This subroutine calculates the distance from a point
501 : // to the plane (of a surface). Used to determine the reveal
502 : // of a heat transfer subsurface.
503 :
504 : // REFERENCE:
505 : // Graphic Gems
506 :
507 : // Return value
508 : Real64 PtDist; // Distance of the point to the plane
509 :
510 : // Locals
511 : // SUBROUTINE ARGUMENT DEFINITIONS:
512 :
513 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
514 : // na
515 :
516 170023 : PtDist = (pleq.x * pt.x) + (pleq.y * pt.y) + (pleq.z * pt.z) + pleq.w;
517 :
518 170023 : return PtDist;
519 : }
520 :
521 74972 : void CreateNewellAreaVector(Array1D<Vector> const &VList, int const NSides, Vector &OutNewellAreaVector)
522 : {
523 :
524 : // SUBROUTINE INFORMATION:
525 : // AUTHOR Linda Lawrie
526 : // DATE WRITTEN May 2004
527 : // MODIFIED na
528 : // RE-ENGINEERED na
529 :
530 : // PURPOSE OF THIS SUBROUTINE:
531 : // This subroutine creates a "Newell" vector from the vector list for a surface
532 : // face. Also the Newell Area vector.
533 :
534 : // METHODOLOGY EMPLOYED:
535 : // na
536 :
537 : // REFERENCES:
538 : // Collaboration with Bill Carroll, LBNL.
539 :
540 : // USE STATEMENTS:
541 : // na
542 :
543 : // Argument array dimensioning
544 :
545 : // Locals
546 : // SUBROUTINE ARGUMENT DEFINITIONS:
547 :
548 : // SUBROUTINE PARAMETER DEFINITIONS:
549 : // na
550 :
551 : // INTERFACE BLOCK SPECIFICATIONS
552 : // na
553 :
554 : // DERIVED TYPE DEFINITIONS
555 : // na
556 :
557 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
558 :
559 : // Object Data
560 149944 : Vector V1;
561 149944 : Vector V2;
562 :
563 74972 : OutNewellAreaVector = 0.0;
564 :
565 74972 : V1 = VList(2) - VList(1);
566 226308 : for (int Vert = 3; Vert <= NSides; ++Vert) {
567 151336 : V2 = VList(Vert) - VList(1);
568 151336 : OutNewellAreaVector += cross(V1, V2);
569 151336 : V1 = V2;
570 : }
571 : // do vert=1,nsides
572 : // write(outputfiledebug,*) vlist(vert)
573 : // enddo
574 :
575 74972 : OutNewellAreaVector /= 2.0;
576 74972 : }
577 :
578 48198 : void CreateNewellSurfaceNormalVector(Array1D<Vector> const &VList, int const NSides, Vector &OutNewellSurfaceNormalVector)
579 : {
580 :
581 : // SUBROUTINE INFORMATION:
582 : // AUTHOR Linda Lawrie
583 : // DATE WRITTEN Jan 2011
584 : // MODIFIED na
585 : // RE-ENGINEERED na
586 :
587 : // PURPOSE OF THIS SUBROUTINE:
588 : // This subroutine creates a "Newell" surface normal vector from the vector list
589 : // for a surface face.
590 :
591 : // METHODOLOGY EMPLOYED:
592 : // na
593 :
594 : // REFERENCES:
595 : // September 2010: from OpenGL.org
596 : // Begin Function CalculateSurfaceNormal (Input Polygon) Returns Vector
597 : // Set Vertex Normal to (0, 0, 0)
598 : // Begin Cycle for Index in [0, Polygon.vertexNumber)
599 : // Set Vertex Current to Polygon.verts[Index]
600 : // Set Vertex Next to Polygon.verts[(Index plus 1) mod Polygon.vertexNumber]
601 : // Set Normal.x to Sum of Normal.x and (multiply (Current.y minus Next.y) by (Current.z plus Next.z)
602 : // Set Normal.y to Sum of Normal.y and (multiply (Current.z minus Next.z) by (Current.x plus Next.x)
603 : // Set Normal.z to Sum of Normal.z and (multiply (Current.x minus Next.x) by (Current.y plus Next.y)
604 : // End Cycle
605 : // Returning Normalize(Normal)
606 : // End Function
607 :
608 : // USE STATEMENTS:
609 : // na
610 :
611 : // Argument array dimensioning
612 :
613 : // Locals
614 : // SUBROUTINE ARGUMENT DEFINITIONS:
615 :
616 : // SUBROUTINE PARAMETER DEFINITIONS:
617 : // na
618 :
619 : // INTERFACE BLOCK SPECIFICATIONS
620 : // na
621 :
622 : // DERIVED TYPE DEFINITIONS
623 : // na
624 :
625 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
626 : // TYPE(Vector) :: U
627 : // TYPE(Vector) :: V
628 : int Side;
629 : int curVert;
630 : int nextVert;
631 : Real64 xvalue;
632 : Real64 yvalue;
633 : Real64 zvalue;
634 :
635 48198 : OutNewellSurfaceNormalVector = 0.0;
636 48198 : xvalue = 0.0;
637 48198 : yvalue = 0.0;
638 48198 : zvalue = 0.0;
639 :
640 : // IF (NSides > 3) THEN
641 241663 : for (Side = 1; Side <= NSides; ++Side) {
642 193465 : curVert = Side;
643 193465 : nextVert = Side + 1;
644 193465 : if (nextVert > NSides) nextVert = 1;
645 193465 : xvalue += (VList(curVert).y - VList(nextVert).y) * (VList(curVert).z + VList(nextVert).z);
646 193465 : yvalue += (VList(curVert).z - VList(nextVert).z) * (VList(curVert).x + VList(nextVert).x);
647 193465 : zvalue += (VList(curVert).x - VList(nextVert).x) * (VList(curVert).y + VList(nextVert).y);
648 : }
649 : // ELSE ! Triangle
650 : // U=VList(2)-VList(1)
651 : // V=VList(3)-VList(1)
652 : // xvalue=(U%y*V%z)-(U%z*V%y)
653 : // yvalue=(U%z*V%x)-(U%x*V%z)
654 : // zvalue=(U%x*V%y)-(U%y*V%x)
655 : // ENDIF
656 :
657 48198 : OutNewellSurfaceNormalVector.x = xvalue;
658 48198 : OutNewellSurfaceNormalVector.y = yvalue;
659 48198 : OutNewellSurfaceNormalVector.z = zvalue;
660 48198 : OutNewellSurfaceNormalVector = VecNormalize(OutNewellSurfaceNormalVector);
661 48198 : }
662 :
663 6442 : void CompareTwoVectors(Vector const &vector1, // standard vector
664 : Vector const &vector2, // standard vector
665 : bool &areSame, // true if the two vectors are the same within specified tolerance
666 : Real64 const tolerance // specified tolerance
667 : )
668 : {
669 :
670 : // SUBROUTINE INFORMATION:
671 : // AUTHOR Linda Lawrie
672 : // DATE WRITTEN February 2012
673 : // MODIFIED na
674 : // RE-ENGINEERED na
675 :
676 : // PURPOSE OF THIS SUBROUTINE:
677 : // This routine will provide the ability to compare two vectors (e.g. surface normals)
678 : // to be the same within a specified tolerance.
679 :
680 : // METHODOLOGY EMPLOYED:
681 : // compare each element (x,y,z)
682 :
683 : // REFERENCES:
684 : // na
685 :
686 : // USE STATEMENTS:
687 : // na
688 :
689 : // Locals
690 : // SUBROUTINE ARGUMENT DEFINITIONS:
691 :
692 : // SUBROUTINE PARAMETER DEFINITIONS:
693 : // na
694 :
695 : // INTERFACE BLOCK SPECIFICATIONS:
696 : // na
697 :
698 : // DERIVED TYPE DEFINITIONS:
699 : // na
700 :
701 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
702 : // na
703 6442 : areSame = true;
704 6442 : if (std::abs(vector1.x - vector2.x) > tolerance) areSame = false;
705 6442 : if (std::abs(vector1.y - vector2.y) > tolerance) areSame = false;
706 6442 : if (std::abs(vector1.z - vector2.z) > tolerance) areSame = false;
707 6442 : }
708 :
709 40707 : void CalcCoPlanarNess(Array1D<Vector> &Surf, int const NSides, bool &IsCoPlanar, Real64 &MaxDist, int &ErrorVertex)
710 : {
711 :
712 : // SUBROUTINE INFORMATION:
713 : // AUTHOR Linda Lawrie
714 : // DATE WRITTEN June 2004
715 : // MODIFIED na
716 : // RE-ENGINEERED na
717 :
718 : // PURPOSE OF THIS SUBROUTINE:
719 : // This subroutine provides the calculation to determine if the
720 : // surface is planar or not.
721 :
722 : // METHODOLOGY EMPLOYED:
723 : // na
724 :
725 : // REFERENCES:
726 : // Eric W. Weisstein. "Coplanar." From MathWorld--A Wolfram Web Resource.
727 : // http://mathworld.wolfram.com/Coplanar.html
728 :
729 : // USE STATEMENTS:
730 : // na
731 :
732 : // Argument array dimensioning
733 40707 : EP_SIZE_CHECK(Surf, NSides);
734 :
735 : // Locals
736 : // SUBROUTINE ARGUMENT DEFINITIONS:
737 :
738 : // SUBROUTINE PARAMETER DEFINITIONS:
739 40707 : Real64 constexpr DistTooSmall(1.e-4);
740 :
741 : // INTERFACE BLOCK SPECIFICATIONS
742 : // na
743 :
744 : // DERIVED TYPE DEFINITIONS
745 : // na
746 :
747 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
748 : bool plerror;
749 : Real64 dist;
750 :
751 : // Object Data
752 40707 : PlaneEq NewellPlane;
753 :
754 40707 : IsCoPlanar = true;
755 40707 : MaxDist = 0.0;
756 40707 : ErrorVertex = 0;
757 :
758 : // Use first three to determine plane
759 40707 : PlaneEquation(Surf, NSides, NewellPlane, plerror);
760 :
761 204230 : for (int vert = 1; vert <= NSides; ++vert) {
762 163523 : dist = Pt2Plane(Surf(vert), NewellPlane);
763 163523 : if (std::abs(dist) > MaxDist) {
764 10052 : MaxDist = std::abs(dist);
765 10052 : ErrorVertex = vert;
766 : }
767 : }
768 :
769 40707 : if (std::abs(MaxDist) > DistTooSmall) IsCoPlanar = false;
770 40707 : }
771 :
772 15 : std::vector<int> PointsInPlane(Array1D<Vector> &BaseSurf, int const BaseSides, Array1D<Vector> &QuerySurf, int const QuerySides, bool &ErrorFound)
773 : {
774 15 : std::vector<int> pointIndices;
775 :
776 15 : PlaneEq NewellPlane;
777 15 : PlaneEquation(BaseSurf, BaseSides, NewellPlane, ErrorFound);
778 :
779 75 : for (int vert = 1; vert <= QuerySides; ++vert) {
780 60 : Real64 dist = Pt2Plane(QuerySurf(vert), NewellPlane);
781 60 : if (std::abs(dist) < 1.e-4) { // point on query surface is co-planar with base surface
782 30 : pointIndices.push_back(vert);
783 : }
784 : }
785 15 : return pointIndices;
786 : }
787 :
788 4795 : Real64 CalcPolyhedronVolume(EnergyPlusData &state, Polyhedron const &Poly)
789 : {
790 :
791 : // SUBROUTINE INFORMATION:
792 : // AUTHOR Linda Lawrie
793 : // DATE WRITTEN June 2004
794 : // MODIFIED na
795 : // RE-ENGINEERED na
796 :
797 : // PURPOSE OF THIS SUBROUTINE:
798 : // This subroutine provides the volume calculation for a polyhedron
799 : // (i.e. Zone).
800 :
801 : // METHODOLOGY EMPLOYED:
802 : // na
803 :
804 : // REFERENCES:
805 : // Conversations with Bill Carroll, LBNL.
806 :
807 : // SUBROUTINE LOCAL VARIABLE DECLARATIONS:
808 : Real64 PyramidVolume;
809 :
810 : // Object Data
811 9590 : Vector p3FaceOrigin;
812 :
813 4795 : Real64 Volume = 0.0;
814 :
815 38121 : for (int NFace = 1; NFace <= Poly.NumSurfaceFaces; ++NFace) {
816 33326 : p3FaceOrigin = Poly.SurfaceFace(NFace).FacePoints(2);
817 33326 : PyramidVolume = dot(Poly.SurfaceFace(NFace).NewellAreaVector, (p3FaceOrigin - state.dataVectors->p0));
818 33326 : Volume += PyramidVolume / 3.0;
819 : }
820 9590 : return Volume;
821 : }
822 :
823 2313 : } // namespace EnergyPlus::Vectors
|
