Struct Mesh

pub struct Mesh {Show 18 fields
    pub id: u32,
    pub geometry: Geometry,
    pub particle: Particle,
    pub imesh: Vec<f64>,
    pub iints: usize,
    pub jmesh: Vec<f64>,
    pub jints: usize,
    pub kmesh: Vec<f64>,
    pub kints: usize,
    pub emesh: Vec<f64>,
    pub eints: usize,
    pub tmesh: Vec<f64>,
    pub tints: usize,
    pub origin: [f64; 3],
    pub axs: [f64; 3],
    pub vec: [f64; 3],
    pub voxels: Vec<Voxel>,
    pub format: Format,
}

Expand description

Common data structure representing a mesh tally

Mesh attributes correspond closely to MCNP input cards for consistency and intuitive use. Units are unchanged from the MCNP defaults.

All parsed output formats are stored as a Mesh to provide a common interface for all post-processing operations. For example: conversion to VTK formats, weight window generation, data extraction, etc…

§Terminology notes

§I, J, K generics

Coordinate systems use different names i.e. (X,Y,Z) and (R,Z,Theta).

The generic (I,J,K) are used to represent all systems, in keeping with MCNP user manuals.

§Groups

A full set of (I,J,K) voxels are repeated for every time and energy bin in the mesh.

It is also possible to have additional ‘Total’ groups if the EMESH or TMESH cards contain multiple bins.

Time and energy bins are therefore often mapped to Group::Value and Group::Total variants to be explicit when absolute values are handled (See Group for details).

§Examples

§Reading meshtal files

Basic reading of files is very simple, regardless of format.

// Extract all meshes from a file into a Vec<Mesh>
let mesh_list = read("/path/to/meshtal.msht").unwrap();
let mesh_list = read("/path/to/runtape.r.h5").unwrap();

// Extract just one target mesh from a file into a single Mesh
let mesh = read_target("/path/to/meshtal.msht", 104).unwrap();
let mesh = read_target("/path/to/runtape.r.h5", 104).unwrap();

All the parsing and interpretation are done for you, and the data are in a common Mesh type. This means that all Mesh methods are available for any format mesh of any geometry type.

Fields§

§id: u32

Mesh tally number e.g fmesh104 => id = 104

§geometry: Geometry

Mesh geometry type, usually rectangular for MCNP default

§particle: Particle

Name of the particle type

§imesh: Vec<f64>

i mesh boundaries

§iints: usize

Number of voxels in i

§jmesh: Vec<f64>

j mesh boundaries

§jints: usize

Number of voxels in j

§kmesh: Vec<f64>

k mesh boundaries

§kints: usize

Number of voxels in j

§emesh: Vec<f64>

Energy bins

§eints: usize

Number of energy bins, EXCLUDING ‘total’ group

§tmesh: Vec<f64>

Time bins [shakes]

§tints: usize

Number of time bins, EXCLUDING ‘total’ group

§origin: [f64; 3]

ORIGIN card, [0.0, 0.0, 0.0] for MCNP default

§axs: [f64; 3]

AXS card, [0.0, 0.0, 1.0] for MCNP default

§vec: [f64; 3]

VEC card, [1.0, 0.0, 0.0] for MCNP default

§voxels: Vec<Voxel>

List of every Voxel in the mesh

§format: Format

Output format from which this mesh was derived

Implementations§

§

impl Mesh

Common methods

pub fn new(id: u32) -> Self

Initialise new mesh with known id

The id is the tally number used on the FMESH card in the input deck. For example, Fmesh204:n => id=204. This will initialise a Mesh with all of the default values.

pub fn scale(&mut self, factor: f64)

Multiply all voxel results by a constant factor

Uncertanties are relative and are therfore unaffected.

// Simple test mesh with three voxels
let mut mesh = Mesh {
    voxels: vec![
        Voxel {result: 12.0, ..Default::default()},
        Voxel {result: 18.0, ..Default::default()},
        Voxel {result:  4.0, ..Default::default()},
    ],
    ..Default::default()
};

