Merge pull request #306 from nhorelik/plot_meshlines

Added capability to plot mesh boundaries on top of slice plots

docs/source/usersguide/input.rst

@@ -1261,7 +1261,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 ``<plot>`` element,
@@ -1301,6 +1301,35 @@ attributes or sub-elements.  These are not used in "voxel" plots:
 
     *Default*: None
 
+  :meshlines:
+    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.  Optionally, a color may be specified for the overlay:
+
+    :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. This element is only required for ``meshtype="tally"``.
+
+    :linewidth:
+      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 whitespace.  This element is optional.
+    
+      *Default*: 0 0 0 (black)
+
+    *Default*: None
+
 .. _usersguide_cmfd:
 
 ------------------------------

examples/lattice/simple/plots.xml

@@ -5,6 +5,7 @@
     <origin>0. 0. 0.</origin>
     <width>4.0 4.0</width>
     <pixels>400 400</pixels>
+    <!-- <meshlines mesh="1" linewidth="2" color="0 255 0"/> -->
   </plot>
 
 </plots>

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)))

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
@@ -517,6 +517,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", &
@@ -548,6 +549,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
@@ -1895,7 +1901,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
 
@@ -2043,6 +2049,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
 
@@ -2767,20 +2776,25 @@ contains
 
   subroutine read_plots_xml()
 
-    integer i, j
-    integer n_cols, col_id, n_comp, n_masks
+    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()
     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"
@@ -3032,6 +3046,144 @@ 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 (0)
+            ! Skip if no meshlines are specified
+          case (1)
+
+            ! Get pointer to meshlines
+            call get_list_item(node_meshline_list, 1, node_meshlines)
+            
+            ! Check mesh type
+            if (check_for_node(node_meshlines, "meshtype")) then
+              call get_node_value(node_meshlines, "meshtype", meshtype)
+            else
+              message = "Must specify a meshtype 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 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
+
+            ! 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
+              
+              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
+              
+              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')
+
+              ! 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
+
+          case default
+            message = "Mutliple meshlines" // &
+                 " specified in plot " // trim(to_str(pl % id))
+            call fatal_error()
+        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)

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 (associated(pl % meshlines_mesh)) call draw_mesh_lines(pl, img)
+
     ! Write out the ppm to a file
     call output_ppm(pl,img)
 
@@ -180,6 +181,111 @@ contains
 
   end subroutine create_ppm
 
+!===============================================================================
+! DRAW_MESH_LINES draws mesh line boundaries on an image
+!===============================================================================
+  subroutine draw_mesh_lines(pl, img)
+  
+    type(ObjectPlot), pointer, intent(in)    :: pl
+    type(Image),               intent(inout) :: img
+    
+    logical :: in_mesh
+    integer :: out_, in_  ! pixel location
+    integer :: r, g, b    ! RGB color for meshlines pixels
+    integer :: outrange(2), inrange(2) ! range of pixel locations
+    integer :: i, j       ! loop indices
+    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
+    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 => pl % meshlines_mesh
+    
+    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)
+        outer = 1
+        inner = 2
+      case(PLOT_BASIS_XZ)
+        outer = 1
+        inner = 3
+      case(PLOT_BASIS_YZ)
+        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
+
 !===============================================================================
 ! OUTPUT_PPM writes out a previously generated image to a PPM file
 !===============================================================================

src/plot_header.F90

@@ -1,6 +1,7 @@
 module plot_header
 
   use constants
+  use mesh_header,  only: StructuredMesh
 
   implicit none
 
@@ -25,6 +26,9 @@ 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_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
   end type ObjectPlot

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+ } })?
     }*
   }*
 }