subroutine unstructuredvtkB64(qunit)
! output for vtu format to paraview, binary version output
! not parallel, uses calc_grid to compute nwauxio variables
! allows renormalizing using convert factors
use mod_forest, only: Morton_start, Morton_stop, sfc_to_igrid
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 :: normconv(0:nw+nwauxio)
integer, allocatable :: intstatus(:,:)
integer :: itag,ipe,igrid,level,icel,ixCmin1,ixCmin2,ixCmin3,ixCmax1,&
ixCmax2,ixCmax3,ixCCmin1,ixCCmin2,ixCCmin3,ixCCmax1,ixCCmax2,ixCCmax3,&
Morton_no,Morton_length
integer :: nx1,nx2,nx3,nxC1,nxC2,nxC3,nc,np,VTK_type,ix1,ix2,ix3,filenr
integer*8 :: offset
integer:: k,iw
integer:: length,lengthcc,length_coords,length_conn,length_offsets
character:: buf
character(len=80):: filename
character(len=name_len) :: wnamei(1:nw+nwauxio),&
xandwnamei(1:ndim+nw+nwauxio)
character(len=1024) :: outfilehead
logical :: fileopen,cell_corner=.false.
logical, allocatable :: Morton_aim(:),Morton_aim_p(:)
normconv=one
Morton_length=Morton_stop(npe-1)-Morton_start(0)+1
allocate(Morton_aim(Morton_start(0):Morton_stop(npe-1)))
allocate(Morton_aim_p(Morton_start(0):Morton_stop(npe-1)))
Morton_aim=.false.
Morton_aim_p=.false.
do Morton_no=Morton_start(mype),Morton_stop(mype)
igrid=sfc_to_igrid(Morton_no)
level=node(plevel_,igrid)
! we can clip parts of the grid away, select variables, levels etc.
if(writelevel(level)) then
! 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
Morton_aim_p(Morton_no)=.true.
end if
end if
end do
call MPI_ALLREDUCE(Morton_aim_p,Morton_aim,Morton_length,MPI_LOGICAL,&
MPI_LOR,icomm,ierrmpi)
select case(convert_type)
case('vtuB64','vtuBmpi64')
cell_corner=.true.
case('vtuBCC64','vtuBCCmpi64')
cell_corner=.false.
end select
if (mype /= 0) then
do Morton_no=Morton_start(mype),Morton_stop(mype)
if(.not. Morton_aim(Morton_no)) cycle
igrid=sfc_to_igrid(Morton_no)
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.)
itag=Morton_no
call mpi_send_wrapper(xC_TMP,1,type_block_xc_io, 0,itag,icomm,ierrmpi)
if(cell_corner) then
call mpi_send_wrapper(wC_TMP,1,type_block_wc_io, 0,itag,icomm,ierrmpi)
else
call mpi_send_wrapper(wCC_TMP,1,type_block_wcc_io, 0,itag,icomm,ierrmpi)
end if
end do
else
! mype==0
offset=0
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='replace')
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,*) real(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
length=np*size_double
lengthcc=nc*size_double
length_coords=3*length
length_conn=2**3*size_int*nc
length_offsets=nc*size_int
! Note: using the w_write, writelevel, writespshift
do Morton_no=Morton_start(0),Morton_stop(0)
if(.not. Morton_aim(Morton_no)) cycle
if(cell_corner) then
! 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,i16,a)')'<DataArray type="Float64" Name="',&
TRIM(wnamei(iw)), '" format="appended" offset="',offset,'">'
write(qunit,'(a)')'</DataArray>'
offset=offset+length+size_int
end do
write(qunit,'(a)')'</PointData>'
write(qunit,'(a)')'<Points>'
write(qunit,'(a,i16,a)')&
'<DataArray type="Float64" NumberOfComponents="3" format="appended" offset="',&
offset,'"/>'
! write cell corner coordinates in a backward dimensional loop, always 3D output
offset=offset+length_coords+size_int
write(qunit,'(a)')'</Points>'
else
! 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,i16,a)')'<DataArray type="Float64" Name="',&
TRIM(wnamei(iw)), '" format="appended" offset="',offset,'">'
write(qunit,'(a)')'</DataArray>'
offset=offset+lengthcc+size_int
end do
write(qunit,'(a)')'</CellData>'
write(qunit,'(a)')'<Points>'
write(qunit,'(a,i16,a)')&
'<DataArray type="Float64" NumberOfComponents="3" format="appended" offset="',&
offset,'"/>'
! write cell corner coordinates in a backward dimensional loop, always 3D output
offset=offset+length_coords+size_int
write(qunit,'(a)')'</Points>'
end if
write(qunit,'(a)')'<Cells>'
! connectivity part
write(qunit,'(a,i16,a)'&
)'<DataArray type="Int32" Name="connectivity" format="appended" offset="',&
offset,'"/>'
offset=offset+length_conn+size_int
! offsets data array
write(qunit,'(a,i16,a)')&
'<DataArray type="Int32" Name="offsets" format="appended" offset="',offset,&
'"/>'
offset=offset+length_offsets+size_int
! VTK cell type data array
write(qunit,'(a,i16,a)')&
'<DataArray type="Int32" Name="types" format="appended" offset="',&
offset,'"/>'
offset=offset+size_int+nc*size_int
write(qunit,'(a)')'</Cells>'
write(qunit,'(a)')'</Piece>'
end do
! write metadata communicated from other processors
if(npe>1)then
do ipe=1, npe-1
do Morton_no=Morton_start(ipe),Morton_stop(ipe)
if(.not. Morton_aim(Morton_no)) cycle
if(cell_corner) then
! 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,i16,a)')'<DataArray type="Float64" Name="',&
TRIM(wnamei(iw)), '" format="appended" offset="',offset,&
'">'
write(qunit,'(a)')'</DataArray>'
offset=offset+length+size_int
end do
write(qunit,'(a)')'</PointData>'
write(qunit,'(a)')'<Points>'
write(qunit,'(a,i16,a)')&
'<DataArray type="Float64" NumberOfComponents="3" format="appended" offset="',&
offset,'"/>'
! write cell corner coordinates in a backward dimensional loop, always 3D output
offset=offset+length_coords+size_int
write(qunit,'(a)')'</Points>'
else
! 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,i16,a)')'<DataArray type="Float64" Name="',&
TRIM(wnamei(iw)), '" format="appended" offset="',offset,&
'">'
write(qunit,'(a)')'</DataArray>'
offset=offset+lengthcc+size_int
end do
write(qunit,'(a)')'</CellData>'
write(qunit,'(a)')'<Points>'
write(qunit,'(a,i16,a)')&
'<DataArray type="Float64" NumberOfComponents="3" format="appended" offset="',&
offset,'"/>'
! write cell corner coordinates in a backward dimensional loop, always 3D output
offset=offset+length_coords+size_int
write(qunit,'(a)')'</Points>'
end if
write(qunit,'(a)')'<Cells>'
! connectivity part
write(qunit,'(a,i16,a)'&
)'<DataArray type="Int32" Name="connectivity" format="appended" offset="',&
offset,'"/>'
offset=offset+length_conn+size_int
! offsets data array
write(qunit,'(a,i16,a)')&
'<DataArray type="Int32" Name="offsets" format="appended" offset="',offset,&
'"/>'
offset=offset+length_offsets+size_int
! VTK cell type data array
write(qunit,'(a,i16,a)')&
'<DataArray type="Int32" Name="types" format="appended" offset="',offset,&
'"/>'
offset=offset+size_int+nc*size_int
write(qunit,'(a)')'</Cells>'
write(qunit,'(a)')'</Piece>'
end do
end do
end if
write(qunit,'(a)')'</UnstructuredGrid>'
write(qunit,'(a)')'<AppendedData encoding="raw">'
close(qunit)
open(qunit,file=filename,access='stream',form='unformatted',&
position='append')
buf='_'
write(qunit) TRIM(buf)
do Morton_no=Morton_start(0),Morton_stop(0)
if(.not. Morton_aim(Morton_no)) cycle
igrid=sfc_to_igrid(Morton_no)
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.)
do iw=1,nw+nwauxio
if(iw<=nw) then
if(.not.w_write(iw)) cycle
end if
if(cell_corner) then
write(qunit) length
write(qunit) (((wC_TMP(ix1,ix2,ix3,iw)*normconv(iw),ix1=ixCmin1,&
ixCmax1),ix2=ixCmin2,ixCmax2),ix3=ixCmin3,ixCmax3)
else
write(qunit) lengthcc
write(qunit) (((wCC_TMP(ix1,ix2,ix3,iw)*normconv(iw),ix1=ixCCmin1,&
ixCCmax1),ix2=ixCCmin2,ixCCmax2),ix3=ixCCmin3,ixCCmax3)
end if
end do
write(qunit) length_coords
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);
do k=1,3
write(qunit) x_VTK(k)
end do
end do
end do
end do
write(qunit) length_conn
do ix3=1,nx3
do ix2=1,nx2
do ix1=1,nx1
write(qunit)(ix3-1)*nxC2*nxC1+(ix2-1)*nxC1+ix1-1,&
(ix3-1)*nxC2*nxC1+(ix2-1)*nxC1+ix1,&
(ix3-1)*nxC2*nxC1+ ix2*nxC1+ix1-1,&
(ix3-1)*nxC2*nxC1+ ix2*nxC1+ix1,&
ix3*nxC2*nxC1+(ix2-1)*nxC1+ix1-1,ix3*nxC2*nxC1+(ix2-1)*nxC1+ix1,&
ix3*nxC2*nxC1+ ix2*nxC1+ix1-1,ix3*nxC2*nxC1+ ix2*nxC1+ix1
end do
end do
end do
write(qunit) length_offsets
do icel=1,nc
write(qunit) icel*(2**3)
end do
VTK_type=11
write(qunit) size_int*nc
do icel=1,nc
write(qunit) VTK_type
end do
end do
allocate(intstatus(MPI_STATUS_SIZE,1))
if(npe>1)then
ixCCmin1=ixMlo1;ixCCmin2=ixMlo2;ixCCmin3=ixMlo3; ixCCmax1=ixMhi1
ixCCmax2=ixMhi2;ixCCmax3=ixMhi3;
ixCmin1=ixMlo1-1;ixCmin2=ixMlo2-1;ixCmin3=ixMlo3-1; ixCmax1=ixMhi1
ixCmax2=ixMhi2;ixCmax3=ixMhi3;
do ipe=1, npe-1
do Morton_no=Morton_start(ipe),Morton_stop(ipe)
if(.not. Morton_aim(Morton_no)) cycle
itag=Morton_no
call mpi_recv_wrapper(xC_TMP,1,type_block_xc_io, ipe,itag,icomm,&
intstatus(:,1),ierrmpi)
if(cell_corner) then
call mpi_recv_wrapper(wC_TMP,1,type_block_wc_io, ipe,itag,icomm,&
intstatus(:,1),ierrmpi)
else
call mpi_recv_wrapper(wCC_TMP,1,type_block_wcc_io, ipe,itag,icomm,&
intstatus(:,1),ierrmpi)
end if
do iw=1,nw+nwauxio
if(iw<=nw) then
if(.not.w_write(iw)) cycle
end if
if(cell_corner) then
write(qunit) length
write(qunit) (((wC_TMP(ix1,ix2,ix3,iw)*normconv(iw),&
ix1=ixCmin1,ixCmax1),ix2=ixCmin2,ixCmax2),ix3=ixCmin3,&
ixCmax3)
else
write(qunit) lengthcc
write(qunit) (((wCC_TMP(ix1,ix2,ix3,iw)*normconv(iw),&
ix1=ixCCmin1,ixCCmax1),ix2=ixCCmin2,ixCCmax2),ix3=ixCCmin3,&
ixCCmax3)
end if
end do
write(qunit) length_coords
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);
do k=1,3
write(qunit) x_VTK(k)
end do
end do
end do
end do
write(qunit) length_conn
do ix3=1,nx3
do ix2=1,nx2
do ix1=1,nx1
write(qunit)(ix3-1)*nxC2*nxC1+(ix2-1)*nxC1+ix1-1,&
(ix3-1)*nxC2*nxC1+(ix2-1)*nxC1+ix1,&
(ix3-1)*nxC2*nxC1+ ix2*nxC1+ix1-1,&
(ix3-1)*nxC2*nxC1+ ix2*nxC1+ix1,&
ix3*nxC2*nxC1+(ix2-1)*nxC1+ix1-1,&
ix3*nxC2*nxC1+(ix2-1)*nxC1+ix1,&
ix3*nxC2*nxC1+ ix2*nxC1+ix1-1,&
ix3*nxC2*nxC1+ ix2*nxC1+ix1
end do
end do
end do
write(qunit) length_offsets
do icel=1,nc
write(qunit) icel*(2**3)
end do
VTK_type=11
write(qunit) size_int*nc
do icel=1,nc
write(qunit) VTK_type
end do
end do
end do
end if
close(qunit)
open(qunit,file=filename,status='unknown',form='formatted',&
position='append')
write(qunit,'(a)')'</AppendedData>'
write(qunit,'(a)')'</VTKFile>'
close(qunit)
deallocate(intstatus)
end if
deallocate(Morton_aim,Morton_aim_p)
if (npe>1) then
call MPI_BARRIER(icomm,ierrmpi)
end if
end subroutine unstructuredvtkB64