mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 06:05:58 -04:00
Added capability to plot meshlines
This commit is contained in:
Nick Horelik
parent
dc47763f66
commit
bb8a95b3f1
3 changed files with 174 additions and 7 deletions
|
|
@ -2683,7 +2683,7 @@ contains
|
|||
subroutine read_plots_xml()
|
||||
|
||||
integer i, j
|
||||
integer n_cols, col_id, n_comp, n_masks
|
||||
integer n_cols, col_id, n_comp, n_masks, n_meshlines
|
||||
integer, allocatable :: iarray(:)
|
||||
logical :: file_exists ! does plots.xml file exist?
|
||||
character(MAX_LINE_LEN) :: filename ! absolute path to plots.xml
|
||||
|
|
@ -2693,9 +2693,11 @@ contains
|
|||
type(Node), pointer :: node_plot => null()
|
||||
type(Node), pointer :: node_col => null()
|
||||
type(Node), pointer :: node_mask => null()
|
||||
type(Node), pointer :: node_meshlines => null()
|
||||
type(NodeList), pointer :: node_plot_list => null()
|
||||
type(NodeList), pointer :: node_col_list => null()
|
||||
type(NodeList), pointer :: node_mask_list => null()
|
||||
type(NodeList), pointer :: node_meshline_list => null()
|
||||
|
||||
! Check if plots.xml exists
|
||||
filename = trim(path_input) // "plots.xml"
|
||||
|
|
@ -2947,6 +2949,67 @@ contains
|
|||
end do
|
||||
end if
|
||||
|
||||
! Deal with meshlines
|
||||
call get_node_list(node_plot, "meshlines", node_meshline_list)
|
||||
n_meshlines = get_list_size(node_meshline_list)
|
||||
if (n_meshlines /= 0) then
|
||||
|
||||
if (pl % type == PLOT_TYPE_VOXEL) then
|
||||
message = "Meshlines ignored in voxel plot " // &
|
||||
trim(to_str(pl % id))
|
||||
call warning()
|
||||
end if
|
||||
|
||||
select case(n_meshlines)
|
||||
case default
|
||||
message = "Mutliple meshlines" // &
|
||||
" specified in plot " // trim(to_str(pl % id))
|
||||
call fatal_error()
|
||||
case (0)
|
||||
case (1)
|
||||
|
||||
! Get pointer to meshlines
|
||||
call get_list_item(node_meshline_list, 1, node_meshlines)
|
||||
|
||||
! Ensure that there is a mesh id for this meshlines specification
|
||||
if (check_for_node(node_meshlines, "mesh")) then
|
||||
call get_node_value(node_meshlines, "mesh", pl % meshlines_id)
|
||||
else
|
||||
message = "Must specify a mesh id for meshlines " // &
|
||||
"specification in plot " // trim(to_str(pl % id))
|
||||
call fatal_error()
|
||||
end if
|
||||
|
||||
! Ensure that there is a linewidth for this meshlines specification
|
||||
if (check_for_node(node_meshlines, "linewidth")) then
|
||||
call get_node_value(node_meshlines, "linewidth", &
|
||||
pl % meshlines_width)
|
||||
else
|
||||
message = "Must specify a linewidth for meshlines " // &
|
||||
"specification in plot " // trim(to_str(pl % id))
|
||||
call fatal_error()
|
||||
end if
|
||||
|
||||
! Check if the specified tally mesh exists
|
||||
if (mesh_dict % has_key(pl % meshlines_id)) then
|
||||
pl % meshlines_id = mesh_dict % get_key(pl % meshlines_id)
|
||||
if (meshes(pl % meshlines_id) % type /= LATTICE_RECT) then
|
||||
message = "Non-rectangular mesh specified in meshlines for" // &
|
||||
" plot " // trim(to_str(pl % id))
|
||||
call fatal_error()
|
||||
end if
|
||||
else
|
||||
message = "Could not find tally mesh " // &
|
||||
trim(to_str(pl % meshlines_id)) // &
|
||||
" specified in meshlines for plot " // &
|
||||
trim(to_str(pl % id))
|
||||
call fatal_error()
|
||||
end if
|
||||
|
||||
end select
|
||||
|
||||
end if
|
||||
|
||||
! Deal with masks
|
||||
call get_node_list(node_plot, "mask", node_mask_list)
|
||||
n_masks = get_list_size(node_mask_list)
|
||||
|
|
|
|||
114
src/plot.F90
114
src/plot.F90
|
|
@ -5,6 +5,7 @@ module plot
|
|||
use geometry, only: find_cell, check_cell_overlap
|
||||
use geometry_header, only: Cell, BASE_UNIVERSE
|
||||
use global
|
||||
use mesh, only: get_mesh_indices
|
||||
use output, only: write_message
|
||||
use particle_header, only: deallocate_coord, Particle
|
||||
use plot_header
|
||||
|
|
@ -118,27 +119,24 @@ contains
|
|||
call init_image(img)
|
||||
call allocate_image(img, pl % pixels(1), pl % pixels(2))
|
||||
|
||||
in_pixel = pl % width(1)/dble(pl % pixels(1))
|
||||
out_pixel = pl % width(2)/dble(pl % pixels(2))
|
||||
|
||||
if (pl % basis == PLOT_BASIS_XY) then
|
||||
in_i = 1
|
||||
out_i = 2
|
||||
in_pixel = pl % width(1)/dble(pl % pixels(1))
|
||||
out_pixel = pl % width(2)/dble(pl % pixels(2))
|
||||
xyz(1) = pl % origin(1) - pl % width(1) / 2.0
|
||||
xyz(2) = pl % origin(2) + pl % width(2) / 2.0
|
||||
xyz(3) = pl % origin(3)
|
||||
else if (pl % basis == PLOT_BASIS_XZ) then
|
||||
in_i = 1
|
||||
out_i = 3
|
||||
in_pixel = pl % width(1)/dble(pl % pixels(1))
|
||||
out_pixel = pl % width(2)/dble(pl % pixels(2))
|
||||
xyz(1) = pl % origin(1) - pl % width(1) / 2.0
|
||||
xyz(2) = pl % origin(2)
|
||||
xyz(3) = pl % origin(3) + pl % width(2) / 2.0
|
||||
else if (pl % basis == PLOT_BASIS_YZ) then
|
||||
in_i = 2
|
||||
out_i = 3
|
||||
in_pixel = pl % width(1)/dble(pl % pixels(1))
|
||||
out_pixel = pl % width(2)/dble(pl % pixels(2))
|
||||
xyz(1) = pl % origin(1)
|
||||
xyz(2) = pl % origin(2) - pl % width(1) / 2.0
|
||||
xyz(3) = pl % origin(3) + pl % width(2) / 2.0
|
||||
|
|
@ -169,6 +167,9 @@ contains
|
|||
p % coord0 % xyz(out_i) = p % coord0 % xyz(out_i) - out_pixel
|
||||
end do
|
||||
|
||||
! Draw tally mesh boundaries on the image if requested
|
||||
if (pl % meshlines_id > 0) call draw_mesh_lines(pl, img)
|
||||
|
||||
! Write out the ppm to a file
|
||||
call output_ppm(pl,img)
|
||||
|
||||
|
|
@ -180,6 +181,107 @@ contains
|
|||
|
||||
end subroutine create_ppm
|
||||
|
||||
!===============================================================================
|
||||
! DRAW_MESH_LINES draws mesh line boundaries on an image
|
||||
!===============================================================================
|
||||
subroutine draw_mesh_lines(pl, img)
|
||||
|
||||
type(ObjectPlot), pointer :: pl
|
||||
type(Image) :: img
|
||||
|
||||
logical :: in_mesh
|
||||
integer :: n
|
||||
integer :: x, y ! pixel location
|
||||
integer :: xrange(2), yrange(2) ! range of pixel locations
|
||||
integer :: i, j ! loop indices
|
||||
integer :: ijk(3)
|
||||
integer :: plus
|
||||
integer :: ijk_ll(3) ! mesh bin ijk indicies of plot lower left
|
||||
integer :: ijk_ur(3) ! mesh bin ijk indicies of plot upper right
|
||||
real(8) :: frac
|
||||
real(8) :: width(3) ! real widths of the plot
|
||||
real(8) :: xyz_ll_plot(3) ! lower left xyz of plot image
|
||||
real(8) :: xyz_ur_plot(3) ! upper right xyz of plot image
|
||||
real(8) :: xyz_ll(3) ! lower left xyz
|
||||
real(8) :: xyz_ur(3) ! upper right xyz
|
||||
type(StructuredMesh), pointer :: m => null()
|
||||
|
||||
m => meshes(pl % meshlines_id)
|
||||
n = m % n_dimension
|
||||
|
||||
select case (pl % basis)
|
||||
case(PLOT_BASIS_XY)
|
||||
xyz_ll_plot(1) = pl % origin(1) - pl % width(1) / 2.0
|
||||
xyz_ll_plot(2) = pl % origin(2) - pl % width(2) / 2.0
|
||||
xyz_ll_plot(3) = pl % origin(3)
|
||||
xyz_ur_plot(1) = pl % origin(1) + pl % width(1) / 2.0
|
||||
xyz_ur_plot(2) = pl % origin(2) + pl % width(2) / 2.0
|
||||
xyz_ur_plot(3) = pl % origin(3)
|
||||
|
||||
width = xyz_ur_plot - xyz_ll_plot
|
||||
|
||||
call get_mesh_indices(m, xyz_ll_plot, ijk_ll(:m % n_dimension), in_mesh)
|
||||
call get_mesh_indices(m, xyz_ur_plot, ijk_ur(:m % n_dimension), in_mesh)
|
||||
|
||||
! sweep through all meshbins on this plane and draw borders
|
||||
ijk(3) = 0
|
||||
do i = ijk_ll(1), ijk_ur(1)
|
||||
do j = ijk_ll(2), ijk_ur(2)
|
||||
! check if we're in the mesh for this ijk
|
||||
if (i > 0 .and. i <= m % dimension(1) .and. &
|
||||
j > 0 .and. j <= m % dimension(2)) then
|
||||
|
||||
! get xyz's of lower left and upper right of this mesh cell
|
||||
ijk(1) = i
|
||||
ijk(2) = j
|
||||
xyz_ll(:m % n_dimension) = &
|
||||
m % lower_left + m % width * (ijk(:m % n_dimension) - 1)
|
||||
xyz_ur(:m % n_dimension) = &
|
||||
m % lower_left + m % width * ijk(:m % n_dimension)
|
||||
|
||||
! map the xyz to pixel locations
|
||||
frac = (xyz_ll(1) - xyz_ll_plot(1)) / width(1)
|
||||
xrange(1) = int(frac * real(img % width,8))
|
||||
frac = (xyz_ur(1) - xyz_ll_plot(1)) / width(1)
|
||||
xrange(2) = int(frac * real(img % width,8))
|
||||
|
||||
frac = (xyz_ll(2) - xyz_ll_plot(2)) / width(2)
|
||||
yrange(1) = img % height - int(frac * real(img % height,8))
|
||||
frac = (xyz_ur(2) - xyz_ll_plot(2)) / width(2)
|
||||
yrange(2) = img % height - int(frac * real(img % height,8))
|
||||
|
||||
! draw black lines
|
||||
do x = xrange(1), xrange(2)
|
||||
do plus = 0, pl % meshlines_width
|
||||
call set_pixel(img, x, yrange(1) + plus, 0, 0, 0)
|
||||
call set_pixel(img, x, yrange(2) + plus, 0, 0, 0)
|
||||
call set_pixel(img, x, yrange(1) - plus, 0, 0, 0)
|
||||
call set_pixel(img, x, yrange(2) - plus, 0, 0, 0)
|
||||
end do
|
||||
end do
|
||||
do y = yrange(2), yrange(1)
|
||||
do plus = 0, pl % meshlines_width
|
||||
call set_pixel(img, xrange(1) + plus, y, 0, 0, 0)
|
||||
call set_pixel(img, xrange(2) + plus, y, 0, 0, 0)
|
||||
call set_pixel(img, xrange(1) - plus, y, 0, 0, 0)
|
||||
call set_pixel(img, xrange(2) - plus, y, 0, 0, 0)
|
||||
end do
|
||||
end do
|
||||
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
|
||||
case(PLOT_BASIS_XZ)
|
||||
message = "Meshline plotting currently only implemented for basis xy"
|
||||
call fatal_error()
|
||||
case(PLOT_BASIS_YZ)
|
||||
message = "Meshline plotting currently only implemented for basis xy"
|
||||
call fatal_error()
|
||||
end select
|
||||
|
||||
end subroutine draw_mesh_lines
|
||||
|
||||
!===============================================================================
|
||||
! OUTPUT_PPM writes out a previously generated image to a PPM file
|
||||
!===============================================================================
|
||||
|
|
|
|||
|
|
@ -25,6 +25,8 @@ module plot_header
|
|||
real(8) :: width(3) ! xyz widths of plot
|
||||
integer :: basis ! direction of plot slice
|
||||
integer :: pixels(3) ! pixel width/height of plot slice
|
||||
integer :: meshlines_id = -1 ! id of the mesh to plot lines of
|
||||
integer :: meshlines_width ! pixel width of meshlines
|
||||
type(ObjectColor) :: not_found ! color for positions where no cell found
|
||||
type(ObjectColor), allocatable :: colors(:) ! colors of cells/mats
|
||||
end type ObjectPlot
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue