Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

FDS Source: TEST_CC_FOR_BLOCKING, add test for centroid location in C… #12213

Merged
merged 1 commit into from
Nov 7, 2023
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
64 changes: 36 additions & 28 deletions Source/geom.f90
Original file line number Diff line number Diff line change
Expand Up @@ -1519,7 +1519,7 @@ SUBROUTINE SET_CUTCELLS_3D
! DO ICC=1,M%N_CUTCELL_MESH+M%N_GCCUTCELL_MESH
! CC=>M%CUT_CELL(ICC)
! DO JCC=1,CC%NCELL
! IF(CC%NOADVANCE(JCC)) WRITE(787,'(5I8)') ICC,JCC,CC%IJK(IAXIS:KAXIS)
! IF(CC%NOADVANCE(JCC)>0) WRITE(LU_ERR,*) 'B',NM,';',ICC,JCC,CC%IJK(IAXIS:KAXIS),CC%NCELL
! ENDDO
! ENDDO
! CLOSE(787)
Expand Down Expand Up @@ -2689,7 +2689,7 @@ SUBROUTINE TEST_CC_FOR_BLOCKING(NM,ICC,NOM,IIO1,JJO1,KKO1,ICC2)
REAL(EB),ALLOCATABLE, DIMENSION(:,:):: XYZVERTIJK,XYZVERTSTN
REAL(EB):: XYZCEN(MAX_DIM),NVEC(MAX_DIM),P0(MAX_DIM),A,B,C,D,XYZ_P(MAX_DIM),PTCEN(IAXIS:JAXIS),X1F,XC2(MAX_DIM),XC3(MAX_DIM),&
XLO,XHI,YLO,YHI,ZLO,ZHI,XLO2,XHI2,YLO2,YHI2,ZLO2,ZHI2,CFCEN(MAX_DIM),XYZC(MAX_DIM,1),N(MAX_DIM,1),S(MAX_DIM,1),&
T(MAX_DIM,1),TBN(MAX_DIM,MAX_DIM),XYZCSTN(MAX_DIM,1),NN(MAX_DIM,1),XN_CEN,XN_INT
T(MAX_DIM,1),TBN(MAX_DIM,MAX_DIM),XYZCSTN(MAX_DIM,1),NN(MAX_DIM,1),XN_CEN,XN_INT,XYZC2(IAXIS:KAXIS,1)
REAL(EB), PARAMETER :: SCALE_FCT=1.E-4_EB
LOGICAL :: IN_CFACE,BLOCK_CELL,FGPOINT
! INTEGER :: LU_CCELL
Expand Down Expand Up @@ -2836,6 +2836,13 @@ SUBROUTINE TEST_CC_FOR_BLOCKING(NM,ICC,NOM,IIO1,JJO1,KKO1,ICC2)
N(:,1) = -NVEC; S(:,1) = XC2/NORM2(XC2); CALL CROSS_PRODUCT(T(:,1),N(:,1),S(:,1))
TBN(1,:)= S(:,1); TBN(2,:)= T(:,1); TBN(3,:)= N(:,1)

! Check that cut-face centroid is within its polygon.
XYZC2(IAXIS:KAXIS,1) = CFCEN(IAXIS:KAXIS); XYZCSTN = MATMUL(TBN,XYZC2)
DO IV = 1,NVERT2; XYZVERTSTN(:,IV) = MATMUL(TBN,XYZVERTIJK(:,IV))-XYZCSTN(:,1); ENDDO
CFELEM2(1:VERT_CUTFACE2) =M%CUT_FACE(IFC1)%CFELEM(1:VERT_CUTFACE2,JFC1)
PTCEN(IAXIS:JAXIS) = 0._EB; CALL POINT_IN_POLYGON(PTCEN,VERT_CUTFACE2,CFELEM2,NVERT2,1,2,XYZVERTSTN,IN_CFACE)
IF(.NOT.IN_CFACE) CYCLE

