FDS Source : address highly distorted cut-cells.

Author: marcosvanella

Source/geom.f90

@@ -38,6 +38,7 @@ MODULE COMPLEX_GEOMETRY
 REAL(EB), PARAMETER :: MIN_VOL_FACTOR   = 5.E-4_EB
 REAL(EB), PARAMETER :: ADIFF_INFO_FACTOR= 1.E-1_EB
 REAL(EB), PARAMETER :: SNAP_DIST_FACTOR = 1.E-5_EB
+REAL(EB), PARAMETER :: MIN_LENGTH_FACTOR= 1.E-2_EB
 
 INTEGER,  SAVE ::      NGUARD = 5        ! Layers of guard-cells.
 INTEGER,  SAVE ::      CCGUARD= 5 - 2    ! Layers of guard cut-cells.
@@ -103,12 +104,12 @@ MODULE COMPLEX_GEOMETRY
 INTEGER, PARAMETER :: CC_ETYPE_CFINB = 2 ! An INBOUNDARY cut-edge.
 
 ! Cut-faces types in FACE_LIST of CUT_CELL:
-INTEGER, PARAMETER :: CC_FTYPE_RGGAS = 0 ! This face of a cut-cell is a regular GASPHASE face.
+INTEGER, PARAMETER :: CC_FTYPE_RGGAS = 0 ! This face is a regular GASPHASE face.
 INTEGER, PARAMETER :: CC_FTYPE_CFGAS = 1 ! A GASPHASE cut-face or cell.
 INTEGER, PARAMETER :: CC_FTYPE_CFINB = 2 ! An INBOUNDARY cut-face.
 INTEGER, PARAMETER :: CC_FTYPE_SVERT = 3 ! A SOLID Vertex.
 ! Extra tagging parameters, for RESTRICT_EP:
-INTEGER, PARAMETER :: CC_FTYPE_RCGAS = 4 ! A Face between a regular cell and a cut-cell.
+INTEGER, PARAMETER :: CC_FTYPE_RCGAS = 4 ! A gas Face between a regular cell and a cut-cell.
 INTEGER, PARAMETER :: CC_FTYPE_CCGAS = 5 ! A regular gas cut-cell.
 INTEGER, PARAMETER :: CC_ETYPE_SCINB =12
 INTEGER, PARAMETER :: CC_ETYPE_RCGAS =14 ! A regular edge next to a cut-face and a regular gasphase face.
@@ -8209,22 +8210,50 @@ SUBROUTINE GET_CELL_LINK_INFO(NM)
 
 ! Local Variables:
 INTEGER :: ICC,JCC,ICC2,JCC2,JCC_LNK,I,J,K,I_LNK,J_LNK,K_LNK,IFC,IFC2,IFACE,IFACE2,IFACE3,IBOD,IWSEL,VAL_UNKZ,&
-           LINK_ITER,INGH,JNGH,KNGH,X1AXIS,ILH,INRM(1:3),DUM,LNK_LEV,ULINK_COUNT,LINK_LEV_UP
-REAL(EB):: AREA,AF,NRML(IAXIS:KAXIS),VAL_CVOL,CCVOL_THRES
-LOGICAL :: QUITLINK_FLG,CRTCELL_FLG,MASK(IAXIS:KAXIS)
+           LINK_ITER,INGH,JNGH,KNGH,X1AXIS,ILH,INRM(1:3),DUM,LNK_LEV,ULINK_COUNT,LINK_LEV_UP,AX_MIN,AX_OTHERS(2)
+REAL(EB):: AREA,AF,NRML(IAXIS:KAXIS),VAL_CVOL,CCVOL_THRES, XYZCELL(IAXIS:KAXIS,LOW_IND:HIGH_IND),&
+           MINMAX_XYZ_CC(IAXIS:KAXIS,LOW_IND:HIGH_IND),CELL_DELTA(IAXIS:KAXIS)
+LOGICAL :: QUITLINK_FLG,CRTCELL_FLG,MASK(IAXIS:KAXIS),BLOCK_SLIM_IF,SOLID_FACES
 CHARACTER(MESSAGE_LENGTH) :: UNLINKED_FILE
 INTEGER, SAVE :: LU_UNLNK
 LOGICAL, SAVE :: UNLINKED_1ST_CALL=.TRUE.
 TYPE (MESH_TYPE), POINTER :: M