// Scale results by a factor of x0.5 (errors relative, thus unchanged)
mesh.scale(0.5);
assert_eq!(mesh.voxels[0].result, 6.0);
assert_eq!(mesh.voxels[1].result, 9.0);
assert_eq!(mesh.voxels[2].result, 2.0);

pub fn translate(&mut self, x: f64, y: f64, z: f64)

Translate all coordinates by (x, y, z)

Simply updates the relevant origin coordiantes and mesh geometry bounds using the cartesian values provided. For cylindrical meshes the voxel bounds will be unaffected.

pub fn n_ebins(&self) -> usize

Returns the number of energy bins

This will include the Total bin in the count for tallies with multiple energy bins defined on the EMESH card.

Card values	#Groups
None	1
0.0 100	1
0.0 20 100	3 (2 + ‘Total’ bin)
0.0 20 100 T	3 (2 + ‘Total’ bin)

pub fn n_tbins(&self) -> usize

Returns the number of time bins

This will include the ‘Total’ bin in the count for tallies with multiple time bins defined on the TMESH card. Defaults to 1 Total bin if none are explicitly defined.

Card values	#Groups
None	1
0.0 1e16	1
0.0 1e16 1e36	3 (2 + ‘Total’ bin)
0.0 1e16 1e36 T	3 (2 + ‘Total’ bin)

pub fn maximum(&self) -> (f64, f64)

Find the maximum (value, error) in the mesh

Will panic! when the mesh contains no voxel data. Use try_maximum() to handle this case explicitly.

For example:

let mesh = Mesh {
    voxels: vec![
        Voxel{index: 0, result: 1.0, error: 0.1},
        Voxel{index: 1, result: 2.0, error: 0.2},
        Voxel{index: 2, result: 3.0, error: 0.3},
    ],
    ..Default::default()
};

assert_eq!(mesh.maximum(), (3.0, 0.3));

To find the maximum of any collection of voxels, the VoxelSliceExt trait is implemented for convenience.

For example:

let voxels = vec![
        Voxel{index: 0, result: 1.0, error: 0.1},
        Voxel{index: 1, result: 2.0, error: 0.2},
        Voxel{index: 2, result: 3.0, error: 0.3} ];

assert_eq!(voxels.maximum_voxel().unwrap()       , &voxels[2]);
assert_eq!(voxels.maximum_result_error().unwrap(), (3.0, 0.3));

pub fn try_maximum(&self) -> Result<(f64, f64), Error>

Find the maximum (value, error) in the mesh

Will fail if mesh.voxels is empty. The maximum() will simply allow this case to panic!.

For example:

let mut mesh = Mesh {
    voxels: vec![
        Voxel{index: 0, result: 1.0, error: 0.1},
        Voxel{index: 1, result: 2.0, error: 0.2},
        Voxel{index: 2, result: 3.0, error: 0.3},
    ],
    ..Default::default()
};

// Successful example
assert_eq!(mesh.try_maximum().unwrap(), (3.0, 0.3));

// Failure example for empty mesh.voxels
mesh.voxels.clear();
assert!(mesh.try_maximum().is_err());

To find the maximum of any collection of voxels, the VoxelSliceExt trait is implemented for convenience.

For example:

let voxels = vec![
        Voxel{index: 0, result: 1.0, error: 0.1},
        Voxel{index: 1, result: 2.0, error: 0.2},
        Voxel{index: 2, result: 3.0, error: 0.3} ];

assert_eq!(voxels.maximum_voxel().unwrap()       , &voxels[2]);
assert_eq!(voxels.maximum_result_error().unwrap(), (3.0, 0.3));

pub fn minimum(&self) -> (f64, f64)

Find the minimum (value, error) in the mesh

Will panic! when the mesh contains no voxel data. Use try_minimum() to handle this case explicitly.

For example:

let mesh = Mesh {
    voxels: vec![
        Voxel{index: 0, result: 1.0, error: 0.1},
        Voxel{index: 1, result: 2.0, error: 0.2},
        Voxel{index: 2, result: 3.0, error: 0.3},
    ],
    ..Default::default()
};

