Merge pull request #12181 from ericvmueller/master

Update cut face normals when split cell blocking is used; add CELL_BLOCK_ORIENTATION input

Manuals/FDS_User_Guide/FDS_User_Guide.tex

@@ -1535,7 +1535,7 @@ \subsection{Reading Geometry Node Locations And Connectivity Data From Binary}
 \subsection{Handling Split Cells and Thin Geometries}
 \label{info:thin_geom}
 
-One current limitation of the unstructured geometry algorithm is that a Cartesian cell cannot be split by a thin geometry or corner, potentially creating two background pressure zones within one Cartesian cell.  By default, if a cell is split by a thin geometry, FDS will fill in the smaller of the two disjoint volumes.  However, it may be desired to leave one side perfectly smooth while allowing the other to absorb the geometry change.  For example, you want the inside wall of a tunnel or ceiling to be smooth but you are not as concerned with the outer portion for your problem.  In this scenario, you can control which side of the geometry gets filled by setting {\ct CELL\_BLOCK\_IOR} on the {\ct GEOM} line.  An example of how this works is shown in Fig.~\ref{figure:thin_geometry}. Currently, this functionality is limited to geometries where the simple {\ct IOR} can define the direction.  For example, this would not work for a hollow sphere.
+One current limitation of the unstructured geometry algorithm is that a Cartesian cell cannot be split by a thin geometry or corner, potentially creating two background pressure zones within one Cartesian cell.  By default, if a cell is split by a thin geometry, FDS will fill in the smaller of the two disjoint volumes.  However, it may be desired to leave one side perfectly smooth while allowing the other to absorb the geometry change.  For example, you want the inside wall of a tunnel or ceiling to be smooth but you are not as concerned with the outer portion for your problem.  In this scenario, you can control which side of the geometry gets filled by setting either {\ct CELL\_BLOCK\_IOR} or {\ct CELL\_BLOCK\_ORIENTATION} on the {\ct GEOM} line. The former is an integer indicating a Cartesian axis while the latter allows the specification of a vector in any direction.  An example of how this works is shown in Fig.~\ref{figure:thin_geometry}. Currently, this functionality is limited to geometries where the simple {\ct IOR} can define the direction.  For example, this would not work for a hollow sphere.
 
 
 \begin{figure}
@@ -11549,6 +11549,7 @@ \section{\texorpdfstring{{\tt GEOM}}{GEOM} (Unstructured Geometry Parameters)}
 {\ct BINARY\_FILE}     & Character              & Section~\ref{subsec:readbin}            &           &  {\ct 'null'}            \\ \hline
 {\ct BNDF\_GEOM}       & Logical                & Section~\ref{info:BNDF}            &           &  {\ct T}                 \\ \hline
 {\ct CELL\_BLOCK\_IOR} & Integer                & Section~\ref{info:thin_geom}       &           &  0                       \\ \hline
+{\ct CELL\_BLOCK\_ORIENTATION} & Real Triplet   & Section~\ref{info:thin_geom}       &           &  0.0,0.0,0.0             \\ \hline
 {\ct COLOR}            & Character              & Section~\ref{info:colors}          &           &  {\ct 'null'}            \\ \hline
 {\ct CYLINDER\_LENGTH} & Real                   & Section~\ref{info:GEOM_cyls}            &   m       &                          \\ \hline
 {\ct CYLINDER\_RADIUS} & Real                   & Section~\ref{info:GEOM_cyls}            &   m       &                          \\ \hline

Source/geom.f90

@@ -785,7 +785,7 @@ SUBROUTINE SET_CUTCELLS_3D
 
 LOGICAL, SAVE :: FIRST_CALL_ARG=.TRUE., FIRST_CALL_ARG2=.TRUE.
 
-REAL(EB):: VERT_AUX(IAXIS:KAXIS)
+REAL(EB):: VERT_AUX(IAXIS:KAXIS),CELL_BLOCK_ORIENTATION(IAXIS:KAXIS)
 INTEGER :: ING,INOD,IWSEL,IEL,FACE_AUX(NOD1:NOD3),VOL_AUX(NOD1:NOD4),N_SPCELLCF_TOT,N_SPCELL_TOT
 CHARACTER(100) :: FILENAME
 
