Merge pull request #1089 from pshriwise/voxel_script

Voxel script

.gitignore

@@ -69,6 +69,7 @@ scripts/G4EMLOW*/
 # Images
 *.ppm
 *.voxel
+*.vti
 
 # PyCharm project configuration files
 .idea

scripts/openmc-voxel-to-silovtk

@@ -1,85 +0,0 @@
-#!/usr/bin/env python3
-
-import struct
-import sys
-from argparse import ArgumentParser
-
-import numpy as np
-import h5py
-
-
-def main():
-    # Process command line arguments
-    parser = ArgumentParser()
-    parser.add_argument('voxel_file', help='Path to voxel file')
-    parser.add_argument('-o', '--output', action='store',
-                        default='plot', help='Path to output SILO or VTK file.')
-    parser.add_argument('-s', '--silo', action='store_true',
-                        default=False, help='Flag to convert to SILO instead of VTK.')
-    args = parser.parse_args()
-
-    # Read data from voxel file
-    fh = h5py.File(args.voxel_file, 'r')
-    dimension = fh.attrs['num_voxels']
-    width = fh.attrs['voxel_width']
-    lower_left = fh.attrs['lower_left']
-    voxel_data = fh['data'].value
-
-    nx, ny, nz = dimension
-    upper_right = lower_left + width*dimension
-
-    if not args.silo:
-        import vtk
-
-        grid = vtk.vtkImageData()
-        grid.SetDimensions(nx+1, ny+1, nz+1)
-        grid.SetOrigin(*lower_left)
-        grid.SetSpacing(*width)
-
-        data = vtk.vtkDoubleArray()
-        data.SetName("id")
-        data.SetNumberOfTuples(nx*ny*nz)
-        for x in range(nx):
-            sys.stdout.write(" {}%\r".format(int(x/nx*100)))
-            sys.stdout.flush()
-            for y in range(ny):
-                for z in range(nz):
-                    i = z*nx*ny + y*nx + x
-                    data.SetValue(i, voxel_data[x, y, z])
-        grid.GetCellData().AddArray(data)
-
-        writer = vtk.vtkXMLImageDataWriter()
-        if vtk.vtkVersion.GetVTKMajorVersion() > 5:
-            writer.SetInputData(grid)
-        else:
-            writer.SetInput(grid)
-        if not args.output.endswith(".vti"):
-            args.output += ".vti"
-        writer.SetFileName(args.output)
-        writer.Write()
-
-    else:
-        import silomesh
-
-        if not args.output.endswith(".silo"):
-            args.output += ".silo"
-        silomesh.init_silo(args.output)
-        meshparams = list(map(int, dimension)) + list(map(float, lower_left)) + \
-                     list(map(float, upper_right))
-        silomesh.init_mesh('plot', *meshparams)
-        silomesh.init_var("id")
-        for x in range(nx):
-            sys.stdout.write(" {}%\r".format(int(x/nx*100)))
-            sys.stdout.flush()
-            for y in range(ny):
-                for z in range(nz):
-                    silomesh.set_value(float(voxel_data[x, y, z]),
-                                       x + 1, y + 1, z + 1)
-        print()
-        silomesh.finalize_var()
-        silomesh.finalize_mesh()
-        silomesh.finalize_silo()
-
-
-if __name__ == '__main__':
-    main()

scripts/openmc-voxel-to-vtk

@@ -0,0 +1,58 @@
+#!/usr/bin/env python3
+
+import struct
+import sys
+from argparse import ArgumentParser
+
+import numpy as np
+import h5py
+import vtk
+
+
+def main():
+    # Process command line arguments
+    parser = ArgumentParser()
+    parser.add_argument('voxel_file', help='Path to voxel file')
+    parser.add_argument('-o', '--output', action='store',
+                        default='plot', help='Path to output VTK file.')
+    args = parser.parse_args()
+
+    # Read data from voxel file
+    fh = h5py.File(args.voxel_file, 'r')
+    dimension = fh.attrs['num_voxels']
+    width = fh.attrs['voxel_width']
+    lower_left = fh.attrs['lower_left']
+    voxel_data = fh['data'].value
+
+    nx, ny, nz = dimension
+    upper_right = lower_left + width*dimension
+    
+    grid = vtk.vtkImageData()
+    grid.SetDimensions(nx+1, ny+1, nz+1)
+    grid.SetOrigin(*lower_left)
+    grid.SetSpacing(*width)
+    
+    data = vtk.vtkDoubleArray()
+    data.SetName("id")
+    data.SetNumberOfTuples(nx*ny*nz)
+    for x in range(nx):
+        sys.stdout.write(" {}%\r".format(int(x/nx*100)))
+        sys.stdout.flush()
+        for y in range(ny):
+            for z in range(nz):
+                i = z*nx*ny + y*nx + x
+                data.SetValue(i, voxel_data[x, y, z])
+    grid.GetCellData().AddArray(data)
+
+    writer = vtk.vtkXMLImageDataWriter()
+    if vtk.vtkVersion.GetVTKMajorVersion() > 5:
+        writer.SetInputData(grid)
+    else:
+        writer.SetInput(grid)
+    if not args.output.endswith(".vti"):
+        args.output += ".vti"
+    writer.SetFileName(args.output)
+    writer.Write()
+
+if __name__ == '__main__':
+    main()

