module vis_vtk3d_m !! 3d Visualisation of mesh: Lifts mesh to 3d (aligned boxes) and outpots to vtk use vtk_fortran, only: vtk_file use precision_m, only: FP, FP_EPS, SP use screen_io_m, only: get_stdout, get_stderr use comm_handling_m, only: is_master use list_operations_m, only: unique_tuples, findindex_2d use equilibrium_m, only: equilibrium_t use mesh_cart_m, only: mesh_cart_t use fieldline_tracer_m, only: trace implicit none public :: write_vtk_mesh private :: build_cells_2d private :: build_cells_3d contains subroutine write_vtk_mesh(equi, mesh, phi, dphi_bwd, dphi_fwd, vtk2d, vtk3d, quadratic_hex, dbgout) !! Writes mesh to vtk files !! Based on VTKFortran library (https://github.com/szaghi/VTKFortran/) !! For description of cell types and their definition of connectivity see https://kitware.github.io/vtk-examples/site/VTKBook/05Chapter5/ !! here used (quad, bi-quadrtatic quad, hexahedron, quadratic hexahedron) !! For their VTK descriptor see https://vtk.org/doc/nightly/html/vtkCellType_8h_source.html class(equilibrium_t), intent(inout) :: equi !! Mesh type(mesh_cart_t), intent(in) :: mesh !! Mesh real(FP), intent(in) :: phi !! Toroidal position at which center 3D mesh is located real(FP), intent(in) :: dphi_bwd !! Toroidal distance in negative phi direction for which 3D aligned mesh is built real(FP), intent(in) :: dphi_fwd !! Toroidal distance in positive phi direction for which 3D aligned mesh is built type(vtk_file), intent(inout) , optional :: vtk2d !! VTK file for 2d vtk mesh (inplane) type(vtk_file), intent(inout) , optional :: vtk3d !! VTK file for 3d vtk mesh integer, intent(in), optional :: dbgout !! Debug output level logical, intent(in), optional :: quadratic_hex !! Indicates if 3d mesh is based on standard hexahedrons or quadratic hexahedrons (default) logical :: quad_hex integer :: outi, verr, k integer :: nnodes_per_cell2d, nnodes2d, ncells2d !! Number of nodes and cells of 2d vtk mesh real(FP), allocatable, dimension(:) :: xnodes2d, ynodes2d, znodes2d_dummy !! Coordinates of 2D vtk nodes integer, allocatable, dimension(:,:) :: connectivity2d !! Connectivity information for 2D cells --> see description of VTKFortran library integer :: nnodes_per_cell3d, nnodes3d, ncells3d !! Number of nodes and cells of 2d vtk mesh real(FP), allocatable, dimension(:) :: xc_nodes3d, yc_nodes3d, zc_nodes3d !! Coordinates of 3D vtk nodes integer, allocatable, dimension(:,:) :: connectivity3d !! Connectivity information for 3D cells --> see description of VTKFortran library integer, allocatable, dimension(:) :: connectivity2d_mod, connectivity3d_mod, offset !! --> see description of VTKFortran library integer(1), allocatable, dimension(:) :: cell_type !! --> see description of VTKFortran library ! Build nodes and cells for 2d mesh ! set output level outi = 0 if (present(dbgout)) then if (is_master()) then outi = dbgout endif if (dbgout >= 3) then outi = dbgout ! every rank writes endif endif if (present(quadratic_hex)) then quad_hex = quadratic_hex else quad_hex = .true. endif if (outi >= 1) then write(get_stdout(), *) '' write(get_stdout(), *) '------------------------------------------------' write(get_stdout(), *) 'Building vtk output for mesh' write(get_stdout(), *) 'quadratic_hex = ', quad_hex endif call build_cells_2d(mesh, quad_hex, nnodes_per_cell2d, nnodes2d, ncells2d, xnodes2d, ynodes2d, connectivity2d) allocate(znodes2d_dummy(nnodes2d),source = 0.0_FP) ! Dummy, for embedding 2d mesh into 3d if (outi >= 1) then write(get_stdout(), 201) nnodes_per_cell2d, nnodes2d, ncells2d 201 FORMAT('...2d nodes and cells built...' / & 3X,'nnodes_per_cell2d = ',I9, / & 3X,'nnodes2d = ',I9, / & 3X,'ncells2d = ',I9 ) endif ! Modify connectivity (1d array, entries refer to 0 based indexing) ! offset and cell_shapes allocate(connectivity2d_mod(nnodes_per_cell2d*ncells2d)) connectivity2d_mod = reshape(connectivity2d, (/nnodes_per_cell2d*ncells2d/) ) connectivity2d_mod = connectivity2d_mod - 1 allocate(offset(ncells2d) ) allocate(cell_type(ncells2d)) !$OMP PARALLEL PRIVATE(k) !$OMP DO do k = 1, ncells2d offset(k) = nnodes_per_cell2d*k if (quad_hex) then cell_type(k) = 28_1 ! Biquadratic quad type else cell_type(k) = 9_1 ! Quad type endif enddo !$OMP END DO !$OMP END PARALLEL ! Write to file if (present(vtk2d)) then if (outi >= 1) then write(get_stdout(), *) '...writing 2d vtk-mesh...' endif verr = vtk2d%xml_writer%write_piece(np=nnodes2d, nc=ncells2d) call vtk_error_stop(verr) verr = vtk2d%xml_writer%write_geo(np=nnodes2d, nc=ncells2d, & x = real(xnodes2d, kind =SP), & y = real(ynodes2d, kind = SP), & z = real(znodes2d_dummy, kind = SP) ) call vtk_error_stop(verr) verr = vtk2d%xml_writer%write_connectivity(nc = ncells2d, & connectivity = connectivity2d_mod, & offset = offset, & cell_type = cell_type ) call vtk_error_stop(verr) endif deallocate(cell_type) deallocate(offset) if (outi >= 1) then write(get_stdout(), *) '...building 3d nodes and cells...' endif call build_cells_3d(equi, mesh, quad_hex, & nnodes_per_cell2d, nnodes2d, ncells2d, xnodes2d, ynodes2d, connectivity2d, & phi, dphi_bwd, dphi_fwd, & nnodes_per_cell3d, nnodes3d, ncells3d, xc_nodes3d, yc_nodes3d, zc_nodes3d, connectivity3d ) if (outi >= 1) then write(get_stdout(), 202) nnodes_per_cell3d, nnodes3d, ncells3d 202 FORMAT('...3d nodes and cells built...' / & 3X,'nnodes_per_cell3d = ',I9, / & 3X,'nnodes3d = ',I9, / & 3X,'ncells3d = ',I9 ) endif ! Modify connectivity (1d array, entries refer to 0 based indexing) ! offset and cell_shapes allocate(connectivity3d_mod(nnodes_per_cell3d*ncells3d)) connectivity3d_mod = reshape(connectivity3d, (/nnodes_per_cell3d*ncells3d/) ) connectivity3d_mod = connectivity3d_mod - 1 allocate(offset(ncells3d) ) allocate(cell_type(ncells3d)) !$OMP PARALLEL PRIVATE(k) !$OMP DO do k = 1, ncells3d offset(k) = nnodes_per_cell3d*k if (quad_hex) then cell_type(k) = 25_1 ! Quadratic Hexahedron type else cell_type(k) = 12_1 ! Hexahedron endif enddo !$OMP END DO !$OMP END PARALLEL if (present(vtk3d)) then if (outi >= 1) then write(get_stdout(), *) '...writing 3d vtk-mesh...' endif verr = vtk3d%xml_writer%write_piece(np=nnodes3d, nc=ncells3d) call vtk_error_stop(verr) verr = vtk3d%xml_writer%write_geo(np=nnodes3d, nc=ncells3d, & x = real(xc_nodes3d, kind = SP), & y = real(yc_nodes3d, kind = SP), & z = real(zc_nodes3d, kind = SP) ) call vtk_error_stop(verr) verr = vtk3d%xml_writer%write_connectivity(nc = ncells3d, & connectivity = connectivity3d_mod, & offset = offset, & cell_type = cell_type ) call vtk_error_stop(verr) endif if (outi >= 1) then write(get_stdout(), *) & '------------------------------------------------' write(get_stdout(), *) '' endif contains subroutine vtk_error_stop(status) ! stop at vtk error and print error number integer, intent(in) ::status if (status /= 0) then write(get_stderr(), *) & 'error(write_vtk_mesh): vtk error', status error stop endif end subroutine end subroutine subroutine build_cells_2d(mesh, quad_hex, nnodes_per_cell2d, nnodes2d, ncells2d, xnodes2d, ynodes2d, connectivity2d) !! Builds inplane 2D cells, i.e. nodes and its connectivity around mesh points type(mesh_cart_t), intent(in) :: mesh !! Mesh logical, intent(in) :: quad_hex !! If quad_hex = false builds standard quads !! If quad_hex = true builds bi-quadratic quads integer, intent(out) :: nnodes_per_cell2d !! Number of nodes per cell integer, intent(out) :: nnodes2d !! Number of nodes integer, intent(out) :: ncells2d !! Number of cells real(FP), allocatable, dimension(:) :: xnodes2d !! x-coordinates of nodes real(FP), allocatable, dimension(:) :: ynodes2d !! y-coordinates of nodes integer, allocatable, dimension(:,:) :: connectivity2d !! Connectivity of cells integer :: np, l, ic, jc, ic_tw, jc_tw, ns, in, jn, i, j, k, lg real(FP) :: ic_f, jc_f integer, allocatable, dimension(:,:) :: ij_nodes ncells2d = mesh%get_n_points() ! Find indizes of nodes (staggered to mesh points) if (quad_hex) then nnodes_per_cell2d = 9 else nnodes_per_cell2d = 4 endif nnodes2d = nnodes_per_cell2d * ncells2d allocate(ij_nodes(2, nnodes2d)) !$OMP PARALLEL PRIVATE(l, ic, jc, ic_tw, jc_tw, ns) !$OMP DO do l = 1, ncells2d ! Discrete coordinates of nodes ic = mesh%get_cart_i(l) jc = mesh%get_cart_j(l) ! Trick: i and j indices are multiplied for two to deal with in plane staggering ic_tw = 2 * ic jc_tw = 2 * jc ns = nnodes_per_cell2d*(l-1)+1 ij_nodes(1, ns) = ic_tw - 1 ij_nodes(2, ns) = jc_tw - 1 ij_nodes(1, ns + 1) = ic_tw + 1 ij_nodes(2, ns + 1) = jc_tw - 1 ij_nodes(1, ns + 2) = ic_tw + 1 ij_nodes(2, ns + 2) = jc_tw + 1 ij_nodes(1, ns + 3) = ic_tw - 1 ij_nodes(2, ns + 3) = jc_tw + 1 if (quad_hex) then ! For buiolding bi-quadratic quads, further nodes are added ij_nodes(1, ns + 4) = ic_tw ij_nodes(2, ns + 4) = jc_tw - 1 ij_nodes(1, ns + 5) = ic_tw + 1 ij_nodes(2, ns + 5) = jc_tw ij_nodes(1, ns + 6) = ic_tw ij_nodes(2, ns + 6) = jc_tw + 1 ij_nodes(1, ns + 7) = ic_tw - 1 ij_nodes(2, ns + 7) = jc_tw ij_nodes(1, ns + 8) = ic_tw ij_nodes(2, ns + 8) = jc_tw endif enddo !$OMP END DO !$OMP END PARALLEL call unique_tuples(2, nnodes2d, ij_nodes) ! Coordinates of nodes allocate(xnodes2d(nnodes2d)) allocate(ynodes2d(nnodes2d)) !$OMP PARALLEL PRIVATE(k, ic_f, jc_f) !$OMP DO do k = 1, nnodes2d ! Undo staggering --> half valued indices ic_f = 0.5_FP * ij_nodes(1,k) jc_f = 0.5_FP * ij_nodes(2,k) xnodes2d(k) = mesh%get_xmin() + (ic_f-1.0_FP) * mesh%get_spacing_f() ynodes2d(k) = mesh%get_ymin() + (jc_f-1.0_FP) * mesh%get_spacing_f() enddo !$OMP END DO !$OMP END PARALLEL ! Connectivity allocate(connectivity2d(nnodes_per_cell2d, ncells2d)) !$OMP PARALLEL PRIVATE(l, ic, jc, ic_tw, jc_tw, in, jn, lg) !$OMP DO do l = 1, ncells2d ic = mesh%get_cart_i(l) jc = mesh%get_cart_j(l) ic_tw = 2 * ic jc_tw = 2 * jc in = ic_tw-1 jn = jc_tw-1 call findindex_2d(1, nnodes2d, ij_nodes(1,:), ij_nodes(2,:) , in, jn, lg) if (lg == 0) then write(get_stderr(), *) & 'error(mesh%vis_inplane) node point not found',l, in, jn error stop endif connectivity2d(1, l) = lg in = ic_tw + 1 jn = jc_tw - 1 call findindex_2d(1, nnodes2d, ij_nodes(1,:), ij_nodes(2,:) , in, jn, lg) if (lg == 0) then write(get_stderr(), *) & 'error(mesh%vis_inplane) node point not found',l, in, jn error stop endif connectivity2d(2, l) = lg in = ic_tw + 1 jn = jc_tw + 1 call findindex_2d(1, nnodes2d, ij_nodes(1,:), ij_nodes(2,:) , in, jn, lg) if (lg == 0) then write(get_stderr(), *) & 'error(mesh%vis_inplane) node point not found',l, in, jn error stop endif connectivity2d(3, l) = lg in = ic_tw - 1 jn = jc_tw + 1 call findindex_2d(1, nnodes2d, ij_nodes(1,:), ij_nodes(2,:) , in, jn, lg) if (lg == 0) then write(get_stderr(), *) & 'error(mesh%vis_inplane) node point not found',l, in, jn error stop endif connectivity2d(4, l) = lg if (quad_hex) then in = ic_tw jn = jc_tw - 1 call findindex_2d(1, nnodes2d, ij_nodes(1,:), ij_nodes(2,:) , in, jn, lg) if (lg == 0) then write(get_stderr(), *) & 'error(mesh%vis_inplane) node point not found',l, in, jn error stop endif connectivity2d(5, l) = lg in = ic_tw + 1 jn = jc_tw call findindex_2d(1, nnodes2d, ij_nodes(1,:), ij_nodes(2,:) , in, jn, lg) if (lg == 0) then write(get_stderr(), *) & 'error(mesh%vis_inplane) node point not found',l, in, jn error stop endif connectivity2d(6, l) = lg in = ic_tw jn = jc_tw + 1 call findindex_2d(1, nnodes2d, ij_nodes(1,:), ij_nodes(2,:) , in, jn, lg) if (lg == 0) then write(get_stderr(), *) & 'error(mesh%vis_inplane) node point not found',l, in, jn error stop endif connectivity2d(7, l) = lg in = ic_tw - 1 jn = jc_tw call findindex_2d(1, nnodes2d, ij_nodes(1,:), ij_nodes(2,:) , in, jn, lg) if (lg == 0) then write(get_stderr(), *) & 'error(mesh%vis_inplane) node point not found',l, in, jn error stop endif connectivity2d(8, l) = lg in = ic_tw jn = jc_tw call findindex_2d(1, nnodes2d, ij_nodes(1,:), ij_nodes(2,:) , in, jn, lg) if (lg == 0) then write(get_stderr(), *) & 'error(mesh%vis_inplane) node point not found',l, in, jn error stop endif connectivity2d(9, l) = lg endif enddo !$OMP END DO !$OMP END PARALLEL end subroutine subroutine build_cells_3d(equi, mesh, quad_hex, & nnodes_per_cell2d, nnodes2d, ncells2d, xnodes2d, ynodes2d, connectivity2d, & phi, dphi_bwd, dphi_fwd, & nnodes_per_cell3d, nnodes3d, ncells3d, xc_nodes3d, yc_nodes3d, zc_nodes3d, connectivity3d ) !! Builds 3D cells, i.e. nodes and its connectivity !! Based on extension of 2D cells class(equilibrium_t), intent(inout) :: equi !! Equilibrium type(mesh_cart_t), intent(in) :: mesh !! Mesh logical, intent(in) :: quad_hex !! If quad_hex = false builds standard hexehedra !! If quad_hex = true builds quadratic hexahedra integer, intent(in) :: nnodes_per_cell2d !! Number of nodes per cell of 2d Cells integer, intent(in) :: nnodes2d !! Number of 2d nodes integer, intent(in) :: ncells2d !! Number of 2d cells real(FP), dimension(nnodes2d), intent(in) :: xnodes2d !! x-coordinates of 2d nodes real(FP), dimension(nnodes2d), intent(in) :: ynodes2d !! y-coordinates of 2d nodes integer, dimension(nnodes_per_cell2d, ncells2d) :: connectivity2d !! Connectivity of 2d cells real(FP), intent(in) :: phi !! Toroidal position at which center 3D mesh is located real(FP), intent(in) :: dphi_bwd !! Toroidal distance in negative phi direction for which 3D aligned mesh is built real(FP), intent(in) :: dphi_fwd !! Toroidal distance in positive phi direction for which 3D aligned mesh is built integer, intent(out) :: nnodes_per_cell3d !! Number of nodes per cell of3d Cells integer, intent(out) :: nnodes3d !! Number of 3d nodes integer, intent(out) :: ncells3d !! Number of 3d cells real(FP), allocatable, dimension(:), intent(out) :: xc_nodes3d !! x-coordinates of 3d nodes real(FP), allocatable, dimension(:), intent(out) :: yc_nodes3d !! y-coordinates of 3d nodes real(FP), allocatable, dimension(:), intent(out) :: zc_nodes3d !! z-coordinates of 3d nodes integer, allocatable, dimension(:,:), intent(out) :: connectivity3d !! Connectivity of cells integer :: k real(FP) :: rmaj, zver, tor if (quad_hex) then nnodes_per_cell3d = 20 nnodes3d = nnodes2d*3 else nnodes_per_cell3d = 8 nnodes3d = nnodes2d*2 endif ncells3d = ncells2d allocate(xc_nodes3d(nnodes3d)) allocate(yc_nodes3d(nnodes3d)) allocate(zc_nodes3d(nnodes3d)) allocate(connectivity3d(nnodes_per_cell3d,ncells3d)) if (quad_hex) then do k = 1, nnodes2d xc_nodes3d(k) = xnodes2d(k) yc_nodes3d(k) = ynodes2d(k) zc_nodes3d(k) = phi call trace(x_init=xnodes2d(k), & y_init=ynodes2d(k), & phi_init=phi, & dphi=-dphi_bwd, & equi=equi, & x_end=xc_nodes3d(k+nnodes2d), & y_end=yc_nodes3d(k+nnodes2d) ) zc_nodes3d(k+nnodes2d) = phi - dphi_bwd call trace(x_init=xnodes2d(k), & y_init=ynodes2d(k), & phi_init=phi, & dphi=dphi_fwd, & equi=equi, & x_end=xc_nodes3d(k+2*nnodes2d), & y_end=yc_nodes3d(k+2*nnodes2d) ) zc_nodes3d(k+2*nnodes2d) = phi + dphi_fwd enddo do k = 1, ncells2d connectivity3d( 9,k) = connectivity2d(1,k) connectivity3d(11,k) = connectivity2d(2,k) connectivity3d(15,k) = connectivity2d(3,k) connectivity3d(13,k) = connectivity2d(4,k) connectivity3d( 2,k) = connectivity2d(1,k) + nnodes2d connectivity3d( 3,k) = connectivity2d(2,k) + nnodes2d connectivity3d( 7,k) = connectivity2d(3,k) + nnodes2d connectivity3d( 6,k) = connectivity2d(4,k) + nnodes2d connectivity3d(10,k) = connectivity2d(5,k) + nnodes2d connectivity3d(19,k) = connectivity2d(6,k) + nnodes2d connectivity3d(14,k) = connectivity2d(7,k) + nnodes2d connectivity3d(18,k) = connectivity2d(8,k) + nnodes2d connectivity3d( 1,k) = connectivity2d(1,k) + 2*nnodes2d connectivity3d( 4,k) = connectivity2d(2,k) + 2*nnodes2d connectivity3d( 8,k) = connectivity2d(3,k) + 2*nnodes2d connectivity3d( 5,k) = connectivity2d(4,k) + 2*nnodes2d connectivity3d(12,k) = connectivity2d(5,k) + 2*nnodes2d connectivity3d(20,k) = connectivity2d(6,k) + 2*nnodes2d connectivity3d(16,k) = connectivity2d(7,k) + 2*nnodes2d connectivity3d(17,k) = connectivity2d(8,k) + 2*nnodes2d enddo else do k = 1, nnodes2d call trace(x_init=xnodes2d(k), & y_init=ynodes2d(k), & phi_init=phi, & dphi=-dphi_bwd, & equi=equi, & x_end=xc_nodes3d(k), & y_end=yc_nodes3d(k) ) zc_nodes3d(k) = phi - dphi_bwd call trace(x_init=xnodes2d(k), & y_init=ynodes2d(k), & phi_init=phi, & dphi=dphi_fwd, & equi=equi, & x_end=xc_nodes3d(k+nnodes2d), & y_end=yc_nodes3d(k+nnodes2d) ) zc_nodes3d(k+nnodes2d) = phi + dphi_fwd enddo do k = 1, ncells2d connectivity3d(1:4,k) = connectivity2d(1:4,k) connectivity3d(5:8,k) = connectivity2d(1:4,k) + nnodes2d enddo endif ! For toroidal geometries, transform data from (R,phi,Z) to (x,y,z) if (equi%jacobian(0.0_FP, 0.0_FP, phi) <= FP_EPS) then !$OMP PARALLEL PRIVATE(k, rmaj, zver, tor) !$OMP DO do k = 1, nnodes3d rmaj = xc_nodes3d(k) zver = yc_nodes3d(k) tor = zc_nodes3d(k) xc_nodes3d(k) = rmaj * cos(tor) yc_nodes3d(k) = rmaj * sin(tor) zc_nodes3d(k) = zver enddo !$OMP END DO !$OMP END PARALLEL endif end subroutine end module