+TYPE (CC_CUTCELL_TYPE), POINTER :: CC
 
 M => MESHES(NM)
 
 ! Initialize UNKZ, used here to define if cell has been noted in linking hierarchy. Assume CCVAR has been allocated:
 M%CCVAR(:,:,:,CC_UNKZ) = CC_UNDEFINED
 DO ICC=1,M%N_CUTCELL_MESH+M%N_GCCUTCELL_MESH
-   M%CUT_CELL(ICC)%UNKZ(:) = CC_UNDEFINED
-   DO JCC=1,M%CUT_CELL(ICC)%NCELL
-      IF (M%CUT_CELL(ICC)%NOADVANCE(JCC)>0) M%CUT_CELL(ICC)%IJK_LINK(1,JCC) = CC_SOLID
+   CC => M%CUT_CELL(ICC); I=CC%IJK(IAXIS); J=CC%IJK(JAXIS); K=CC%IJK(KAXIS)
+   ! Test for sliver trapped cells blocking:
+   XYZCELL(IAXIS,LOW_IND) = XFACE(I-1); XYZCELL(IAXIS,HIGH_IND) = XFACE(I);
+   XYZCELL(JAXIS,LOW_IND) = YFACE(J-1); XYZCELL(JAXIS,HIGH_IND) = YFACE(J);
+   XYZCELL(KAXIS,LOW_IND) = ZFACE(K-1); XYZCELL(KAXIS,HIGH_IND) = ZFACE(K);
+   MINMAX_XYZ_CC(IAXIS:KAXIS,LOW_IND) =  HUGE(EB)
+   MINMAX_XYZ_CC(IAXIS:KAXIS,HIGH_IND)= -HUGE(EB)
+   DO JCC=1,CC%NCELL
+      ! Get cut-cell bounding box:
+      CALL CUT_CELL_BOUNDING_BOX(NM,ICC,JCC,XYZCELL,MINMAX_XYZ_CC)            
+      ! Perform Tests:
+      DO DUM=IAXIS,KAXIS
+         CELL_DELTA(DUM) = ABS(MINMAX_XYZ_CC(DUM,HIGH_IND)-MINMAX_XYZ_CC(DUM,LOW_IND))
+      ENDDO
+      ! Axis with minimum width:
+      AX_MIN = MINLOC(CELL_DELTA(IAXIS:KAXIS),DIM=1)
+      SELECT CASE(AX_MIN)
+      CASE(IAXIS); AX_OTHERS(1:2) = (/ JAXIS, KAXIS /); SOLID_FACES = ALL(M%FCVAR(I-1:I,J,K,CC_FGSC,IAXIS)==CC_SOLID)
+      CASE(JAXIS); AX_OTHERS(1:2) = (/ IAXIS, KAXIS /); SOLID_FACES = ALL(M%FCVAR(I,J-1:J,K,CC_FGSC,JAXIS)==CC_SOLID)
+      CASE(KAXIS); AX_OTHERS(1:2) = (/ IAXIS, JAXIS /); SOLID_FACES = ALL(M%FCVAR(I,J,K-1:K,CC_FGSC,KAXIS)==CC_SOLID)
+      END SELECT   
+      ! Perform Test:
+      BLOCK_SLIM_IF = (CELL_DELTA(AX_MIN)<10._EB*MIN_LENGTH_FACTOR*CELL_DELTA(AX_OTHERS(1))) .AND. &
+                      (CELL_DELTA(AX_MIN)<10._EB*MIN_LENGTH_FACTOR*CELL_DELTA(AX_OTHERS(2))) 
+      IF(BLOCK_SLIM_IF .AND. SOLID_FACES) CC%NOADVANCE(JCC) = BLOCKED_SMALL_CELL
+   ENDDO
+   CC%UNKZ(:) = CC_UNDEFINED
+   DO JCC=1,CC%NCELL
+      IF (CC%NOADVANCE(JCC)>0) CC%IJK_LINK(1,JCC) = CC_SOLID
    ENDDO
 ENDDO
 
@@ -20339,11 +20368,13 @@ SUBROUTINE GET_CARTCELL_CUTCELLS(NM)
 REAL(EB),ALLOCATABLE, DIMENSION(:)   :: VOL ! Cut-cell volumes.
 INTEGER, ALLOCATABLE, DIMENSION(:)   :: NOADVANCE
 
-INTEGER :: IFACE, IEDGE, ISEG, SEG(NOD1:NOD2), ICELL, NFACEI
+REAL(EB) :: XYZCELL(IAXIS:KAXIS,LOW_IND:HIGH_IND),MINMAX_XYZ_CC(IAXIS:KAXIS,LOW_IND:HIGH_IND),CELL_DELTA(IAXIS:KAXIS)
+
+INTEGER :: IFACE, IEDGE, ISEG, SEG(NOD1:NOD2), ICELL, NFACEI, JCC, AX_MIN, AX_OTHERS(2) 
 LOGICAL :: INLIST, TEST1, TEST2, NEWFACE
 INTEGER :: NIEDGE, NEF, LOCSEG, JFACE, KFACE, NFACEK, NUM_FACE, NCUTCELL, NCFACE_CUTCELL
 INTEGER :: DFCT, CFELEM(5), CTVAL, CTVAL2, IBOD, ITRI, IDCF, MAXSEG, N_GAS_CFACES, NIBFACE, THRES, NSPCELL_LIST
-LOGICAL :: CYCLE_CELL
+LOGICAL :: CYCLE_CELL, BLOCK_SLIM_IF
 
 INTEGER :: IBNDINT
 LOGICAL, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: IJK_COUNT
@@ -21063,6 +21094,32 @@ SUBROUTINE GET_CARTCELL_CUTCELLS(NM)
          MESHES(NM)%CUT_CELL(NCUTCELL)%XYZCEN(IAXIS:KAXIS,1:NCELL)      = XYZCEN(IAXIS:KAXIS,1:NCELL)
          MESHES(NM)%CUT_CELL(NCUTCELL)%NOADVANCE(1:NCELL)               = NOADVANCE(1:NCELL)
 
+         ! Test for sliver cells blocking:
+         XYZCELL(IAXIS,LOW_IND) = XFACE(I-1); XYZCELL(IAXIS,HIGH_IND) = XFACE(I);
+         XYZCELL(JAXIS,LOW_IND) = YFACE(J-1); XYZCELL(JAXIS,HIGH_IND) = YFACE(J);
+         XYZCELL(KAXIS,LOW_IND) = ZFACE(K-1); XYZCELL(KAXIS,HIGH_IND) = ZFACE(K);
+         MINMAX_XYZ_CC(IAXIS:KAXIS,LOW_IND) =  HUGE(EB)
+         MINMAX_XYZ_CC(IAXIS:KAXIS,HIGH_IND)= -HUGE(EB)
+         DO JCC=1,NCELL
+            ! Get cut-cell bounding box:
+            CALL CUT_CELL_BOUNDING_BOX(NM,NCUTCELL,JCC,XYZCELL,MINMAX_XYZ_CC)            
+            ! Perform Tests:
+            DO MYAXIS=IAXIS,KAXIS
+               CELL_DELTA(MYAXIS) = ABS(MINMAX_XYZ_CC(MYAXIS,HIGH_IND)-MINMAX_XYZ_CC(MYAXIS,LOW_IND))
+            ENDDO
+            ! Axis with minimum width:
+            AX_MIN = MINLOC(CELL_DELTA(IAXIS:KAXIS),DIM=1)
+            SELECT CASE(AX_MIN)
+            CASE(IAXIS); AX_OTHERS(1:2) = (/ JAXIS, KAXIS /);
+            CASE(JAXIS); AX_OTHERS(1:2) = (/ IAXIS, KAXIS /);
+            CASE(KAXIS); AX_OTHERS(1:2) = (/ IAXIS, JAXIS /);
+            END SELECT   
+            ! Perform Test:
+            BLOCK_SLIM_IF = (CELL_DELTA(AX_MIN)<MIN_LENGTH_FACTOR*CELL_DELTA(AX_OTHERS(1))) .AND. &
+                            (CELL_DELTA(AX_MIN)<MIN_LENGTH_FACTOR*CELL_DELTA(AX_OTHERS(2))) 
+            IF(BLOCK_SLIM_IF) MESHES(NM)%CUT_CELL(NCUTCELL)%NOADVANCE(JCC) = BLOCKED_SMALL_CELL
+         ENDDO
+
          ! IF((NM==3 .AND. I==4 .AND. J==6 .AND. K==36)) THEN
          !    WRITE(LU_ERR,*) 'Found LARGE CUTCELL=',&
          !    MESHES(NM)%N_CUTCELL_MESH+MESHES(NM)%N_GCCUTCELL_MESH+1,VOL(1),SIZE(XYZVERT,DIM=2)
@@ -21167,6 +21224,60 @@ END SUBROUTINE REALLOCATE_FACE_CELL_VERTS
 END SUBROUTINE GET_CARTCELL_CUTCELLS
 
 
+! ------------------------ CUT_CELL_BOUNDING_BOX ------------------------------------
+
+SUBROUTINE CUT_CELL_BOUNDING_BOX(NM,ICC,JCC,XYZCELL,MINMAX_XYZ_JCC)
+
+! Computes bounding box for cut-cell (ICC,JCC) in mesh NM.
+! Underlaying cartesian cell bounds XYZCELL(IAXIS:KAXIS,LOW_IND:HIGH_IND) has to be provided.
+
+INTEGER, INTENT(IN) :: NM,ICC,JCC
+REAL(EB),INTENT(IN) :: XYZCELL(IAXIS:KAXIS,LOW_IND:HIGH_IND)
+REAL(EB),INTENT(OUT):: MINMAX_XYZ_JCC(IAXIS:KAXIS,LOW_IND:HIGH_IND)
+
+! Local Variables:
+INTEGER :: IFC,IFACE,LOHI,HILO,X1AXIS,IFCX,JFCX,IVERT,AXIS
+REAL(EB):: XYZFACE(IAXIS:KAXIS,LOW_IND:HIGH_IND),XYZ(IAXIS:KAXIS)
+TYPE(CC_CUTCELL_TYPE), POINTER :: CC
+TYPE(CC_CUTFACE_TYPE), POINTER :: CF
+
+CC => MESHES(NM)%CUT_CELL(ICC)
+
+! Get cut-cell bounding box:
+MINMAX_XYZ_JCC(IAXIS:KAXIS,LOW_IND) =  HUGE(EB)
+MINMAX_XYZ_JCC(IAXIS:KAXIS,HIGH_IND)= -HUGE(EB)
+DO IFC=1,CC%CCELEM(1,JCC) ! Loop over cut-faces boundary of this cell. 
+   IFACE=CC%CCELEM(IFC+1,JCC)
+   LOHI    =   CC%FACE_LIST(2,IFACE)
+   HILO    =   3-LOHI   !  2 for LOW_IND, 1 for HIGH_IND
+   X1AXIS  =   CC%FACE_LIST(3,IFACE)
+   IFCX    =   CC%FACE_LIST(4,IFACE)
+   JFCX    =   CC%FACE_LIST(5,IFACE)
+
+   SELECT CASE(CC%FACE_LIST(1,IFACE))
+   CASE(CC_FTYPE_RCGAS) ! Regular Gas face with a regular cell on one side and a cut-cell on the other.
+      XYZFACE = XYZCELL; XYZFACE(X1AXIS,HILO) = XYZFACE(X1AXIS,LOHI) ! Same location in X1AXIS for both sides of face.
+      DO AXIS=IAXIS,KAXIS
+         MINMAX_XYZ_JCC(AXIS,LOW_IND) = MIN(MINMAX_XYZ_JCC(AXIS,LOW_IND) ,XYZFACE(AXIS,LOW_IND))
+         MINMAX_XYZ_JCC(AXIS,HIGH_IND)= MAX(MINMAX_XYZ_JCC(AXIS,HIGH_IND),XYZFACE(AXIS,HIGH_IND))
+      ENDDO
+
+   CASE(CC_FTYPE_CFGAS,CC_FTYPE_CFINB) ! GAS or Boundary cut-face:
+      CF => MESHES(NM)%CUT_FACE(IFCX)
+      DO IVERT=1,CF%CFELEM(1,JFCX)
+         XYZ(IAXIS:KAXIS) = CF%XYZVERT(IAXIS:KAXIS,CF%CFELEM(IVERT+1,JFCX))
+         DO AXIS=IAXIS,KAXIS
+            MINMAX_XYZ_JCC(AXIS,LOW_IND) = MIN(MINMAX_XYZ_JCC(AXIS,LOW_IND) ,XYZ(AXIS))
+            MINMAX_XYZ_JCC(AXIS,HIGH_IND)= MAX(MINMAX_XYZ_JCC(AXIS,HIGH_IND),XYZ(AXIS))
+         ENDDO
+      ENDDO
+
+   END SELECT
+ENDDO
+
+END SUBROUTINE CUT_CELL_BOUNDING_BOX
+
+
 ! -------------------------CUT_CELL_ARRAY_REALLOC------------------------------------
 
 SUBROUTINE CUT_CELL_ARRAY_REALLOC(NM,ICC)