FDS Source: use trapezoidal integration for cone cut by vertical mesh face in CONE_MESH_INTERSECTION_VOLUME

Author: ericvmueller

Source/func.f90

@@ -5490,7 +5490,7 @@ REAL(EB) FUNCTION CIRCLE_CELL_INTERSECTION_AREA(X0,Y0,RAD,X1,X2,Y1,Y2)
 
 ! No overlap
 
-IF ((X2 <= X0-RAD) .OR. (X1 >= X0+RAD) .OR. (Y2 <= Y0-RAD) .OR. (Y1 >= Y0+RAD)) THEN
+IF ((X2 <= X0-RAD) .OR. (X1 >= X0+RAD) .OR. (Y2 <= Y0-RAD) .OR. (Y1 >= Y0+RAD) .OR. RAD<TWO_EPSILON_EB) THEN
    CIRCLE_CELL_INTERSECTION_AREA = 0._EB
    RETURN
 ENDIF
@@ -5527,7 +5527,7 @@ REAL(EB) FUNCTION CIRCLE_CELL_INTERSECTION_AREA(X0,Y0,RAD,X1,X2,Y1,Y2)
          RETURN
       ENDIF
       IF (Y1 >= Y0 .AND. Y1 < Y0+RAD) THEN
-         THETA = 2._EB*ACOS((X1-X0)/RAD)         
+         THETA = 2._EB*ACOS((Y1-Y0)/RAD)         
          CIRCLE_CELL_INTERSECTION_AREA = 0.5_EB*R2*(THETA-SIN(THETA))
          RETURN
       ENDIF
@@ -5646,9 +5646,10 @@ REAL(EB) FUNCTION CONE_MESH_INTERSECTION_VOLUME(NM,X0,Y0,Z0,RR0,RRI,HH0,G_FACTOR
 
 INTEGER, INTENT(IN) :: NM,G_FACTOR
 REAL(EB), INTENT(IN) :: X0,Y0,Z0,RR0,RRI,HH0
-INTEGER :: I,J,K,NX,NY,NZ
-REAL(EB) :: XX,YY,ZZ,R2,X_MIN,X_MAX,Y_MIN,Y_MAX,Z_MIN,Z_MAX,DX,DY,DZ,SHORT_DIMENSION,&
-            OUTER_BASE_AREA,INNER_BASE_AREA,HH_IN,RR0_BASE,RRI_BASE,RR0_TOP,RRI_TOP
+INTEGER :: K,NZ
+REAL(EB) :: X_MIN,X_MAX,Y_MIN,Y_MAX,Z_MIN,Z_MAX,DZ,Z1,Z2,&
+            OUTER_BASE_AREA,INNER_BASE_AREA,OUTER_TOP_AREA,INNER_TOP_AREA,HH_IN,&
+            RR0_BASE,RRI_BASE,RR0_TOP,RRI_TOP
 TYPE (MESH_TYPE), POINTER :: M
 
 CONE_MESH_INTERSECTION_VOLUME = 0._EB
@@ -5691,29 +5692,24 @@ REAL(EB) FUNCTION CONE_MESH_INTERSECTION_VOLUME(NM,X0,Y0,Z0,RR0,RRI,HH0,G_FACTOR
    RETURN
 ENDIF
 
-! Case 4: partial hollow cone intersection volume (numerical integration)
-SHORT_DIMENSION = MIN(X_MAX-X_MIN,Y_MAX-Y_MIN,Z_MAX-Z_MIN)
-! Limit resolution length scale 
-SHORT_DIMENSION = MAX(SHORT_DIMENSION, 0.1_EB*MIN(M%DXMIN,M%DYMIN,M%DZMIN))
-NX = CEILING(50*(X_MAX-X_MIN)/SHORT_DIMENSION)
-NY = CEILING(50*(Y_MAX-Y_MIN)/SHORT_DIMENSION)
-NZ = CEILING(50*(Z_MAX-Z_MIN)/SHORT_DIMENSION)
-DX = (X_MAX-X_MIN)/REAL(NX,EB)
-DY = (Y_MAX-Y_MIN)/REAL(NY,EB)
-DZ = (Z_MAX-Z_MIN)/REAL(NZ,EB)
+! Case 4: partial hollow cone intersection volume (trapezoidal integration)
+NZ = 100
+DZ = HH_IN/REAL(NZ,EB)
 
 DO K=1,NZ
-   ZZ = Z_MIN + (K-0.5_EB)*DZ
-   IF (ZZ<Z0 .OR. ZZ>Z0+HH0) CYCLE
-   DO J=1,NY
-      YY = Y_MIN + (J-0.5_EB)*DY
-      DO I=1,NX
-         XX = X_MIN + (I-0.5_EB)*DX
-         R2 = (XX-X0)**2+(YY-Y0)**2
-         IF ((R2<(RR0*(1._EB-G_FACTOR*(ZZ-Z0)/HH0))**2) .AND. (R2>(RRI*(1._EB-G_FACTOR*(ZZ-Z0)/HH0))**2)) &
-            CONE_MESH_INTERSECTION_VOLUME = CONE_MESH_INTERSECTION_VOLUME + DX*DY*DZ
-      ENDDO
-   ENDDO
+   Z1 = Z_MIN + (K-1)*DZ
+   Z2 = Z1 + DZ
+   RR0_BASE = RR0*(1-(Z1-Z0)/HH0)
+   RRI_BASE = RRI*(1-(Z1-Z0)/HH0)
+   RR0_TOP = RR0*(1-(Z2-Z0)/HH0)
+   RRI_TOP = RRI*(1-(Z2-Z0)/HH0)
+   OUTER_BASE_AREA = CIRCLE_CELL_INTERSECTION_AREA(X0,Y0,RR0_BASE,M%XS,M%XF,M%YS,M%YF)
+   INNER_BASE_AREA = CIRCLE_CELL_INTERSECTION_AREA(X0,Y0,RRI_BASE,M%XS,M%XF,M%YS,M%YF)
+   OUTER_TOP_AREA = CIRCLE_CELL_INTERSECTION_AREA(X0,Y0,RR0_TOP,M%XS,M%XF,M%YS,M%YF)
+   INNER_TOP_AREA = CIRCLE_CELL_INTERSECTION_AREA(X0,Y0,RRI_TOP,M%XS,M%XF,M%YS,M%YF)
+
+   CONE_MESH_INTERSECTION_VOLUME = CONE_MESH_INTERSECTION_VOLUME + &
+      DZ*0.5_EB*(OUTER_BASE_AREA+OUTER_TOP_AREA-INNER_BASE_AREA-INNER_TOP_AREA)
 ENDDO
 
 END FUNCTION CONE_MESH_INTERSECTION_VOLUME