Merge pull request #16238 from marcosvanella/master

FDS Source : Block special cut-cells in mesh corner edges.

Author: marcosvanella

Source/geom.f90

@@ -731,7 +731,7 @@ SUBROUTINE SET_CUTCELLS_3D
 REAL(EB), DIMENSION(MAX_DIM) :: PLNORMAL
 INTEGER,  DIMENSION(MAX_DIM) :: INDX1
 REAL(EB) :: X1PLN, X3RAY
-LOGICAL :: TRI_ONPLANE_ONLY, RAYTRACE_X2_ONLY
+LOGICAL :: TRI_ONPLANE_ONLY, RAYTRACE_X2_ONLY, MESH_BOUNDARY_FACE
 LOGICAL, SAVE :: FIRST_CALL = .TRUE.
 INTEGER :: NCUTFACE_IAXIS, NCUTFACE_JAXIS, NCUTFACE_KAXIS, ICE1, ICF1, NFACE, IERR, &
            NCUTEDGE_IBCC, NCUTEDGE_IBCF
@@ -1401,7 +1401,6 @@ SUBROUTINE SET_CUTCELLS_3D
       ENDIF
    ENDDO
 
-
    IF (ONE_CC_PER_CARTESIAN_CELL) THEN
       ! Here Block all cells that have volume less (or equal) than the first largest cell found.
       DO ICC1=1,MESHES(NM)%N_CUTCELL_MESH+MESHES(NM)%N_GCCUTCELL_MESH
@@ -1540,6 +1539,51 @@ SUBROUTINE SET_CUTCELLS_3D
          K=K+1
       ENDDO
    ENDDO
+
+   ! Block cells trapped against a mesh boundary whose only same-mesh gas opening is a small cut-face.
+   DO ICC1=1,MESHES(NM)%N_CUTCELL_MESH+MESHES(NM)%N_GCCUTCELL_MESH
+      CC=>MESHES(NM)%CUT_CELL(ICC1)
+      DO J=1,CC%NCELL
+         IF(CC%NOADVANCE(J)/=NOT_BLOCKED) CYCLE
+         SUM_FACE=0; SUM_CCELL=0; ICF1=0
+         DO I=1,CC%CCELEM(1,J)
+            IFACE  = CC%CCELEM(I+1,J)
+            SELECT CASE(CC%FACE_LIST(1,IFACE))
+            CASE(CC_FTYPE_CFGAS,CC_FTYPE_RCGAS)
+               SIDE   = CC%FACE_LIST(2,IFACE)
+               X1AXIS = CC%FACE_LIST(3,IFACE)
+               IIF = CC%IJK(IAXIS) + ADDI(SIDE,X1AXIS)
+               JJF = CC%IJK(JAXIS) + ADDJ(SIDE,X1AXIS)
+               KKF = CC%IJK(KAXIS) + ADDK(SIDE,X1AXIS)
+               MESH_BOUNDARY_FACE = .FALSE.
+               SELECT CASE(X1AXIS)
+               CASE(IAXIS); MESH_BOUNDARY_FACE = IIF==0 .OR. IIF==M%IBAR
+               CASE(JAXIS); MESH_BOUNDARY_FACE = JJF==0 .OR. JJF==M%JBAR
+               CASE(KAXIS); MESH_BOUNDARY_FACE = KKF==0 .OR. KKF==M%KBAR
+               END SELECT
+               IF(MESH_BOUNDARY_FACE) CYCLE
+               SELECT CASE(CC%FACE_LIST(1,IFACE))
+               CASE(CC_FTYPE_CFGAS)
+                  SUM_FACE = SUM_FACE+1
+                  ICF1     = CC%FACE_LIST(4,IFACE)
+               CASE(CC_FTYPE_RCGAS)
+                  SUM_CCELL = SUM_CCELL+1
+               END SELECT
+            END SELECT
+         ENDDO
+         IF(SUM_FACE/=1 .OR. SUM_CCELL>0) CYCLE
+         CF => M%CUT_FACE(ICF1)
+         IIF=CF%IJK(IAXIS); JJF=CF%IJK(JAXIS); KKF=CF%IJK(KAXIS)
+         SELECT CASE(CF%IJK(KAXIS+1))
+         CASE(IAXIS); ACRT = DY(JJF)*DZ(KKF)
+         CASE(JAXIS); ACRT = DZ(KKF)*DX(IIF)
+         CASE(KAXIS); ACRT = DX(IIF)*DY(JJF)
+         END SELECT
+         IF(SUM(CF%AREA(1:CF%NFACE))/ACRT>=CCVOL_LINK) CYCLE
+         CC%NOADVANCE(J)=BLOCKED_CAVITY_CELL
+         K=K+1
+      ENDDO
+   ENDDO
    IF (K>0) THEN
       CALL BLOCK_SMALL_UNLINKED_CUTCELLS(NM,SUM_CCELL)
       IF(SUM_CCELL>0) THEN ! Rebuild incidences and cell linking information: