subroutine unstructuredvtk(qunit)
! output for vtu format to paraview
! not parallel, uses calc_grid to compute nwauxio variables
! allows renormalizing using convert factors
use mod_global_parameters
use mod_calculate_xw
integer, intent(in) :: qunit
double precision :: x_VTK(1:3)
double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,&
ixMlo3-1:ixMhi3,ndim) :: xC_TMP
double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo3:ixMhi3,&
ndim) :: xCC_TMP
double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,&
ixMlo3-1:ixMhi3,ndim) :: xC
double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo3:ixMhi3,&
ndim) :: xCC
double precision, dimension(ixMlo1-1:ixMhi1,ixMlo2-1:ixMhi2,&
ixMlo3-1:ixMhi3,nw+nwauxio) :: wC_TMP
double precision, dimension(ixMlo1:ixMhi1,ixMlo2:ixMhi2,ixMlo3:ixMhi3,&
nw+nwauxio) :: wCC_TMP
double precision, dimension(0:nw+nwauxio) :: normconv
integer:: igrid,iigrid,level,igonlevel,icel,ixCmin1,ixCmin2,&
ixCmin3,ixCmax1,ixCmax2,ixCmax3,ixCCmin1,ixCCmin2,ixCCmin3,ixCCmax1,&
ixCCmax2,ixCCmax3,iw
integer:: NumGridsOnLevel(1:nlevelshi)
integer :: nx1,nx2,nx3,nxC1,nxC2,nxC3,nodesonlevel,elemsonlevel,nc,np,&
VTK_type,ix1,ix2,ix3
character(len=80):: filename
integer :: filenr
character(len=name_len) :: wnamei(1:nw+nwauxio),&
xandwnamei(1:ndim+nw+nwauxio)
character(len=1024) :: outfilehead
logical :: fileopen
if(npe>1)then
if(mype==0) PRINT *,'unstructuredvtk not parallel, use vtumpi'
call mpistop('npe>1, unstructuredvtk')
end if
inquire(qunit,opened=fileopen)
if(.not.fileopen)then
! generate filename
filenr=snapshotini
if (autoconvert) filenr=snapshotnext
write(filename,'(a,i4.4,a)') TRIM(base_filename),filenr,".vtu"
! Open the file for the header part
open(qunit,file=filename,status='unknown')
end if
call getheadernames(wnamei,xandwnamei,outfilehead)
! generate xml header
write(qunit,'(a)')'<?xml version="1.0"?>'
write(qunit,'(a)',advance='no') '<VTKFile type="UnstructuredGrid"'
write(qunit,'(a)')' version="0.1" byte_order="LittleEndian">'
write(qunit,'(a)')'<UnstructuredGrid>'
write(qunit,'(a)')'<FieldData>'
write(qunit,'(2a)')'<DataArray type="Float32" Name="TIME" ',&
'NumberOfTuples="1" format="ascii">'
write(qunit,*) dble(dble(global_time)*time_convert_factor)
write(qunit,'(a)')'</DataArray>'
write(qunit,'(a)')'</FieldData>'
! number of cells, number of corner points, per grid.
nx1=ixMhi1-ixMlo1+1;nx2=ixMhi2-ixMlo2+1;nx3=ixMhi3-ixMlo3+1;
nxC1=nx1+1;nxC2=nx2+1;nxC3=nx3+1;
nc=nx1*nx2*nx3
np=nxC1*nxC2*nxC3
! Note: using the w_write, writelevel, writespshift
! we can clip parts of the grid away, select variables, levels etc.
do level=levmin,levmax
if (writelevel(level)) then
do iigrid=1,igridstail; igrid=igrids(iigrid);
if (node(plevel_,igrid)/=level) cycle
block=>ps(igrid)
! only output a grid when fully within clipped region selected
! by writespshift array
if((rnode(rpxmin1_,igrid)>=xprobmin1+&
(xprobmax1-xprobmin1)*writespshift(1,1).and.rnode(rpxmin2_,&
igrid)>=xprobmin2+(xprobmax2-xprobmin2)*writespshift(2,&
1).and.rnode(rpxmin3_,igrid)>=&
xprobmin3+(xprobmax3-xprobmin3)*writespshift(3,&
1)).and.(rnode(rpxmax1_,igrid)<=&
xprobmax1-(xprobmax1-xprobmin1)*writespshift(1,&
2).and.rnode(rpxmax2_,igrid)<=&
xprobmax2-(xprobmax2-xprobmin2)*writespshift(2,&
2).and.rnode(rpxmax3_,igrid)<=&
xprobmax3-(xprobmax3-xprobmin3)*writespshift(3,2))) then
call calc_x(igrid,xC,xCC)
call calc_grid(qunit,igrid,xC,xCC,xC_TMP,xCC_TMP,wC_TMP,wCC_TMP,&
normconv,ixCmin1,ixCmin2,ixCmin3,ixCmax1,ixCmax2,ixCmax3,&
ixCCmin1,ixCCmin2,ixCCmin3,ixCCmax1,ixCCmax2,ixCCmax3,.true.)
select case(convert_type)
case('vtu')
! we write out every grid as one VTK PIECE
write(qunit,'(a,i7,a,i7,a)') '<Piece NumberOfPoints="',np,&
'" NumberOfCells="',nc,'">'
write(qunit,'(a)')'<PointData>'
do iw=1,nw+nwauxio
if(iw<=nw) then
if(.not.w_write(iw)) cycle
end if
write(qunit,'(a,a,a)')'<DataArray type="Float64" Name="',&
TRIM(wnamei(iw)),'" format="ascii">'
write(qunit,'(200(1pe14.6))') (((wC_TMP(ix1,ix2,ix3,&
iw)*normconv(iw),ix1=ixCmin1,ixCmax1),ix2=ixCmin2,ixCmax2),&
ix3=ixCmin3,ixCmax3)
write(qunit,'(a)')'</DataArray>'
end do
write(qunit,'(a)')'</PointData>'
write(qunit,'(a)')'<Points>'
write(qunit,'(a)'&
)'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
! write cell corner coordinates in a backward dimensional loop, always 3D output
do ix3=ixCmin3,ixCmax3
do ix2=ixCmin2,ixCmax2
do ix1=ixCmin1,ixCmax1
x_VTK(1:3)=zero;
x_VTK(1:ndim)=xC_TMP(ix1,ix2,ix3,1:ndim)*normconv(0);
write(qunit,'(3(1pe14.6))') x_VTK
end do
end do
end do
write(qunit,'(a)')'</DataArray>'
write(qunit,'(a)')'</Points>'
case('vtuCC')
! we write out every grid as one VTK PIECE
write(qunit,'(a,i7,a,i7,a)') '<Piece NumberOfPoints="',np,&
'" NumberOfCells="',nc,'">'
write(qunit,'(a)')'<CellData>'
do iw=1,nw+nwauxio
if(iw<=nw) then
if(.not.w_write(iw)) cycle
end if
write(qunit,'(a,a,a)')'<DataArray type="Float64" Name="',&
TRIM(wnamei(iw)),'" format="ascii">'
write(qunit,'(200(1pe14.6))') (((wCC_TMP(ix1,ix2,ix3,&
iw)*normconv(iw),ix1=ixCCmin1,ixCCmax1),ix2=ixCCmin2,&
ixCCmax2),ix3=ixCCmin3,ixCCmax3)
write(qunit,'(a)')'</DataArray>'
end do
write(qunit,'(a)')'</CellData>'
write(qunit,'(a)')'<Points>'
write(qunit,'(a)'&
)'<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
! write cell corner coordinates in a backward dimensional loop, always 3D output
do ix3=ixCmin3,ixCmax3
do ix2=ixCmin2,ixCmax2
do ix1=ixCmin1,ixCmax1
x_VTK(1:3)=zero;
x_VTK(1:ndim)=xC_TMP(ix1,ix2,ix3,1:ndim)*normconv(0);
write(qunit,'(3(1pe14.6))') x_VTK
end do
end do
end do
write(qunit,'(a)')'</DataArray>'
write(qunit,'(a)')'</Points>'
end select
write(qunit,'(a)')'<Cells>'
! connectivity part
write(qunit,'(a)'&
)'<DataArray type="Int32" Name="connectivity" format="ascii">'
call save_connvtk(qunit,igrid)
write(qunit,'(a)')'</DataArray>'
! offsets data array
write(qunit,'(a)'&
)'<DataArray type="Int32" Name="offsets" format="ascii">'
do icel=1,nc
write(qunit,'(i7)') icel*(2**3)
end do
write(qunit,'(a)')'</DataArray>'
! VTK cell type data array
write(qunit,'(a)'&
)'<DataArray type="Int32" Name="types" format="ascii">'
! VTK_LINE=3; VTK_PIXEL=8; VTK_VOXEL=11 -> vtk-syntax
VTK_type=11
do icel=1,nc
write(qunit,'(i2)') VTK_type
end do
write(qunit,'(a)')'</DataArray>'
write(qunit,'(a)')'</Cells>'
write(qunit,'(a)')'</Piece>'
end if
end do
end if
end do
write(qunit,'(a)')'</UnstructuredGrid>'
write(qunit,'(a)')'</VTKFile>'
close(qunit)
end subroutine unstructuredvtk