Xdmf

From Core Physics Wiki
Jump to: navigation, search

I was recently looking to use XDMF to visualize some HDF5 data and I had trouble finding some examples in Fortran.

The other source of examples I found is to download the XDMF library from Kitware XDMF Page


2D curvilinear mesh

This is a Fortran example which generatse an HDF file with 2D curvilinear mesh data, and an XMF file to tell Visit how to plot the HDF data. This example is a translation of the C++ Program at VisitUsers

!=======================================================================
!
!  Sample Fortran program to demonstrate using XDMF to read HDF files
!
!  This program is a Fortran version of the C++ program given at:
!  http://visitusers.org/index.php?title=Using_XDMF_to_read_HDF5
!
!  This program must be compiled with the HDF5 libraries
!
!=======================================================================
 
      program makexdmf
      implicit none
 
      integer :: NX=30
      integer :: NY=20
 
!  write sample HDF file
 
      call write_hdf5_data(nx,ny)
 
!  write corresponding XDMF file
 
      call write_xdmf_xml(nx,ny)
 
      end
 
!=======================================================================
 
      subroutine write_hdf5_data(nx,ny)
      use hdf5
      implicit none
 
      integer, intent(in) :: nx, ny
 
      integer :: i, j
      integer :: ndx
 
      double precision :: xcoord((nx+1)*(ny+1))      ! automatic arrays
      double precision :: ycoord((nx+1)*(ny+1))      ! automatic arrays
      double precision :: pressure(nx*ny)            ! automatic arrays
      double precision :: velocityx((nx+1)*(ny+1))   ! automatic arrays
 
      double precision :: m_pi=3.1415926535897932384d0
 
      double precision :: xt, yt, angle, r
 
      integer(hsize_t) :: idims(3)
 
      integer(hid_t)   :: file_id
      integer(hid_t)   :: dataset_id
      integer(hid_t)   :: dataspace_id
      integer(hid_t)   :: dtype_id      ! data type to be written
      integer          :: ierror
 
      character(len=1) :: xname='X'
      character(len=1) :: yname='Y'
 
!--- Initialize fortran interface.
 
      call h5open_f(ierror)
      if (ierror.ne.0) stop 'ierror: h5open_f'
 
!--- create file
 
      call h5fcreate_f("xdmf2d.h5", H5F_ACC_TRUNC_F, file_id, ierror)
      if (ierror.ne.0) stop 'could not create file'
 
!--- Create the coordinate data.
 
      ndx=0
      do j=0, ny
        yt = dble(j) / dble(NY)
        angle = yt * M_PI
 
        do i=0, nx
 
           xt = dble(i) / dble(NX)
           R = (1.0d0-xt)*2.0d0 + xt*5.0d0
 
           ndx=ndx+1
           xcoord(ndx) = R * cos(angle)
           ycoord(ndx) = R * sin(angle)
!!         write (*,*) ndx, i, j, r, angle, xcoord(ndx), ycoord(ndx)
        enddo
      enddo
 
!--- Create the scalar data.
 
      do j=0, ny-1
        do i=0, nx-1
          ndx = j*nx + i + 1
          pressure(ndx) = dble(j)
        enddo
      enddo
 
      do j=0, ny
        do i=0, nx
           ndx = j * (NX+1) + i + 1
           velocityx(ndx) = dble(i)
        enddo
      enddo
 
!--- Write the data file
 
      dtype_id=H5T_NATIVE_DOUBLE  !  FLOAT
 
!  write x coordinates
 
        write (0,*) 'writing x'
 
        idims(1) = (NY + 1)
        idims(2) = (NX + 1)
        idims(3) = 0
 
        call h5screate_simple_f(2, idims, dataspace_id, ierror)
 
        call h5dcreate_f(file_id, xname, dtype_id, dataspace_id, dataset_id, ierror)
 
        call h5dwrite_f(dataset_id, dtype_id, xcoord, idims, ierror)
 
        call h5dclose_f(dataset_id, ierror)
        call h5sclose_f(dataspace_id, ierror)
 