@@ -1380,47 +1380,52 @@ SUBROUTINE SET_CUTCELLS_3D
       DO ICC1=1,MESHES(NM)%N_CUTCELL_MESH+MESHES(NM)%N_GCCUTCELL_MESH
          CC=>MESHES(NM)%CUT_CELL(ICC1); IF(CC%NCELL<2) CYCLE
          ! Find if any GEOMETRY related to CC_INBOUNDARY faces has CELL_BLOCK_IOR>0:
-         CELL_BLOCK_IOR=0
+         CELL_BLOCK_IOR=0; CELL_BLOCK_ORIENTATION = 0._EB
          NCELL_LOOP_1 : DO J=1,CC%NCELL
             DO I=2,CC%CCELEM(1,J)+1
                IF(CC%FACE_LIST(1,CC%CCELEM(I,J))==CC_FTYPE_CFINB) THEN
                   ICF=CC%FACE_LIST(4,CC%CCELEM(I,J)); JCF=CC%FACE_LIST(5,CC%CCELEM(I,J))
                   IG=MESHES(NM)%CUT_FACE(ICF)%BODTRI(1,JCF)
                   IF(IG>0) THEN
-                     IF (ABS(GEOMETRY(IG)%CELL_BLOCK_IOR)>0) THEN
+                     IF (NORM2(GEOMETRY(IG)%CELL_BLOCK_ORIENTATION(IAXIS:KAXIS))>TWO_EPSILON_EB) THEN
+                        CELL_BLOCK_ORIENTATION = GEOMETRY(IG)%CELL_BLOCK_ORIENTATION
+                     ELSEIF (ABS(GEOMETRY(IG)%CELL_BLOCK_IOR)>0) THEN
                         CELL_BLOCK_IOR = GEOMETRY(IG)%CELL_BLOCK_IOR
                         EXIT NCELL_LOOP_1
                      ENDIF
                   ENDIF
                ENDIF
             ENDDO
          ENDDO NCELL_LOOP_1
-         ALLOCATE(VOLUME(1:CC%NCELL))
-         ! Make search for double precision min/max unambiguous.
-         SELECT CASE(CELL_BLOCK_IOR)
-         CASE(-IAXIS,IAXIS); VOLUME(1:CC%NCELL) = CC%XYZCEN(IAXIS,1:CC%NCELL)
-         CASE(-JAXIS,JAXIS); VOLUME(1:CC%NCELL) = CC%XYZCEN(JAXIS,1:CC%NCELL)
-         CASE(-KAXIS,KAXIS); VOLUME(1:CC%NCELL) = CC%XYZCEN(KAXIS,1:CC%NCELL)
-         CASE DEFAULT;       VOLUME(1:CC%NCELL) = CC%VOLUME(1:CC%NCELL)
-         END SELECT
-         ! IF(NORM2(CELL_BLOCK_NVEC)>TWO_EPSILON_EB) THEN ! Block in the NVEC direction:
-         !    DO J=1,CC%NCELL; VOLUME(J) = DOT_PRODUCT(CELL_BLOCK_NVEC,CC%XYZCEN(IAXIS:KAXIS,J)); ENDDO
-         ! ENDIF
-
+         ALLOCATE(VOLUME(1:CC%NCELL)); VOLUME(1:CC%NCELL) = CC%VOLUME(1:CC%NCELL)
          DO J=1,CC%NCELL; VOLUME(J) = VOLUME(J)*(1._EB+REAL(J-1,EB)*GEOMEPS); ENDDO
