Skip to content
Snippets Groups Projects
Commit d4f6b0d1 authored by Behrens, Prof. Dr. Jörn's avatar Behrens, Prof. Dr. Jörn
Browse files

added plotting capability for stacked spherical grid

parent 5f8f57ee
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
flash2d/src/options-sphere/IO_vtu_stacked.F90
+ 118
− 24
View file @ d4f6b0d1
......@@ -34,10 +34,11 @@ PUBLIC :: plot_vtu, t_vtu_data
! This defines a tetrahedron in VTK format
INTEGER(KIND=GRID_SI), PARAMETER :: i_vtucelltype = 10
#elif defined VTU_OUTPUTSTC
INTEGER(KIND=GRID_SI), PARAMETER :: i_nodespercell = GRID_elementnodes*2
INTEGER(KIND=GRID_SI), PARAMETER :: i_nodespercell = GRID_elementnodes
INTEGER(KIND=GRID_SI), PARAMETER :: i_nodesperface = GRID_edgenodes
! This defines a wedge/prism in VTK format
INTEGER(KIND=GRID_SI), PARAMETER :: i_nodespervtucell = GRID_elementnodes*2
INTEGER(KIND=GRID_SI), PARAMETER :: i_vtucelltype = 13
#endif
......@@ -111,7 +112,7 @@ SUBROUTINE plot_vtu(p_mesh, c_filename, &
LOGICAL, INTENT(IN), OPTIONAL :: l_continuous
INTEGER( KIND = GRID_SI) :: i_alct
INTEGER( KIND = GRID_SI) :: i_cnt
INTEGER( KIND = GRID_SI) :: i_cnt, i_hcnt
INTEGER( KIND = GRID_SI) :: i_fst
INTEGER( KIND = GRID_SI) :: i_ncnt
INTEGER( KIND = GRID_SI) :: i_lay
......@@ -122,6 +123,16 @@ SUBROUTINE plot_vtu(p_mesh, c_filename, &
INTEGER( KIND = GRID_SI) :: i_numberofcells
INTEGER( KIND = GRID_SI) :: i_numberoffaces
INTEGER( KIND = GRID_SI) :: i_numberofpoints
#ifdef VTU_OUTPUTSTC
INTEGER( KIND = GRID_SI) :: i_numberofvtupoints
INTEGER( KIND = GRID_SI) :: i_numberofvtufaces
INTEGER( KIND = GRID_SI) :: i_numberofvtucells
REAL ( KIND = GRID_SR ) :: r_heightfactor
TYPE(t_vtu_data), DIMENSION(6) :: p_stcnodedata
REAL (KIND = GRID_SR), DIMENSION(:,:), POINTER :: r_stcval
INTEGER( KIND = GRID_SI) :: i_stcbeg
INTEGER( KIND = GRID_SI) :: i_stcend
#endif
INTEGER (KIND = GRID_SI), DIMENSION(:,:), ALLOCATABLE :: i_cellnodes
REAL (KIND = GRID_SR), DIMENSION(:,:), ALLOCATABLE :: r_nodecoor
......@@ -165,15 +176,17 @@ SUBROUTINE plot_vtu(p_mesh, c_filename, &
END IF
! Set dimension depending types
#if defined VTU_OUTPUT2D || defined VTU_OUTPUTSPH
#if defined VTU_OUTPUT2D || defined VTU_OUTPUTSPH || defined VTU_OUTPUTSTC
i_numberofcells = p_mesh%i_enumfine
i_numberoffaces = p_mesh%i_gnumfine
#elif defined VTU_OUTPUT3D
i_numberofcells = p_mesh%i_tnumfine
i_numberoffaces = p_mesh%i_enumfine
#elif defined VTU_OUTPUTSTC
i_numberofcells = p_mesh%i_enumfine* (i_lay-1)
i_numberoffaces = p_mesh%i_gnumfine* i_lay + p_mesh%i_nnumber* (i_lay-1)
#endif
#if defined VTU_OUTPUTSTC
i_numberofvtucells = p_mesh%i_enumfine* (i_lay-1)
i_numberoffvtuaces = p_mesh%i_gnumfine* i_lay + p_mesh%i_nnumber* (i_lay-1)
#endif
IF(l_continuous_data) THEN
......@@ -181,6 +194,13 @@ SUBROUTINE plot_vtu(p_mesh, c_filename, &
ELSE
i_numberofpoints = i_numberofcells * i_nodespercell
END IF
#ifdef VTU_OUTPUTSTC
IF(l_continuous_data) THEN
i_numberofvtupoints = p_mesh%i_nnumber* i_lay
ELSE
i_numberofvtupoints = i_numberofvtucells * i_nodespercell
END IF
#endif
! We only need the node coordinates and the nodes per tetra
! as additional information.
......@@ -193,7 +213,7 @@ SUBROUTINE plot_vtu(p_mesh, c_filename, &
END IF
CALL GRID_getinfo(p_mesh, r_nodecoordinates = r_nodecoor, &
#if defined VTU_OUTPUT2D || defined VTU_OUTPUTSPH
#if defined VTU_OUTPUT2D || defined VTU_OUTPUTSPH || defined VTU_OUTPUTSTC
i_elementnodes = i_cellnodes)
#elif defined VTU_OUTPUT3D
i_tetranodes = i_cellnodes)
......@@ -208,11 +228,19 @@ SUBROUTINE plot_vtu(p_mesh, c_filename, &
WRITE(i_fhandle, "(A)") '<?xml version="1.0"?>'
WRITE(i_fhandle, "(A)") '<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">'
WRITE(i_fhandle, *) '<UnstructuredGrid>'
#if defined VTU_OUTPUT2D || defined VTU_OUTPUTSPH || defined VTU_OUTPUT3D
WRITE(i_fhandle, "(A,I8,A,I8,A)") '<Piece NumberOfPoints="', &
i_numberofpoints, &
'" NumberOfCells="', &
i_numberofcells, &
'">'
#elif defined VTU_OUTPUTSTC
WRITE(i_fhandle, "(A,I8,A,I8,A)") '<Piece NumberOfPoints="', &
i_numberofvtupoints, &
'" NumberOfCells="', &
i_numberofvtucells, &
'">'
#endif
! write the node coordinates to VTU file
WRITE(i_fhandle, *) '<Points>'
......@@ -222,10 +250,12 @@ SUBROUTINE plot_vtu(p_mesh, c_filename, &
! Distunguish between continuous and discontinuous case
IF(l_continuous_data) THEN
!> The continuous case is represented by amatos' grid topology.
DO i_cnt = 1, p_mesh%i_nnumber
#if defined VTU_OUTPUT3D || defined VTU_OUTPUTSPH
DO i_cnt = 1, p_mesh%i_nnumber
WRITE(i_fhandle, *) r_nodecoor(:, i_cnt)
END DO
#elif defined VTU_OUTPUT2D
DO i_cnt = 1, p_mesh%i_nnumber
IF(PRESENT(r_zcoordinate)) THEN
WRITE(i_fhandle, *) r_nodecoor(:, i_cnt), r_zcoordinate(1, i_cnt)
ELSE IF(PRESENT(i_zcoordinate)) THEN
......@@ -234,17 +264,28 @@ SUBROUTINE plot_vtu(p_mesh, c_filename, &
ELSE
WRITE(i_fhandle, *) r_nodecoor(:, i_cnt), 0.0
END IF
#endif
END DO
#elif defined VTU_OUTPUTSTC
DO i_hcnt = 1,i_lay
r_heightfactor = 1._GRID_SR + (1000. - r_height(i_hcnt))* 0.001
DO i_cnt = 1, p_mesh%i_nnumber
WRITE(i_fhandle, *) r_nodecoor(:, i_cnt)*r_heightfactor
END DO
END DO
#endif
ELSE
! Write the node coordinates. In this discontinuous case one node in
! paraview belongs only to one cell.
#if defined VTU_OUTPUT2D || defined VTU_OUTPUTSPH
DO i_cnt = 1, i_numberofcells
DO i_ncnt = 1, i_nodespercell
#if defined VTU_OUTPUT2D || defined VTU_OUTPUTSPH
WRITE(i_fhandle, *) r_nodecoor(:, i_cellnodes(i_ncnt, i_cnt))
END DO
END DO
#elif defined VTU_OUTPUT2D
DO i_cnt = 1, i_numberofcells
DO i_ncnt = 1, i_nodespercell
IF(PRESENT(r_zcoordinate)) THEN
WRITE(i_fhandle, *) r_nodecoor(:, i_cellnodes(i_ncnt, i_cnt)), &
r_zcoordinate(1, i_cellnodes(i_ncnt, i_cnt))
......@@ -254,9 +295,11 @@ SUBROUTINE plot_vtu(p_mesh, c_filename, &
ELSE
WRITE(i_fhandle, *) r_nodecoor(:, i_cellnodes(i_ncnt, i_cnt)), 0.0
END IF
#endif
END DO
END DO
#elif defined VTU_OUTPUTSTC
CALL GRID_error(c_error='[plot_vtu] non-consecutive plotting not supported in stacked mode')
#endif
END IF
WRITE(i_fhandle, *) '</DataArray>'
......@@ -270,14 +313,20 @@ SUBROUTINE plot_vtu(p_mesh, c_filename, &
IF(l_continuous_data) THEN
! The indexing in paraview starts with 0. Therefore we subtract one from
! the node connectivity information.
#if defined VTU_OUTPUT2D || defined VTU_OUTPUTSPH
DO i_cnt = 1, p_mesh%i_enumfine
#elif defined VTU_OUTPUT3D
DO i_cnt = 1, p_mesh%i_tnumfine
#endif
#if defined VTU_OUTPUT2D || defined VTU_OUTPUTSPH || defined VTU_OUTPUT2D
DO i_cnt = 1, i_numberofcells
WRITE(i_fhandle, *) i_cellnodes(:, i_cnt) - 1
END DO
#elif defined VTU_OUTPUTSTC
DO i_cnt = 1, i_numberofcells
DO i_hcnt = 1, i_lay-1
WRITE(i_fhandle, *) (i_cellnodes(:, i_cnt)+(p_mesh%i_nnumber*(i_hcnt-1)) - 1), &
(i_cellnodes(:, i_cnt)+(p_mesh%i_nnumber*i_hcnt) - 1)
END DO
END DO
#endif
ELSE
#if defined VTU_OUTPUT2D || defined VTU_OUTPUTSPH || defined VTU_OUTPUT2D
! Write element nodes (indexing starts by 0!)
! In the discontinuous case each tetrahedron has its own nodes
ALLOCATE(i_nodes(i_nodespercell), stat = i_alct)
......@@ -296,6 +345,9 @@ SUBROUTINE plot_vtu(p_mesh, c_filename, &
END DO
DEALLOCATE(i_nodes)
#elif defined VTU_OUTPUTSTC
CALL GRID_error(c_error='[plot_vtu] non-consecutive plotting not supported in stacked mode')
#endif
END IF
WRITE(i_fhandle, *) '</DataArray>'
......@@ -331,9 +383,28 @@ SUBROUTINE plot_vtu(p_mesh, c_filename, &
! When variables for the node data are present, save them.
IF(PRESENT(i_nodedata) .AND. PRESENT(p_nodedata)) THEN
#if defined VTU_OUTPUT2D || defined VTU_OUTPUTSPH || defined VTU_OUTPUT2D
DO i_cnt = 1, i_nodedata
CALL write_vtu_data(i_fhandle, p_nodedata(i_cnt))
END DO
#elif defined VTU_OUTPUTSTC
ALLOCATE(r_stcval(1,i_numberofvtupoints), stat=i_alct)
IF(i_alct /= 0) THEN
CALL grid_error(c_error='[plot_vtu]: could not allocate aux. stc values arrays')
END IF
DO i_cnt = 1, i_nodedata
IF (p_nodedata(i_cnt)%c_name(1:5) == 'layer') THEN
READ(i_hcnt,'I8') p_nodedata(i_cnt)%c_name(6:32)
i_stcbeg = p_mesh%i_nnumber*(i_hcnt-1)+1
i_stcend = p_mesh%i_nnumber*i_hcnt
r_stcval(1:1,i_stcbeg:i_stcend) = p_nodedata(i_cnt)%p_vdata(1:1,:)
END IF
END DO
p_stcnodedata(1)%c_name = 'stc-tracer'
p_stcnodedata(1)%i_size = 1
p_stcnodedata(1)%p_vdata=>r_stcval
CALL write_vtu_data(i_fhandle, p_stcnodedata(1))
#endif
END IF
WRITE(i_fhandle, *) '</PointData>'
......@@ -355,11 +426,31 @@ SUBROUTINE plot_vtu(p_mesh, c_filename, &
END IF
! Write the variable belonging to the tetradra.
#if defined VTU_OUTPUT2D || defined VTU_OUTPUTSPH || defined VTU_OUTPUT2D
IF(PRESENT(i_celldata) .AND. PRESENT(p_celldata)) THEN
DO i_cnt = 1, i_celldata
CALL write_vtu_data(i_fhandle, p_celldata(i_cnt))
END DO
END IF
#elif defined VTU_OUTPUTSTC
! at this time we assign the same cell values to the whole column of stacked cells
IF(ALLOCATED(r_stcval)) DEALLOCATE(r_stcval)
ALLOCATE(r_stcval(1,i_numberofvtucells), stat=i_alct)
IF(i_alct /= 0) THEN
CALL grid_error(c_error='[plot_vtu]: could not allocate aux. stc values arrays')
END IF
DO i_cnt = 1, i_celldata
DO i_hcnt=1,i_numberofcells
i_stcbeg = (i_hcnt-1)*(i_lay-1)+1
i_stcend = i_hcnt*(i_lay-1)
r_stcval(i_stcbeg:i_stcend) = p_celldata(i_cnt)%p_vdata(i_hct)
END DO
p_stcnodedata(1)%c_name = p_celldata(i_cnt)%c_name
p_stcnodedata(1)%i_size = p_celldata(i_cnt)%i_size
p_stcnodedata(1)%p_vdata=>r_stcval
CALL write_vtu_data(i_fhandle, p_celldata(i_cnt))
END DO
#endif
! write the end of cell data and the footer.
WRITE(i_fhandle, *) '</CellData>'
......@@ -370,6 +461,9 @@ SUBROUTINE plot_vtu(p_mesh, c_filename, &
! tidy up
CLOSE(i_fhandle)
DEALLOCATE(r_nodecoor, i_cellnodes)
#if defined VTU_OUTPUTSTC
IF(ALLOCATED(r_stcval)) DEALLOCATE(r_stcval)
#endif
END SUBROUTINE plot_vtu
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment