static_data_vtk_s.f90 Source File


Contents

Source Code


Source Code

submodule (static_data_m) static_data_vtk_s
    !! Handles built of VTK meshes or visualisation
    implicit none

    logical :: quadratic_hex = .false.
    !! If true a higher order polyhedra are used for visualisation
    !! (Maybe very slow)
    logical :: build_vtk_cano = .true.
    !! If true, builds vtk mesh for canonical mesh
    logical :: build_vtk_stag = .false.
    !! If true, builds vtk mesh for staggered mesh

    namelist / vtk_visualisation / quadratic_hex, build_vtk_cano, build_vtk_stag

contains

    subroutine read_params_vtk(par_filepath)
        character(len=*), intent(in), optional :: par_filepath

        character(len=:), allocatable :: fpath

        integer :: io_unit, io_error
        character(len=256) :: io_errmsg

        if (present(par_filepath)) then
            fpath = par_filepath
        else
            fpath = static_data_parpath
        endif

        open(newunit=io_unit, file=fpath, status='old', action='read', iostat=io_error, iomsg=io_errmsg)
        if (io_error /= 0) then
            call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__)
        endif
        read(io_unit, nml=vtk_visualisation, iostat=io_error, iomsg=io_errmsg)
        if (io_error /= 0) then
            call handle_error(io_errmsg, GRILLIX_ERR_NAMELIST, __LINE__, __FILE__)
        endif
        close(io_unit)

    end subroutine

    subroutine display_params_vtk()

        if (is_rank_info_writer) then
            write(get_stdout(), nml=vtk_visualisation)
        endif

    end subroutine

#ifdef ENABLE_VTK
    ! VTK build routine only available if compiled with -DENABLE_VTK=ON

    module subroutine build_vtk_mesh(comm_handler, par_filepath, dirpath)
        use vtk_fortran, only : vtk_file
        use vis_vtk3d_m, only : write_vtk_mesh

        type(comm_handler_t), intent(in) :: comm_handler
        character(len=*), intent(in), optional :: par_filepath
        character(len=*), intent(in), optional :: dirpath

        logical, dimension(2) :: mask_build
        character(len=5) :: plane_c
        character(len=:), allocatable :: dpath
        character(len=PATHLEN_MAX), dimension(2,2) :: dpath_cr
        character(len=:), allocatable :: fpath
        integer :: cmd_err, cmd_exit, ierr, k, dm, vtk_stat
        real(GP) :: phi_cano, phi_stag, dphi
        type(vtk_file), allocatable :: vtk2d_fid
        type(vtk_file), allocatable :: vtk3d_fid

        if (cntrl_build_vtkmesh == NONE) return

        if ((cntrl_build_vtkmesh == READ) .or. (cntrl_build_vtkmesh == BUILD_NOWRITE))  then
            call handle_error("READ and BUILD_NOWRITE functionality not available for VTK mesh build", &
                              GRILLIX_ERR_STATIC_DATA, __LINE__, __FILE__)
        endif

        if (is_rank_info_writer) then
            write(get_stdout(),*)''
            write(get_stdout(),*)repeat('-',80)
            write(get_stdout(),*)'Building VTK meshes for 3D visualisation'
        endif

        if (present(par_filepath)) then
            call read_params_vtk(par_filepath)
        else
            call read_params_vtk()
        endif

        call display_params_vtk()

        phi_cano = mesh_cano%get_phi()
        phi_stag = mesh_stag%get_phi()
        dphi = map%get_dphi()

        ! Setup directory structure
        if (present(dirpath)) then
            dpath = trim(dirpath)
        else
            dpath = 'vtk'
        endif

        dpath_cr(1,1) = trim(dpath) // '/cano_2d'
        dpath_cr(1,2) = trim(dpath) // '/cano_3d'
        dpath_cr(2,1) = trim(dpath) // '/stag_2d'
        dpath_cr(2,2) = trim(dpath) // '/stag_3d'

        write(plane_c,'(I5.5)')comm_handler%get_plane()

        ! k=1: cano, k=2: stag
        mask_build = [build_vtk_cano, build_vtk_stag]
        do k = 1, 2
            if (.not.mask_build(k)) then
                cycle
            endif

            do dm = 1, 2
                ! Create folder structure (dm=1: 2D, dm=2: 3D)
                if (is_rank_info_writer) then
                    call execute_command_line('mkdir -p ' // trim(dpath_cr(k, dm)), &
                                            exitstat=cmd_exit, cmdstat=cmd_err)
                    if (cmd_err /= 0) then
                        call handle_error("Command execution failed: mkdir -p " // trim(dpath_cr(k, dm)), &
                                        GRILLIX_ERR_CMD, __LINE__, __FILE__, &
                                        error_info_t('cmd_err: ',[cmd_err]))
                    endif
                    if (cmd_exit /= 0) then
                        call handle_error("Command execution returned with warning: mkdir -p" &
                                        // trim(dpath_cr(k, dm)), &
                                        GRILLIX_ERR_CMD, __LINE__, __FILE__, &
                                        error_info_t('cmd_exit: ',[cmd_exit]))
                    endif
                endif
                call MPI_barrier(comm_handler%get_comm_cart(), ierr)
            enddo

            fpath = trim(dpath_cr(k, 1))//'/vtk2d_plane'//plane_c//'.vtu'
            allocate(vtk2d_fid)
            vtk_stat = vtk2d_fid%initialize(format='binary', filename=fpath, mesh_topology='UnstructuredGrid')

            fpath = trim(dpath_cr(k, 2))//'/vtk3d_plane'//plane_c//'.vtu'
            allocate(vtk3d_fid)
            vtk_stat = vtk3d_fid%initialize(format='binary', filename=fpath, mesh_topology='UnstructuredGrid')

            if (k == 1) then
                call write_vtk_mesh(equi, mesh_cano, phi_cano, dphi/2.0_GP, dphi/2.0_GP, &
                                    vtk2d_fid, vtk3d_fid, quadratic_hex=quadratic_hex, dbgout=1)
            elseif (k==2) then
                call write_vtk_mesh(equi, mesh_stag, phi_stag, dphi/2.0_GP, dphi/2.0_GP, &
                                    vtk2d_fid, vtk3d_fid, quadratic_hex=quadratic_hex, dbgout=1)
            endif

            vtk_stat = vtk2d_fid%xml_writer%write_piece()
            vtk_stat = vtk2d_fid%finalize()
            
            vtk_stat = vtk3d_fid%xml_writer%write_piece()
            vtk_stat = vtk3d_fid%finalize()

            deallocate(vtk3d_fid)
            deallocate(vtk2d_fid)
        enddo

        if (is_rank_info_writer) then
            write(get_stdout(),*)repeat('-',80)
            write(get_stdout(),*)''
        endif

    end subroutine

#else
    ! If not compiled with -DENABLE_VTK=ON, fall back to empty placeholder routine

    module subroutine build_vtk_mesh(comm_handler, par_filepath, dirpath)
        type(comm_handler_t), intent(in) :: comm_handler
        character(len=*), intent(in), optional :: par_filepath
        character(len=*), intent(in), optional :: dirpath

        if (cntrl_build_vtkmesh == NONE) then
            return
        else
            call handle_error("You are requesting a VTK mesh build, but code is not compiled with this feature, &
                               Recompile with -DENABLE_VTK=ON ", &
                              GRILLIX_ERR_STATIC_DATA, __LINE__, __FILE__)
        endif

    end subroutine

#endif

end submodule