-         SELECT CASE(CELL_BLOCK_IOR)
-         CASE(-IAXIS); I=MAXLOC(VOLUME(1:CC%NCELL),DIM=1)
-         CASE( IAXIS); I=MINLOC(VOLUME(1:CC%NCELL),DIM=1)
-         CASE(-JAXIS); I=MAXLOC(VOLUME(1:CC%NCELL),DIM=1)
-         CASE( JAXIS); I=MINLOC(VOLUME(1:CC%NCELL),DIM=1)
-         CASE(-KAXIS); I=MAXLOC(VOLUME(1:CC%NCELL),DIM=1)
-         CASE( KAXIS); I=MINLOC(VOLUME(1:CC%NCELL),DIM=1)
-         CASE DEFAULT; I=MAXLOC(VOLUME(1:CC%NCELL),DIM=1)
-         END SELECT
+         I=MAXLOC(VOLUME(1:CC%NCELL),DIM=1)
+         IF(NORM2(CELL_BLOCK_ORIENTATION)>TWO_EPSILON_EB) THEN 
+            ! Cell Block Orientation:
+            DO J=1,CC%NCELL; VOLUME(J) = DOT_PRODUCT(CELL_BLOCK_ORIENTATION,CC%XYZCEN(IAXIS:KAXIS,J)); ENDDO
+            DO J=1,CC%NCELL; VOLUME(J) = VOLUME(J)*(1._EB+REAL(J-1,EB)*GEOMEPS); ENDDO
+            I=MINLOC(VOLUME(1:CC%NCELL),DIM=1)
+         ELSEIF(ABS(CELL_BLOCK_IOR)>0) THEN
+            ! Make search for double precision min/max unambiguous.
+            SELECT CASE(CELL_BLOCK_IOR)
+            CASE(-IAXIS,IAXIS); VOLUME(1:CC%NCELL) = CC%XYZCEN(IAXIS,1:CC%NCELL)
+            CASE(-JAXIS,JAXIS); VOLUME(1:CC%NCELL) = CC%XYZCEN(JAXIS,1:CC%NCELL)
+            CASE(-KAXIS,KAXIS); VOLUME(1:CC%NCELL) = CC%XYZCEN(KAXIS,1:CC%NCELL)
+            END SELECT
+            DO J=1,CC%NCELL; VOLUME(J) = VOLUME(J)*(1._EB+REAL(J-1,EB)*GEOMEPS); ENDDO
+            SELECT CASE(CELL_BLOCK_IOR)
+            CASE(-IAXIS); I=MAXLOC(VOLUME(1:CC%NCELL),DIM=1)
+            CASE( IAXIS); I=MINLOC(VOLUME(1:CC%NCELL),DIM=1)
+            CASE(-JAXIS); I=MAXLOC(VOLUME(1:CC%NCELL),DIM=1)
+            CASE( JAXIS); I=MINLOC(VOLUME(1:CC%NCELL),DIM=1)
+            CASE(-KAXIS); I=MAXLOC(VOLUME(1:CC%NCELL),DIM=1)
+            CASE( KAXIS); I=MINLOC(VOLUME(1:CC%NCELL),DIM=1)
+            END SELECT
+         ENDIF
          DEALLOCATE(VOLUME)
          NCELL_LOOP_2 : DO J=1,CC%NCELL
             IF(J==I) CYCLE NCELL_LOOP_2
-            CC%NOADVANCE(J) = BLOCKED_SPLIT_CELL
+            IF(CC%NOADVANCE(J)==NOT_BLOCKED) CC%NOADVANCE(J) = BLOCKED_SPLIT_CELL
          ENDDO NCELL_LOOP_2
       ENDDO
    ENDIF
@@ -1473,7 +1478,7 @@ SUBROUTINE SET_CUTCELLS_3D
             END SELECT
          ENDDO
          IF(SUM_FACE>1 .OR. SUM_CCELL>0) CYCLE
-         CC%NOADVANCE(J)=BLOCKED_CAVITY_CELL
+         IF(CC%NOADVANCE(J)==NOT_BLOCKED) CC%NOADVANCE(J)=BLOCKED_CAVITY_CELL
          K=K+1
       ENDDO
    ENDDO
@@ -2455,7 +2460,7 @@ SUBROUTINE TAG_CC_BLOCKING_REFINEMENT
                         IFACE  = CC%CCELEM(IFC,JCC)
                         IF( (CC%FACE_LIST(1,IFACE)==CC_FTYPE_RCGAS .OR. CC%FACE_LIST(1,IFACE)==CC_FTYPE_CFGAS) .AND. &
                            CC%FACE_LIST(2,IFACE)==LOHIF .AND. CC%FACE_LIST(3,IFACE)==X1AXIS ) THEN
-                           CC%NOADVANCE(JCC)=BLOCKED_REFI_INTER
+                           IF(CC%NOADVANCE(JCC)==NOT_BLOCKED) CC%NOADVANCE(JCC)=BLOCKED_REFI_INTER
                            CYCLE JCC_LOOP_1
                         ENDIF
                      ENDDO
@@ -2528,7 +2533,7 @@ SUBROUTINE TAG_CC_BLOCKING_REFINEMENT
                                  IFACE  = CC%CCELEM(IFC,JCC)
                                  IF( (CC%FACE_LIST(1,IFACE)==CC_FTYPE_RCGAS .OR. CC%FACE_LIST(1,IFACE)==CC_FTYPE_CFGAS) .AND. &
                                     CC%FACE_LIST(2,IFACE)==LOHIF .AND. CC%FACE_LIST(3,IFACE)==X1AXIS ) THEN
-                                    CC%NOADVANCE(JCC)=BLOCKED_REFI_INTER
+                                    IF(CC%NOADVANCE(JCC)==NOT_BLOCKED) CC%NOADVANCE(JCC)=BLOCKED_REFI_INTER
                                     CYCLE JCC_LOOP_3
                                  ENDIF
                               ENDDO
@@ -2602,7 +2607,8 @@ SUBROUTINE TAG_BLOCK_CELL(NM,II1,JJ1,KK1,NOM,IIO1,JJO1,KKO1,FINE_CELL)
       ! Here Test if cut-cells in II,KK,KK are blocked or not in IIO,JJO,KKO:
       IF(ICC>0) THEN
          IF(M2%CCVAR(IIO1,JJO1,KKO1,CC_CGSC)==CC_SOLID) THEN
-            M%CUT_CELL(ICC)%NOADVANCE(1:M%CUT_CELL(ICC)%NCELL) = BLOCKED_REFI_INTER
+            WHERE(M%CUT_CELL(ICC)%NOADVANCE(1:M%CUT_CELL(ICC)%NCELL)==NOT_BLOCKED) &
+                  M%CUT_CELL(ICC)%NOADVANCE(1:M%CUT_CELL(ICC)%NCELL) = BLOCKED_REFI_INTER
          ELSE; CALL TEST_CC_FOR_BLOCKING(NM,ICC,NOM,IIO1,JJO1,KKO1,ICC2)
          ENDIF
       ENDIF
@@ -2921,7 +2927,7 @@ SUBROUTINE TEST_CC_FOR_BLOCKING(NM,ICC,NOM,IIO1,JJO1,KKO1,ICC2)
       ! ENDIF
    ENDDO INBFC_LOC_LOOP
    ! Here set no ADVANCE if BLOCK_CELL=T:
-   IF(BLOCK_CELL) CC%NOADVANCE(JCC) = BLOCKED_REFI_INTER
+   IF(BLOCK_CELL .AND. CC%NOADVANCE(JCC)==NOT_BLOCKED) CC%NOADVANCE(JCC) = BLOCKED_REFI_INTER
 ENDDO JCC_LOOP
 
 ! IF(NM==1 .AND. ICC<30) CLOSE(LU_CCELL)
@@ -4148,7 +4154,7 @@ SUBROUTINE GET_EXT_INB_CUTFACES_TO_CFACE
    MESHES(NM)%INTERNAL_CFACE_CELLS_LB = MESHES(NM)%N_EXTERNAL_CFACE_CELLS + MESHES(NM)%N_INTWALL_CFACE_CELLS
    ! Define pointers among CC_INBOUNDARY CUT_FACE and CFACE (N_INTERNAL_CFACE_CELLS):
    DO ICF=1,MESHES(NM)%N_CUTFACE_MESH