!  write y coordinates
 
        write (0,*) 'writing y'
 
        idims(1) = (NY + 1)
        idims(2) = (NX + 1)
        idims(3) = 0
 
        call h5screate_simple_f(2, idims, dataspace_id, ierror)
 
        call h5dcreate_f(file_id, yname, dtype_id, dataspace_id, dataset_id, ierror)
 
        call h5dwrite_f(dataset_id, dtype_id, ycoord, idims, ierror)
 
        call h5dclose_f(dataset_id, ierror)
        call h5sclose_f(dataspace_id, ierror)
 
 
!---  Write the scalar data.
 
! write presssure
 
      write (0,*) 'writing pressure'
      idims(1) = NY
      idims(2) = NX
      idims(3) = 0
 
      call h5screate_simple_f(2, idims, dataspace_id, ierror)
 
      call h5dcreate_f(file_id, "/Pressure", dtype_id, dataspace_id, dataset_id, ierror)
 
      call h5dwrite_f(dataset_id, dtype_id, pressure, idims, ierror)
 
      call h5dclose_f(dataset_id, ierror)
      call h5sclose_f(dataspace_id, ierror)
 
! write velocity
 
      write (0,*) 'writing velocity'
      idims(1) = NY + 1
      idims(2) = NX + 1
      idims(3) = 0
      call h5screate_simple_f(2, idims, dataspace_id, ierror)
 
      call h5dcreate_f(file_id, "/VelocityX", dtype_id, dataspace_id, dataset_id, ierror)
 
      call h5dwrite_f(dataset_id, dtype_id, velocityx, idims, ierror)
 
      call h5dclose_f(dataset_id, ierror)
      call h5sclose_f(dataspace_id, ierror)
 
!--- close HDF file
 
      write (0,*) 'closing HDF file'
 
      call h5fclose_f(file_id, ierror)
 
      return
      end subroutine
 
!=======================================================================
      subroutine write_xdmf_xml(nx,ny)
      implicit none
 
      integer, intent(in) :: nx, ny
 
      integer :: xmf=12
 
!     Open the file and write the XML description of the mesh..
 
      write (0,*) 'writing XMF file'
 
      open (xmf, file="xdmf2d.xmf")
      write (xmf,20) '<?xml version="1.0" ?>'
      write (xmf,20) '<!DOCTYPE Xdmf SYSTEM "Xdmf.dtd" []>'
      write (xmf,20) '<Xdmf Version="2.0">'
      write (xmf,20) ' <Domain>'
      write (xmf,20) '   <Grid Name="mesh1" GridType="Uniform">'
      write (xmf,33) '     <Topology TopologyType="2DSMesh" NumberOfElements="',ny+1,nx+1,'"/>'
      write (xmf,20) '     <Geometry GeometryType="X_Y">'
      write (xmf,33) '       <DataItem Dimensions="',ny+1, nx+1,'" NumberType="Float" Precision="4" Format="HDF">'
      write (xmf,20) '        xdmf2d.h5:/X'
      write (xmf,20) '       </DataItem>'
      write (xmf,33) '       <DataItem Dimensions="',ny+1, nx+1,'" NumberType="Float" Precision="4" Format="HDF">'
      write (xmf,20) '        xdmf2d.h5:/Y'
      write (xmf,20) '       </DataItem>'
      write (xmf,20) '     </Geometry>'
      write (xmf,20) '     <Attribute Name="Pressure" AttributeType="Scalar" Center="Cell">'
      write (xmf,33) '       <DataItem Dimensions="',ny,nx,'" NumberType="Float" Precision="4" Format="HDF">'
      write (xmf,20) '        xdmf2d.h5:/Pressure'
      write (xmf,20) '       </DataItem>'
      write (xmf,20) '     </Attribute>'
      write (xmf,20) '     <Attribute Name="VelocityX" AttributeType="Scalar" Center="Node">'
      write (xmf,33) '       <DataItem Dimensions="',ny+1,nx+1,'" NumberType="Float" Precision="4" Format="HDF">'
      write (xmf,20) '        xdmf2d.h5:/VelocityX'
      write (xmf,20) '       </DataItem>'
      write (xmf,20) '     </Attribute>'
      write (xmf,20) '   </Grid>'
      write (xmf,20) ' </Domain>'
      write (xmf,20) '</Xdmf>'
  20  format (a)
  33  format (a,i0,1x,i0,a)
      close(xmf)
 
      return
      end subroutine write_xdmf_xml
!=======================================================================