diff --git a/Source/func.f90 b/Source/func.f90 index db19ad581f7..3be5d082e00 100644 --- a/Source/func.f90 +++ b/Source/func.f90 @@ -9436,7 +9436,8 @@ 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 +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 TYPE (MESH_TYPE), POINTER :: M CONE_MESH_INTERSECTION_VOLUME = 0._EB @@ -9448,18 +9449,38 @@ REAL(EB) FUNCTION CONE_MESH_INTERSECTION_VOLUME(NM,X0,Y0,Z0,RR0,RRI,HH0,G_FACTOR Y_MAX = MIN(M%YF,Y0+RR0) Z_MIN = MAX(M%ZS,Z0) Z_MAX = MIN(M%ZF,Z0+HH0) +HH_IN = Z_MAX-Z_MIN + +! Case 1: volume fully outside the mesh IF (X_MAX<=X_MIN .OR. Y_MAX<=Y_MIN .OR. Z_MAX<=Z_MIN) RETURN + +! Case 2: base area fully inside the mesh IF ((X_MIN-M%XS)>TWO_EPSILON_EB .AND. (M%XF-X_MAX)>TWO_EPSILON_EB & - .AND. (Y_MIN-M%YS)>TWO_EPSILON_EB .AND. (M%YF-Y_MAX)>TWO_EPSILON_EB & - .AND. (Z_MIN-M%ZS)>TWO_EPSILON_EB .AND. (M%ZF-Z_MAX)>TWO_EPSILON_EB) THEN + .AND. (Y_MIN-M%YS)>TWO_EPSILON_EB .AND. (M%YF-Y_MAX)>TWO_EPSILON_EB) THEN IF (G_FACTOR==0) THEN - CONE_MESH_INTERSECTION_VOLUME = PI*(RR0**2._EB-RRI**2._EB)*HH0 + ! Subtract inner from outer cylinder volume + CONE_MESH_INTERSECTION_VOLUME = PI*(RR0**2._EB-RRI**2._EB)*HH_IN ELSEIF (G_FACTOR==1) THEN - CONE_MESH_INTERSECTION_VOLUME = PI*(RR0**2._EB-RRI**2._EB)*HH0/3._EB + RR0_BASE = RR0*(1-(Z_MIN-Z0)/HH0) + RRI_BASE = RRI*(1-(Z_MIN-Z0)/HH0) + RR0_TOP = RR0*(1-(Z_MAX-Z0)/HH0) + RRI_TOP = RRI*(1-(Z_MAX-Z0)/HH0) + ! Subtract inner from outer truncated cone volume + CONE_MESH_INTERSECTION_VOLUME = PI*HH_IN/3._EB*& + (RR0_BASE**2._EB+RR0_BASE*RR0_TOP+RR0_TOP**2._EB-RRI_BASE**2._EB-RRI_BASE*RRI_TOP-RRI_TOP**2._EB) ENDIF RETURN ENDIF +! Case 3: partial hollow cylinder intersection volume +IF (G_FACTOR==0) THEN + OUTER_BASE_AREA = CIRCLE_CELL_INTERSECTION_AREA(X0,Y0,RR0,M%XS,M%XF,M%YS,M%YF) + INNER_BASE_AREA = CIRCLE_CELL_INTERSECTION_AREA(X0,Y0,RRI,M%XS,M%XF,M%YS,M%YF) + CONE_MESH_INTERSECTION_VOLUME = (OUTER_BASE_AREA-INNER_BASE_AREA)*HH_IN + 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))