-      CF => CUT_FACE(ICF); IF(CF%STATUS /= CC_INBOUNDARY) CYCLE
+      CF => MESHES(NM)%CUT_FACE(ICF); IF(CF%STATUS /= CC_INBOUNDARY) CYCLE
       I = CF%IJK(IAXIS); J = CF%IJK(JAXIS); K = CF%IJK(KAXIS)
       ! Don't count INB cut-faces inside an OBST:
       IF (CELL(CELL_INDEX(I,J,K))%SOLID) CYCLE
@@ -7715,6 +7721,7 @@ SUBROUTINE DROP_CUTFACE(NM,FTYPE,I,J,K,ILHF,X1AXIS,IFC,JFC)
           CF%BODTRI(     :,IND(DUM)) = CF%BODTRI(     :,DUM)
           CF%SURF_INDEX(   IND(DUM)) = CF%SURF_INDEX(   DUM)
           CF%BLK_TAG(      IND(DUM)) = CF%BLK_TAG(      DUM)
+          CF%CFACE_ORIGIN( IND(DUM)) = CF%CFACE_ORIGIN( DUM)
           CF%AREA_ADJUST(  IND(DUM)) = CF%AREA_ADJUST(  DUM)
       ENDIF
       DO ILH=LOW_IND,CT
@@ -7743,6 +7750,7 @@ SUBROUTINE DROP_CUTFACE(NM,FTYPE,I,J,K,ILHF,X1AXIS,IFC,JFC)
    IF(FTYPE==CC_FTYPE_CFINB) THEN
       CF%BODTRI(     :,CF%NFACE) = CC_UNDEFINED
       CF%SURF_INDEX(   CF%NFACE) = CC_UNDEFINED
+      CF%CFACE_ORIGIN( CF%NFACE) = CC_UNDEFINED
    ENDIF
    DEALLOCATE(IND)
 ENDIF NFACE_IF_1
@@ -8330,7 +8338,7 @@ SUBROUTINE GET_CELL_LINK_INFO(NM)
 DO ICC=1,M%N_CUTCELL_MESH
    DO JCC=1,M%CUT_CELL(ICC)%NCELL
       IF(M%CUT_CELL(ICC)%IJK_LINK(1,JCC)==CC_SOLID) THEN
-         M%CUT_CELL(ICC)%NOADVANCE(JCC) = BLOCKED_UNLINK_CELL
+         IF(M%CUT_CELL(ICC)%NOADVANCE(JCC)==NOT_BLOCKED) M%CUT_CELL(ICC)%NOADVANCE(JCC) = BLOCKED_UNLINK_CELL
          M%CUT_CELL(ICC)%LINK_LEV(JCC)  = CC_UNDEFINED
          M%CUT_CELL(ICC)%IJK_LINK(2:5,JCC)= CC_UNDEFINED
       ELSEIF(M%CUT_CELL(ICC)%LINK_LEV(JCC) < M%FINEST_LINK_LEV) THEN
@@ -8433,6 +8441,7 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB,
 
 USE GEOMETRY_FUNCTIONS, ONLY : SEARCH_OTHER_MESHES
 USE MEMORY_FUNCTIONS, ONLY: ALLOCATE_STORAGE
+USE MATH_FUNCTIONS, ONLY : CROSS_PRODUCT
 
 ! Routine that initializes new CFACE with index CFACE_INDEX.
 ! Geometry information for CFACE is loaded from MESHES(NM)%CUT_FACE(ICF)%AREA(IFACE), etc.
@@ -8453,9 +8462,12 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB,
 TYPE (WALL_TYPE), POINTER :: WC=>NULL()
 TYPE (MESH_TYPE), POINTER :: M
 TYPE (CFACE_TYPE), POINTER :: CFA
+TYPE (CC_CUTFACE_TYPE), POINTER :: CF
+
 
 M => MESHES(NM)
 SF=> SURFACE(SURF_INDEX)
+CF=> CUT_FACE(ICF)
 
 STAGE_FLG_BRANCH : SELECT CASE(STAGE_FLG)
 
@@ -8469,17 +8481,17 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB,
 
    CFA%SURF_INDEX = SURF_INDEX
 
-   BC%X = CUT_FACE(ICF)%XYZCEN(IAXIS,IFACE)
-   BC%Y = CUT_FACE(ICF)%XYZCEN(JAXIS,IFACE)
-   BC%Z = CUT_FACE(ICF)%XYZCEN(KAXIS,IFACE)
-   CFA%AREA = CUT_FACE(ICF)%AREA(IFACE)
+   BC%X = CF%XYZCEN(IAXIS,IFACE)
+   BC%Y = CF%XYZCEN(JAXIS,IFACE)
+   BC%Z = CF%XYZCEN(KAXIS,IFACE)
+   CFA%AREA = CF%AREA(IFACE)
 
    ! Now populate cut-face information:
    CFA%CUT_FACE_IND1 = ICF
    CFA%CUT_FACE_IND2 = IFACE
 
    INS_INB_COND_1 : IF (IS_INB) THEN
-      CFA%VEL_ERR_NEW=CUT_FACE(ICF)%VEL(IFACE) - 0._EB ! Assumes zero veloc of solid.
+      CFA%VEL_ERR_NEW=CF%VEL(IFACE) - 0._EB ! Assumes zero veloc of solid.
 
       ! Check if fire spreads radially over this surface type
       IF (SF%FIRE_SPREAD_RATE>0._EB) THEN
@@ -8489,21 +8501,28 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB,
       ENDIF
 
       ! Normal to cut-face:
-      IBOD =CUT_FACE(ICF)%BODTRI(1,IFACE)
-      IWSEL=CUT_FACE(ICF)%BODTRI(2,IFACE)
-      CFA%NVEC(IAXIS:KAXIS) = GEOMETRY(IBOD)%FACES_NORMAL(IAXIS:KAXIS,IWSEL)
+      V2(IAXIS:KAXIS) = CF%XYZVERT(IAXIS:KAXIS,CF%CFELEM(2,IFACE))-CF%XYZCEN(IAXIS:KAXIS,IFACE)
+      V3(IAXIS:KAXIS) = CF%XYZVERT(IAXIS:KAXIS,CF%CFELEM(3,IFACE))-CF%XYZCEN(IAXIS:KAXIS,IFACE)
+      CALL CROSS_PRODUCT(CFA%NVEC(IAXIS:KAXIS),V2,V3)
+      IF(NORM2(CFA%NVEC)>TWO_EPSILON_EB .AND. CF%CFACE_ORIGIN(IFACE)==BLOCKED_SPLIT_CELL) THEN
+         CFA%NVEC(IAXIS:KAXIS) = CFA%NVEC(IAXIS:KAXIS)/NORM2(CFA%NVEC)
+      ELSE
+         IBOD =CF%BODTRI(1,IFACE)
+         IWSEL=CF%BODTRI(2,IFACE)
+         CFA%NVEC(IAXIS:KAXIS) = GEOMETRY(IBOD)%FACES_NORMAL(IAXIS:KAXIS,IWSEL)
+      ENDIF
 
       ! Boundary CFACES processed are defined of type SOLID_BOUNDARY
       CFA%BOUNDARY_TYPE = SOLID_BOUNDARY
 
       ! Might need to rethink this, but for the time being...
-      BC%II = CUT_FACE(ICF)%IJK(IAXIS)
-      BC%JJ = CUT_FACE(ICF)%IJK(JAXIS)
-      BC%KK = CUT_FACE(ICF)%IJK(KAXIS)
+      BC%II = CF%IJK(IAXIS)
+      BC%JJ = CF%IJK(JAXIS)
+      BC%KK = CF%IJK(KAXIS)
 
-      BC%IIG = CUT_FACE(ICF)%IJK(IAXIS)
-      BC%JJG = CUT_FACE(ICF)%IJK(JAXIS)
-      BC%KKG = CUT_FACE(ICF)%IJK(KAXIS)
+      BC%IIG = CF%IJK(IAXIS)
+      BC%JJG = CF%IJK(JAXIS)
+      BC%KKG = CF%IJK(KAXIS)
 
    ELSE INS_INB_COND_1 ! External mesh boundary CFACE
 