setup.py

@@ -64,7 +64,7 @@ kwargs = {
     # Optional dependencies
     'extras_require': {
         'test': ['pytest', 'pytest-cov'],
-        'vtk': ['vtk', 'silomesh'],
+        'vtk': ['vtk'],
     },
 }
 

tests/regression_tests/plot_voxel/geometry.xml

@@ -0,0 +1,13 @@
+<?xml version="1.0"?>
+<geometry>
+
+  <surface id="1" type="z-cylinder" coeffs="0 0 2" />
+  <surface id="2" type="z-cylinder" coeffs="0 0 5" />
+  <surface id="3" type="z-cylinder" coeffs="0 0 10" boundary="vacuum" />
+  <surface id="4" type="y-plane" coeffs="5" boundary="vacuum" />
+
+  <cell id="1" material="1" region="  -1 -4" />
+  <cell id="2" material="3" region="1 -2 -4" />
+  <cell id="3" material="2" region="2 -3 -4" />
+
+</geometry>

tests/regression_tests/plot_voxel/materials.xml

@@ -0,0 +1,19 @@
+<?xml version="1.0"?>
+<materials>
+
+  <material id="1">
+    <density value="4.5" units="g/cc" />
+    <nuclide name="U235" ao="1.0" />
+  </material>
+
+  <material id="2">
+    <density value="4.5" units="g/cc" />
+    <nuclide name="U235" ao="1.0" />
+  </material>
+
+  <material id="3">
+    <density value="4.5" units="g/cc" />
+    <nuclide name="U235" ao="1.0" />
+  </material>
+
+</materials>

tests/regression_tests/plot_voxel/plots.xml

@@ -0,0 +1,10 @@
+<?xml version="1.0"?>
+<plots>
+
+  <plot id="4" type="voxel">
+    <pixels>50 50 10</pixels>
+    <origin>0. 0. 0.</origin>
+    <width>20 20 10</width>
+  </plot>
+
+</plots>

tests/regression_tests/plot_voxel/results_true.dat

@@ -0,0 +1 @@
+20a0c3598bb698efa4bba65c5e55d71f4385e5b1bcf0067964e45001c12f5c0a729935a96409c760dbd3ec439ae6c676fce77d6e578b41f8f3e0ac68c9277acb

tests/regression_tests/plot_voxel/settings.xml

@@ -0,0 +1,13 @@
+<?xml version="1.0"?>
+<settings>
+
+  <run_mode>plot</run_mode>
+
+  <mesh id="1">
+    <dimension>5 4 3</dimension>
+    <lower_left>-10 -10 -10</lower_left>
+    <upper_right>10 10 10</upper_right>
+  </mesh>
+  <entropy_mesh>1</entropy_mesh>
+
+</settings>

tests/regression_tests/plot_voxel/test.py

@@ -0,0 +1,62 @@
+import glob
+import hashlib
+import os
+from subprocess import check_call
+
+import h5py
+import openmc
+import pytest
+
+from tests.testing_harness import TestHarness
+from tests.regression_tests import config
+
+
+vtk = pytest.importorskip('vtk')
+class PlotVoxelTestHarness(TestHarness):
+    """Specialized TestHarness for running OpenMC voxel plot tests."""
+    def __init__(self, plot_names):
+        super().__init__(None)
+        self._plot_names = plot_names
+
+    def _run_openmc(self):
+        openmc.plot_geometry(openmc_exec=config['exe'])
+
+        check_call(['../../../scripts/openmc-voxel-to-vtk'] +
+                   glob.glob('plot_4.h5'))
+                
+    def _test_output_created(self):
+        """Make sure *.ppm has been created."""
+        for fname in self._plot_names:
+            assert os.path.exists(fname), 'Plot output file does not exist.'
+
+    def _cleanup(self):
+        super()._cleanup()
+        for fname in self._plot_names:
+            if os.path.exists(fname):
+                os.remove(fname)
+
+    def _get_results(self):
+        """Return a string hash of the plot files."""
+        outstr = bytes()
+
+        for fname in self._plot_names:
+            if fname.endswith('.h5'):
+                # Add voxel data to results
+                with h5py.File(fname, 'r') as fh:
+                    outstr += fh.attrs['filetype']
+                    outstr += fh.attrs['num_voxels'].tostring()
+                    outstr += fh.attrs['lower_left'].tostring()
+                    outstr += fh.attrs['voxel_width'].tostring()
+                    outstr += fh['data'].value.tostring()
+
+        # Hash the information and return.
+        sha512 = hashlib.sha512()
+        sha512.update(outstr)
+        outstr = sha512.hexdigest()
+
+        return outstr
+
+
+def test_plot_voxel():
+    harness = PlotVoxelTestHarness(('plot_4.h5', 'plot.vti'))
+    harness.main()

tools/ci/travis-install.sh

@@ -23,8 +23,13 @@ fi
 # Build and install OpenMC executable
 python tools/ci/travis-install.py
 
-# Install Python API in editable mode
-pip install -e .[test]
+if [[ $TRAVIS_PYTHON_VERSION == "3.7" ]]; then
+    # Install Python API in editable mode
+    pip install -e .[test]
+else
+    # Conditionally install vtk
+    pip install -e .[test,vtk]
+fi
 
 # For uploading to coveralls
 pip install coveralls