assert_eq!(mesh.minimum(), (1.0, 0.1));

To find the minimum of any collection of voxels, the VoxelSliceExt trait is implemented for convenience.

For example:

let voxels = vec![
        Voxel{index: 0, result: 1.0, error: 0.1},
        Voxel{index: 1, result: 2.0, error: 0.2},
        Voxel{index: 2, result: 3.0, error: 0.3} ];

assert_eq!(voxels.minimum_voxel().unwrap()       , &voxels[0]);
assert_eq!(voxels.minimum_result_error().unwrap(), (1.0, 0.1));

pub fn try_minimum(&self) -> Result<(f64, f64), Error>

Find the minimum (value, error) in the mesh

Will fail if mesh.voxels is empty. The minimum() will simply allow this case to panic!.

For example:

let mut mesh = Mesh {
    voxels: vec![
        Voxel{index: 0, result: 1.0, error: 0.1},
        Voxel{index: 1, result: 2.0, error: 0.2},
        Voxel{index: 2, result: 3.0, error: 0.3},
    ],
    ..Default::default()
};

// Successful example
assert_eq!(mesh.try_minimum().unwrap(), (1.0, 0.1));

// Failure example for empty mesh.voxels
mesh.voxels.clear();
assert!(mesh.try_minimum().is_err());

To find the minimum of any collection of voxels, the VoxelSliceExt trait is implemented for convenience.

For example:

let voxels = vec![
        Voxel{index: 0, result: 1.0, error: 0.1},
        Voxel{index: 1, result: 2.0, error: 0.2},
        Voxel{index: 2, result: 3.0, error: 0.3} ];

assert_eq!(voxels.minimum_voxel().unwrap()       , &voxels[0]);
assert_eq!(voxels.minimum_result_error().unwrap(), (1.0, 0.1));

pub fn average(&self) -> (f64, f64)

Find the averege (value, error) in the mesh

Will panic! when the mesh contains no voxel data. Use try_average() to handle this case explicitly.

For example:

let mesh = Mesh {
    voxels: vec![
        Voxel{index: 0, result: 1.0, error: 0.1},
        Voxel{index: 1, result: 2.0, error: 0.2},
        Voxel{index: 2, result: 3.0, error: 0.3},
    ],
    ..Default::default()
};

assert_eq!(mesh.average(), (2.0, 0.49497474683058323));

To find the average of any collection of voxels, the VoxelSliceExt trait is implemented for convenience.

For example:

let voxels = vec![
        Voxel{index: 0, result: 1.0, error: 0.1},
        Voxel{index: 1, result: 2.0, error: 0.2},
        Voxel{index: 2, result: 3.0, error: 0.3} ];

assert_eq!(voxels.average_result_error().unwrap(), (2.0, 0.49497474683058323));

pub fn try_average(&self) -> Result<(f64, f64), Error>

Find the average (value, error) in the mesh

Will fail if mesh.voxels is empty. The average() will simply allow this case to panic!.

For example:

let mut mesh = Mesh {
    voxels: vec![
        Voxel{index: 0, result: 1.0, error: 0.1},
        Voxel{index: 1, result: 2.0, error: 0.2},
        Voxel{index: 2, result: 3.0, error: 0.3},
    ],
    ..Default::default()
};

// Successful example
assert_eq!(mesh.try_average().unwrap(), (2.0, 0.49497474683058323));

// Failure example for empty mesh.voxels
mesh.voxels.clear();
assert!(mesh.try_average().is_err());

To find the average of any collection of voxels, the VoxelSliceExt trait is implemented for convenience.

For example:

let voxels = vec![
        Voxel{index: 0, result: 1.0, error: 0.1},
        Voxel{index: 1, result: 2.0, error: 0.2},
        Voxel{index: 2, result: 3.0, error: 0.3} ];

assert_eq!(voxels.average_result_error().unwrap(), (2.0, 0.49497474683058323));

§

impl Mesh

Point method implementations for the Mesh type

pub fn find_point_data( &self, point: Point, boundary: BoundaryTreatment, ) -> Option<(f64, f64)>

Find the result at a Point

Results are averaged between adjacent voxels when the point is on a boundary. Points outside the mesh return None.

A small tolerance is granted for detecting points on boundaries. By default this is set to within 0.01% of the voxel length in each axis.

For example, for a voxel spanning 0.0 - 1.0 in the x-axis, a Point with x = 0.999 is considered to be on the boundary. The result will therefore be the avaerage of this and the appropriate adjacent voxel.

pub fn find_points_data( &self, points: &[Point], boundary: BoundaryTreatment, ) -> Vec<Option<(f64, f64)>>

Find the results for a list of Points

Equivalent to looping over find_point_data() in a loop for multiple points, collecting the results to a vector.

Results are averaged between adjacent voxels when the point is on a boundary. Points outside the mesh return None.

See find_point_data() for details.

pub fn find_point_voxels( &self, point: Point, boundary: BoundaryTreatment, ) -> Result<Vec<Voxel>, Error>

Find all relevant voxels corresponding to a Point

Try to find all adjacent voxels that a Point touches.

A small tolerance may be granted for detecting points on boundaries. For example, tol=0.01 would consider points to be on a boudary if within 1% of the total voxel length for each axis.

§

impl Mesh

Voxels and voxel slicing

pub fn n_voxels(&self) -> usize

Returns the number of voxels

pub fn n_voxels_expected(&self) -> usize

Calculate how many voxels there should be

Useful for common sense checking against the value returned by a mesh.voxels.len() call. This is the product of the number of energy bins, time bins, and fine mesh bins in each dimension.

// Generate a simple test mesh
let mesh = Mesh {
    eints: 1,   // 1 energy group (total only)
    tints: 2,   // 3 time groups (2+total)
    iints: 2,
    jints: 3,
    kints: 4,
    ..Default::default()
};
// Calculate the coordinates for the voxel
assert_eq!(mesh.n_voxels_expected(), 72);

pub fn n_voxels_per_group(&self) -> usize

Expected number of voxels per group set

Returns the product of the number of fine mesh bins in each dimension. This is equivalent to the total number of voxels in the mesh geometry itself. For example, a 5x5x5 rectangular mesh should have 100 voxels for every unique energy/time group set.

// Generate a simple test mesh
let mesh = Mesh {
    iints: 2,
    jints: 3,
    kints: 4,
    ..Default::default()
};
// Calculate the coordinates for the voxel
assert_eq!(mesh.n_voxels_per_group(), 24);

pub fn voxel_coordinates(&self, index: usize) -> Result<VoxelCoordinate, Error>

Get a full set of coordinates for a voxel

// Get some voxel from an example file
let mesh = read_target("./data/meshes/fmesh_114.msht", 114).unwrap();
// Calculate the coordinates for the voxel
assert_eq!( mesh.voxel_coordinates(63).unwrap(),
            VoxelCoordinate {
                energy: Group::Value(1.0),
                time: Group::Value(1.0E+30),
                i: 9.375,
                j: 4.500,
                k: 2.500} );

pub fn voxel_data(&self) -> Vec<(f64, f64)>

Collect (value, error) pairs for all voxels in the mesh

For example:

let mesh = Mesh {
    voxels: vec![
        Voxel{index: 0, result: 1.0, error: 0.1},
        Voxel{index: 1, result: 2.0, error: 0.2},
        Voxel{index: 2, result: 3.0, error: 0.3},
    ],
    ..Default::default()
};

assert_eq!(mesh.voxel_data(), vec![(1.0, 0.1), (2.0, 0.2), (3.0, 0.3)]);

pub fn voxels_by_group_index( &self, e_idx: usize, t_idx: usize, ) -> Result<&[Voxel], Error>

Slice the full list of mesh Voxels by energy/time index

Very fast, but operates on indicies and therefore relies on the voxels being sorted. For an arbitrary list of Voxels, use the filter methods which explicitly check the groups of every voxel provided.

pub fn voxels_by_group_value( &self, energy_group: Group, time_group: Group, ) -> Result<&[Voxel], Error>

Slice the full list of mesh Voxels by both energy/time groups

Very fast, but operates on indicies and therefore relies on the voxels being sorted. For an arbitrary list of Voxels, use the filter methods which explicitly check the groups of every voxel provided.

§

impl Mesh

pub fn energy_bins_upper(&self) -> &[f64]

Returns slice of emesh for upper energy bin edges

let mesh = Mesh {
    emesh: vec![0.0, 1.0, 2.0, 3.0, 4.0],
    ..Default::default()
};
// Get a slice of the upper edges of energy bins
assert_eq!(mesh.energy_bins_upper(), vec![1.0, 2.0, 3.0, 4.0]);

pub fn energy_bins_lower(&self) -> &[f64]

Returns slice of emesh for lower energy bin edges

let mesh = Mesh {
    emesh: vec![0.0, 1.0, 2.0, 3.0, 4.0],
    ..Default::default()
};
// Get a slice of the lower edges of energy bins
assert_eq!(mesh.energy_bins_lower(), vec![0.0, 1.0, 2.0, 3.0]);

pub fn energy_groups(&self) -> Vec<Group>

Returns a collection of all energy groups, including total

Builds a full list of the energy groups in the mesh, which can include both Group variants.

let mesh = Mesh {
    eints: 2,
    emesh: vec![0.0, 1.0, 2.0],
    ..Default::default()
};
// See what energy groups the voxels may be split into
assert_eq!(mesh.energy_groups(), vec![Group::Value(1.0), // bin 0.0 - 1.0
                                      Group::Value(2.0), // bin 1.0 - 2.0
                                      Group::Total]);

pub fn time_bins_upper(&self) -> &[f64]

Returns slice of tmesh for upper time bin edges

let mesh = Mesh {
    tmesh: vec![0.0, 1e12, 1e16, 1e20],
    ..Default::default()
};
// Get a slice of the upper edges of time bins
assert_eq!(mesh.time_bins_upper(), vec![1e12, 1e16, 1e20]);

pub fn time_bins_lower(&self) -> &[f64]

Returns slice of emesh for lower time bin edges

let mesh = Mesh {
    tmesh: vec![0.0, 1e12, 1e16, 1e20],
    ..Default::default()
};
// Get a slice of the lower edges of time bins
assert_eq!(mesh.time_bins_lower(), vec![0.0, 1e12, 1e16]);

pub fn time_groups(&self) -> Vec<Group>

Returns a collection of all time groups, including total

Builds a full list of the time groups in the mesh, which can include both Group variants.

let mesh = Mesh {
    tints: 2,
    tmesh: vec![0.0, 1e12, 1e16],
    ..Default::default()
};
// See what time groups the voxels may be split into
assert_eq!(mesh.time_groups(), vec![Group::Value(1e12), // bin 0.0  - 1e12
                                    Group::Value(1e16), // bin 1e12 - 1e16
                                    Group::Total]);     // 'total' group

§

impl Mesh

Indexing and conversion helpers

pub fn voxel_index_from_etijk( &self, e_idx: usize, t_idx: usize, i_idx: usize, j_idx: usize, k_idx: usize, ) -> usize

Find the global voxel index from (e,t,i,j,k) indicies

The voxel index corresponds to the order seen in the column format output files. All formats are coerced into this sorting for consistency.

MCNP writes voxels as (k,j,i,t,e) loops, i.e:

for e in energy_groups
    for t in time_groups
        for i in imesh
            for j in jmesh
                for k in kmesh
                    // ... write voxel

Frustratingly, MCNP internally uses the “cell index” ordering which loops in the reverse for (i,j,k). See cell_index_from_etijk() for details.

pub fn cell_index_from_etijk( &self, e_idx: usize, t_idx: usize, i_idx: usize, j_idx: usize, k_idx: usize, ) -> usize

Find the global cell index from (e,t,i,j,k) indicies