@@ -8532,14 +8551,14 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB,
       ENDIF
    ENDIF INS_INB_COND_1
 
-   B1%AREA = CUT_FACE(ICF)%AREA(IFACE) ! Init to CFACE AREA.
+   B1%AREA = CF%AREA(IFACE) ! Init to CFACE AREA.
 
 CASE(INTEGER_TWO) ! Assign AREA_ADJUST for CFACE, BCs information for CFACE.
 
    CFA  => M%CFACE(CFACE_INDEX)
    B1   => M%BOUNDARY_PROP1(CFA%B1_INDEX)
    ! First: Assign AREA_ADJUST for CFACEs.
-   B1%AREA_ADJUST = CUT_FACE(ICF)%AREA_ADJUST(IFACE)
+   B1%AREA_ADJUST = CF%AREA_ADJUST(IFACE)
 
 CASE(INTEGER_THREE)
 
@@ -8551,11 +8570,11 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB,
 
       ! Associated cut-cell location in CUT_CELL array.
       ! This CFACE initialization assumes TMP,RHO,ZZ have been initialized in cut-cell ICC,JCC.
-      ICC = CUT_FACE(ICF)%CELL_LIST(2,LOW_IND,IFACE)
-      JCC = CUT_FACE(ICF)%CELL_LIST(3,LOW_IND,IFACE)
+      ICC = CF%CELL_LIST(2,LOW_IND,IFACE)
+      JCC = CF%CELL_LIST(3,LOW_IND,IFACE)
 
       ! Set TMP_F to Surface value and rest to ambient in underlying cartesian cell.
-      CFA%TMP_G = TMP_0(CUT_FACE(ICF)%IJK(KAXIS))
+      CFA%TMP_G = TMP_0(CF%IJK(KAXIS))
       IF (SF%TMP_FRONT > 0._EB) THEN
          B1%TMP_F = SF%TMP_FRONT
       ELSE
@@ -8574,7 +8593,7 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB,
       ! Assign normal velocity to CFACE from SURF input:
       B1%U_NORMAL_0 = SF%VEL
       ! Assign normal velocity from VOLUME_FLOW :
-      IBOD =CUT_FACE(ICF)%BODTRI(1,IFACE)
+      IBOD =CF%BODTRI(1,IFACE)
       IF(IBOD>0 .AND. ABS(SF%VOLUME_FLOW)>=TWO_EPSILON_EB) B1%U_NORMAL_0 = SF%VOLUME_FLOW / FDS_AREA_GEOM(SURF_INDEX,IBOD)
       ! Assign normal velocity from MASS_FLUX_TOTAL :
       IF(ABS(SF%MASS_FLUX_TOTAL)>=TWO_EPSILON_EB) B1%U_NORMAL_0 = SF%MASS_FLUX_TOTAL / RHOA * B1%AREA_ADJUST
@@ -8583,8 +8602,8 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB,
 
       ! Case of exposed Backing we need to find CFACE_INDEX of BACK CFACE.
       IF (SF%BACKING==EXPOSED .AND. SF%THERMAL_BC_INDEX==THERMALLY_THICK) THEN
-         IG  = CUT_FACE(ICF)%BODTRI(1,IFACE)
-         TRI = CUT_FACE(ICF)%BODTRI(2,IFACE)
+         IG  = CF%BODTRI(1,IFACE)
+         TRI = CF%BODTRI(2,IFACE)
          XP(IAXIS:KAXIS)  = (/ BC%X, BC%Y, BC%Z /) ! CFACE centroid location.
          RDIR(IAXIS:KAXIS)= - GEOMETRY(IG)%FACES_NORMAL(IAXIS:KAXIS,TRI) ! Normal into the body.
          TRI_LOOP : DO IWSEL=1,GEOMETRY(IG)%N_FACES