! Run again over all CFACES of the JCC cut-cell (except IFC) and check for other intersections within their polygons:
! 1. First of all compute XYZCENSTN, allocate XYZVERTSTN and populate it. Compute XYZVERTSTN-XYZCENSTN.
XYZCSTN = MATMUL(TBN,XYZC)
Expand Down Expand Up @@ -3139,10 +3146,10 @@ SUBROUTINE GET_CC_FACE_CELL_LIST_INFO(NM,PHASE)
ICF1=M%FCVAR(I,J,K,CC_IDCF,X1AXIS); CF=>M%CUT_FACE(ICF1)
WRITE(33,'(I8,I8,I8,I8,I8)') CF%IJK(1:4),CF%NFACE,CF%STATUS
DO ICF2=1,CF%NFACE
WRITE(33,'(I8,3F12.8,F12.8)') ICF2,CF%XYZCEN(:,ICF2),CF%AREA(ICF2)
WRITE(33,'(I8,3F16.8,F16.8)') ICF2,CF%XYZCEN(IAXIS:KAXIS,ICF2),CF%AREA(ICF2)
ICC=CF%CELL_LIST(2,LOW_IND,ICF2); JCC=CF%CELL_LIST(3,LOW_IND,ICF2)
WRITE(33,'(3I8,3F12.8,F12.8)') M%CUT_CELL(ICC)%IJK(1:3),&
M%CUT_CELL(ICC)%VOLUME(JCC),M%CUT_CELL(ICC)%XYZCEN(:,JCC)
WRITE(33,'(3I8,F16.8,3F16.8)') M%CUT_CELL(ICC)%IJK(1:3),&
M%CUT_CELL(ICC)%VOLUME(JCC),M%CUT_CELL(ICC)%XYZCEN(IAXIS:KAXIS,JCC)
CC=>M%CUT_CELL(ICC)
IFACE = CC%CCELEM(CF%CELL_LIST(4,LOW_IND,ICF2)+1,JCC)
IF(ICF1/=CC%FACE_LIST(4,IFACE) .OR. ICF2/=CC%FACE_LIST(5,IFACE)) THEN
Expand All @@ -3152,8 +3159,8 @@ SUBROUTINE GET_CC_FACE_CELL_LIST_INFO(NM,PHASE)

