From ad9276997e949a431ad373afa0bb362835f6401b Mon Sep 17 00:00:00 2001 From: marcosvanella Date: Tue, 7 Nov 2023 15:31:57 -0500 Subject: [PATCH] FDS Source: TEST_CC_FOR_BLOCKING, add test for centroid location in CFACES. --- Source/geom.f90 | 64 +++++++++++++++++++++++++++---------------------- 1 file changed, 36 insertions(+), 28 deletions(-) diff --git a/Source/geom.f90 b/Source/geom.f90 index 793f85099c8..3fed07ec5da 100644 --- a/Source/geom.f90 +++ b/Source/geom.f90 @@ -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) @@ -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 @@ -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) @@ -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 @@ -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 @@ -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 @@ -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 @@ -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),& @@ -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) @@ -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),& @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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) @@ -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)