@@ -8720,8 +8739,8 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB,
          B1%U_NORMAL_0 = WC_B1%U_NORMAL_0
 
          ! Here downscale velocity:
-         IF (IFACE==CUT_FACE(ICF)%NFACE) WC_B1%U_NORMAL_0 = &
-         WC_B1%U_NORMAL_0 * SUM(CUT_FACE(ICF)%AREA(1:CUT_FACE(ICF)%NFACE))/WC_B1%AREA
+         IF (IFACE==CF%NFACE) WC_B1%U_NORMAL_0 = &
+         WC_B1%U_NORMAL_0 * SUM(CF%AREA(1:CF%NFACE))/WC_B1%AREA
 
          ! Vegetation T_IGN setup:
          B1%T_IGN      = WC_B1%T_IGN
@@ -21646,7 +21665,7 @@ SUBROUTINE READ_GEOM
 REAL(EB) :: SPHERE_ORIGIN(3),SPHERE_RADIUS,TEXTURE_ORIGIN(3),TEXTURE_SCALE(2),XB(6),DX,BOX_XYZ(3),&
             ZMIN,VOLUME,TXMIN,TXMAX,TYMIN,TYMAX,TX,TY,DV1(MAX_DIM),DV2(MAX_DIM),&
             NVECI(MAX_DIM),DXCEN(MAX_DIM),DOTI,TRANSPARENCY,CYLINDER_ORIGIN(3),CYLINDER_AXIS(3),&
-            CYLINDER_RADIUS,CYLINDER_LENGTH,EXTRUDE
+            CYLINDER_RADIUS,CYLINDER_LENGTH,EXTRUDE,CELL_BLOCK_ORIENTATION(3)
 
 INTEGER :: MAX_IDS=0,MAX_SURF_IDS=0,MAX_ZVALS=0,MAX_VERTS=0,MAX_FACES=0,MAX_VOLUS=0,MAX_POLY_VERTS=0,&
            N_VERTS,N_FACES,N_FACES_TEMP,N_VOLUS,N_ZVALS,N_SURF_ID,N_SURF_ID2,N_POLY_VERTS,&
@@ -21693,7 +21712,7 @@ SUBROUTINE READ_GEOM
 INTEGER, PARAMETER :: DELTA_GEOM_LINE=1000
 INTEGER :: GEOM_LINE_SIZE
 
-NAMELIST /GEOM/ BNDF_GEOM,BINARY_FILE,CELL_BLOCK_IOR,COLOR,CYLINDER_ORIGIN,CYLINDER_AXIS,&
+NAMELIST /GEOM/ BNDF_GEOM,BINARY_FILE,CELL_BLOCK_IOR,CELL_BLOCK_ORIENTATION,COLOR,CYLINDER_ORIGIN,CYLINDER_AXIS,&
                 CYLINDER_RADIUS,CYLINDER_LENGTH,CYLINDER_NSEG_THETA,CYLINDER_NSEG_AXIS,&
                 EXTRUDE,EXTEND_TERRAIN,FACES,ID,IJK,IS_TERRAIN,MOVE_ID,N_LAT,N_LEVELS,N_LONG,POLY,&
                 RGB,SPHERE_ORIGIN,SPHERE_RADIUS,SPHERE_TYPE,SURF_ID,SURF_IDS,SURF_ID6,&
@@ -21770,6 +21789,7 @@ SUBROUTINE READ_GEOM
       CALL COLOR2RGB(RGB,COLOR)
    ENDIF
    G%CELL_BLOCK_IOR = CELL_BLOCK_IOR
+   G%CELL_BLOCK_ORIENTATION = CELL_BLOCK_ORIENTATION
    G%RGB = RGB
    G%TRANSPARENCY = TRANSPARENCY
    N_VERTS=0
@@ -23792,6 +23812,7 @@ SUBROUTINE SET_GEOM_DEFAULTS
    BINARY_FILE = 'null'
    RGB=-1
    CELL_BLOCK_IOR=0
+   CELL_BLOCK_ORIENTATION = 0._EB
 
 END SUBROUTINE SET_GEOM_DEFAULTS