diff --git a/flash2d/src/options-sphere/IO_vtu_stacked.F90 b/flash2d/src/options-sphere/IO_vtu_stacked.F90
index f2d0481aac47e8267a7dee53a174feaebea8a231..7fdcb7e5f912648b063d88466b52f888327e6be9 100644
--- a/flash2d/src/options-sphere/IO_vtu_stacked.F90
+++ b/flash2d/src/options-sphere/IO_vtu_stacked.F90
@@ -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,29 +250,42 @@ 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
- IF(PRESENT(r_zcoordinate)) THEN
- WRITE(i_fhandle, *) r_nodecoor(:, i_cnt), r_zcoordinate(1, i_cnt)
- ELSE IF(PRESENT(i_zcoordinate)) THEN
- WRITE(i_fhandle, *) r_nodecoor(:, i_cnt), &
- p_nodedata(i_zcoordinate)%p_vdata(1, i_cnt)
- ELSE
- WRITE(i_fhandle, *) r_nodecoor(:, i_cnt), 0.0
- END IF
-#endif
+ 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
+ WRITE(i_fhandle, *) r_nodecoor(:, i_cnt), &
+ p_nodedata(i_zcoordinate)%p_vdata(1, i_cnt)
+ ELSE
+ WRITE(i_fhandle, *) r_nodecoor(:, i_cnt), 0.0
+ END IF
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