From bb8a95b3f1cca5311e20c6094bda296be9a1e118 Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Wed, 20 Aug 2014 10:22:00 -0400 Subject: [PATCH 01/17] Added capability to plot meshlines --- src/input_xml.F90 | 65 ++++++++++++++++++++++++- src/plot.F90 | 114 +++++++++++++++++++++++++++++++++++++++++--- src/plot_header.F90 | 2 + 3 files changed, 174 insertions(+), 7 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 1869b9709d..3f61f8d764 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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) diff --git a/src/plot.F90 b/src/plot.F90 index e7642314f9..842b1c089f 100644 --- a/src/plot.F90 +++ b/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 !=============================================================================== diff --git a/src/plot_header.F90 b/src/plot_header.F90 index 6956bfaf93..5a301f8809 100644 --- a/src/plot_header.F90 +++ b/src/plot_header.F90 @@ -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 From f77bbb2336547b49fc9272c812ac06df8495fca2 Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Wed, 20 Aug 2014 10:28:02 -0400 Subject: [PATCH 02/17] Removed unused variable --- src/plot.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/plot.F90 b/src/plot.F90 index 842b1c089f..25c0c7123a 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -190,7 +190,6 @@ contains 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 @@ -207,7 +206,6 @@ contains type(StructuredMesh), pointer :: m => null() m => meshes(pl % meshlines_id) - n = m % n_dimension select case (pl % basis) case(PLOT_BASIS_XY) From 583f73764a5664db3e628fe02de11eaadb315c76 Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Wed, 20 Aug 2014 10:31:16 -0400 Subject: [PATCH 03/17] Updated users guide for meshline plotting --- docs/source/usersguide/input.rst | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index c9aa8cc866..290713c35a 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1229,6 +1229,25 @@ attributes or sub-elements. These are not used in "voxel" plots: *Default*: None + :meshlines: + The special ``meshlines`` sub-element allows for plotting the boundaries of + a tally mesh on top of a plot. Only one ``meshlines`` element is allowed per + ``plot`` element, and it must contain as attributes or sub-elements a mesh + id and a linewidth: + + :mesh: + A single integer id number for the mesh specified on ``tallies.xml`` that + should be plotted. + + :linewidth: + A single integer number of pixels to specify the linewidth of the mesh + boundaries. + + *Default*: None + + .. warning:: Meshline plotting is currently only implemented for plots with + an "xy" basis. + ------------------------------ CMFD Specification -- cmfd.xml ------------------------------ From 2ae8fe4e8ae0f00f43e2079a3a66085a00fb5f08 Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Wed, 20 Aug 2014 10:45:42 -0400 Subject: [PATCH 04/17] Added commented-out example of meshline plotting --- examples/lattice/simple/plots.xml | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/lattice/simple/plots.xml b/examples/lattice/simple/plots.xml index a8c6fbb7c5..a25a5a8280 100644 --- a/examples/lattice/simple/plots.xml +++ b/examples/lattice/simple/plots.xml @@ -5,6 +5,7 @@ 0. 0. 0. 4.0 4.0 400 400 + From cdcbd277ea8ebea9d441fbd0a24831f76f9143fc Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Wed, 20 Aug 2014 10:46:51 -0400 Subject: [PATCH 05/17] Added support for setting the color of a meshlines plot overlay --- src/input_xml.F90 | 18 ++++++++++++++++++ src/plot.F90 | 21 +++++++++++++-------- src/plot_header.F90 | 1 + 3 files changed, 32 insertions(+), 8 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 3f61f8d764..5d9d54a976 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2990,6 +2990,24 @@ contains call fatal_error() end if + ! Check for color + if (check_for_node(node_meshlines, "color")) then + + ! Check and make sure 3 values are specified for RGB + if (get_arraysize_double(node_meshlines, "color") /= 3) then + message = "Bad RGB for meshlines color " & + // "in plot " // trim(to_str(pl % id)) + call fatal_error() + end if + + call get_node_array(node_meshlines, "color", & + pl % meshlines_color % rgb) + else + + pl % meshlines_color % rgb = (/ 0, 0, 0 /) + + 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) diff --git a/src/plot.F90 b/src/plot.F90 index 25c0c7123a..888b50dff7 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -191,6 +191,7 @@ contains logical :: in_mesh integer :: x, y ! pixel location + integer :: r, g, b ! RGB color for meshlines pixels integer :: xrange(2), yrange(2) ! range of pixel locations integer :: i, j ! loop indices integer :: ijk(3) @@ -207,6 +208,10 @@ contains m => meshes(pl % meshlines_id) + r = pl % meshlines_color % rgb(1) + g = pl % meshlines_color % rgb(2) + b = pl % meshlines_color % rgb(3) + select case (pl % basis) case(PLOT_BASIS_XY) xyz_ll_plot(1) = pl % origin(1) - pl % width(1) / 2.0 @@ -251,18 +256,18 @@ contains ! 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) + call set_pixel(img, x, yrange(1) + plus, r, g, b) + call set_pixel(img, x, yrange(2) + plus, r, g, b) + call set_pixel(img, x, yrange(1) - plus, r, g, b) + call set_pixel(img, x, yrange(2) - plus, r, g, b) 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) + call set_pixel(img, xrange(1) + plus, y, r, g, b) + call set_pixel(img, xrange(2) + plus, y, r, g, b) + call set_pixel(img, xrange(1) - plus, y, r, g, b) + call set_pixel(img, xrange(2) - plus, y, r, g, b) end do end do diff --git a/src/plot_header.F90 b/src/plot_header.F90 index 5a301f8809..e66bf9077e 100644 --- a/src/plot_header.F90 +++ b/src/plot_header.F90 @@ -27,6 +27,7 @@ module plot_header 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) :: meshlines_color ! Color for meshlines type(ObjectColor) :: not_found ! color for positions where no cell found type(ObjectColor), allocatable :: colors(:) ! colors of cells/mats end type ObjectPlot From cdcfa22c51679acc1aa7bc3aa724b7471a9215ad Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Wed, 20 Aug 2014 10:51:44 -0400 Subject: [PATCH 06/17] Updated plot documentation for colored meshlines --- docs/source/usersguide/input.rst | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 290713c35a..fcb979614e 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1189,7 +1189,7 @@ attributes or sub-elements. These are not used in "voxel" plots: Specifies the RGB color of the regions where no OpenMC cell can be found. Should be three integers separated by spaces. - *Default*: 0 0 0 (white) + *Default*: 0 0 0 (black) :col_spec: Any number of this optional tag may be included in each ```` element, @@ -1233,15 +1233,23 @@ attributes or sub-elements. These are not used in "voxel" plots: The special ``meshlines`` sub-element allows for plotting the boundaries of a tally mesh on top of a plot. Only one ``meshlines`` element is allowed per ``plot`` element, and it must contain as attributes or sub-elements a mesh - id and a linewidth: + id and a linewidth, as well as an optional color: :mesh: A single integer id number for the mesh specified on ``tallies.xml`` that should be plotted. :linewidth: - A single integer number of pixels to specify the linewidth of the mesh - boundaries. + A single integer number of pixels of linewidth to specify for the mesh + boundaries. Specifying this as 0 indicates that lines will be 1 pixel + thick, specifying 1 indicates 3 pixels thick, specifying 2 indicates + 5 pixels thick, etc. + + :color: + Specifies the custom color for the meshlines boundaries. Should be 3 + integers separated by spaces. This element is optional. + + *Default*: 0 0 0 (black) *Default*: None From 20ad38e52832a8b6b0f31703f216a4d7ee9eaec3 Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Sat, 30 Aug 2014 02:53:08 -0400 Subject: [PATCH 07/17] Added support for other slice orientations --- src/plot.F90 | 137 ++++++++++++++++++++++++++------------------------- 1 file changed, 69 insertions(+), 68 deletions(-) diff --git a/src/plot.F90 b/src/plot.F90 index 888b50dff7..97ead21cb1 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -190,14 +190,14 @@ contains type(Image) :: img logical :: in_mesh - integer :: x, y ! pixel location + integer :: out_, in_ ! pixel location integer :: r, g, b ! RGB color for meshlines pixels - integer :: xrange(2), yrange(2) ! range of pixel locations + integer :: outrange(2), inrange(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 + integer :: outer, inner real(8) :: frac real(8) :: width(3) ! real widths of the plot real(8) :: xyz_ll_plot(3) ! lower left xyz of plot image @@ -214,74 +214,75 @@ contains 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, r, g, b) - call set_pixel(img, x, yrange(2) + plus, r, g, b) - call set_pixel(img, x, yrange(1) - plus, r, g, b) - call set_pixel(img, x, yrange(2) - plus, r, g, b) - end do - end do - do y = yrange(2), yrange(1) - do plus = 0, pl % meshlines_width - call set_pixel(img, xrange(1) + plus, y, r, g, b) - call set_pixel(img, xrange(2) + plus, y, r, g, b) - call set_pixel(img, xrange(1) - plus, y, r, g, b) - call set_pixel(img, xrange(2) - plus, y, r, g, b) - end do - end do - - end if - end do - end do - + outer = 1 + inner = 2 case(PLOT_BASIS_XZ) - message = "Meshline plotting currently only implemented for basis xy" - call fatal_error() + outer = 1 + inner = 3 case(PLOT_BASIS_YZ) - message = "Meshline plotting currently only implemented for basis xy" - call fatal_error() + outer = 2 + inner = 3 end select + + xyz_ll_plot = pl % origin + xyz_ur_plot = pl % origin + + xyz_ll_plot(outer) = pl % origin(1) - pl % width(1) / 2.0 + xyz_ll_plot(inner) = pl % origin(2) - pl % width(2) / 2.0 + xyz_ur_plot(outer) = pl % origin(1) + pl % width(1) / 2.0 + xyz_ur_plot(inner) = pl % origin(2) + pl % width(2) / 2.0 + + 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 + do i = ijk_ll(outer), ijk_ur(outer) + do j = ijk_ll(inner), ijk_ur(inner) + ! check if we're in the mesh for this ijk + if (i > 0 .and. i <= m % dimension(outer) .and. & + j > 0 .and. j <= m % dimension(inner)) then + + ! get xyz's of lower left and upper right of this mesh cell + xyz_ll(outer) = m % lower_left(outer) + m % width(outer) * (i - 1) + xyz_ll(inner) = m % lower_left(inner) + m % width(inner) * (j - 1) + xyz_ur(outer) = m % lower_left(outer) + m % width(outer) * i + xyz_ur(inner) = m % lower_left(inner) + m % width(inner) * j + + ! map the xyz ranges to pixel ranges + + frac = (xyz_ll(outer) - xyz_ll_plot(outer)) / width(outer) + outrange(1) = int(frac * real(img % width, 8)) + frac = (xyz_ur(outer) - xyz_ll_plot(outer)) / width(outer) + outrange(2) = int(frac * real(img % width, 8)) + + frac = (xyz_ll(inner) - xyz_ll_plot(inner)) / width(inner) + inrange(1) = int(frac * real(img % height, 8)) + frac = (xyz_ur(inner) - xyz_ll_plot(inner)) / width(inner) + inrange(2) = int(frac * real(img % height, 8)) + + ! draw lines + do out_ = outrange(1), outrange(2) + do plus = 0, pl % meshlines_width + call set_pixel(img, out_, inrange(1) + plus, r, g, b) + call set_pixel(img, out_, inrange(2) + plus, r, g, b) + call set_pixel(img, out_, inrange(1) - plus, r, g, b) + call set_pixel(img, out_, inrange(2) - plus, r, g, b) + end do + end do + do in_ = inrange(1), inrange(2) + do plus = 0, pl % meshlines_width + call set_pixel(img, outrange(1) + plus, in_, r, g, b) + call set_pixel(img, outrange(2) + plus, in_, r, g, b) + call set_pixel(img, outrange(1) - plus, in_, r, g, b) + call set_pixel(img, outrange(2) - plus, in_, r, g, b) + end do + end do + + end if + end do + end do end subroutine draw_mesh_lines From f108e75be42ead7ba2f210826278df2090eef684 Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Sat, 30 Aug 2014 03:01:47 -0400 Subject: [PATCH 08/17] Switched to pointer for plot mesh --- src/input_xml.F90 | 17 +++++++++-------- src/plot.F90 | 4 ++-- src/plot_header.F90 | 3 ++- 3 files changed, 13 insertions(+), 11 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 5d9d54a976..8e43ab3973 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2682,8 +2682,9 @@ contains subroutine read_plots_xml() - integer i, j - integer n_cols, col_id, n_comp, n_masks, n_meshlines + integer :: i, j + integer :: n_cols, col_id, n_comp, n_masks, n_meshlines + integer :: meshid integer, allocatable :: iarray(:) logical :: file_exists ! does plots.xml file exist? character(MAX_LINE_LEN) :: filename ! absolute path to plots.xml @@ -2973,7 +2974,7 @@ contains ! 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) + call get_node_value(node_meshlines, "mesh", meshid) else message = "Must specify a mesh id for meshlines " // & "specification in plot " // trim(to_str(pl % id)) @@ -3009,16 +3010,16 @@ contains 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 + if (mesh_dict % has_key(meshid)) then + pl % meshlines_mesh => meshes(mesh_dict % get_key(meshid)) + if (meshes(meshid) % 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)) // & + message = "Could not find mesh " // & + trim(to_str(meshid)) // & " specified in meshlines for plot " // & trim(to_str(pl % id)) call fatal_error() diff --git a/src/plot.F90 b/src/plot.F90 index 97ead21cb1..f6a41c9513 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -168,7 +168,7 @@ contains end do ! Draw tally mesh boundaries on the image if requested - if (pl % meshlines_id > 0) call draw_mesh_lines(pl, img) + if (associated(pl % meshlines_mesh)) call draw_mesh_lines(pl, img) ! Write out the ppm to a file call output_ppm(pl,img) @@ -206,7 +206,7 @@ contains real(8) :: xyz_ur(3) ! upper right xyz type(StructuredMesh), pointer :: m => null() - m => meshes(pl % meshlines_id) + m => pl % meshlines_mesh r = pl % meshlines_color % rgb(1) g = pl % meshlines_color % rgb(2) diff --git a/src/plot_header.F90 b/src/plot_header.F90 index e66bf9077e..68eb89a11a 100644 --- a/src/plot_header.F90 +++ b/src/plot_header.F90 @@ -1,6 +1,7 @@ module plot_header use constants + use mesh_header, only: StructuredMesh implicit none @@ -25,8 +26,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(StructuredMesh), pointer :: meshlines_mesh => null() ! mesh to plot type(ObjectColor) :: meshlines_color ! Color for meshlines type(ObjectColor) :: not_found ! color for positions where no cell found type(ObjectColor), allocatable :: colors(:) ! colors of cells/mats From 9ef26276c82c54162b22927a9ac6ed73b7a8b813 Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Fri, 12 Sep 2014 10:55:19 -0400 Subject: [PATCH 09/17] Added support for plotting ufs, entropy, and cmfd meshes --- src/input_xml.F90 | 99 ++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 80 insertions(+), 19 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 6f1aa5fd01..6ff284152c 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -512,6 +512,7 @@ contains allocate(entropy_mesh) allocate(entropy_mesh % lower_left(3)) allocate(entropy_mesh % upper_right(3)) + allocate(entropy_mesh % width(3)) ! Copy values call get_node_array(node_entropy, "lower_left", & @@ -543,6 +544,11 @@ contains ! Copy dimensions call get_node_array(node_entropy, "dimension", entropy_mesh % dimension) + + ! Calculate width + entropy_mesh % width = (entropy_mesh % upper_right - & + entropy_mesh % lower_left) / entropy_mesh % dimension + end if ! Turn on Shannon entropy calculation @@ -2773,10 +2779,12 @@ contains integer :: i, j integer :: n_cols, col_id, n_comp, n_masks, n_meshlines integer :: meshid + integer :: i_mesh integer, allocatable :: iarray(:) logical :: file_exists ! does plots.xml file exist? character(MAX_LINE_LEN) :: filename ! absolute path to plots.xml character(MAX_LINE_LEN) :: temp_str + character(MAX_WORD_LEN) :: meshtype type(ObjectPlot), pointer :: pl => null() type(Node), pointer :: doc => null() type(Node), pointer :: node_plot => null() @@ -3060,11 +3068,11 @@ contains ! 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", meshid) + ! Check mesh type + if (check_for_node(node_meshlines, "meshtype")) then + call get_node_value(node_meshlines, "meshtype", meshtype) else - message = "Must specify a mesh id for meshlines " // & + message = "Must specify a meshtype for meshlines " // & "specification in plot " // trim(to_str(pl % id)) call fatal_error() end if @@ -3084,8 +3092,8 @@ contains ! Check and make sure 3 values are specified for RGB if (get_arraysize_double(node_meshlines, "color") /= 3) then - message = "Bad RGB for meshlines color " & - // "in plot " // trim(to_str(pl % id)) + message = "Bad RGB for meshlines color " // & + "in plot " // trim(to_str(pl % id)) call fatal_error() end if @@ -3097,22 +3105,75 @@ contains end if - ! Check if the specified tally mesh exists - if (mesh_dict % has_key(meshid)) then - pl % meshlines_mesh => meshes(mesh_dict % get_key(meshid)) - if (meshes(meshid) % type /= LATTICE_RECT) then - message = "Non-rectangular mesh specified in meshlines for" // & - " plot " // trim(to_str(pl % id)) + ! Set mesh based on type + select case (trim(meshtype)) + case ('ufs') + + if (.not. associated(ufs_mesh)) then + message = "No UFS mesh for meshlines on plot " // & + trim(to_str(pl % id)) call fatal_error() end if - else - message = "Could not find mesh " // & - trim(to_str(meshid)) // & - " specified in meshlines for plot " // & - trim(to_str(pl % id)) - call fatal_error() - end if + + pl % meshlines_mesh => ufs_mesh + + case ('cmfd') + + if (.not. cmfd_run) then + message = "Need CMFD run to plot CMFD mesh for meshlines " // & + "on plot " // trim(to_str(pl % id)) + call fatal_error() + end if + + i_mesh = cmfd_tallies(1) % & + filters(cmfd_tallies(1) % find_filter(FILTER_MESH)) % & + int_bins(1) + pl % meshlines_mesh => meshes(i_mesh) + + case ('entropy') + if (.not. associated(entropy_mesh)) then + message = "No entropy mesh for meshlines on plot " // & + trim(to_str(pl % id)) + call fatal_error() + end if + + pl % meshlines_mesh => entropy_mesh + + case ('tally') + + ! Ensure that there is a mesh id if the type is tally + if (check_for_node(node_meshlines, "id")) then + call get_node_value(node_meshlines, "id", meshid) + else + message = "Must specify a mesh id for meshlines tally mesh" // & + "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(meshid)) then + pl % meshlines_mesh => meshes(mesh_dict % get_key(meshid)) + if (meshes(meshid) % 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 mesh " // & + trim(to_str(meshid)) // & + " specified in meshlines for plot " // & + trim(to_str(pl % id)) + call fatal_error() + end if + + case default + message = "Invalid type for meshlines on plot " // & + trim(to_str(pl % id)) // ": " // trim(meshtype) + call fatal_error() + end select + + end select end if From 378aa92e668e5a0f064b09e96f22fb0c2fe3f28c Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Fri, 12 Sep 2014 10:58:15 -0400 Subject: [PATCH 10/17] Allowed read_tallies_xml to run in plotting mode, only to read mesh info --- src/input_xml.F90 | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 6ff284152c..42b684abf6 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -36,7 +36,7 @@ contains if (run_mode /= MODE_PLOTTING) call read_cross_sections_xml() call read_geometry_xml() call read_materials_xml() - if (run_mode /= MODE_PLOTTING) call read_tallies_xml() + call read_tallies_xml() if (cmfd_run) call configure_cmfd() end subroutine read_input_xml @@ -1909,7 +1909,7 @@ contains end if ! Allocate tally array - if (n_user_tallies > 0) then + if (n_user_tallies > 0 .and. run_mode /= MODE_PLOTTING) then call add_tallies("user", n_user_tallies) end if @@ -2058,6 +2058,9 @@ contains call mesh_dict % add_key(m % id, i) end do + ! We only need the mesh info for plotting + if (run_mode == MODE_PLOTTING) return + ! ========================================================================== ! READ TALLY DATA From d5ce05cf363b654c66fc6ef6944116c57fcfea57 Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Fri, 12 Sep 2014 11:04:05 -0400 Subject: [PATCH 11/17] Updates documentation for meshline plotting --- docs/source/usersguide/input.rst | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 5cc595ed50..44f67e62c1 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1298,11 +1298,16 @@ attributes or sub-elements. These are not used in "voxel" plots: The special ``meshlines`` sub-element allows for plotting the boundaries of a tally mesh on top of a plot. Only one ``meshlines`` element is allowed per ``plot`` element, and it must contain as attributes or sub-elements a mesh - id and a linewidth, as well as an optional color: + type and a linewidth, as well as an optional color: - :mesh: + :meshtype: + The type of the mesh to be plotted. Valid options are "tally", "entropy", + "ufs", and "cmfd". If plotting "tally" meshes, the id of the mesh to plot + must be specified with the ``id`` sub-element. + + :id: A single integer id number for the mesh specified on ``tallies.xml`` that - should be plotted. + should be plotted. This element is only required for ``meshtype="tally"``. :linewidth: A single integer number of pixels of linewidth to specify for the mesh @@ -1318,9 +1323,6 @@ attributes or sub-elements. These are not used in "voxel" plots: *Default*: None - .. warning:: Meshline plotting is currently only implemented for plots with - an "xy" basis. - ------------------------------ CMFD Specification -- cmfd.xml ------------------------------ From e2c9d7550793ae01cb732987e2fd40dac5ddd4f8 Mon Sep 17 00:00:00 2001 From: Nicholas Horelik Date: Tue, 16 Sep 2014 22:42:41 -0400 Subject: [PATCH 12/17] Clarified documentation --- docs/source/usersguide/input.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 44f67e62c1..c97ab8ac1e 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1295,10 +1295,10 @@ attributes or sub-elements. These are not used in "voxel" plots: *Default*: None :meshlines: - The special ``meshlines`` sub-element allows for plotting the boundaries of + The ``meshlines`` sub-element allows for plotting the boundaries of a tally mesh on top of a plot. Only one ``meshlines`` element is allowed per ``plot`` element, and it must contain as attributes or sub-elements a mesh - type and a linewidth, as well as an optional color: + type and a linewidth. Optionally, a color may be specified for the overlay: :meshtype: The type of the mesh to be plotted. Valid options are "tally", "entropy", @@ -1317,7 +1317,7 @@ attributes or sub-elements. These are not used in "voxel" plots: :color: Specifies the custom color for the meshlines boundaries. Should be 3 - integers separated by spaces. This element is optional. + integers separated by whitespace. This element is optional. *Default*: 0 0 0 (black) From 224dcb48c878df10934ddf47c033f1688e2891f1 Mon Sep 17 00:00:00 2001 From: Nicholas Horelik Date: Tue, 16 Sep 2014 22:55:28 -0400 Subject: [PATCH 13/17] Tidied up select case --- src/input_xml.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 42b684abf6..f852d2243d 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3061,11 +3061,8 @@ contains end if select case(n_meshlines) - case default - message = "Mutliple meshlines" // & - " specified in plot " // trim(to_str(pl % id)) - call fatal_error() case (0) + ! Skip if no meshlines are specified case (1) ! Get pointer to meshlines @@ -3176,7 +3173,10 @@ contains call fatal_error() end select - + case default + message = "Mutliple meshlines" // & + " specified in plot " // trim(to_str(pl % id)) + call fatal_error() end select end if From 1c939cb2403d429b2af0f8467ff6b6347ee50b07 Mon Sep 17 00:00:00 2001 From: Nicholas Horelik Date: Tue, 16 Sep 2014 22:57:22 -0400 Subject: [PATCH 14/17] Added intent --- src/plot.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/plot.F90 b/src/plot.F90 index f6a41c9513..91b2520eb4 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -186,8 +186,8 @@ contains !=============================================================================== subroutine draw_mesh_lines(pl, img) - type(ObjectPlot), pointer :: pl - type(Image) :: img + type(ObjectPlot), pointer, intent(in) :: pl + type(Image), intent(inout) :: img logical :: in_mesh integer :: out_, in_ ! pixel location From 2c433d3061a7cd99fc33ba181ab95b634fb1ef4a Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Mon, 6 Oct 2014 10:33:22 -0400 Subject: [PATCH 15/17] Fixed entropy mesh width calculation --- src/eigenvalue.F90 | 8 ++++---- src/input_xml.F90 | 6 ++++++ 2 files changed, 10 insertions(+), 4 deletions(-) diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 37edb922c4..cd95568849 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -536,12 +536,12 @@ contains m % n_dimension = 3 allocate(m % dimension(3)) m % dimension = n + + ! determine width + m % width = (m % upper_right - m % lower_left) / m % dimension + end if - ! allocate and determine width - allocate(m % width(3)) - m % width = (m % upper_right - m % lower_left) / m % dimension - ! allocate p allocate(entropy_p(1, m % dimension(1), m % dimension(2), & m % dimension(3))) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index a55fc363f4..522427807d 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3135,6 +3135,12 @@ contains call fatal_error() end if + if (.not. allocated(entropy_mesh % dimension)) then + message = "No dimension specified on entropy mesh for " // + "meshlines on plot " // trim(to_str(pl % id)) + call fatal_error() + end if + pl % meshlines_mesh => entropy_mesh case ('tally') From e67211b09cfe994ae51bb88135f07a5bd58f4c74 Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Tue, 7 Oct 2014 00:56:16 -0400 Subject: [PATCH 16/17] Added missing ampersand --- src/input_xml.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 522427807d..e84f072ea8 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3136,7 +3136,7 @@ contains end if if (.not. allocated(entropy_mesh % dimension)) then - message = "No dimension specified on entropy mesh for " // + message = "No dimension specified on entropy mesh for " // & "meshlines on plot " // trim(to_str(pl % id)) call fatal_error() end if From 2a382d70908c153098595f128d21b13ea55baa6d Mon Sep 17 00:00:00 2001 From: Nick Horelik Date: Tue, 7 Oct 2014 00:59:39 -0400 Subject: [PATCH 17/17] Updated rnc file for meshlines plot xml changes --- src/relaxng/plots.rnc | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/relaxng/plots.rnc b/src/relaxng/plots.rnc index 3c29ba4e0b..a7c657864d 100644 --- a/src/relaxng/plots.rnc +++ b/src/relaxng/plots.rnc @@ -26,6 +26,14 @@ element plots { attribute components { list { xsd:int+ } }) & (element background { list { xsd:int+ } } | attribute background { list { xsd:int+ } }) + }* & + element meshlines { + (element meshtype { ( "tally" | "entropy" | "ufs" | "cmfd" ) } | + attribute meshtype { ( "tally" | "entropy" | "ufs" | "cmfd" ) }) & + (element id { xsd:int } | attribute id { xsd:int })? & + (element linewidth { xsd:int } | attribute linewidth { xsd:int }) & + (element color { list { xsd:int+ } } | + attribute color { list { xsd:int+ } })? }* }* }