IF(CF%STATUS==CC_GASPHASE) THEN
ICC=CF%CELL_LIST(2,HIGH_IND,ICF2); JCC=CF%CELL_LIST(3,HIGH_IND,ICF2)
WRITE(33,'(3I8,3F12.8,F12.8)') M%CUT_CELL(ICC)%IJK(1:3),&
M%CUT_CELL(ICC)%VOLUME(JCC),M%CUT_CELL(ICC)%XYZCEN(:,JCC)
WRITE(33,'(3I8,F16.8,3F16.8)') M%CUT_CELL(ICC)%IJK(1:3),&
M%CUT_CELL(ICC)%VOLUME(JCC),M%CUT_CELL(ICC)%XYZCEN(IAXIS:KAXIS,JCC)
CC=>M%CUT_CELL(ICC)
IFACE = CC%CCELEM(CF%CELL_LIST(4,HIGH_IND,ICF2)+1,JCC)
IF(ICF1/=CC%FACE_LIST(4,IFACE) .OR. ICF2/=CC%FACE_LIST(5,IFACE)) THEN
Expand All @@ -3170,10 +3177,10 @@ SUBROUTINE GET_CC_FACE_CELL_LIST_INFO(NM,PHASE)
ICF1=M%CCVAR(I,J,K,CC_IDCF); CF=>M%CUT_FACE(ICF1)
WRITE(33,'(I8,I8,I8,I8,I8)') CF%IJK(1:4),CF%NFACE,CF%STATUS
DO ICF2=1,CF%NFACE
WRITE(33,'(I8,3F12.8,F12.8)') ICF2,CF%XYZCEN(:,ICF2),CF%AREA(ICF2)
WRITE(33,'(I8,3F16.8,F16.8)') ICF2,CF%XYZCEN(IAXIS:KAXIS,ICF2),CF%AREA(ICF2)
ICC=CF%CELL_LIST(2,LOW_IND,ICF2); JCC=CF%CELL_LIST(3,LOW_IND,ICF2)
WRITE(33,'(3I8,3F12.8,F12.8)') M%CUT_CELL(ICC)%IJK(1:3),&
M%CUT_CELL(ICC)%VOLUME(JCC),M%CUT_CELL(ICC)%XYZCEN(:,JCC)
WRITE(33,'(3I8,F16.8,3F16.8)') M%CUT_CELL(ICC)%IJK(1:3),&
M%CUT_CELL(ICC)%VOLUME(JCC),M%CUT_CELL(ICC)%XYZCEN(IAXIS:KAXIS,JCC)
CC=>M%CUT_CELL(ICC)
IFACE = CC%CCELEM(CF%CELL_LIST(4,LOW_IND,ICF2)+1,JCC)
IF(ICF1/=CC%FACE_LIST(4,IFACE) .OR. ICF2/=CC%FACE_LIST(5,IFACE)) THEN
Expand Down Expand Up @@ -3206,7 +3213,7 @@ SUBROUTINE GET_CC_FACE_CELL_LIST_INFO(NM,PHASE)
! Face JCF:
ICF1=CE%FACE_LIST(1,JCF,JCE); ICF2=CE%FACE_LIST(2,JCF,JCE)
CF=>M%CUT_FACE(ICF1)
WRITE(33,'(4I8,I8,3F12.8,F12.8)') CF%IJK(1:4),ICF2,CF%XYZCEN(:,ICF2),CF%AREA(ICF2)
WRITE(33,'(4I8,I8,3F16.8,F16.8)') CF%IJK(1:4),ICF2,CF%XYZCEN(IAXIS:KAXIS,ICF2),CF%AREA(ICF2)
ENDDO
ENDDO
ENDIF
Expand Down Expand Up @@ -4838,13 +4845,13 @@ SUBROUTINE BLOCK_SMALL_UNLINKED_CUTCELLS(NM,NBLKCELLS)
CE=>MESHES(NM)%CUT_EDGE(IEC)
WRITE(33,'(I6,I6,I6,I6,4I6)') IEC,CE%STATUS,CE%NVERT,CE%NEDGE,CE%IJK(1:4)
DO IVR=1,CE%NVERT
WRITE(33,'(3F18.14)') CE%XYZVERT(IAXIS:KAXIS,IVR)
WRITE(33,'(3F18.10)') CE%XYZVERT(IAXIS:KAXIS,IVR)
ENDDO
DO IVR=1,CE%NVERT
WRITE(33,'(4I6)') CE%VERT_LIST(1:4,IVR)
ENDDO
DO JEC=1,CE%NEDGE
WRITE(33,'(2I6,6I6)') CE%CEELEM(NOD1:NOD2,JEC),CE%INDSEG(1:4,JEC)
WRITE(33,'(2I8,6I8)') CE%CEELEM(NOD1:NOD2,JEC),CE%INDSEG(1:4,JEC)
ENDDO
DO JEC=1,CE%NEDGE
WRITE(33,'(2I6,2I6,2I6,2I6)') CE%FACE_LIST(1:2,-2,JEC),CE%FACE_LIST(1:2,-1,JEC),&
Expand All @@ -4861,7 +4868,7 @@ SUBROUTINE BLOCK_SMALL_UNLINKED_CUTCELLS(NM,NBLKCELLS)
IF(ALLOCATED(CF%EDGE_LIST)) NSEG=SIZE(CF%EDGE_LIST,DIM=2); IF(CF%STATUS==CC_GASPHASE) NSEG=NSEG-1
WRITE(33,'(I6,I6,I6,I6,I6,4I6)') IFC,CF%STATUS,CF%NVERT,CF%NFACE,NSEG,CF%IJK(1:4)
DO IVR=1,CF%NVERT
WRITE(33,'(3F18.14)') CF%XYZVERT(IAXIS:KAXIS,IVR)
WRITE(33,'(3F18.10)') CF%XYZVERT(IAXIS:KAXIS,IVR)
ENDDO
DO JFC=1,CF%NFACE
WRITE(33,'(I6,I6)') CF%CFELEM(1,JFC),CF%CEDGES(1,JFC)
Expand Down Expand Up @@ -4943,13 +4950,13 @@ SUBROUTINE BLOCK_SMALL_UNLINKED_CUTCELLS(NM,NBLKCELLS)
CE=>MESHES(NM)%CUT_EDGE(IEC)
WRITE(33,'(I6,I6,I6,I6,4I6)') IEC,CE%STATUS,CE%NVERT,CE%NEDGE,CE%IJK(1:4)
DO IVR=1,CE%NVERT
WRITE(33,'(3F18.14)') CE%XYZVERT(IAXIS:KAXIS,IVR)
WRITE(33,'(3F18.10)') CE%XYZVERT(IAXIS:KAXIS,IVR)
ENDDO
DO IVR=1,CE%NVERT
WRITE(33,'(4I6)') CE%VERT_LIST(1:4,IVR)
ENDDO
DO JEC=1,CE%NEDGE
WRITE(33,'(2I6,6I6)') CE%CEELEM(NOD1:NOD2,JEC),CE%INDSEG(1:4,JEC)
WRITE(33,'(2I8,6I8)') CE%CEELEM(NOD1:NOD2,JEC),CE%INDSEG(1:4,JEC)
ENDDO
DO JEC=1,CE%NEDGE
WRITE(33,'(2I6,2I6,2I6,2I6)') CE%FACE_LIST(1:2,-2,JEC),CE%FACE_LIST(1:2,-1,JEC),&
Expand All @@ -4966,10 +4973,10 @@ SUBROUTINE BLOCK_SMALL_UNLINKED_CUTCELLS(NM,NBLKCELLS)
IF(ALLOCATED(CF%EDGE_LIST)) NSEG=SIZE(CF%EDGE_LIST,DIM=2); IF(CF%STATUS==CC_GASPHASE) NSEG=NSEG-1
WRITE(33,'(I6,I6,I6,I6,I6,4I6)') IFC,CF%STATUS,CF%NVERT,CF%NFACE,NSEG,CF%IJK(1:4)
DO IVR=1,CF%NVERT
WRITE(33,'(3F18.14)') CF%XYZVERT(IAXIS:KAXIS,IVR)
WRITE(33,'(3F18.10)') CF%XYZVERT(IAXIS:KAXIS,IVR)
ENDDO
DO JFC=1,CF%NFACE
WRITE(33,'(I6,I6)') CF%CFELEM(1,JFC),CF%CEDGES(1,JFC)
WRITE(33,'(I8,I8)') CF%CFELEM(1,JFC),CF%CEDGES(1,JFC)
DO DUM=1,CF%CFELEM(1,JFC)
WRITE(33,'(I6)') CF%CFELEM(DUM+1,JFC)
ENDDO
Expand Down Expand Up @@ -5003,12 +5010,13 @@ SUBROUTINE BLOCK_SMALL_UNLINKED_CUTCELLS(NM,NBLKCELLS)
DO JEC=1,CE%NEDGE
INOD1=CE%CEELEM(NOD1,JEC)
INOD2=CE%CEELEM(NOD2,JEC)
WRITE(33,'(I8,3F12.8,3F12.8)') CE%IJK(4),CE%XYZVERT(:,INOD1),CE%XYZVERT(:,INOD2)
WRITE(33,'(I8,3F16.8,3F16.8)') CE%IJK(4),CE%XYZVERT(:,INOD1),CE%XYZVERT(:,INOD2)
WRITE(33,'(I8,I8,A,4I8,4I8)') CE%IJK(4),JEC,';',CE%VERT_LIST(1:4,INOD1),CE%VERT_LIST(1:4,INOD2)
IF(CE%VERT_LIST(1,INOD1)==CE%VERT_LIST(1,INOD2) .AND. &
CE%VERT_LIST(2,INOD1)==CE%VERT_LIST(2,INOD2) .AND. &
CE%VERT_LIST(3,INOD1)==CE%VERT_LIST(3,INOD2) .AND. &
CE%VERT_LIST(4,INOD1)==CE%VERT_LIST(4,INOD2)) THEN
IF(CE%VERT_LIST(1,INOD1)/=CC_VTYPE_NINB) &
WRITE(LU_ERR,*) 'Edge with same node types=',IEC,JEC,CE%NEDGE,CE%XYZVERT(:,INOD1),&
CE%XYZVERT(:,INOD2),CE%VERT_LIST(1:4,INOD1)
ENDIF
Expand All @@ -5035,7 +5043,7 @@ SUBROUTINE BLOCK_SMALL_UNLINKED_CUTCELLS(NM,NBLKCELLS)
DO JEC=1,CE%NEDGE
INOD1=CE%CEELEM(NOD1,JEC)
INOD2=CE%CEELEM(NOD2,JEC)
WRITE(33,'(I8,3F12.8,3F12.8)') CE%IJK(4),CE%XYZVERT(:,INOD1),CE%XYZVERT(:,INOD2)
WRITE(33,'(I8,3F16.8,3F16.8)') CE%IJK(4),CE%XYZVERT(:,INOD1),CE%XYZVERT(:,INOD2)
WRITE(33,'(I8,I8,A,4I8,4I8)') CE%IJK(4),JEC,';',CE%VERT_LIST(1:4,INOD1),CE%VERT_LIST(1:4,INOD2)
ENDDO
ENDIF
Expand All @@ -5061,7 +5069,7 @@ SUBROUTINE BLOCK_SMALL_UNLINKED_CUTCELLS(NM,NBLKCELLS)
WRITE(LU_ERR,*) 'CUT FACE does not match FCVAR',I,J,K,X1AXIS,':',CF%IJK(IAXIS:KAXIS+1)
WRITE(33,'(I8,I8,I8,I8,I8)') CF%IJK(1:4),CF%NFACE
DO JEC=1,CF%NFACE
WRITE(33,'(I8,3F12.8,F12.8)') CF%IJK(4),CF%XYZCEN(:,JEC),CF%AREA(JEC)
WRITE(33,'(I8,3F16.8,F16.8)') CF%IJK(4),CF%XYZCEN(:,JEC),CF%AREA(JEC)
ENDDO
ENDIF
ENDDO
Expand All @@ -5083,7 +5091,7 @@ SUBROUTINE BLOCK_SMALL_UNLINKED_CUTCELLS(NM,NBLKCELLS)
WRITE(LU_ERR,*) 'CUT CELL does not match CCVAR',I,J,K,':',CC%IJK(IAXIS:KAXIS)
WRITE(33,'(I8,I8,I8,I8,I8)') CC%IJK(1:3),CC%NCELL
DO JEC=1,CC%NCELL
WRITE(33,'(I8,3F12.8,F12.8)') JEC,CC%XYZCEN(:,JEC),CC%VOLUME(JEC)
WRITE(33,'(I8,3F16.8,F16.8)') JEC,CC%XYZCEN(:,JEC),CC%VOLUME(JEC)
ENDDO
ENDIF
ENDDO
Expand Down Expand Up @@ -9412,11 +9420,11 @@ SUBROUTINE GET_REGULAR_CUT_EDGES_BC(NM)
IE=M%CC_RCEDGE(ECOUNT)%IJK(KAXIS+1)
SELECT CASE(IE)
CASE(IAXIS)
WRITE(LU_DB_SETCC,'(4I4,4F13.8)') I,J,K,IE,DX(I),XC(I),Y(J),Z(K)
WRITE(LU_DB_SETCC,'(4I8,4F16.8)') I,J,K,IE,DX(I),XC(I),Y(J),Z(K)
CASE(JAXIS)
WRITE(LU_DB_SETCC,'(4I4,4F13.8)') I,J,K,IE,DY(J),X(I),YC(J),Z(K)
WRITE(LU_DB_SETCC,'(4I8,4F16.8)') I,J,K,IE,DY(J),X(I),YC(J),Z(K)
CASE(KAXIS)
WRITE(LU_DB_SETCC,'(4I4,4F13.8)') I,J,K,IE,DZ(K),X(I),Y(J),ZC(K)
WRITE(LU_DB_SETCC,'(4I8,4F16.8)') I,J,K,IE,DZ(K),X(I),Y(J),ZC(K)
END SELECT
ENDDO
CLOSE(LU_DB_SETCC)
Expand Down Expand Up @@ -10103,11 +10111,11 @@ SUBROUTINE GET_SOLID_CUTCELL_EDGES_BC(NM)
IE=MESHES(NM)%CC_IBEDGE(ECOUNT)%IJK(KAXIS+1)
SELECT CASE(IE)
CASE(IAXIS)
WRITE(LU_DB_SETCC,'(4I4,4F13.8)') I,J,K,IE,DX(I),XC(I),Y(J),Z(K)
WRITE(LU_DB_SETCC,'(4I8,4F16.8)') I,J,K,IE,DX(I),XC(I),Y(J),Z(K)
CASE(JAXIS)
WRITE(LU_DB_SETCC,'(4I4,4F13.8)') I,J,K,IE,DY(J),X(I),YC(J),Z(K)
WRITE(LU_DB_SETCC,'(4I8,4F16.8)') I,J,K,IE,DY(J),X(I),YC(J),Z(K)
CASE(KAXIS)
WRITE(LU_DB_SETCC,'(4I4,4F13.8)') I,J,K,IE,DZ(K),X(I),Y(J),ZC(K)
WRITE(LU_DB_SETCC,'(4I8,4F16.8)') I,J,K,IE,DZ(K),X(I),Y(J),ZC(K)
END SELECT
ENDDO
CLOSE(LU_DB_SETCC)
Expand Down