The cell index corresponds to the internal ordering used by MCNP, not what is seen in output files. The cell index number loops energy and time groups as expected, but (i,j,k) indicies are calculated as:

    # FORTAN indexing
    cell index number = i + (j − 1) x I + (k − 1) x I x J

where I and J are the total number of imesh and jmesh bins.

In other words, this loops as:

for e in energy_groups
    for t in time_groups
        for k in kmesh
            for j in jmesh
                for i in imesh
                    // ...some voxel

This is the order needed by VTK formats and weight window files.

pub fn etijk_from_voxel_index( &self, idx: usize, ) -> (usize, usize, usize, usize, usize)

Find the (e,t,i,j,k) indicies for a given voxel index

The reverse of voxel_index_to_etijk().

pub fn cell_index_from_voxel_index(&self, idx: usize) -> usize

Convert voxel index to a cell index

Besides the obvious convenience, this is needed for UKAEA CuV formats. For some reason it was thought hilarious to order the Number_of_material_cells_per_voxel array by cell index while ordering the output column data by voxel index.

pub fn voxel_index_from_cell_index(&self, idx: usize) -> usize

Convert a cell index to a voxel index

pub fn etijk_from_cell_index( &self, idx: usize, ) -> (usize, usize, usize, usize, usize)

Find the (e,t,i,j,k) indicies for a given cell index

pub fn energy_group_from_value(&self, energy: f64) -> Result<Group, Error>

For a given energy, find what group the results are under

Error returned for handling if the energy is outside of the emesh bounds and should be handled by the caller.

Important: the group an energy is under is determined by find_bin_inclusive() rules for MCNP output, even though internally it follows the find_bin_exclusive() rules.

pub fn energy_group_from_index(&self, e_idx: usize) -> Result<Group, Error>

Get the energy Group for the energy index e_idx

For example:

// Generate a simple test mesh
let mesh = Mesh {
    eints: 2,
    emesh: vec![0.0, 1.0, 20.0],
    ..Default::default()
};
// Find the group for an energy bin index
assert_eq!(mesh.energy_group_from_index(0).unwrap(), Group::Value(1.0));
assert_eq!(mesh.energy_group_from_index(1).unwrap(), Group::Value(20.0));
assert_eq!(mesh.energy_group_from_index(2).unwrap(), Group::Total);
assert!(mesh.energy_group_from_index(3).is_err());

pub fn energy_index_from_group( &self, energy_group: Group, ) -> Result<usize, Error>

For a given energy group, find the corresponding emesh bin index

Important: the group an energy is under is determined by find_bin_inclusive() rules for MCNP output, even though internally it follows the find_bin_exclusive() rules.

pub fn energy_index_from_value(&self, energy: f64) -> Result<usize, Error>

For a given energy, find the corresponding emesh bin index

pub fn time_group_from_value(&self, time: f64) -> Result<Group, Error>

For a given time, find what group the results are under

Error returned for handling if the time is outside of the emesh bounds and should be handled by the caller.

Important: the group a time is under is determined by find_bin_inclusive() rules for MCNP output, even though internally it follows the find_bin_exclusive() rules.

pub fn time_group_from_index(&self, t_idx: usize) -> Result<Group, Error>

Get the time Group for the time index t_idx

For example:

// Generate a simple test mesh
let mesh = Mesh {
    tints: 2,
    tmesh: vec![0.0, 1e16, 1e30],
    ..Default::default()
};
// Find the group for an energy bin index
assert_eq!(mesh.time_group_from_index(0).unwrap(), Group::Value(1e16));
assert_eq!(mesh.time_group_from_index(1).unwrap(), Group::Value(1e30));
assert_eq!(mesh.time_group_from_index(2).unwrap(), Group::Total);
assert!(mesh.time_group_from_index(3).is_err());

pub fn time_index_from_group(&self, time_group: Group) -> Result<usize, Error>

For a given time group, find the corresponding tmesh bin index

Important: the group a time is under is determined by find_bin_inclusive() rules for MCNP output, even though internally it follows the find_bin_exclusive() rules.