Skip to content

Commit

Permalink
FDS Source: TEST_CC_FOR_BLOCKING, add test for centroid location in C…
Browse files Browse the repository at this point in the history 
…FACES.
  • Loading branch information
@marcosvanella
marcosvanella committed Nov 7, 2023
1 parent 5a42cc6 commit ad92769
Showing 1 changed file with 36 additions and 28 deletions.
64 changes: 36 additions & 28 deletions Source/geom.f90
View file
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

0 comments on commit ad92769

Please sign in to comment.