Merge pull request #1066 from paulromano/cpp-mesh

Convert mesh type, Shannon entropy, UFS to C++

CMakeLists.txt

@@ -312,7 +312,6 @@ add_library(libopenmc SHARED
   src/material_header.F90
   src/math.F90
   src/matrix_header.F90
-  src/mesh.F90
   src/mesh_header.F90
   src/message_passing.F90
   src/mgxs_data.F90
@@ -380,11 +379,13 @@ add_library(libopenmc SHARED
   src/tallies/trigger.F90
   src/tallies/trigger_header.F90
   src/cell.cpp
+  src/cmfd_execute.cpp
   src/distribution.cpp
   src/distribution_angle.cpp
   src/distribution_energy.cpp
   src/distribution_multi.cpp
   src/distribution_spatial.cpp
+  src/eigenvalue.cpp
   src/endf.cpp
   src/initialize.cpp
   src/finalize.cpp
@@ -394,6 +395,7 @@ add_library(libopenmc SHARED
   src/lattice.cpp
   src/material.cpp
   src/math_functions.cpp
+  src/mesh.cpp
   src/message_passing.cpp
   src/mgxs.cpp
   src/mgxs_interface.cpp

include/openmc/capi.h

@@ -37,6 +37,7 @@ extern "C" {
   int openmc_filter_set_type(int32_t index, const char* type);
   int openmc_finalize();
   int openmc_find_cell(double* xyz, int32_t* index, int32_t* instance);
+  int openmc_fission_bank(struct Bank** ptr, int64_t* n);
   int openmc_get_cell_index(int32_t id, int32_t* index);
   int openmc_get_filter_index(int32_t id, int32_t* index);
   void openmc_get_filter_next_id(int32_t* id);
@@ -71,7 +72,7 @@ extern "C" {
   int openmc_mesh_get_params(int32_t index, double** ll, double** ur, double** width, int* n);
   int openmc_mesh_set_id(int32_t index, int32_t id);
   int openmc_mesh_set_dimension(int32_t index, int n, const int* dims);
-  int openmc_mesh_set_params(int32_t index, const double* ll, const double* ur, const double* width, int n);
+  int openmc_mesh_set_params(int32_t index, int n, const double* ll, const double* ur, const double* width);
   int openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh);
   int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh);
   int openmc_next_batch(int* status);
@@ -139,7 +140,6 @@ extern "C" {
   extern int32_t n_filters;
   extern int32_t n_lattices;
   extern int32_t n_materials;
-  extern int32_t n_meshes;
   extern int n_nuclides;
   extern int32_t n_plots;
   extern int32_t n_realizations;

include/openmc/eigenvalue.h

@@ -0,0 +1,41 @@
+#ifndef OPENMC_EIGENVALUE_H
+#define OPENMC_EIGENVALUE_H
+
+#include <cstdint> // for int64_t
+#include <vector>
+
+#include "xtensor/xtensor.hpp"
+
+#include "openmc/particle.h"
+
+namespace openmc {
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+extern std::vector<double> entropy; //!< Shannon entropy at each generation
+extern xt::xtensor<double, 1> source_frac; //!< Source fraction for UFS
+
+extern "C" int64_t n_bank;
+#pragma omp threadprivate(n_bank)
+
+//==============================================================================
+// Non-member functions
+//==============================================================================
+
+//! Calculates the Shannon entropy of the fission source distribution to assess
+//! source convergence
+extern "C" void shannon_entropy();
+
+//! Determines the source fraction in each UFS mesh cell and reweights the
+//! source bank so that the sum of the weights is equal to n_particles. The
+//! 'source_frac' variable is used later to bias the production of fission sites
+extern "C" void ufs_count_sites();
+
+//! Get UFS weight corresponding to particle's location
+extern "C" double ufs_get_weight(const Particle* p);
+
+} // namespace openmc
+
+#endif // OPENMC_EIGENVALUE_H

include/openmc/hdf5_interface.h

@@ -317,6 +317,14 @@ write_attribute(hid_t obj_id, const char* name, const std::array<T, N>& buffer)
   write_attr(obj_id, 1, dims, name, H5TypeMap<T>::type_id, buffer.data());
 }
 
+template<typename T> inline void
+write_attribute(hid_t obj_id, const char* name, const std::vector<T>& buffer)
+{
+  hsize_t dims[] {buffer.size()};
+  write_attr(obj_id, 1, dims, name, H5TypeMap<T>::type_id, buffer.data());
+}
+
+
 //==============================================================================
 // Templates/overloads for write_dataset
 //==============================================================================
@@ -347,6 +355,15 @@ write_dataset(hid_t obj_id, const char* name, const std::vector<T>& buffer)
   write_dataset(obj_id, 1, dims, name, H5TypeMap<T>::type_id, buffer.data(), false);
 }
 
+template<typename T> inline void
+write_dataset(hid_t obj_id, const char* name, const xt::xarray<T>& arr)
+{
+  auto s = arr.shape();
+  std::vector<hsize_t> dims {s.cbegin(), s.cend()};
+  write_dataset(obj_id, dims.size(), dims.data(), name, H5TypeMap<T>::type_id,
+                arr.data(), false);
+}
+
 inline void
 write_dataset(hid_t obj_id, const char* name, Position r)
 {

include/openmc/mesh.h

@@ -0,0 +1,131 @@
+//! \file mesh.h
+//! \brief Mesh types used for tallies, Shannon entropy, CMFD, etc.
+
+#ifndef OPENMC_MESH_H
+#define OPENMC_MESH_H
+
+#include <memory> // for unique_ptr
+#include <vector>
+#include <unordered_map>
+
+#include "hdf5.h"
+#include "pugixml.hpp"
+#include "xtensor/xarray.hpp"
+
+#include "openmc/particle.h"
+#include "openmc/position.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Tessellation of n-dimensional Euclidean space by congruent squares or cubes
+//==============================================================================
+
+class RegularMesh {
+public:
+  // Constructors
+  RegularMesh() = default;
+  RegularMesh(pugi::xml_node node);
+
+  // Methods
+
+  //! Determine which bins were crossed by a particle
+  //!
+  //! \param[in] p Particle to check
+  //! \param[out] bins Bins that were crossed
+  //! \param[out] lengths Fraction of tracklength in each bin
+  void bins_crossed(const Particle* p, std::vector<int>& bins,
+                    std::vector<double>& lengths) const;
+
+  //! Determine which surface bins were crossed by a particle
+  //!
+  //! \param[in] p Particle to check
+  //! \param[out] bins Surface bins that were crossed
+  void surface_bins_crossed(const Particle* p, std::vector<int>& bins) const;
+
+  //! Get bin at a given position in space
+  //!
+  //! \param[in] r Position to get bin for
+  //! \return Mesh bin
+  int get_bin(Position r) const;
+
+  //! Get bin given mesh indices
+  //!
+  //! \param[in] Array of mesh indices
+  //! \return Mesh bin
+  int get_bin_from_indices(const int* ijk) const;
+
+  //! Get mesh indices given a position
+  //!
+  //! \param[in] r Position to get indices for
+  //! \param[out] ijk Array of mesh indices
+  //! \param[out] in_mesh Whether position is in mesh
+  void get_indices(Position r, int* ijk, bool* in_mesh) const;
+
+  //! Get mesh indices corresponding to a mesh bin
+  //!
+  //! \param[in] bin Mesh bin
+  //! \param[out] ijk Mesh indices
+  void get_indices_from_bin(int bin, int* ijk) const;
+
+  //! Check if a line connected by two points intersects the mesh
+  //!
+  //! \param[in] r0 Starting position
+  //! \param[in] r1 Ending position
+  //! \return Whether line connecting r0 and r1 intersects mesh
+  bool intersects(Position r0, Position r1) const;
+
+  //! Write mesh data to an HDF5 group
+  //!
+  //! \param[in] group HDF5 group
+  void to_hdf5(hid_t group) const;
+
+  //! Count number of bank sites in each mesh bin / energy bin
+  //!
+  //! \param[in] n Number of bank sites
+  //! \param[in] bank Array of bank sites
+  //! \param[in] n_energy Number of energies
+  //! \param[in] energies Array of energies
+  //! \param[out] Whether any bank sites are outside the mesh
+  //! \return Array indicating number of sites in each mesh/energy bin
+  xt::xarray<double> count_sites(int64_t n, const Bank* bank,
+    int n_energy, const double* energies, bool* outside) const;
+
+  int id_ {-1};  //!< User-specified ID
+  int n_dimension_; //!< Number of dimensions
+  double volume_frac_; //!< Volume fraction of each mesh element
+  xt::xarray<int> shape_; //!< Number of mesh elements in each dimension
+  xt::xarray<double> lower_left_; //!< Lower-left coordinates of mesh
+  xt::xarray<double> upper_right_; //!< Upper-right coordinates of mesh
+  xt::xarray<double> width_; //!< Width of each mesh element
+
+private:
+  bool intersects_1d(Position r0, Position r1) const;
+  bool intersects_2d(Position r0, Position r1) const;
+  bool intersects_3d(Position r0, Position r1) const;
+};
+
+//==============================================================================
+// Non-member functions
+//==============================================================================
+
+//! Read meshes from either settings/tallies
+//! \param[in] root XML node
+extern "C" void read_meshes(pugi::xml_node* root);
+
+//! Write mesh data to an HDF5 group
+//! \param[in] group HDF5 group
+extern "C" void meshes_to_hdf5(hid_t group);
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+extern std::vector<std::unique_ptr<RegularMesh>> meshes;
+
+extern std::unordered_map<int32_t, int32_t> mesh_map;
+
+
+} // namespace openmc
+
+#endif // OPENMC_MESH_H

include/openmc/message_passing.h

@@ -10,6 +10,7 @@ namespace mpi {
 
   extern int rank;
   extern int n_procs;
+  extern bool master;
 
 #ifdef OPENMC_MPI
   extern MPI_Datatype bank;

include/openmc/output.h

@@ -22,5 +22,7 @@ void header(const char* msg, int level);
 
 extern "C" void print_overlap_check();
 
+extern "C" void title();
+
 } // namespace openmc
 #endif // OPENMC_OUTPUT_H

include/openmc/particle.h

@@ -161,7 +161,7 @@ extern "C" {
     {mark_as_lost(message.str());}
 
     //! create a particle restart HDF5 file
-    void write_restart();
+    void write_restart() const;
   };
 
 

include/openmc/position.h

@@ -1,6 +1,7 @@
 #ifndef OPENMC_POSITION_H
 #define OPENMC_POSITION_H
 
+#include <cmath>
 #include <vector>
 
 namespace openmc {
@@ -46,6 +47,9 @@ struct Position {
   inline double dot(Position other) {
     return x*other.x + y*other.y + z*other.z;
   }
+  inline double norm() {
+    return std::sqrt(x*x + y*y + z*z);
+  }
 
   // Data members
   double x = 0.;

include/openmc/settings.h

@@ -89,13 +89,12 @@ extern "C" double weight_survive;        //!< Survival weight after Russian roul
 
 } // namespace settings
 
-//==============================================================================
 //! Read settings from XML file
 //! \param[in] root XML node for <settings>
-//==============================================================================
-
 extern "C" void read_settings_xml();
 
+extern "C" void read_settings_xml_f(pugi::xml_node_struct* root_ptr);
+
 } // namespace openmc
 
 #endif // OPENMC_SETTINGS_H

include/openmc/xml_interface.h

@@ -1,11 +1,14 @@
 #ifndef OPENMC_XML_INTERFACE_H
 #define OPENMC_XML_INTERFACE_H
 
+#include <cstddef> // for size_t
 #include <sstream> // for stringstream
 #include <string>
 #include <vector>
 
 #include "pugixml.hpp"
+#include "xtensor/xarray.hpp"
+#include "xtensor/xadapt.hpp"
 
 
 namespace openmc {
@@ -38,5 +41,14 @@ std::vector<T> get_node_array(pugi::xml_node node, const char* name,
   return values;
 }
 
+template <typename T>
+xt::xarray<T> get_node_xarray(pugi::xml_node node, const char* name,
+                              bool lowercase=false)
+{
+  std::vector<T> v = get_node_array<T>(node, name, lowercase);
+  std::vector<std::size_t> shape = {v.size()};
+  return xt::adapt(v, shape);
+}
+
 } // namespace openmc
 #endif // OPENMC_XML_INTERFACE_H

openmc/capi/mesh.py

@@ -41,6 +41,8 @@ _dll.openmc_mesh_set_params.errcheck = _error_handler
 _dll.openmc_get_mesh_index.argtypes = [c_int32, POINTER(c_int32)]
 _dll.openmc_get_mesh_index.restype = c_int
 _dll.openmc_get_mesh_index.errcheck = _error_handler
+_dll.n_meshes.argtypes = []
+_dll.n_meshes.restype = c_int
 
 
 class Mesh(_FortranObjectWithID):
@@ -172,10 +174,10 @@ class _MeshMapping(Mapping):
 
     def __iter__(self):
         for i in range(len(self)):
-            yield Mesh(index=i + 1).id
+            yield Mesh(index=i).id
 
     def __len__(self):
-        return c_int32.in_dll(_dll, 'n_meshes').value
+        return _dll.n_meshes()
 
     def __repr__(self):
         return repr(dict(self))

src/api.F90

@@ -11,7 +11,6 @@ module openmc_api
   use hdf5_interface
   use material_header
   use math
-  use mesh_header
   use message_passing
   use nuclide_header
   use initialize,      only: openmc_init_f
@@ -305,6 +304,9 @@ contains
     interface
       subroutine free_memory_source() bind(C)
       end subroutine
+
+      subroutine free_memory_mesh() bind(C)
+      end subroutine free_memory_mesh
     end interface
 
     call free_memory_geometry()

src/bank_header.F90

@@ -28,7 +28,7 @@ module bank_header
   type(Bank), allocatable, target :: master_fission_bank(:)
 #endif
 
-  integer(8) :: n_bank       ! # of sites in fission bank
+  integer(C_INT64_T), bind(C) :: n_bank       ! # of sites in fission bank
 
 !$omp threadprivate(fission_bank, n_bank)
 
@@ -71,4 +71,20 @@ contains
     end if
   end function openmc_source_bank
 
+  function openmc_fission_bank(ptr, n) result(err) bind(C)
+    ! Return a pointer to the source bank
+    type(C_PTR), intent(out) :: ptr
+    integer(C_INT64_T), intent(out) :: n
+    integer(C_INT) :: err
+
+    if (.not. allocated(fission_bank)) then
+      err = E_ALLOCATE
+      call set_errmsg("Fission bank has not been allocated.")
+    else
+      err = 0
+      ptr = C_LOC(fission_bank)
+      n = size(fission_bank)
+    end if
+  end function openmc_fission_bank
+
 end module bank_header

src/cmfd_data.F90

@@ -80,7 +80,7 @@ contains
     integer :: i_mesh        ! flattend index for mesh
     logical :: energy_filters! energy filters present
     real(8) :: flux          ! temp variable for flux
-    type(RegularMesh), pointer :: m ! pointer for mesh object
+    type(RegularMesh) :: m ! pointer for mesh object
 
     ! Extract spatial and energy indices from object
     nx = cmfd % indices(1)
@@ -99,7 +99,7 @@ contains
 
     select type(filt => filters(i_filter_mesh) % obj)
     type is (MeshFilter)
-      m => meshes(filt % mesh)
+      m = meshes(filt % mesh)
     end select
 
     ! Set mesh widths
@@ -354,9 +354,6 @@ contains
     ! Normalize openmc source distribution
     cmfd % openmc_src = cmfd % openmc_src/sum(cmfd % openmc_src)*cmfd%norm
 
-    ! Nullify all pointers
-    if (associated(m)) nullify(m)
-
   end subroutine compute_xs
 
 !===============================================================================

src/cmfd_execute.F90

@@ -214,8 +214,6 @@ contains
     use bank_header, only: source_bank
     use constants,   only: ZERO, ONE
     use error,       only: warning, fatal_error
-    use mesh_header, only: RegularMesh
-    use mesh,        only: count_bank_sites
     use message_passing
     use string,      only: to_str
 
@@ -224,7 +222,7 @@ contains
     integer :: nx       ! maximum number of cells in x direction
     integer :: ny       ! maximum number of cells in y direction
     integer :: nz       ! maximum number of cells in z direction
-    integer :: ng       ! maximum number of energy groups
+    integer(C_INT) :: ng       ! maximum number of energy groups
     integer :: i        ! iteration counter
     integer :: g        ! index for group
     integer :: ijk(3)   ! spatial bin location
@@ -232,12 +230,22 @@ contains
     integer :: mesh_bin ! mesh bin of soruce particle
     integer :: n_groups ! number of energy groups
     real(8) :: norm     ! normalization factor
-    logical :: outside  ! any source sites outside mesh
+    logical(C_BOOL) :: outside  ! any source sites outside mesh
     logical :: in_mesh  ! source site is inside mesh
 #ifdef OPENMC_MPI
     integer :: mpi_err
 #endif
 
+    interface
+      subroutine cmfd_populate_sourcecounts(ng, energies, source_counts, outside) bind(C)
+        import C_INT, C_DOUBLE, C_BOOL
+        integer(C_INT), value :: ng
+        real(C_DOUBLE), intent(in) :: energies
+        real(C_DOUBLE), intent(out) :: source_counts
+        logical(C_BOOL), intent(out) :: outside
+      end subroutine
+    end interface
+
     ! Get maximum of spatial and group indices
     nx = cmfd % indices(1)
     ny = cmfd % indices(2)
@@ -261,8 +269,8 @@ contains
       cmfd%weightfactors = ONE
 
       ! Count bank sites in mesh and reverse due to egrid structure
-      call count_bank_sites(cmfd_mesh, source_bank, cmfd%sourcecounts, &
-           cmfd % egrid, sites_outside=outside, size_bank=work)
+      call cmfd_populate_sourcecounts(ng + 1, cmfd % egrid(1), &
+           cmfd % sourcecounts(1,1), outside)
 
       ! Check for sites outside of the mesh
       if (master .and. outside) then

src/cmfd_execute.cpp

@@ -0,0 +1,34 @@
+#include <algorithm> // for copy
+#include <cstdint>
+#include <iostream>
+
+#include "xtensor/xarray.hpp"
+#include "xtensor/xio.hpp"
+
+#include "openmc/capi.h"
+#include "openmc/mesh.h"
+
+namespace openmc {
+
+extern "C" int index_cmfd_mesh;
+
+extern "C" void
+cmfd_populate_sourcecounts(int n_energy, const double* energies,
+  double* source_counts, bool* outside)
+{
+  // Get pointer to source bank
+  Bank* source_bank;
+  int64_t n;
+  openmc_source_bank(&source_bank, &n);
+
+  // Get source counts in each mesh bin / energy bin
+  auto& m = meshes.at(index_cmfd_mesh);
+  xt::xarray<double> counts = m->count_sites(openmc_work, source_bank, n_energy, energies, outside);
+
+  std::cout << counts << "\n";
+
+  // Copy data from the xarray into the source counts array
+  std::copy(counts.begin(), counts.end(), source_counts);
+}
+
+} // namespace openmc

src/cmfd_header.F90

@@ -1,5 +1,7 @@
 module cmfd_header
 
+  use, intrinsic :: ISO_C_BINDING
+
   use constants,  only: CMFD_NOACCEL, ZERO, ONE
   use mesh_header, only: RegularMesh
   use set_header, only: SetInt
@@ -24,7 +26,7 @@ module cmfd_header
     integer :: mat_dim = CMFD_NOACCEL
 
     ! Energy grid
-    real(8), allocatable :: egrid(:)
+    real(C_DOUBLE), allocatable :: egrid(:)
 
     ! Cross sections
     real(8), allocatable :: totalxs(:,:,:,:)
@@ -53,7 +55,7 @@ module cmfd_header
     real(8), allocatable :: openmc_src(:,:,:,:)
 
     ! Source sites in each mesh box
-    real(8), allocatable :: sourcecounts(:,:)
+    real(C_DOUBLE), allocatable :: sourcecounts(:,:)
 
     ! Weight adjustment factors
     real(8), allocatable :: weightfactors(:,:,:,:)
@@ -95,7 +97,8 @@ module cmfd_header
   ! Main object
   type(cmfd_type), public :: cmfd
 
-  type(RegularMesh), public, pointer :: cmfd_mesh => null()
+  integer(C_INT), public, bind(C) :: index_cmfd_mesh
+  type(RegularMesh), public :: cmfd_mesh
 
   ! Pointers for different tallies
   type(TallyContainer), public, pointer :: cmfd_tallies(:) => null()

src/cmfd_input.F90

@@ -3,7 +3,7 @@ module cmfd_input
   use, intrinsic :: ISO_C_BINDING
 
   use cmfd_header
-  use mesh_header,    only: mesh_dict
+  use mesh_header
   use mgxs_interface, only: energy_bins, num_energy_groups
   use tally
   use tally_header
@@ -241,7 +241,7 @@ contains
 
     use constants,        only: MAX_LINE_LEN
     use error,            only: fatal_error, warning
-    use mesh_header,      only: RegularMesh, openmc_extend_meshes
+    use mesh_header
     use string
     use tally,            only: openmc_tally_allocate
     use tally_header,     only: openmc_extend_tallies
@@ -263,115 +263,23 @@ contains
     integer :: i_filt     ! index in filters array
     integer :: filt_id
     integer :: tally_id
-    integer :: iarray3(3) ! temp integer array
-    real(8) :: rarray3(3) ! temp double array
     real(C_DOUBLE), allocatable :: energies(:)
-    type(RegularMesh), pointer :: m
     type(XMLNode) :: node_mesh
 
-    err = openmc_extend_meshes(1, i_start)
+    ! Read CMFD mesh
+    call read_meshes(root % ptr)
 
-    ! Allocate mesh
-    cmfd_mesh => meshes(i_start)
-    m => meshes(i_start)
+    ! Get index of cmfd mesh and set ID
+    i_start = n_meshes() - 1
+    err = openmc_mesh_set_id(i_start, i_start)
 
-    ! Set mesh id
-    m % id = i_start
-
-    ! Set mesh type to rectangular
-    m % type = MESH_REGULAR
+    ! Save reference to CMFD mesh
+    index_cmfd_mesh = i_start
+    cmfd_mesh = meshes(i_start)
 
     ! Get pointer to mesh XML node
     node_mesh = root % child("mesh")
 
-    ! Determine number of dimensions for mesh
-    n = node_word_count(node_mesh, "dimension")
-    if (n /= 2 .and. n /= 3) then
-      call fatal_error("Mesh must be two or three dimensions.")
-    end if
-    m % n_dimension = n
-
-    ! Allocate attribute arrays
-    allocate(m % dimension(n))
-    allocate(m % lower_left(n))
-    allocate(m % width(n))
-    allocate(m % upper_right(n))
-
-    ! Check that dimensions are all greater than zero
-    call get_node_array(node_mesh, "dimension", iarray3(1:n))
-    if (any(iarray3(1:n) <= 0)) then
-      call fatal_error("All entries on the <dimension> element for a tally mesh&
-           & must be positive.")
-    end if
-
-    ! Read dimensions in each direction
-    m % dimension = iarray3(1:n)
-
-    ! Read mesh lower-left corner location
-    if (m % n_dimension /= node_word_count(node_mesh, "lower_left")) then
-      call fatal_error("Number of entries on <lower_left> must be the same as &
-           &the number of entries on <dimension>.")
-    end if
-    call get_node_array(node_mesh, "lower_left", m % lower_left)
-
-    ! Make sure both upper-right or width were specified
-    if (check_for_node(node_mesh, "upper_right") .and. &
-         check_for_node(node_mesh, "width")) then
-      call fatal_error("Cannot specify both <upper_right> and <width> on a &
-           &tally mesh.")
-    end if
-
-    ! Make sure either upper-right or width was specified
-    if (.not.check_for_node(node_mesh, "upper_right") .and. &
-         .not.check_for_node(node_mesh, "width")) then
-      call fatal_error("Must specify either <upper_right> and <width> on a &
-           &tally mesh.")
-    end if
-
-    if (check_for_node(node_mesh, "width")) then
-      ! Check to ensure width has same dimensions
-      if (node_word_count(node_mesh, "width") /= &
-           node_word_count(node_mesh, "lower_left")) then
-        call fatal_error("Number of entries on <width> must be the same as the &
-             &number of entries on <lower_left>.")
-      end if
-
-      ! Check for negative widths
-      call get_node_array(node_mesh, "width", rarray3(1:n))
-      if (any(rarray3(1:n) < ZERO)) then
-        call fatal_error("Cannot have a negative <width> on a tally mesh.")
-      end if
-
-      ! Set width and upper right coordinate
-      m % width = rarray3(1:n)
-      m % upper_right = m % lower_left + m % dimension * m % width
-
-    elseif (check_for_node(node_mesh, "upper_right")) then
-      ! Check to ensure width has same dimensions
-      if (node_word_count(node_mesh, "upper_right") /= &
-           node_word_count(node_mesh, "lower_left")) then
-        call fatal_error("Number of entries on <upper_right> must be the same &
-             &as the number of entries on <lower_left>.")
-      end if
-
-      ! Check that upper-right is above lower-left
-      call get_node_array(node_mesh, "upper_right", rarray3(1:n))
-      if (any(rarray3(1:n) < m % lower_left)) then
-        call fatal_error("The <upper_right> coordinates must be greater than &
-             &the <lower_left> coordinates on a tally mesh.")
-      end if
-
-      ! Set upper right coordinate and width
-      m % upper_right = rarray3(1:n)
-      m % width = (m % upper_right - m % lower_left) / real(m % dimension, 8)
-    end if
-
-    ! Set volume fraction
-    m % volume_frac = ONE/real(product(m % dimension),8)
-
-    ! Add mesh to dictionary
-    call mesh_dict % set(m % id, i_start)
-
     ! Determine number of filters
     energy_filters = check_for_node(node_mesh, "energy")
     n = merge(5, 3, energy_filters)

src/eigenvalue.F90

@@ -6,8 +6,6 @@ module eigenvalue
   use constants,   only: ZERO
   use error,       only: fatal_error, warning
   use math,        only: t_percentile
-  use mesh,        only: count_bank_sites
-  use mesh_header, only: RegularMesh, meshes
   use message_passing
   use random_lcg,  only: prn, set_particle_seed, advance_prn_seed
   use settings
@@ -294,46 +292,6 @@ contains
 
   end subroutine synchronize_bank
 
-!===============================================================================
-! SHANNON_ENTROPY calculates the Shannon entropy of the fission source
-! distribution to assess source convergence
-!===============================================================================
-
-  subroutine shannon_entropy()
-
-    integer :: i              ! index for mesh elements
-    real(8) :: entropy_gen    ! entropy at this generation
-    logical :: sites_outside  ! were there sites outside entropy box?
-
-    associate (m => meshes(index_entropy_mesh))
-      ! count number of fission sites over mesh
-      call count_bank_sites(m, fission_bank, entropy_p, &
-           size_bank=n_bank, sites_outside=sites_outside)
-
-      ! display warning message if there were sites outside entropy box
-      if (sites_outside) then
-        if (master) call warning("Fission source site(s) outside of entropy box.")
-      end if
-
-      ! sum values to obtain shannon entropy
-      if (master) then
-        ! Normalize to total weight of bank sites
-        entropy_p = entropy_p / sum(entropy_p)
-
-        entropy_gen = ZERO
-        do i = 1, size(entropy_p, 2)
-          if (entropy_p(1,i) > ZERO) then
-            entropy_gen = entropy_gen - &
-                 entropy_p(1,i) * log(entropy_p(1,i))/log(TWO)
-          end if
-        end do
-
-        ! Add value to vector
-        call entropy % push_back(entropy_gen)
-      end if
-    end associate
-  end subroutine shannon_entropy
-
 !===============================================================================
 ! CALCULATE_GENERATION_KEFF collects the single-processor tracklength k's onto
 ! the master processor and normalizes them. This should work whether or not the
@@ -577,61 +535,6 @@ contains
 
   end function openmc_get_keff
 
-!===============================================================================
-! COUNT_SOURCE_FOR_UFS determines the source fraction in each UFS mesh cell and
-! reweights the source bank so that the sum of the weights is equal to
-! n_particles. The 'source_frac' variable is used later to bias the production
-! of fission sites
-!===============================================================================
-
-  subroutine count_source_for_ufs()
-
-    real(8) :: total         ! total weight in source bank
-    logical :: sites_outside ! were there sites outside the ufs mesh?
-#ifdef OPENMC_MPI
-    integer :: n             ! total number of ufs mesh cells
-    integer :: mpi_err       ! MPI error code
-#endif
-
-    associate (m => meshes(index_ufs_mesh))
-
-    if (current_batch == 1 .and. current_gen == 1) then
-      ! On the first generation, just assume that the source is already evenly
-      ! distributed so that effectively the production of fission sites is not
-      ! biased
-
-      source_frac = m % volume_frac
-
-    else
-      ! count number of source sites in each ufs mesh cell
-      call count_bank_sites(m, source_bank, source_frac, &
-           sites_outside=sites_outside, size_bank=work)
-
-      ! Check for sites outside of the mesh
-      if (master .and. sites_outside) then
-        call fatal_error("Source sites outside of the UFS mesh!")
-      end if
-
-#ifdef OPENMC_MPI
-      ! Send source fraction to all processors
-      n = product(m % dimension)
-      call MPI_BCAST(source_frac, n, MPI_REAL8, 0, mpi_intracomm, mpi_err)
-#endif
-
-      ! Normalize to total weight to get fraction of source in each cell
-      total = sum(source_frac)
-      source_frac = source_frac / total
-
-      ! Since the total starting weight is not equal to n_particles, we need to
-      ! renormalize the weight of the source sites
-
-      source_bank % wgt = source_bank % wgt * n_particles / total
-    end if
-
-    end associate
-
-  end subroutine count_source_for_ufs
-
 #ifdef _OPENMP
 !===============================================================================
 ! JOIN_BANK_FROM_THREADS joins threadprivate fission banks into a single fission

src/eigenvalue.cpp

@@ -0,0 +1,155 @@
+#include "openmc/eigenvalue.h"
+
+#include "xtensor/xmath.hpp"
+#include "xtensor/xtensor.hpp"
+#include "xtensor/xview.hpp"
+
+#include "openmc/capi.h"
+#include "openmc/error.h"
+#include "openmc/hdf5_interface.h"
+#include "openmc/mesh.h"
+#include "openmc/message_passing.h"
+#include "openmc/settings.h"
+#include "openmc/simulation.h"
+
+namespace openmc {
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+std::vector<double> entropy;
+xt::xtensor<double, 1> source_frac;
+
+//==============================================================================
+// Non-member functions
+//==============================================================================
+
+void shannon_entropy()
+{
+  // Get pointer to entropy mesh
+  auto& m = meshes[settings::index_entropy_mesh];
+
+  // Get pointer to fission bank
+  Bank* fission_bank;
+  int64_t n;
+  openmc_fission_bank(&fission_bank, &n);
+
+  // Get source weight in each mesh bin
+  bool sites_outside;
+  xt::xtensor<double, 1> p = m->count_sites(
+    n_bank, fission_bank, 0, nullptr, &sites_outside);
+
+  // display warning message if there were sites outside entropy box
+  if (sites_outside) {
+    if (mpi::master) warning("Fission source site(s) outside of entropy box.");
+  }
+
+  // sum values to obtain shannon entropy
+  if (mpi::master) {
+    // Normalize to total weight of bank sites
+    p /= xt::sum(p);
+
+    double H = 0.0;
+    for (auto p_i : p) {
+      if (p_i > 0.0) {
+        H -= p_i * std::log(p_i)/std::log(2.0);
+      }
+    }
+
+    // Add value to vector
+    entropy.push_back(H);
+  }
+}
+
+void ufs_count_sites()
+{
+  auto &m = meshes[settings::index_ufs_mesh];
+
+  if (openmc_current_batch == 1 && openmc_current_gen == 1) {
+    // On the first generation, just assume that the source is already evenly
+    // distributed so that effectively the production of fission sites is not
+    // biased
+
+    auto s = xt::view(source_frac, xt::all());
+    s = m->volume_frac_;
+
+  } else {
+    // Get pointer to source bank
+    Bank* source_bank;
+    int64_t n;
+    openmc_source_bank(&source_bank, &n);
+
+    // count number of source sites in each ufs mesh cell
+    bool sites_outside;
+    source_frac = m->count_sites(openmc_work, source_bank, 0, nullptr,
+      &sites_outside);
+
+    // Check for sites outside of the mesh
+    if (mpi::master && sites_outside) {
+      fatal_error("Source sites outside of the UFS mesh!");
+    }
+
+#ifdef OPENMC_MPI
+    // Send source fraction to all processors
+    int n_bins = xt::prod(m->shape_)();
+    MPI_Bcast(source_frac.data(), n_bins, MPI_DOUBLE, 0, mpi::intracomm);
+#endif
+
+    // Normalize to total weight to get fraction of source in each cell
+    double total = xt::sum(source_frac)();
+    source_frac /= total;
+
+    // Since the total starting weight is not equal to n_particles, we need to
+    // renormalize the weight of the source sites
+    for (int i = 0; i < openmc_work; ++i) {
+      source_bank[i].wgt *= settings::n_particles / total;
+    }
+  }
+}
+
+double ufs_get_weight(const Particle* p)
+{
+  auto& m = meshes[settings::index_ufs_mesh];
+
+  // Determine indices on ufs mesh for current location
+  // TODO: off by one
+  int mesh_bin = m->get_bin({p->coord[0].xyz}) - 1;
+  if (mesh_bin < 0) {
+    p->write_restart();
+    fatal_error("Source site outside UFS mesh!");
+  }
+
+  if (source_frac(mesh_bin) != 0.0) {
+    return m->volume_frac_ / source_frac(mesh_bin);
+  } else {
+    return 1.0;
+  }
+}
+
+extern "C" void entropy_to_hdf5(hid_t group)
+{
+  if (settings::entropy_on) {
+    write_dataset(group, "entropy", entropy);
+  }
+}
+
+extern "C" void entropy_from_hdf5(hid_t group)
+{
+  if (settings::entropy_on) {
+    read_dataset(group, "entropy", entropy);
+  }
+}
+
+extern "C" double entropy_c(int i)
+{
+  return entropy.at(i - 1);
+}
+
+extern "C" double entropy_clear()
+{
+  entropy.clear();
+}
+
+
+} // namespace openmc

src/initialize.cpp

@@ -61,7 +61,7 @@ namespace openmc {
 #ifdef OPENMC_MPI
 void initialize_mpi(MPI_Comm intracomm)
 {
-  openmc::mpi::intracomm = intracomm;
+  mpi::intracomm = intracomm;
 
   // Initialize MPI
   int flag;
@@ -69,13 +69,13 @@ void initialize_mpi(MPI_Comm intracomm)
   if (!flag) MPI_Init(nullptr, nullptr);
 
   // Determine number of processes and rank for each
-  MPI_Comm_size(intracomm, &openmc::mpi::n_procs);
-  MPI_Comm_rank(intracomm, &openmc::mpi::rank);
+  MPI_Comm_size(intracomm, &mpi::n_procs);
+  MPI_Comm_rank(intracomm, &mpi::rank);
 
   // Set variable for Fortran side
-  openmc_n_procs = openmc::mpi::n_procs;
-  openmc_rank = openmc::mpi::rank;
-  openmc_master = (openmc::mpi::rank == 0);
+  openmc_n_procs = mpi::n_procs;
+  openmc_rank = mpi::rank;
+  openmc_master = mpi::master = (mpi::rank == 0);
 
   // Create bank datatype
   Bank b;
@@ -88,8 +88,8 @@ void initialize_mpi(MPI_Comm intracomm)
   };
   int blocks[] {1, 3, 3, 1, 1};
   MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_INT};
-  MPI_Type_create_struct(5, blocks, disp, types, &openmc::mpi::bank);
-  MPI_Type_commit(&openmc::mpi::bank);
+  MPI_Type_create_struct(5, blocks, disp, types, &mpi::bank);
+  MPI_Type_commit(&mpi::bank);
 }
 #endif // OPENMC_MPI
 

src/input_xml.F90

@@ -4,7 +4,7 @@ module input_xml
 
   use algorithm,        only: find
   use cmfd_input,       only: configure_cmfd
-  use cmfd_header,      only: cmfd_mesh
+  use cmfd_header,      only: index_cmfd_mesh
   use constants
   use dict_header,      only: DictIntInt, DictCharInt, DictEntryCI
   use endf,             only: reaction_name
@@ -203,22 +203,14 @@ contains
   subroutine read_settings_xml_f(root_ptr) bind(C)
     type(C_PTR), value :: root_ptr
 
-    character(MAX_LINE_LEN) :: temp_str
     integer :: i
     integer :: n
-    integer :: temp_int
-    integer :: temp_int_array3(3)
-    integer(C_INT32_T) :: i_start, i_end
-    integer(C_INT) :: err
     integer, allocatable :: temp_int_array(:)
     integer :: n_tracks
     type(XMLNode) :: root
-    type(XMLNode) :: node_entropy
-    type(XMLNode) :: node_ufs
     type(XMLNode) :: node_sp
     type(XMLNode) :: node_res_scat
     type(XMLNode) :: node_vol
-    type(XMLNode), allocatable :: node_mesh_list(:)
     type(XMLNode), allocatable :: node_vol_list(:)
 
     ! Get proper XMLNode type given pointer
@@ -227,7 +219,6 @@ contains
     if (run_mode == MODE_EIGENVALUE) then
       ! Preallocate space for keff and entropy by generation
       call k_generation % reserve(n_max_batches*gen_per_batch)
-      call entropy % reserve(n_max_batches*gen_per_batch)
     end if
 
     ! Particle tracks
@@ -249,110 +240,6 @@ contains
       track_identifiers = reshape(temp_int_array, [3, n_tracks/3])
     end if
 
-    ! Read meshes
-    call get_node_list(root, "mesh", node_mesh_list)
-
-    ! Check for user meshes and allocate
-    n = size(node_mesh_list)
-    if (n > 0) then
-      err = openmc_extend_meshes(n, i_start, i_end)
-    end if
-
-    do i = 1, n
-      associate (m => meshes(i_start + i - 1))
-        ! Instantiate mesh from XML node
-        call m % from_xml(node_mesh_list(i))
-
-        ! Add mesh to dictionary
-        call mesh_dict % set(m % id, i_start + i - 1)
-      end associate
-    end do
-
-    ! Shannon Entropy mesh
-    if (check_for_node(root, "entropy_mesh")) then
-      call get_node_value(root, "entropy_mesh", temp_int)
-      if (mesh_dict % has(temp_int)) then
-        index_entropy_mesh = mesh_dict % get(temp_int)
-      else
-        call fatal_error("Mesh " // to_str(temp_int) // " specified for &
-             &Shannon entropy does not exist.")
-      end if
-    elseif (check_for_node(root, "entropy")) then
-      call warning("Specifying a Shannon entropy mesh via the <entropy> element &
-           &is deprecated. Please create a mesh using <mesh> and then reference &
-           &it by specifying its ID in an <entropy_mesh> element.")
-
-      ! Get pointer to entropy node
-      node_entropy = root % child("entropy")
-
-      err = openmc_extend_meshes(1, index_entropy_mesh)
-
-      associate (m => meshes(index_entropy_mesh))
-        ! Assign ID
-        m % id = 10000
-
-        call m % from_xml(node_entropy)
-      end associate
-    end if
-
-    if (index_entropy_mesh > 0) then
-      associate(m => meshes(index_entropy_mesh))
-        if (.not. allocated(m % dimension)) then
-          ! If the user did not specify how many mesh cells are to be used in
-          ! each direction, we automatically determine an appropriate number of
-          ! cells
-          m % n_dimension = 3
-          allocate(m % dimension(3))
-          m % dimension = ceiling((n_particles/20)**(ONE/THREE))
-
-          ! Calculate width
-          m % width = (m % upper_right - m % lower_left) / m % dimension
-        end if
-
-        ! Allocate space for storing number of fission sites in each mesh cell
-        allocate(entropy_p(1, product(m % dimension)))
-      end associate
-
-      ! Turn on Shannon entropy calculation
-      entropy_on = .true.
-    end if
-
-    ! Uniform fission source weighting mesh
-    if (check_for_node(root, "ufs_mesh")) then
-      call get_node_value(root, "ufs_mesh", temp_int)
-      if (mesh_dict % has(temp_int)) then
-        index_ufs_mesh = mesh_dict % get(temp_int)
-      else
-        call fatal_error("Mesh " // to_str(temp_int) // " specified for &
-             &uniform fission site method does not exist.")
-      end if
-    elseif (check_for_node(root, "uniform_fs")) then
-      call warning("Specifying a UFS mesh via the <uniform_fs> element &
-           &is deprecated. Please create a mesh using <mesh> and then reference &
-           &it by specifying its ID in a <ufs_mesh> element.")
-
-      ! Get pointer to ufs node
-      node_ufs = root % child("uniform_fs")
-
-      err = openmc_extend_meshes(1, index_ufs_mesh)
-
-      ! Allocate mesh object and coordinates on mesh
-      associate (m => meshes(index_ufs_mesh))
-        ! Assign ID
-        m % id = 10001
-
-        call m % from_xml(node_ufs)
-      end associate
-    end if
-
-    if (index_ufs_mesh > 0) then
-      ! Allocate array to store source fraction for UFS
-      allocate(source_frac(1, product(meshes(index_ufs_mesh) % dimension)))
-
-      ! Turn on uniform fission source weighting
-      ufs = .true.
-    end if
-
     ! Check if the user has specified to write state points
     if (check_for_node(root, "state_point")) then
 
@@ -1269,13 +1156,11 @@ contains
     character(MAX_WORD_LEN), allocatable :: sarray(:)
     type(DictCharInt) :: trigger_scores
     type(TallyFilterContainer), pointer :: f
-    type(RegularMesh), pointer :: m
     type(XMLDocument) :: doc
     type(XMLNode) :: root
     type(XMLNode) :: node_tal
     type(XMLNode) :: node_filt
     type(XMLNode) :: node_trigger
-    type(XMLNode), allocatable :: node_mesh_list(:)
     type(XMLNode), allocatable :: node_tal_list(:)
     type(XMLNode), allocatable :: node_filt_list(:)
     type(XMLNode), allocatable :: node_trigger_list(:)
@@ -1306,9 +1191,6 @@ contains
     ! ==========================================================================
     ! DETERMINE SIZE OF ARRAYS AND ALLOCATE
 
-    ! Get pointer list to XML <mesh>
-    call get_node_list(root, "mesh", node_mesh_list)
-
     ! Get pointer list to XML <filter>
     call get_node_list(root, "filter", node_filt_list)
 
@@ -1324,20 +1206,7 @@ contains
     ! READ MESH DATA
 
     ! Check for user meshes and allocate
-    n = size(node_mesh_list)
-    if (n > 0) then
-      err = openmc_extend_meshes(n, i_start, i_end)
-    end if
-
-    do i = 1, n
-      m => meshes(i_start + i - 1)
-
-      ! Instantiate mesh from XML node
-      call m % from_xml(node_mesh_list(i))
-
-      ! Add mesh to dictionary
-      call mesh_dict % set(m % id, i_start + i - 1)
-    end do
+    call read_meshes(root % ptr)
 
     ! We only need the mesh info for plotting
     if (run_mode == MODE_PLOTTING) then
@@ -2179,6 +2048,7 @@ contains
     integer :: i, j
     integer :: n_cols, col_id, n_comp, n_masks, n_meshlines
     integer :: meshid
+    integer(C_INT) :: err, idx
     integer, allocatable :: iarray(:)
     logical :: file_exists              ! does plots.xml file exist?
     character(MAX_LINE_LEN) :: filename ! absolute path to plots.xml
@@ -2501,7 +2371,7 @@ contains
                      // trim(to_str(pl % id)))
               end if
 
-              pl % meshlines_mesh => meshes(index_ufs_mesh)
+              pl % index_meshlines_mesh = index_ufs_mesh
 
             case ('cmfd')
 
@@ -2510,7 +2380,7 @@ contains
                      &meshlines on plot " // trim(to_str(pl % id)))
               end if
 
-              pl % meshlines_mesh => cmfd_mesh
+              pl % index_meshlines_mesh = index_cmfd_mesh
 
             case ('entropy')
 
@@ -2519,7 +2389,7 @@ contains
                      // trim(to_str(pl % id)))
               end if
 
-              pl % meshlines_mesh => meshes(index_entropy_mesh)
+              pl % index_meshlines_mesh = index_entropy_mesh
 
             case ('tally')
 
@@ -2532,17 +2402,13 @@ contains
               end if
 
               ! Check if the specified tally mesh exists
-              if (mesh_dict % has(meshid)) then
-                pl % meshlines_mesh => meshes(mesh_dict % get(meshid))
-                if (meshes(meshid) % type /= MESH_REGULAR) then
-                  call fatal_error("Non-rectangular mesh specified in &
-                       &meshlines for plot " // trim(to_str(pl % id)))
-                end if
-              else
+              err = openmc_get_mesh_index(meshid, idx)
+              if (err /= 0) then
                 call fatal_error("Could not find mesh " &
                      // trim(to_str(meshid)) // " specified in meshlines for &
                      &plot " // trim(to_str(pl % id)))
               end if
+              pl % index_meshlines_mesh = idx
 
             case default
               call fatal_error("Invalid type for meshlines on plot " &

src/mesh.F90

@@ -1,106 +0,0 @@
-module mesh
-
-  use algorithm,  only: binary_search
-  use bank_header, only: bank
-  use constants
-  use mesh_header
-  use message_passing
-
-  implicit none
-
-contains
-
-!===============================================================================
-! COUNT_BANK_SITES determines the number of fission bank sites in each cell of a
-! given mesh as well as an optional energy group structure. This can be used for
-! a variety of purposes (Shannon entropy, CMFD, uniform fission source
-! weighting)
-!===============================================================================
-
-  subroutine count_bank_sites(m, bank_array, cnt, energies, size_bank, &
-       sites_outside)
-
-    type(RegularMesh), intent(in) :: m             ! mesh to count sites
-    type(Bank), intent(in)     :: bank_array(:) ! fission or source bank
-    real(8),    intent(out)    :: cnt(:,:)      ! weight of sites in each
-    ! cell and energy group
-    real(8), intent(in),    optional :: energies(:)   ! energy grid to search
-    integer(8), intent(in), optional :: size_bank     ! # of bank sites (on each proc)
-    logical, intent(inout), optional :: sites_outside ! were there sites outside mesh?
-    real(8), allocatable :: cnt_(:,:)
-
-    integer :: i        ! loop index for local fission sites
-    integer :: n_sites  ! size of bank array
-    integer :: n        ! number of energy groups / size
-    integer :: mesh_bin ! mesh bin
-    integer :: e_bin    ! energy bin
-#ifdef OPENMC_MPI
-    integer :: mpi_err  ! MPI error code
-#endif
-    logical :: outside  ! was any site outside mesh?
-
-    ! initialize variables
-    allocate(cnt_(size(cnt,1), size(cnt,2)))
-    cnt_ = ZERO
-    outside = .false.
-
-    ! Set size of bank
-    if (present(size_bank)) then
-      n_sites = int(size_bank,4)
-    else
-      n_sites = size(bank_array)
-    end if
-
-    ! Determine number of energies in group structure
-    if (present(energies)) then
-      n = size(energies) - 1
-    else
-      n = 1
-    end if
-
-    ! loop over fission sites and count how many are in each mesh box
-    FISSION_SITES: do i = 1, n_sites
-      ! determine scoring bin for entropy mesh
-      call m % get_bin(bank_array(i) % xyz, mesh_bin)
-
-      ! if outside mesh, skip particle
-      if (mesh_bin == NO_BIN_FOUND) then
-        outside = .true.
-        cycle
-      end if
-
-      ! determine energy bin
-      if (present(energies)) then
-        if (bank_array(i) % E < energies(1)) then
-          e_bin = 1
-        elseif (bank_array(i) % E > energies(n + 1)) then
-          e_bin = n
-        else
-          e_bin = binary_search(energies, n + 1, bank_array(i) % E)
-        end if
-      else
-        e_bin = 1
-      end if
-
-      ! add to appropriate mesh box
-      cnt_(e_bin, mesh_bin) = cnt_(e_bin, mesh_bin) + bank_array(i) % wgt
-    end do FISSION_SITES
-
-#ifdef OPENMC_MPI
-    ! collect values from all processors
-    n = size(cnt_)
-    call MPI_REDUCE(cnt_, cnt, n, MPI_REAL8, MPI_SUM, 0, mpi_intracomm, mpi_err)
-
-    ! Check if there were sites outside the mesh for any processor
-    if (present(sites_outside)) then
-      call MPI_REDUCE(outside, sites_outside, 1, MPI_LOGICAL, MPI_LOR, 0, &
-           mpi_intracomm, mpi_err)
-    end if
-#else
-    sites_outside = outside
-    cnt = cnt_
-#endif
-
-  end subroutine count_bank_sites
-
-end module mesh

src/mesh.cpp

@@ -0,0 +1,976 @@
+#include "openmc/mesh.h"
+
+#include <algorithm> // for copy, min
+#include <cstddef> // for size_t
+#include <cmath>  // for ceil
+#include <string>
+
+#ifdef OPENMC_MPI
+#include "mpi.h"
+#endif
+#include "xtensor/xbuilder.hpp"
+#include "xtensor/xeval.hpp"
+#include "xtensor/xmath.hpp"
+#include "xtensor/xsort.hpp"
+#include "xtensor/xtensor.hpp"
+
+#include "openmc/capi.h"
+#include "openmc/constants.h"
+#include "openmc/error.h"
+#include "openmc/hdf5_interface.h"
+#include "openmc/message_passing.h"
+#include "openmc/search.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+std::vector<std::unique_ptr<RegularMesh>> meshes;
+
+std::unordered_map<int32_t, int32_t> mesh_map;
+
+//==============================================================================
+// RegularMesh implementation
+//==============================================================================
+
+RegularMesh::RegularMesh(pugi::xml_node node)
+{
+  // Copy mesh id
+  if (check_for_node(node, "id")) {
+    id_ = std::stoi(get_node_value(node, "id"));
+
+    // Check to make sure 'id' hasn't been used
+    if (mesh_map.find(id_) != mesh_map.end()) {
+      fatal_error("Two or more meshes use the same unique ID: " +
+        std::to_string(id_));
+    }
+  }
+
+  // Read mesh type
+  if (check_for_node(node, "type")) {
+    auto temp = get_node_value(node, "type", true, true);
+    if (temp == "regular") {
+      // TODO: move elsewhere
+    } else {
+      fatal_error("Invalid mesh type: " + temp);
+    }
+  }
+
+  // Determine number of dimensions for mesh
+  if (check_for_node(node, "dimension")) {
+    shape_ = get_node_xarray<int>(node, "dimension");
+    int n = n_dimension_ = shape_.size();
+    if (n != 1 && n != 2 && n != 3) {
+      fatal_error("Mesh must be one, two, or three dimensions.");
+    }
+
+    // Check that dimensions are all greater than zero
+    if (xt::any(shape_ <= 0)) {
+      fatal_error("All entries on the <dimension> element for a tally "
+        "mesh must be positive.");
+    }
+  }
+
+  // Check for lower-left coordinates
+  if (check_for_node(node, "lower_left")) {
+    // Read mesh lower-left corner location
+    lower_left_ = get_node_xarray<double>(node, "lower_left");
+  } else {
+    fatal_error("Must specify <lower_left> on a mesh.");
+  }
+
+  if (check_for_node(node, "width")) {
+    // Make sure both upper-right or width were specified
+    if (check_for_node(node, "upper_right")) {
+      fatal_error("Cannot specify both <upper_right> and <width> on a mesh.");
+    }
+
+    width_ = get_node_xarray<double>(node, "width");
+
+    // Check to ensure width has same dimensions
+    auto n = width_.size();
+    if (n != lower_left_.size()) {
+      fatal_error("Number of entries on <width> must be the same as "
+        "the number of entries on <lower_left>.");
+    }
+
+    // Check for negative widths
+    if (xt::any(width_ < 0.0)) {
+      fatal_error("Cannot have a negative <width> on a tally mesh.");
+    }
+
+    // Set width and upper right coordinate
+    upper_right_ = xt::eval(lower_left_ + shape_ * width_);
+
+  } else if (check_for_node(node, "upper_right")) {
+    upper_right_ = get_node_xarray<double>(node, "upper_right");
+
+    // Check to ensure width has same dimensions
+    auto n = upper_right_.size();
+    if (n != lower_left_.size()) {
+      fatal_error("Number of entries on <upper_right> must be the "
+        "same as the number of entries on <lower_left>.");
+    }
+
+    // Check that upper-right is above lower-left
+    if (xt::any(upper_right_ < lower_left_)) {
+      fatal_error("The <upper_right> coordinates must be greater than "
+        "the <lower_left> coordinates on a tally mesh.");
+    }
+
+    // Set width and upper right coordinate
+    width_ = xt::eval((upper_right_ - lower_left_) / shape_);
+  } else {
+    fatal_error("Must specify either <upper_right> and <width> on a mesh.");
+  }
+
+  if (shape_.dimension() > 0) {
+    if (shape_.size() != lower_left_.size()) {
+      fatal_error("Number of entries on <lower_left> must be the same "
+        "as the number of entries on <dimension>.");
+    }
+
+    // Set volume fraction
+    volume_frac_ = 1.0/xt::prod(shape_)();
+  }
+}
+
+int RegularMesh::get_bin(Position r) const
+{
+  // Loop over the dimensions of the mesh
+  for (int i = 0; i < n_dimension_; ++i) {
+    // Check for cases where particle is outside of mesh
+    if (r[i] < lower_left_[i]) {
+      return -1;
+    } else if (r[i] > upper_right_[i]) {
+      return -1;
+    }
+  }
+
+  // Determine indices
+  int ijk[n_dimension_];
+  bool in_mesh;
+  get_indices(r, ijk, &in_mesh);
+  if (!in_mesh) return -1;
+
+  // Convert indices to bin
+  return get_bin_from_indices(ijk);
+}
+
+int RegularMesh::get_bin_from_indices(const int* ijk) const
+{
+  if (n_dimension_ == 1) {
+    return ijk[0];
+  } else if (n_dimension_ == 2) {
+    return (ijk[1] - 1)*shape_[0] + ijk[0];
+  } else if (n_dimension_ == 3) {
+    return ((ijk[2] - 1)*shape_[1] + (ijk[1] - 1))*shape_[0] + ijk[0];
+  }
+}
+
+void RegularMesh::get_indices(Position r, int* ijk, bool* in_mesh) const
+{
+  // Find particle in mesh
+  *in_mesh = true;
+  for (int i = 0; i < n_dimension_; ++i) {
+    ijk[i] = std::ceil((r[i] - lower_left_[i]) / width_[i]);
+
+    // Check if indices are within bounds
+    if (ijk[i] < 1 || ijk[i] > shape_[i]) *in_mesh = false;
+  }
+}
+
+void RegularMesh::get_indices_from_bin(int bin, int* ijk) const
+{
+  if (n_dimension_ == 1) {
+    ijk[0] = bin;
+  } else if (n_dimension_ == 2) {
+    ijk[0] = (bin - 1) % shape_[0] + 1;
+    ijk[1] = (bin - 1) / shape_[0] + 1;
+  } else if (n_dimension_ == 3) {
+    ijk[0] = (bin - 1) % shape_[0] + 1;
+    ijk[1] = ((bin - 1) % (shape_[0] * shape_[1])) / shape_[0] + 1;
+    ijk[2] = (bin - 1) / (shape_[0] * shape_[1]) + 1;
+  }
+}
+
+bool RegularMesh::intersects(Position r0, Position r1) const
+{
+  switch(n_dimension_) {
+    case 1:
+      return intersects_1d(r0, r1);
+    case 2:
+      return intersects_2d(r0, r1);
+    case 3:
+      return intersects_3d(r0, r1);
+  }
+}
+
+bool RegularMesh::intersects_1d(Position r0, Position r1) const
+{
+  // Copy coordinates of mesh lower_left and upper_right
+  double left = lower_left_[0];
+  double right = upper_right_[0];
+
+  // Check if line intersects either left or right surface
+  if (r0.x < left) {
+    return r1.x > left;
+  } else if (r0.x < right) {
+    return r1.x < left || r1.x > right;
+  } else {
+    return r1.x < right;
+  }
+}
+
+bool RegularMesh::intersects_2d(Position r0, Position r1) const
+{
+  // Copy coordinates of starting point
+  double x0 = r0.x;
+  double y0 = r0.y;
+
+  // Copy coordinates of ending point
+  double x1 = r1.x;
+  double y1 = r1.y;
+
+  // Copy coordinates of mesh lower_left
+  double xm0 = lower_left_[0];
+  double ym0 = lower_left_[1];
+
+  // Copy coordinates of mesh upper_right
+  double xm1 = upper_right_[0];
+  double ym1 = upper_right_[1];
+
+  // Check if line intersects left surface -- calculate the intersection point y
+  if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) {
+    double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0);
+    if (yi >= ym0 && yi < ym1) {
+      return true;
+    }
+  }
+
+  // Check if line intersects back surface -- calculate the intersection point
+  // x
+  if ((y0 < ym0 && y1 > ym0) || (y0 > ym0 && y1 < ym0)) {
+    double xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0);
+    if (xi >= xm0 && xi < xm1) {
+      return true;
+    }
+  }
+
+  // Check if line intersects right surface -- calculate the intersection
+  // point y
+  if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) {
+    double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0);
+    if (yi >= ym0 && yi < ym1) {
+      return true;
+    }
+  }
+
+  // Check if line intersects front surface -- calculate the intersection point
+  // x
+  if ((y0 < ym1 && y1 > ym1) || (y0 > ym1 && y1 < ym1)) {
+    double xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0);
+    if (xi >= xm0 && xi < xm1) {
+      return true;
+    }
+  }
+  return false;
+}
+
+bool RegularMesh::intersects_3d(Position r0, Position r1) const
+{
+  // Copy coordinates of starting point
+  double x0 = r0.x;
+  double y0 = r0.y;
+  double z0 = r0.z;
+
+  // Copy coordinates of ending point
+  double x1 = r1.x;
+  double y1 = r1.y;
+  double z1 = r1.z;
+
+  // Copy coordinates of mesh lower_left
+  double xm0 = lower_left_[0];
+  double ym0 = lower_left_[1];
+  double zm0 = lower_left_[2];
+
+  // Copy coordinates of mesh upper_right
+  double xm1 = upper_right_[0];
+  double ym1 = upper_right_[1];
+  double zm1 = upper_right_[2];
+
+  // Check if line intersects left surface -- calculate the intersection point
+  // (y,z)
+  if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) {
+    double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0);
+    double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0);
+    if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) {
+      return true;
+    }
+  }
+
+  // Check if line intersects back surface -- calculate the intersection point
+  // (x,z)
+  if ((y0 < ym0 && y1 > ym0) || (y0 > ym0 && y1 < ym0)) {
+    double xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0);
+    double zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0);
+    if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) {
+      return true;
+    }
+  }
+
+  // Check if line intersects bottom surface -- calculate the intersection
+  // point (x,y)
+  if ((z0 < zm0 && z1 > zm0) || (z0 > zm0 && z1 < zm0)) {
+    double xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0);
+    double yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0);
+    if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) {
+      return true;
+    }
+  }
+
+  // Check if line intersects right surface -- calculate the intersection point
+  // (y,z)
+  if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) {
+    double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0);
+    double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0);
+    if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) {
+      return true;
+    }
+  }
+
+  // Check if line intersects front surface -- calculate the intersection point
+  // (x,z)
+  if ((y0 < ym1 && y1 > ym1) || (y0 > ym1 && y1 < ym1)) {
+    double xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0);
+    double zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0);
+    if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) {
+      return true;
+    }
+  }
+
+  // Check if line intersects top surface -- calculate the intersection point
+  // (x,y)
+  if ((z0 < zm1 && z1 > zm1) || (z0 > zm1 && z1 < zm1)) {
+    double xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0);
+    double yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0);
+    if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) {
+      return true;
+    }
+  }
+  return false;
+}
+
+void RegularMesh::bins_crossed(const Particle* p, std::vector<int>& bins,
+                               std::vector<double>& lengths) const
+{
+  constexpr int MAX_SEARCH_ITER = 100;
+
+  // ========================================================================
+  // Determine if the track intersects the tally mesh.
+
+  // Copy the starting and ending coordinates of the particle.  Offset these
+  // just a bit for the purposes of determining if there was an intersection
+  // in case the mesh surfaces coincide with lattice/geometric surfaces which
+  // might produce finite-precision errors.
+  Position last_r {p->last_xyz};
+  Position r {p->coord[0].xyz};
+  Direction u {p->coord[0].uvw};
+
+  Position r0 = last_r + TINY_BIT*u;
+  Position r1 = r - TINY_BIT*u;
+
+  // Determine indices for starting and ending location.
+  int n = n_dimension_;
+  xt::xtensor<int, 1> ijk0 = xt::empty<int>({n});
+  bool start_in_mesh;
+  get_indices(r0, ijk0.data(), &start_in_mesh);
+  xt::xtensor<int, 1> ijk1 = xt::empty<int>({n});
+  bool end_in_mesh;
+  get_indices(r1, ijk1.data(), &end_in_mesh);
+
+  // Check if the track intersects any part of the mesh.
+  if (!start_in_mesh && !end_in_mesh) {
+    if (!intersects(r0, r1)) return;
+  }
+
+  // ========================================================================
+  // Figure out which mesh cell to tally.
+
+  // Copy the un-modified coordinates the particle direction.
+  r0 = last_r;
+  r1 = r;
+
+  // Compute the length of the entire track.
+  double total_distance = (r1 - r0).norm();
+
+  // We are looking for the first valid mesh bin.  Check to see if the
+  // particle starts inside the mesh.
+  if (!start_in_mesh) {
+    xt::xtensor<double, 1> d = xt::zeros<double>({n});
+
+    // The particle does not start in the mesh.  Note that we nudged the
+    // start and end coordinates by a TINY_BIT each so we will have
+    // difficulty resolving tracks that are less than 2*TINY_BIT in length.
+    // If the track is that short, it is also insignificant so we can
+    // safely ignore it in the tallies.
+    if (total_distance < 2*TINY_BIT) return;
+
+    // The particle does not start in the mesh so keep iterating the ijk0
+    // indices to cross the nearest mesh surface until we've found a valid
+    // bin.  MAX_SEARCH_ITER prevents an infinite loop.
+    int search_iter = 0;
+    int j;
+    while (xt::any(ijk0 < 1) || xt::any(ijk0 > shape_)) {
+      if (search_iter == MAX_SEARCH_ITER) {
+        warning("Failed to find a mesh intersection on a tally mesh filter.");
+        return;
+      }
+
+      for (j = 0; j < n; ++j) {
+        if (std::fabs(u[j]) < FP_PRECISION) {
+          d(j) = INFTY;
+        } else if (u[j] > 0.0) {
+          double xyz_cross = lower_left_[j] + ijk0(j) * width_[j];
+          d(j) = (xyz_cross - r0[j]) / u[j];
+        } else {
+          double xyz_cross = lower_left_[j] + (ijk0(j) - 1) * width_[j];
+          d(j) = (xyz_cross - r0[j]) / u[j];
+        }
+      }
+
+      j = xt::argmin(d)(0);
+      if (u[j] > 0.0) {
+        ++ijk0(j);
+      } else {
+        --ijk0(j);
+      }
+
+      ++search_iter;
+    }
+
+    // Advance position
+    r0 += d(j) * u;
+  }
+
+  while (true) {
+    // ========================================================================
+    // Compute the length of the track segment in the each mesh cell and return
+
+    double distance;
+    int j;
+    if (ijk0 == ijk1) {
+      // The track ends in this cell.  Use the particle end location rather
+      // than the mesh surface.
+      distance = (r1 - r0).norm();
+    } else {
+      // The track exits this cell.  Determine the distance to the closest mesh
+      // surface.
+      xt::xtensor<double, 1> d = xt::zeros<double>({n});
+      for (int j = 0; j < n; ++j) {
+        if (std::fabs(u[j]) < FP_PRECISION) {
+          d(j) = INFTY;
+        } else if (u[j] > 0) {
+          double xyz_cross = lower_left_[j] + ijk0(j) * width_[j];
+          d(j) = (xyz_cross - r0[j]) / u[j];
+        } else {
+          double xyz_cross = lower_left_[j] + (ijk0(j) - 1) * width_[j];
+          d(j) = (xyz_cross - r0[j]) / u[j];
+        }
+      }
+      j = xt::argmin(d)(0);
+      distance = d(j);
+    }
+
+    // Assign the next tally bin and the score.
+    int bin = get_bin_from_indices(ijk0.data());
+    bins.push_back(bin);
+    lengths.push_back(distance / total_distance);
+
+    // If the particle track ends in that bin, then we are done.
+    if (ijk0 == ijk1) break;
+
+    // Translate the starting coordintes by the distance to that face. This
+    // should be the xyz that we computed the distance to in the last
+    // iteration of the filter loop.
+    r0 += distance * u;
+
+    // Increment the indices into the next mesh cell.
+    if (u[j] > 0.0) {
+      ++ijk0(j);
+    } else {
+      --ijk0(j);
+    }
+
+    // If the next indices are invalid, then the track has left the mesh and
+    // we are done.
+    if (xt::any(ijk0 < 1) || xt::any(ijk0 > shape_)) break;
+  }
+}
+
+void RegularMesh::surface_bins_crossed(const Particle* p, std::vector<int>& bins) const
+{
+  // ========================================================================
+  // Determine if the track intersects the tally mesh.
+
+  // Copy the starting and ending coordinates of the particle.
+  Position r0 {p->last_xyz_current};
+  Position r1 {p->coord[0].xyz};
+  Direction u {p->coord[0].uvw};
+
+  // Determine indices for starting and ending location.
+  int n = n_dimension_;
+  xt::xtensor<int, 1> ijk0 = xt::empty<int>({n});
+  bool start_in_mesh;
+  get_indices(r0, ijk0.data(), &start_in_mesh);
+  xt::xtensor<int, 1> ijk1 = xt::empty<int>({n});
+  bool end_in_mesh;
+  get_indices(r1, ijk1.data(), &end_in_mesh);
+
+  // Check if the track intersects any part of the mesh.
+  if (!start_in_mesh && !end_in_mesh) {
+    if (!intersects(r0, r1)) return;
+  }
+
+  // ========================================================================
+  // Figure out which mesh cell to tally.
+
+  // Calculate number of surface crossings
+  int n_cross = xt::sum(xt::abs(ijk1 - ijk0))();
+  if (n_cross == 0) return;
+
+  // Bounding coordinates
+  Position xyz_cross;
+  for (int i = 0; i < n; ++i) {
+    if (u[i] > 0.0) {
+      xyz_cross[i] = lower_left_[i] + ijk0[i] * width_[i];
+    } else {
+      xyz_cross[i] = lower_left_[i] + (ijk0[i] - 1) * width_[i];
+    }
+  }
+
+  for (int j = 0; j < n_cross; ++j) {
+    // Set the distances to infinity
+    Position d {INFTY, INFTY, INFTY};
+
+    // Determine closest bounding surface. We need to treat
+    // special case where the cosine of the angle is zero since this would
+    // result in a divide-by-zero.
+    double distance = INFTY;
+    for (int i = 0; i < n; ++i) {
+      if (u[i] == 0) {
+        d[i] = INFINITY;
+      } else {
+        d[i] = (xyz_cross[i] - r0[i])/u[i];
+      }
+      distance = std::min(distance, d[i]);
+    }
+
+    // Loop over the dimensions
+    for (int i = 0; i < n; ++i) {
+      // Check whether distance is the shortest distance
+      if (distance == d[i]) {
+
+        // Check whether particle is moving in positive i direction
+        if (u[i] > 0) {
+
+          // Outward current on i max surface
+          if (xt::all(ijk0 >= 1) && xt::all(ijk0 <= shape_)) {
+            int i_surf = 4*i + 3;
+            int i_mesh = get_bin_from_indices(ijk0.data());
+            int i_bin = 4*n*(i_mesh - 1) + i_surf;
+
+            bins.push_back(i_bin);
+          }
+
+          // Advance position
+          ++ijk0[i];
+          xyz_cross[i] += width_[i];
+
+          // If the particle crossed the surface, tally the inward current on
+          // i min surface
+          if (xt::all(ijk0 >= 1) && xt::all(ijk0 <= shape_)) {
+            int i_surf = 4*i + 2;
+            int i_mesh = get_bin_from_indices(ijk0.data());
+            int i_bin = 4*n*(i_mesh - 1) + i_surf;
+
+            bins.push_back(i_bin);
+          }
+
+        } else {
+          // The particle is moving in the negative i direction
+
+          // Outward current on i min surface
+          if (xt::all(ijk0 >= 1) && xt::all(ijk0 <= shape_) ){
+            int i_surf = 4*i + 1;
+            int i_mesh = get_bin_from_indices(ijk0.data());
+            int i_bin = 4*n*(i_mesh - 1) + i_surf;
+
+            bins.push_back(i_bin);
+          }
+
+          // Advance position
+          --ijk0[i];
+          xyz_cross[i] -= width_[i];
+
+          // If the particle crossed the surface, tally the inward current on
+          // i max surface
+          if (xt::all(ijk0 >= 1) && xt::all(ijk0 <= shape_)) {
+            int i_surf = 4*i + 4;
+            int i_mesh = get_bin_from_indices(ijk0.data());
+            int i_bin = 4*n*(i_mesh - 1) + i_surf;
+
+            bins.push_back(i_bin);
+          }
+        }
+      }
+    }
+
+    // Calculate new coordinates
+    r0 += distance * u;
+  }
+}
+
+void RegularMesh::to_hdf5(hid_t group) const
+{
+  hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_));
+
+  write_dataset(mesh_group, "type", "regular");
+  write_dataset(mesh_group, "dimension", shape_);
+  write_dataset(mesh_group, "lower_left", lower_left_);
+  write_dataset(mesh_group, "upper_right", upper_right_);
+  write_dataset(mesh_group, "width", width_);
+
+  close_group(mesh_group);
+}
+
+xt::xarray<double> RegularMesh::count_sites(int64_t n, const Bank* bank,
+  int n_energy, const double* energies, bool* outside) const
+{
+  // Determine shape of array for counts
+  std::size_t m = xt::prod(shape_)();
+  std::vector<std::size_t> shape;
+  if (n_energy > 0) {
+    shape = {m, static_cast<std::size_t>(n_energy - 1)};
+  } else {
+    shape = {m};
+  }
+
+  // Create array of zeros
+  xt::xarray<double> cnt {shape, 0.0};
+  bool outside_ = false;
+
+  for (int64_t i = 0; i < n; ++i) {
+    // determine scoring bin for entropy mesh
+    // TODO: off-by-one
+    int mesh_bin = get_bin({bank[i].xyz}) - 1;
+
+    // if outside mesh, skip particle
+    if (mesh_bin < 0) {
+      outside_ = true;
+      continue;
+    }
+
+    if (n_energy > 0) {
+      double E = bank[i].E;
+      if (E >= energies[0] && E <= energies[n_energy - 1]) {
+        // determine energy bin
+        int e_bin = lower_bound_index(energies, energies + n_energy, E);
+
+        // Add to appropriate bin
+        cnt(mesh_bin, e_bin) += bank[i].wgt;
+      }
+    } else {
+      // Add to appropriate bin
+      cnt(mesh_bin) += bank[i].wgt;
+    }
+  }
+
+  // Create copy of count data
+  int total = cnt.size();
+  double* cnt_reduced = new double[total];
+
+#ifdef OPENMC_MPI
+  // collect values from all processors
+  MPI_Reduce(cnt.data(), cnt_reduced, total, MPI_DOUBLE, MPI_SUM, 0,
+    mpi::intracomm);
+
+  // Check if there were sites outside the mesh for any processor
+  if (outside) {
+    MPI_Reduce(&outside_, outside, 1, MPI_C_BOOL, MPI_LOR, 0, mpi::intracomm);
+  }
+#else
+  std::copy(cnt.data(), cnt.data() + total, cnt_reduced);
+  if (outside) *outside = outside_;
+#endif
+
+  // Adapt reduced values in array back into an xarray
+  auto arr = xt::adapt(cnt_reduced, total, xt::acquire_ownership(), shape);
+  xt::xarray<double> counts = arr;
+
+  return counts;
+}
+
+//==============================================================================
+// C API functions
+//==============================================================================
+
+//! Extend the meshes array by n elements
+extern "C" int
+openmc_extend_meshes(int32_t n, int32_t* index_start, int32_t* index_end)
+{
+  if (index_start) *index_start = meshes.size();
+  for (int i = 0; i < n; ++i) {
+    meshes.emplace_back(new RegularMesh{});
+  }
+  if (index_end) *index_end = meshes.size() - 1;
+
+  return 0;
+}
+
+//! Return the index in the meshes array of a mesh with a given ID
+extern "C" int
+openmc_get_mesh_index(int32_t id, int32_t* index)
+{
+  auto pair = mesh_map.find(id);
+  if (pair == mesh_map.end()) {
+    set_errmsg("No mesh exists with ID=" + std::to_string(id) + ".");
+    return OPENMC_E_INVALID_ID;
+  }
+  *index = pair->second;
+  return 0;
+}
+
+// Return the ID of a mesh
+extern "C" int
+openmc_mesh_get_id(int32_t index, int32_t* id)
+{
+  if (index < 0 || index >= meshes.size()) {
+    set_errmsg("Index in meshes array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+  *id = meshes[index]->id_;
+  return 0;
+}
+
+//! Set the ID of a mesh
+extern "C" int
+openmc_mesh_set_id(int32_t index, int32_t id)
+{
+  if (index < 0 || index >= meshes.size()) {
+    set_errmsg("Index in meshes array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+  meshes[index]->id_ = id;
+  mesh_map[id] = index;
+  return 0;
+}
+
+//! Get the dimension of a mesh
+extern "C" int
+openmc_mesh_get_dimension(int32_t index, int** dims, int* n)
+{
+  if (index < 0 || index >= meshes.size()) {
+    set_errmsg("Index in meshes array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+  *dims = meshes[index]->shape_.data();
+  *n = meshes[index]->n_dimension_;
+  return 0;
+}
+
+//! Set the dimension of a mesh
+extern "C" int
+openmc_mesh_set_dimension(int32_t index, int n, const int* dims)
+{
+  if (index < 0 || index >= meshes.size()) {
+    set_errmsg("Index in meshes array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+
+  // Copy dimension
+  std::vector<std::size_t> shape = {static_cast<std::size_t>(n)};
+  auto& m = meshes[index];
+  m->shape_ = xt::adapt(dims, n, xt::no_ownership(), shape);
+  m->n_dimension_ = m->shape_.size();
+
+  return 0;
+}
+
+//! Get the mesh parameters
+extern "C" int
+openmc_mesh_get_params(int32_t index, double** ll, double** ur, double** width, int* n)
+{
+  if (index < 0 || index >= meshes.size()) {
+    set_errmsg("Index in meshes array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+
+  auto& m = meshes[index];
+  if (m->lower_left_.dimension() == 0) {
+    set_errmsg("Mesh parameters have not been set.");
+    return OPENMC_E_ALLOCATE;
+  }
+
+  *ll = m->lower_left_.data();
+  *ur = m->upper_right_.data();
+  *width = m->width_.data();
+  *n = m->n_dimension_;
+  return 0;
+}
+
+//! Set the mesh parameters
+extern "C" int
+openmc_mesh_set_params(int32_t index, int n, const double* ll, const double* ur,
+                       const double* width)
+{
+  if (index < 0 || index >= meshes.size()) {
+    set_errmsg("Index in meshes array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+
+  auto& m = meshes[index];
+  std::vector<std::size_t> shape = {static_cast<std::size_t>(n)};
+  if (ll && ur) {
+    m->lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape);
+    m->upper_right_ = xt::adapt(ur, n, xt::no_ownership(), shape);
+    m->width_ = (m->upper_right_ - m->lower_left_) / m->shape_;
+  } else if (ll && width) {
+    m->lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape);
+    m->width_ = xt::adapt(width, n, xt::no_ownership(), shape);
+    m->upper_right_ = m->lower_left_ + m->shape_ * m->width_;
+  } else if (ur && width) {
+    m->upper_right_ = xt::adapt(ur, n, xt::no_ownership(), shape);
+    m->width_ = xt::adapt(width, n, xt::no_ownership(), shape);
+    m->lower_left_ = m->upper_right_ - m->shape_ * m->width_;
+  } else {
+    set_errmsg("At least two parameters must be specified.");
+    return OPENMC_E_INVALID_ARGUMENT;
+  }
+
+  return 0;
+}
+
+//==============================================================================
+// Non-member functions
+//==============================================================================
+
+void read_meshes(pugi::xml_node* root)
+{
+  for (auto node : root->children("mesh")) {
+    // Read mesh and add to vector
+    meshes.emplace_back(new RegularMesh{node});
+
+    // Map ID to position in vector
+    mesh_map[meshes.back()->id_] = meshes.size() - 1;
+  }
+}
+
+void meshes_to_hdf5(hid_t group)
+{
+  // Write number of meshes
+  hid_t meshes_group = create_group(group, "meshes");
+  int32_t n_meshes = meshes.size();
+  write_attribute(meshes_group, "n_meshes", n_meshes);
+
+  if (n_meshes > 0) {
+    // Write IDs of meshes
+    std::vector<int> ids;
+    for (const auto& m : meshes) {
+      m->to_hdf5(meshes_group);
+      ids.push_back(m->id_);
+    }
+    write_attribute(meshes_group, "ids", ids);
+  }
+
+  close_group(meshes_group);
+}
+
+//==============================================================================
+// Fortran compatibility
+//==============================================================================
+
+extern "C" {
+  // Declaration of Fortran procedures
+  void vector_int_push_back(void* ptr, int value);
+  void vector_real_push_back(void* ptr, double value);
+
+  int n_meshes() { return meshes.size(); }
+
+  RegularMesh* mesh_ptr(int i) { return meshes.at(i).get(); }
+
+  int32_t mesh_id(RegularMesh* m) { return m->id_; }
+
+  double mesh_volume_frac(RegularMesh* m) { return m->volume_frac_; }
+
+  int mesh_n_dimension(RegularMesh* m) { return m->n_dimension_; }
+
+  int mesh_dimension(RegularMesh* m, int i) { return m->shape_(i - 1); }
+
+  double mesh_lower_left(RegularMesh* m, int i) { return m->lower_left_(i - 1); }
+
+  double mesh_upper_right(RegularMesh* m, int i) { return m->upper_right_(i - 1); }
+
+  double mesh_width(RegularMesh* m, int i) { return m->width_(i - 1); }
+
+  int mesh_get_bin(RegularMesh* m, const double* xyz)
+  {
+    return m->get_bin({xyz});
+  }
+
+  int mesh_get_bin_from_indices(RegularMesh* m, const int* ijk)
+  {
+    return m->get_bin_from_indices(ijk);
+  }
+
+  void mesh_get_indices(RegularMesh* m, const double* xyz, int* ijk, bool* in_mesh)
+  {
+    m->get_indices({xyz}, ijk, in_mesh);
+  }
+
+  void mesh_get_indices_from_bin(RegularMesh* m, int bin, int* ijk)
+  {
+    m->get_indices_from_bin(bin, ijk);
+  }
+
+  void mesh_bins_crossed(RegularMesh* m, const Particle* p, void* match_bins,
+    void* match_weights)
+  {
+    // Get bins crossed
+    std::vector<int> bins;
+    std::vector<double> lengths;
+    m->bins_crossed(p, bins, lengths);
+
+    // Call bindings for VectorInt and VectorReal on Fortran side
+    for (int i = 0; i < bins.size(); ++i) {
+      vector_int_push_back(match_bins, bins[i]);
+      vector_real_push_back(match_weights, lengths[i]);
+    }
+  }
+
+  void mesh_surface_bins_crossed(RegularMesh* m, const Particle* p,
+    void* match_bins, void* match_weights)
+  {
+    // Get surface bins crossed
+    std::vector<int> bins;
+    m->surface_bins_crossed(p, bins);
+
+    // Call bindings for VectorInt and VectorReal
+    for (auto b : bins) {
+      vector_int_push_back(match_bins, b);
+      vector_real_push_back(match_weights, 1.0);
+    }
+  }
+
+  void free_memory_mesh()
+  {
+    meshes.clear();
+    mesh_map.clear();
+  }
+}
+
+
+} // namespace openmc

src/mesh_header.F90

@@ -2,185 +2,242 @@ module mesh_header
 
   use, intrinsic :: ISO_C_BINDING
 
-  use constants
-  use dict_header, only: DictIntInt
-  use error
-  use hdf5_interface
-  use string, only: to_str, to_lower
-  use xml_interface
-
   implicit none
-  private
-  public :: free_memory_mesh
-  public :: openmc_extend_meshes
-  public :: openmc_get_mesh_index
-  public :: openmc_mesh_get_id
-  public :: openmc_mesh_get_dimension
-  public :: openmc_mesh_get_params
-  public :: openmc_mesh_set_id
-  public :: openmc_mesh_set_dimension
-  public :: openmc_mesh_set_params
 
 !===============================================================================
 ! STRUCTUREDMESH represents a tessellation of n-dimensional Euclidean space by
 ! congruent squares or cubes
 !===============================================================================
 
-  type, public :: RegularMesh
-    integer :: id = -1                     ! user-specified id
-    integer :: type = MESH_REGULAR         ! rectangular, hexagonal
-    integer(C_INT) :: n_dimension          ! rank of mesh
-    real(8) :: volume_frac                 ! volume fraction of each cell
-    integer(C_INT), allocatable :: dimension(:)   ! number of cells in each direction
-    real(C_DOUBLE), allocatable :: lower_left(:)  ! lower-left corner of mesh
-    real(C_DOUBLE), allocatable :: upper_right(:) ! upper-right corner of mesh
-    real(C_DOUBLE), allocatable :: width(:)       ! width of each mesh cell
+  type :: RegularMesh
+    type(C_PTR) :: ptr
   contains
-    procedure :: from_xml => regular_from_xml
+    procedure :: id => regular_id
+    procedure :: volume_frac => regular_volume_frac
+    procedure :: n_dimension => regular_n_dimension
+    procedure :: dimension => regular_dimension
+    procedure :: lower_left => regular_lower_left
+    procedure :: upper_right => regular_upper_right
+    procedure :: width => regular_width
+
     procedure :: get_bin => regular_get_bin
     procedure :: get_indices => regular_get_indices
     procedure :: get_bin_from_indices => regular_get_bin_from_indices
     procedure :: get_indices_from_bin => regular_get_indices_from_bin
-    procedure :: intersects => regular_intersects
-    procedure :: to_hdf5 => regular_to_hdf5
   end type RegularMesh
 
-  integer(C_INT32_T), public, bind(C) :: n_meshes = 0 ! # of structured meshes
+  interface
+    function openmc_extend_meshes(n, index_start, index_end) result(err) bind(C)
+      import C_INT32_T, C_INT
+      integer(C_INT32_T), value, intent(in) :: n
+      integer(C_INT32_T), optional, intent(out) :: index_start
+      integer(C_INT32_T), optional, intent(out) :: index_end
+      integer(C_INT) :: err
+    end function openmc_extend_meshes
 
-  type(RegularMesh), public, allocatable, target :: meshes(:)
+    function openmc_get_mesh_index(id, index) result(err) bind(C)
+      import C_INT32_T, C_INT
+      integer(C_INT32_T), value :: id
+      integer(C_INT32_T), intent(out) :: index
+      integer(C_INT) :: err
+    end function openmc_get_mesh_index
 
-  ! Dictionary that maps user IDs to indices in 'meshes'
-  type(DictIntInt), public :: mesh_dict
+    function openmc_mesh_get_id(index, id) result(err) bind(C)
+      import C_INT32_T, C_INT
+      integer(C_INT32_T), value       :: index
+      integer(C_INT32_T), intent(out) :: id
+      integer(C_INT) :: err
+    end function openmc_mesh_get_id
+
+    function openmc_mesh_set_id(index, id) result(err) bind(C)
+      import C_INT32_T, C_INT
+      integer(C_INT32_T), value, intent(in) :: index
+      integer(C_INT32_T), value, intent(in) :: id
+      integer(C_INT) :: err
+    end function openmc_mesh_set_id
+
+    function openmc_mesh_get_dimension(index, dims, n) result(err) bind(C)
+      import C_INT32_T, C_PTR, C_INT
+      integer(C_INT32_T), value, intent(in) :: index
+      type(C_PTR),               intent(out) :: dims
+      integer(C_INT),            intent(out) :: n
+      integer(C_INT) :: err
+    end function openmc_mesh_get_dimension
+
+    function openmc_mesh_set_dimension(index, n, dims) result(err) bind(C)
+      import C_INT32_T, C_INT
+      integer(C_INT32_T), value, intent(in) :: index
+      integer(C_INT),     value, intent(in) :: n
+      integer(C_INT),            intent(in) :: dims(n)
+      integer(C_INT) :: err
+    end function openmc_mesh_set_dimension
+
+    function openmc_mesh_get_params(index, ll, ur, width, n) result(err) bind(C)
+      import C_INT32_T, C_PTR, C_INT
+      integer(C_INT32_T), value, intent(in) :: index
+      type(C_PTR), intent(out) :: ll
+      type(C_PTR), intent(out) :: ur
+      type(C_PTR), intent(out) :: width
+      integer(C_INT), intent(out) :: n
+      integer(C_INT) :: err
+    end function openmc_mesh_get_params
+
+    function openmc_mesh_set_params(index, n, ll, ur, width) result(err) bind(C)
+      import C_INT32_T, C_INT, C_DOUBLE
+      integer(C_INT32_T), value, intent(in) :: index
+      integer(C_INT),     value, intent(in) :: n
+      real(C_DOUBLE), intent(in), optional :: ll(n)
+      real(C_DOUBLE), intent(in), optional :: ur(n)
+      real(C_DOUBLE), intent(in), optional :: width(n)
+      integer(C_INT) :: err
+    end function openmc_mesh_set_params
+
+    function mesh_id(ptr) result(id) bind(C)
+      import C_PTR, C_INT32_T
+      type(C_PTR), value :: ptr
+      integer(C_INT32_T) :: id
+    end function
+
+    function mesh_volume_frac(ptr) result(volume_frac) bind(C)
+      import C_PTR, C_DOUBLE
+      type(C_PTR), value :: ptr
+      real(C_DOUBLE) :: volume_frac
+    end function
+
+    function mesh_n_dimension(ptr) result(n) bind(C)
+      import C_PTR, C_INT
+      type(C_PTR), value :: ptr
+      integer(C_INT) :: n
+    end function
+
+    function mesh_dimension(ptr, i) result(d) bind(C)
+      import C_PTR, C_INT
+      type(C_PTR), value :: ptr
+      integer(C_INT), value :: i
+      integer(C_INT) :: d
+    end function
+
+    function mesh_lower_left(ptr, i) result(ll) bind(C)
+      import C_PTR, C_INT, C_DOUBLE
+      type(C_PTR), value :: ptr
+      integer(C_INT), value :: i
+      real(C_DOUBLE) :: ll
+    end function
+
+    function mesh_upper_right(ptr, i) result(ur) bind(C)
+      import C_PTR, C_INT, C_DOUBLE
+      type(C_PTR), value :: ptr
+      integer(C_INT), value :: i
+      real(C_DOUBLE) :: ur
+    end function
+
+    function mesh_width(ptr, i) result(w) bind(C)
+      import C_PTR, C_INT, C_DOUBLE
+      type(C_PTR), value :: ptr
+      integer(C_INT), value :: i
+      real(C_DOUBLE) :: w
+    end function
+
+    pure function mesh_get_bin(ptr, xyz) result(bin) bind(C)
+      import C_PTR, C_DOUBLE, C_INT
+      type(C_PTR), value :: ptr
+      real(C_DOUBLE), intent(in) :: xyz(*)
+      integer(C_INT) :: bin
+    end function
+
+    pure function mesh_get_bin_from_indices(ptr, ijk) result(bin) bind(C)
+      import C_PTR, C_INT
+      type(C_PTR), value :: ptr
+      integer(C_INT), intent(in) :: ijk(*)
+      integer(C_INT) :: bin
+    end function
+
+    pure subroutine mesh_get_indices(ptr, xyz, ijk, in_mesh) bind(C)
+      import C_PTR, C_DOUBLE, C_INT, C_BOOL
+      type(C_PTR), value :: ptr
+      real(C_DOUBLE), intent(in) :: xyz(*)
+      integer(C_INT), intent(out) :: ijk(*)
+      logical(C_BOOL), intent(out) :: in_mesh
+    end subroutine
+
+    pure subroutine mesh_get_indices_from_bin(ptr, bin, ijk) bind(C)
+      import C_PTR, C_INT
+      type(C_PTR), value :: ptr
+      integer(C_INT), value :: bin
+      integer(C_INT), intent(out) :: ijk(*)
+    end subroutine
+
+    function mesh_ptr(i) result(ptr) bind(C)
+      import C_INT, C_PTR
+      integer(C_INT), value :: i
+      type(C_PTR) :: ptr
+    end function
+
+    subroutine read_meshes(node_ptr) bind(C)
+      import C_PTR
+      type(C_PTR) :: node_ptr
+    end subroutine
+
+    function n_meshes() result(n) bind(C)
+      import C_INT
+      integer(C_INT) :: n
+    end function
+  end interface
 
 contains
 
-  subroutine regular_from_xml(this, node)
-    class(RegularMesh), intent(inout) :: this
-    type(XMLNode), intent(in) :: node
+  function meshes(i) result(m)
+    integer, intent(in) :: i
+    type(RegularMesh) :: m
 
-    integer :: n
-    character(MAX_LINE_LEN) :: temp_str
+    m % ptr = mesh_ptr(i)
+  end function
 
-    ! Copy mesh id
-    if (check_for_node(node, "id")) then
-      call get_node_value(node, "id", this % id)
+  function regular_id(this) result(id)
+    class(RegularMesh), intent(in) :: this
+    integer(C_INT32_T) :: id
+    id = mesh_id(this % ptr)
+  end function
 
-      ! Check to make sure 'id' hasn't been used
-      if (mesh_dict % has(this % id)) then
-        call fatal_error("Two or more meshes use the same unique ID: " &
-             // to_str(this % id))
-      end if
-    end if
+  function regular_volume_frac(this) result(volume_frac)
+    class(RegularMesh), intent(in) :: this
+    real(C_DOUBLE) :: volume_frac
+    volume_frac = mesh_volume_frac(this % ptr)
+  end function
 
-    ! Read mesh type
-    if (check_for_node(node, "type")) then
-      call get_node_value(node, "type", temp_str)
-      select case (to_lower(temp_str))
-      case ('rect', 'rectangle', 'rectangular')
-        call warning("Mesh type '" // trim(temp_str) // "' is deprecated. &
-             &Please use 'regular' instead.")
-        this % type = MESH_REGULAR
-      case ('regular')
-        this % type = MESH_REGULAR
-      case default
-        call fatal_error("Invalid mesh type: " // trim(temp_str))
-      end select
-    else
-      this % type = MESH_REGULAR
-    end if
+  function regular_n_dimension(this) result(n)
+    class(RegularMesh), intent(in) :: this
+    integer(C_INT) :: n
+    n = mesh_n_dimension(this % ptr)
+  end function
 
-    ! Determine number of dimensions for mesh
-    if (check_for_node(node, "dimension")) then
-      n = node_word_count(node, "dimension")
-      if (n /= 1 .and. n /= 2 .and. n /= 3) then
-        call fatal_error("Mesh must be one, two, or three dimensions.")
-      end if
-      this % n_dimension = n
+  function regular_dimension(this, i) result(d)
+    class(RegularMesh), intent(in) :: this
+    integer(C_INT), intent(in) :: i
+    integer(C_INT) :: d
+    d = mesh_dimension(this % ptr, i)
+  end function
 
-      ! Allocate attribute arrays
-      allocate(this % dimension(n))
+  function regular_lower_left(this, i) result(ll)
+    class(RegularMesh), intent(in) :: this
+    integer(C_INT), intent(in) :: i
+    real(C_DOUBLE) :: ll
+    ll = mesh_lower_left(this % ptr, i)
+  end function
 
-      ! Check that dimensions are all greater than zero
-      call get_node_array(node, "dimension", this % dimension)
-      if (any(this % dimension <= 0)) then
-        call fatal_error("All entries on the <dimension> element for a tally &
-             &mesh must be positive.")
-      end if
-    end if
+  function regular_upper_right(this, i) result(ur)
+    class(RegularMesh), intent(in) :: this
+    integer(C_INT), intent(in) :: i
+    real(C_DOUBLE) :: ur
+    ur = mesh_upper_right(this % ptr, i)
+  end function
 
-    ! Check for lower-left coordinates
-    if (check_for_node(node, "lower_left")) then
-      n = node_word_count(node, "lower_left")
-      allocate(this % lower_left(n))
-
-      ! Read mesh lower-left corner location
-      call get_node_array(node, "lower_left", this % lower_left)
-    else
-      call fatal_error("Must specify <lower_left> on a mesh.")
-    end if
-
-    if (check_for_node(node, "width")) then
-      ! Make sure both upper-right or width were specified
-      if (check_for_node(node, "upper_right")) then
-        call fatal_error("Cannot specify both <upper_right> and <width> on a &
-             &mesh.")
-      end if
-
-      n = node_word_count(node, "width")
-      allocate(this % width(n))
-      allocate(this % upper_right(n))
-
-      ! Check to ensure width has same dimensions
-      if (n /= size(this % lower_left)) then
-        call fatal_error("Number of entries on <width> must be the same as &
-             &the number of entries on <lower_left>.")
-      end if
-
-      ! Check for negative widths
-      call get_node_array(node, "width", this % width)
-      if (any(this % width < ZERO)) then
-        call fatal_error("Cannot have a negative <width> on a tally mesh.")
-      end if
-
-      ! Set width and upper right coordinate
-      this % upper_right = this % lower_left + this % dimension * this % width
-
-    elseif (check_for_node(node, "upper_right")) then
-      n = node_word_count(node, "upper_right")
-      allocate(this % upper_right(n))
-      allocate(this % width(n))
-
-      ! Check to ensure width has same dimensions
-      if (n /= size(this % lower_left)) then
-        call fatal_error("Number of entries on <upper_right> must be the &
-             &same as the number of entries on <lower_left>.")
-      end if
-
-      ! Check that upper-right is above lower-left
-      call get_node_array(node, "upper_right", this % upper_right)
-      if (any(this % upper_right < this % lower_left)) then
-        call fatal_error("The <upper_right> coordinates must be greater than &
-             &the <lower_left> coordinates on a tally mesh.")
-      end if
-
-      ! Set width and upper right coordinate
-      this % width = (this % upper_right - this % lower_left) / this % dimension
-    else
-      call fatal_error("Must specify either <upper_right> and <width> on a &
-           &mesh.")
-    end if
-
-    if (allocated(this % dimension)) then
-      if (size(this % dimension) /= size(this % lower_left)) then
-        call fatal_error("Number of entries on <lower_left> must be the same &
-             &as the number of entries on <dimension>.")
-      end if
-
-      ! Set volume fraction
-      this % volume_frac = ONE/real(product(this % dimension),8)
-    end if
-
-  end subroutine regular_from_xml
+  function regular_width(this, i) result(w)
+    class(RegularMesh), intent(in) :: this
+    integer(C_INT), intent(in) :: i
+    real(C_DOUBLE) :: w
+    w = mesh_width(this % ptr, i)
+  end function
 
 !===============================================================================
 ! GET_MESH_BIN determines the tally bin for a particle in a structured mesh
@@ -191,38 +248,8 @@ contains
     real(8), intent(in)           :: xyz(:) ! coordinates
     integer, intent(out)          :: bin    ! tally bin
 
-    integer :: n       ! size of mesh
-    integer :: d       ! mesh dimension index
-    integer :: ijk(3)  ! indices in mesh
-    logical :: in_mesh ! was given coordinate in mesh at all?
-
-    ! Get number of dimensions
-    n = this % n_dimension
-
-    ! Loop over the dimensions of the mesh
-    do d = 1, n
-
-      ! Check for cases where particle is outside of mesh
-      if (xyz(d) < this % lower_left(d)) then
-        bin = NO_BIN_FOUND
-        return
-      elseif (xyz(d) > this % upper_right(d)) then
-        bin = NO_BIN_FOUND
-        return
-      end if
-    end do
-
-    ! Determine indices
-    call this % get_indices(xyz, ijk, in_mesh)
-
-    ! Convert indices to bin
-    if (in_mesh) then
-      bin = this % get_bin_from_indices(ijk)
-    else
-      bin = NO_BIN_FOUND
-    end if
-
-  end subroutine regular_get_bin
+    bin = mesh_get_bin(this % ptr, xyz)
+  end subroutine
 
 !===============================================================================
 ! GET_MESH_INDICES determines the indices of a particle in a structured mesh
@@ -234,18 +261,10 @@ contains
     integer, intent(out)          :: ijk(:)  ! indices in mesh
     logical, intent(out)          :: in_mesh ! were given coords in mesh?
 
-    ! Find particle in mesh
-    ijk(:this % n_dimension) = ceiling((xyz(:this % n_dimension) - &
-         this % lower_left)/this % width)
-
-    ! Determine if particle is in mesh
-    if (any(ijk(:this % n_dimension) < 1) .or. &
-         any(ijk(:this % n_dimension) > this % dimension)) then
-      in_mesh = .false.
-    else
-      in_mesh = .true.
-    end if
+    logical(C_BOOL) :: in_mesh_
 
+    call mesh_get_indices(this % ptr, xyz, ijk, in_mesh_)
+    in_mesh = in_mesh_
   end subroutine regular_get_indices
 
 !===============================================================================
@@ -258,15 +277,7 @@ contains
     integer, intent(in)           :: ijk(:)
     integer                       :: bin
 
-    if (this % n_dimension == 1) then
-      bin = ijk(1)
-    elseif (this % n_dimension == 2) then
-      bin = (ijk(2) - 1) * this % dimension(1) + ijk(1)
-    elseif (this % n_dimension == 3) then
-      bin = ((ijk(3) - 1) * this % dimension(2) + (ijk(2) - 1)) &
-           * this % dimension(1) + ijk(1)
-    end if
-
+    bin = mesh_get_bin_from_indices(this % ptr, ijk)
   end function regular_get_bin_from_indices
 
 !===============================================================================
@@ -279,498 +290,7 @@ contains
     integer, intent(in)           :: bin
     integer, intent(out)          :: ijk(:)
 
-    if (this % n_dimension == 1) then
-      ijk(1) = bin
-    else if (this % n_dimension == 2) then
-      ijk(1) = mod(bin - 1, this % dimension(1)) + 1
-      ijk(2) = (bin - 1)/this % dimension(1) + 1
-    else if (this % n_dimension == 3) then
-      ijk(1) = mod(bin - 1, this % dimension(1)) + 1
-      ijk(2) = mod(bin - 1, this % dimension(1) * this % dimension(2)) &
-           / this % dimension(1) + 1
-      ijk(3) = (bin - 1)/(this % dimension(1) * this % dimension(2)) + 1
-    end if
-
+    call mesh_get_indices_from_bin(this % ptr, bin, ijk)
   end subroutine regular_get_indices_from_bin
 
-!===============================================================================
-! MESH_INTERSECTS determines if a line between xyz0 and xyz1 intersects the
-! outer boundary of the given mesh. This is important for determining whether a
-! track will score to a mesh tally.
-!===============================================================================
-
-  pure function regular_intersects(this, xyz0, xyz1) result(intersects)
-    class(RegularMesh), intent(in) :: this
-    real(8), intent(in) :: xyz0(:)
-    real(8), intent(in) :: xyz1(:)
-    logical :: intersects
-
-    select case(this % n_dimension)
-    case (1)
-      intersects = mesh_intersects_1d(this, xyz0, xyz1)
-    case (2)
-      intersects = mesh_intersects_2d(this, xyz0, xyz1)
-    case (3)
-      intersects = mesh_intersects_3d(this, xyz0, xyz1)
-    end select
-  end function regular_intersects
-
-  pure function mesh_intersects_1d(m, xyz0, xyz1) result(intersects)
-    type(RegularMesh), intent(in) :: m
-    real(8), intent(in) :: xyz0(:)
-    real(8), intent(in) :: xyz1(:)
-    logical :: intersects
-
-    real(8) :: x0   ! track start point
-    real(8) :: x1   ! track end point
-    real(8) :: xm0 ! lower-left coordinates of mesh
-    real(8) :: xm1 ! upper-right coordinates of mesh
-
-    ! Copy coordinates of starting point
-    x0 = xyz0(1)
-
-    ! Copy coordinates of ending point
-    x1 = xyz1(1)
-
-    ! Copy coordinates of mesh lower_left
-    xm0 = m % lower_left(1)
-
-    ! Copy coordinates of mesh upper_right
-    xm1 = m % upper_right(1)
-
-    ! Set default value for intersects
-    intersects = .false.
-
-    ! Check if line intersects left surface
-    if ((x0 < xm0 .and. x1 > xm0) .or. (x0 > xm0 .and. x1 < xm0)) then
-      intersects = .true.
-      return
-    end if
-
-    ! Check if line intersects right surface
-    if ((x0 < xm1 .and. x1 > xm1) .or. (x0 > xm1 .and. x1 < xm1)) then
-      intersects = .true.
-      return
-    end if
-
-  end function mesh_intersects_1d
-
-  pure function mesh_intersects_2d(m, xyz0, xyz1) result(intersects)
-    type(RegularMesh), intent(in) :: m
-    real(8), intent(in) :: xyz0(:)
-    real(8), intent(in) :: xyz1(:)
-    logical :: intersects
-
-    real(8) :: x0, y0   ! track start point
-    real(8) :: x1, y1   ! track end point
-    real(8) :: xi, yi   ! track intersection point with mesh
-    real(8) :: xm0, ym0 ! lower-left coordinates of mesh
-    real(8) :: xm1, ym1 ! upper-right coordinates of mesh
-
-    ! Copy coordinates of starting point
-    x0 = xyz0(1)
-    y0 = xyz0(2)
-
-    ! Copy coordinates of ending point
-    x1 = xyz1(1)
-    y1 = xyz1(2)
-
-    ! Copy coordinates of mesh lower_left
-    xm0 = m % lower_left(1)
-    ym0 = m % lower_left(2)
-
-    ! Copy coordinates of mesh upper_right
-    xm1 = m % upper_right(1)
-    ym1 = m % upper_right(2)
-
-    ! Set default value for intersects
-    intersects = .false.
-
-    ! Check if line intersects left surface -- calculate the intersection point
-    ! y
-    if ((x0 < xm0 .and. x1 > xm0) .or. (x0 > xm0 .and. x1 < xm0)) then
-      yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0)
-      if (yi >= ym0 .and. yi < ym1) then
-        intersects = .true.
-        return
-      end if
-    end if
-
-    ! Check if line intersects back surface -- calculate the intersection point
-    ! x
-    if ((y0 < ym0 .and. y1 > ym0) .or. (y0 > ym0 .and. y1 < ym0)) then
-      xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0)
-      if (xi >= xm0 .and. xi < xm1) then
-        intersects = .true.
-        return
-      end if
-    end if
-
-    ! Check if line intersects right surface -- calculate the intersection
-    ! point y
-    if ((x0 < xm1 .and. x1 > xm1) .or. (x0 > xm1 .and. x1 < xm1)) then
-      yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0)
-      if (yi >= ym0 .and. yi < ym1) then
-        intersects = .true.
-        return
-      end if
-    end if
-
-    ! Check if line intersects front surface -- calculate the intersection point
-    ! x
-    if ((y0 < ym1 .and. y1 > ym1) .or. (y0 > ym1 .and. y1 < ym1)) then
-      xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0)
-      if (xi >= xm0 .and. xi < xm1) then
-        intersects = .true.
-        return
-      end if
-    end if
-
-  end function mesh_intersects_2d
-
-  pure function mesh_intersects_3d(m, xyz0, xyz1) result(intersects)
-    type(RegularMesh), intent(in) :: m
-    real(8), intent(in) :: xyz0(:)
-    real(8), intent(in) :: xyz1(:)
-    logical :: intersects
-
-    real(8) :: x0, y0, z0    ! track start point
-    real(8) :: x1, y1, z1    ! track end point
-    real(8) :: xi, yi, zi    ! track intersection point with mesh
-    real(8) :: xm0, ym0, zm0 ! lower-left coordinates of mesh
-    real(8) :: xm1, ym1, zm1 ! upper-right coordinates of mesh
-
-    ! Copy coordinates of starting point
-    x0 = xyz0(1)
-    y0 = xyz0(2)
-    z0 = xyz0(3)
-
-    ! Copy coordinates of ending point
-    x1 = xyz1(1)
-    y1 = xyz1(2)
-    z1 = xyz1(3)
-
-    ! Copy coordinates of mesh lower_left
-    xm0 = m % lower_left(1)
-    ym0 = m % lower_left(2)
-    zm0 = m % lower_left(3)
-
-    ! Copy coordinates of mesh upper_right
-    xm1 = m % upper_right(1)
-    ym1 = m % upper_right(2)
-    zm1 = m % upper_right(3)
-
-    ! Set default value for intersects
-    intersects = .false.
-
-    ! Check if line intersects left surface -- calculate the intersection point
-    ! (y,z)
-    if ((x0 < xm0 .and. x1 > xm0) .or. (x0 > xm0 .and. x1 < xm0)) then
-      yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0)
-      zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0)
-      if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then
-        intersects = .true.
-        return
-      end if
-    end if
-
-    ! Check if line intersects back surface -- calculate the intersection point
-    ! (x,z)
-    if ((y0 < ym0 .and. y1 > ym0) .or. (y0 > ym0 .and. y1 < ym0)) then
-      xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0)
-      zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0)
-      if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then
-        intersects = .true.
-        return
-      end if
-    end if
-
-    ! Check if line intersects bottom surface -- calculate the intersection
-    ! point (x,y)
-    if ((z0 < zm0 .and. z1 > zm0) .or. (z0 > zm0 .and. z1 < zm0)) then
-      xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0)
-      yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0)
-      if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then
-        intersects = .true.
-        return
-      end if
-    end if
-
-    ! Check if line intersects right surface -- calculate the intersection point
-    ! (y,z)
-    if ((x0 < xm1 .and. x1 > xm1) .or. (x0 > xm1 .and. x1 < xm1)) then
-      yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0)
-      zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0)
-      if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then
-        intersects = .true.
-        return
-      end if
-    end if
-
-    ! Check if line intersects front surface -- calculate the intersection point
-    ! (x,z)
-    if ((y0 < ym1 .and. y1 > ym1) .or. (y0 > ym1 .and. y1 < ym1)) then
-      xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0)
-      zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0)
-      if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then
-        intersects = .true.
-        return
-      end if
-    end if
-
-    ! Check if line intersects top surface -- calculate the intersection point
-    ! (x,y)
-    if ((z0 < zm1 .and. z1 > zm1) .or. (z0 > zm1 .and. z1 < zm1)) then
-      xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0)
-      yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0)
-      if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then
-        intersects = .true.
-        return
-      end if
-    end if
-
-  end function mesh_intersects_3d
-
-!===============================================================================
-! TO_HDF5 writes the mesh data to an HDF5 group
-!===============================================================================
-
-  subroutine regular_to_hdf5(this, group)
-    class(RegularMesh), intent(in) :: this
-    integer(HID_T), intent(in) :: group
-
-    integer(HID_T) :: mesh_group
-
-    mesh_group = create_group(group, "mesh " // trim(to_str(this % id)))
-
-    call write_dataset(mesh_group, "type", "regular")
-    call write_dataset(mesh_group, "dimension", this % dimension)
-    call write_dataset(mesh_group, "lower_left", this % lower_left)
-    call write_dataset(mesh_group, "upper_right", this % upper_right)
-    call write_dataset(mesh_group, "width", this % width)
-
-    call close_group(mesh_group)
-  end subroutine regular_to_hdf5
-
-!===============================================================================
-! FREE_MEMORY_MESH deallocates global arrays defined in this module
-!===============================================================================
-
-  subroutine free_memory_mesh()
-    n_meshes = 0
-    if (allocated(meshes)) deallocate(meshes)
-    call mesh_dict % clear()
-  end subroutine free_memory_mesh
-
-!===============================================================================
-!                               C API FUNCTIONS
-!===============================================================================
-
-  function openmc_extend_meshes(n, index_start, index_end) result(err) bind(C)
-    ! Extend the meshes array by n elements
-    integer(C_INT32_T), value, intent(in) :: n
-    integer(C_INT32_T), optional, intent(out) :: index_start
-    integer(C_INT32_T), optional, intent(out) :: index_end
-    integer(C_INT) :: err
-
-    type(RegularMesh), allocatable :: temp(:) ! temporary meshes array
-
-    if (n_meshes == 0) then
-      ! Allocate meshes array
-      allocate(meshes(n))
-    else
-      ! Allocate meshes array with increased size
-      allocate(temp(n_meshes + n))
-
-      ! Copy original meshes to temporary array
-      temp(1:n_meshes) = meshes
-
-      ! Move allocation from temporary array
-      call move_alloc(FROM=temp, TO=meshes)
-    end if
-
-    ! Return indices in meshes array
-    if (present(index_start)) index_start = n_meshes + 1
-    if (present(index_end)) index_end = n_meshes + n
-    n_meshes = n_meshes + n
-
-    err = 0
-  end function openmc_extend_meshes
-
-
-  function openmc_get_mesh_index(id, index) result(err) bind(C)
-    ! Return the index in the meshes array of a mesh with a given ID
-    integer(C_INT32_T), value :: id
-    integer(C_INT32_T), intent(out) :: index
-    integer(C_INT) :: err
-
-    if (allocated(meshes)) then
-      if (mesh_dict % has(id)) then
-        index = mesh_dict % get(id)
-        err = 0
-      else
-        err = E_INVALID_ID
-        call set_errmsg("No mesh exists with ID=" // trim(to_str(id)) // ".")
-      end if
-    else
-      err = E_ALLOCATE
-      call set_errmsg("Memory has not been allocated for meshes.")
-    end if
-  end function openmc_get_mesh_index
-
-
-  function openmc_mesh_get_id(index, id) result(err) bind(C)
-    ! Return the ID of a mesh
-    integer(C_INT32_T), value       :: index
-    integer(C_INT32_T), intent(out) :: id
-    integer(C_INT) :: err
-
-    if (index >= 1 .and. index <= size(meshes)) then
-      id = meshes(index) % id
-      err = 0
-    else
-      err = E_OUT_OF_BOUNDS
-      call set_errmsg("Index in meshes array is out of bounds.")
-    end if
-  end function openmc_mesh_get_id
-
-
-  function openmc_mesh_set_id(index, id) result(err) bind(C)
-    ! Set the ID of a mesh
-    integer(C_INT32_T), value, intent(in) :: index
-    integer(C_INT32_T), value, intent(in) :: id
-    integer(C_INT) :: err
-
-    if (index >= 1 .and. index <= n_meshes) then
-      meshes(index) % id = id
-      call mesh_dict % set(id, index)
-      err = 0
-    else
-      err = E_OUT_OF_BOUNDS
-      call set_errmsg("Index in meshes array is out of bounds.")
-    end if
-  end function openmc_mesh_set_id
-
-
-  function openmc_mesh_get_dimension(index, dims, n) result(err) bind(C)
-    ! Get the dimension of a mesh
-    integer(C_INT32_T), value, intent(in) :: index
-    type(C_PTR),               intent(out) :: dims
-    integer(C_INT),            intent(out) :: n
-    integer(C_INT) :: err
-
-    if (index >= 1 .and. index <= n_meshes) then
-      dims = C_LOC(meshes(index) % dimension)
-      n = meshes(index) % n_dimension
-      err = 0
-    else
-      err = E_OUT_OF_BOUNDS
-      call set_errmsg("Index in meshes array is out of bounds.")
-    end if
-  end function openmc_mesh_get_dimension
-
-
-  function openmc_mesh_set_dimension(index, n, dims) result(err) bind(C)
-    ! Set the dimension of a mesh
-    integer(C_INT32_T), value, intent(in) :: index
-    integer(C_INT),     value, intent(in) :: n
-    integer(C_INT),            intent(in) :: dims(n)
-    integer(C_INT) :: err
-
-    if (index >= 1 .and. index <= n_meshes) then
-      associate (m => meshes(index))
-        if (allocated(m % dimension)) deallocate (m % dimension)
-        if (allocated(m % lower_left)) deallocate (m % lower_left)
-        if (allocated(m % upper_right)) deallocate (m % upper_right)
-        if (allocated(m % width)) deallocate (m % width)
-
-        m % n_dimension = n
-        allocate(m % dimension(n))
-        allocate(m % lower_left(n))
-        allocate(m % upper_right(n))
-        allocate(m % width(n))
-
-        ! Copy dimension
-        m % dimension(:) = dims
-      end associate
-      err = 0
-    else
-      err = E_OUT_OF_BOUNDS
-      call set_errmsg("Index in meshes array is out of bounds.")
-    end if
-  end function openmc_mesh_set_dimension
-
-
-  function openmc_mesh_get_params(index, ll, ur, width, n) result(err) bind(C)
-    ! Get the mesh parameters
-    integer(C_INT32_T), value, intent(in) :: index
-    type(C_PTR), intent(out) :: ll
-    type(C_PTR), intent(out) :: ur
-    type(C_PTR), intent(out) :: width
-    integer(C_INT), intent(out) :: n
-    integer(C_INT) :: err
-
-    err = 0
-    if (index >= 1 .and. index <= n_meshes) then
-      associate (m => meshes(index))
-        if (allocated(m % lower_left)) then
-          ll = C_LOC(m % lower_left(1))
-          ur = C_LOC(m % upper_right(1))
-          width = C_LOC(m % width(1))
-          n = m % n_dimension
-        else
-          err = E_ALLOCATE
-          call set_errmsg("Mesh parameters have not been set.")
-        end if
-      end associate
-    else
-      err = E_OUT_OF_BOUNDS
-      call set_errmsg("Index in meshes array is out of bounds.")
-    end if
-  end function openmc_mesh_get_params
-
-
-  function openmc_mesh_set_params(index, n, ll, ur, width) result(err) bind(C)
-    ! Set the mesh parameters
-    integer(C_INT32_T), value, intent(in) :: index
-    integer(C_INT),     value, intent(in) :: n
-    real(C_DOUBLE), intent(in), optional :: ll(n)
-    real(C_DOUBLE), intent(in), optional :: ur(n)
-    real(C_DOUBLE), intent(in), optional :: width(n)
-    integer(C_INT) :: err
-
-    err = 0
-    if (index >= 1 .and. index <= n_meshes) then
-      associate (m => meshes(index))
-        if (allocated(m % lower_left)) deallocate (m % lower_left)
-        if (allocated(m % upper_right)) deallocate (m % upper_right)
-        if (allocated(m % width)) deallocate (m % width)
-
-        allocate(m % lower_left(n))
-        allocate(m % upper_right(n))
-        allocate(m % width(n))
-
-        if (present(ll) .and. present(ur)) then
-          m % lower_left(:) = ll
-          m % upper_right(:) = ur
-          m % width(:) = (ur - ll) / m % dimension
-        elseif (present(ll) .and. present(width)) then
-          m % lower_left(:) = ll
-          m % width(:) = width
-          m % upper_right(:) = ll + width * m % dimension
-        elseif (present(ur) .and. present(width)) then
-          m % upper_right(:) = ur
-          m % width(:) = width
-          m % lower_left(:) = ur - width * m % dimension
-        else
-          err = E_INVALID_ARGUMENT
-          call set_errmsg("At least two parameters must be specified.")
-        end if
-      end associate
-    else
-      err = E_OUT_OF_BOUNDS
-      call set_errmsg("Index in meshes array is out of bounds.")
-    end if
-  end function openmc_mesh_set_params
-
 end module mesh_header

src/message_passing.cpp

@@ -5,6 +5,7 @@ namespace mpi {
 
 int rank {0};
 int n_procs {1};
+bool master {true};
 
 #ifdef OPENMC_MPI
 MPI_Comm intracomm;

src/multipole_header.F90

@@ -1,7 +1,6 @@
 module multipole_header
 
   use constants
-  use dict_header,      only: DictIntInt
   use error,            only: fatal_error
   use hdf5_interface
 
@@ -74,12 +73,11 @@ contains
     character(len=*),      intent(in)    :: filename
 
     character(len=10) :: version
-    integer :: i, n_poles, n_residues, n_windows
+    integer :: n_poles, n_residues, n_windows
     integer(HSIZE_T) :: dims_1d(1), dims_2d(2), dims_3d(3)
     integer(HID_T) :: file_id
     integer(HID_T) :: group_id
     integer(HID_T) :: dset
-    type(DictIntInt) :: l_val_dict
 
     ! Open file for reading and move into the /isotope group
     file_id = file_open(filename, 'r', parallel=.true.)

src/nuclide_header.F90

@@ -600,7 +600,6 @@ contains
 
     integer :: i, j, k, l
     integer :: t
-    integer :: m
     integer :: n
     integer :: n_grid
     integer :: i_fission

src/output.F90

@@ -10,7 +10,6 @@ module output
   use error,           only: fatal_error, warning
   use geometry_header
   use math,            only: t_percentile
-  use mesh_header,     only: RegularMesh, meshes
   use message_passing, only: master, n_procs
   use mgxs_interface
   use nuclide_header
@@ -34,6 +33,14 @@ module output
   integer :: ou = OUTPUT_UNIT
   integer :: eu = ERROR_UNIT
 
+  interface
+    function entropy(i) result(h) bind(C, name='entropy_c')
+      import C_INT, C_DOUBLE
+      integer(C_INT), value :: i
+      real(C_DOUBLE) :: h
+    end function
+  end interface
+
 contains
 
 !===============================================================================
@@ -336,7 +343,7 @@ contains
 
     ! write out entropy info
     if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
-         entropy % data(i)
+         entropy(i)
 
     if (n > 1) then
       write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') &
@@ -370,7 +377,7 @@ contains
 
     ! write out entropy info
     if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
-         entropy % data(i)
+         entropy(i)
 
     ! write out accumulated k-effective if after first active batch
     if (n > 1) then
@@ -653,6 +660,10 @@ contains
     if (n_tallies == 0) return
 
     allocate(matches(n_filters))
+    do i = 1, n_filters
+      allocate(matches(i) % bins)
+      allocate(matches(i) % weights)
+    end do
 
     ! Initialize names for scores
     score_names(abs(SCORE_FLUX))               = "Flux"
@@ -848,6 +859,11 @@ contains
 
     close(UNIT=unit_tally)
 
+    do i = 1, n_filters
+      deallocate(matches(i) % bins)
+      deallocate(matches(i) % weights)
+    end do
+
   end subroutine write_tallies
 
 !===============================================================================

src/particle.cpp

@@ -141,7 +141,7 @@ Particle::mark_as_lost(const char* message)
 }
 
 void
-Particle::write_restart()
+Particle::write_restart() const
 {
   // Dont write another restart file if in particle restart mode
   if (settings::run_mode == RUN_MODE_PARTICLE) return;

src/physics.F90

@@ -6,7 +6,6 @@ module physics
   use error,                  only: fatal_error, warning, write_message
   use material_header,        only: Material, materials
   use math
-  use mesh_header,            only: meshes
   use message_passing
   use nuclide_header
   use particle_header
@@ -1182,11 +1181,18 @@ contains
     integer :: nu_d(MAX_DELAYED_GROUPS) ! number of delayed neutrons born
     integer :: i                        ! loop index
     integer :: nu                       ! actual number of neutrons produced
-    integer :: mesh_bin                 ! mesh bin for source site
     real(8) :: nu_t                     ! total nu
     real(8) :: weight                   ! weight adjustment for ufs method
     type(Nuclide),  pointer :: nuc
 
+    interface
+      function ufs_get_weight(p) result(weight) bind(C)
+        import Particle, C_DOUBLE
+        type(Particle), intent(in) :: p
+        real(C_DOUBLE) :: WEIGHT
+      end function
+    end interface
+
     ! Get pointers
     nuc => nuclides(i_nuclide)
 
@@ -1196,20 +1202,7 @@ contains
     ! the expected number of fission sites produced
 
     if (ufs) then
-      associate (m => meshes(index_ufs_mesh))
-        ! Determine indices on ufs mesh for current location
-        call m % get_bin(p % coord(1) % xyz, mesh_bin)
-        if (mesh_bin == NO_BIN_FOUND) then
-          call particle_write_restart(p)
-          call fatal_error("Source site outside UFS mesh!")
-        end if
-
-        if (source_frac(1, mesh_bin) /= ZERO) then
-          weight = m % volume_frac / source_frac(1, mesh_bin)
-        else
-          weight = ONE
-        end if
-      end associate
+      weight = ufs_get_weight(p)
     else
       weight = ONE
     end if

src/physics_mg.F90

@@ -7,7 +7,6 @@ module physics_mg
   use error,                  only: fatal_error, warning, write_message
   use material_header,        only: Material, materials
   use math,                   only: rotate_angle
-  use mesh_header,            only: meshes
   use mgxs_interface
   use message_passing
   use nuclide_header,         only: material_xs
@@ -168,33 +167,26 @@ contains
     integer :: dg                       ! delayed group
     integer :: gout                     ! group out
     integer :: nu                       ! actual number of neutrons produced
-    integer :: mesh_bin                 ! mesh bin for source site
     real(8) :: nu_t                     ! total nu
     real(8) :: mu                       ! fission neutron angular cosine
     real(8) :: phi                      ! fission neutron azimuthal angle
     real(8) :: weight                   ! weight adjustment for ufs method
 
+    interface
+      function ufs_get_weight(p) result(weight) bind(C)
+        import Particle, C_DOUBLE
+        type(Particle), intent(in) :: p
+        real(C_DOUBLE) :: WEIGHT
+      end function
+    end interface
+
     ! TODO: Heat generation from fission
 
     ! If uniform fission source weighting is turned on, we increase of decrease
     ! the expected number of fission sites produced
 
     if (ufs) then
-      associate (m => meshes(index_ufs_mesh))
-        ! Determine indices on ufs mesh for current location
-        call m % get_bin(p % coord(1) % xyz, mesh_bin)
-
-        if (mesh_bin == NO_BIN_FOUND) then
-          call particle_write_restart(p)
-          call fatal_error("Source site outside UFS mesh!")
-        end if
-
-        if (source_frac(1, mesh_bin) /= ZERO) then
-          weight = m % volume_frac / source_frac(1, mesh_bin)
-        else
-          weight = ONE
-        end if
-      end associate
+      weight = ufs_get_weight(p)
     else
       weight = ONE
     end if

src/plot.F90

@@ -9,6 +9,7 @@ module plot
   use hdf5_interface
   use output,          only: time_stamp
   use material_header, only: materials
+  use mesh_header,     only: meshes, RegularMesh
   use particle_header
   use plot_header
   use progress_header, only: ProgressBar
@@ -184,7 +185,7 @@ contains
 !$omp end parallel do
 
     ! Draw tally mesh boundaries on the image if requested
-    if (associated(pl % meshlines_mesh)) call draw_mesh_lines(pl, data)
+    if (pl % index_meshlines_mesh >= 0) call draw_mesh_lines(pl, data)
 
     ! Write out the ppm to a file
     call output_ppm(pl, data)
@@ -214,6 +215,7 @@ contains
     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(RegularMesh) :: m
 
     rgb(:) = pl % meshlines_color % rgb
 
@@ -239,57 +241,56 @@ contains
 
     width = xyz_ur_plot - xyz_ll_plot
 
-    associate (m => pl % meshlines_mesh)
-      call m % get_indices(xyz_ll_plot, ijk_ll(:m % n_dimension), in_mesh)
-      call m % get_indices(xyz_ur_plot, ijk_ur(:m % n_dimension), in_mesh)
+    m = meshes(pl % index_meshlines_mesh)
+    call m % get_indices(xyz_ll_plot, ijk_ll(:m % n_dimension()), in_mesh)
+    call m % get_indices(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
+    ! 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
+          ! 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
+          ! map the xyz ranges to pixel ranges
 
-            frac = (xyz_ll(outer) - xyz_ll_plot(outer)) / width(outer)
-            outrange(1) = int(frac * real(pl % pixels(1), 8))
-            frac = (xyz_ur(outer) - xyz_ll_plot(outer)) / width(outer)
-            outrange(2) = int(frac * real(pl % pixels(1), 8))
+          frac = (xyz_ll(outer) - xyz_ll_plot(outer)) / width(outer)
+          outrange(1) = int(frac * real(pl % pixels(1), 8))
+          frac = (xyz_ur(outer) - xyz_ll_plot(outer)) / width(outer)
+          outrange(2) = int(frac * real(pl % pixels(1), 8))
 
-            frac = (xyz_ur(inner) - xyz_ll_plot(inner)) / width(inner)
-            inrange(1) = int((ONE - frac) * real(pl % pixels(2), 8))
-            frac = (xyz_ll(inner) - xyz_ll_plot(inner)) / width(inner)
-            inrange(2) = int((ONE - frac) * real(pl % pixels(2), 8))
+          frac = (xyz_ur(inner) - xyz_ll_plot(inner)) / width(inner)
+          inrange(1) = int((ONE - frac) * real(pl % pixels(2), 8))
+          frac = (xyz_ll(inner) - xyz_ll_plot(inner)) / width(inner)
+          inrange(2) = int((ONE - frac) * real(pl % pixels(2), 8))
 
-            ! draw lines
-            do out_ = outrange(1), outrange(2)
-              do plus = 0, pl % meshlines_width
-                data(:, out_ + 1, inrange(1) + plus + 1) = rgb
-                data(:, out_ + 1, inrange(2) + plus + 1) = rgb
-                data(:, out_ + 1, inrange(1) - plus + 1) = rgb
-                data(:, out_ + 1, inrange(2) - plus + 1) = rgb
-              end do
+          ! draw lines
+          do out_ = outrange(1), outrange(2)
+            do plus = 0, pl % meshlines_width
+              data(:, out_ + 1, inrange(1) + plus + 1) = rgb
+              data(:, out_ + 1, inrange(2) + plus + 1) = rgb
+              data(:, out_ + 1, inrange(1) - plus + 1) = rgb
+              data(:, out_ + 1, inrange(2) - plus + 1) = rgb
             end do
-            do in_ = inrange(1), inrange(2)
-              do plus = 0, pl % meshlines_width
-                data(:, outrange(1) + plus + 1, in_ + 1) = rgb
-                data(:, outrange(2) + plus + 1, in_ + 1) = rgb
-                data(:, outrange(1) - plus + 1, in_ + 1) = rgb
-                data(:, outrange(2) - plus + 1, in_ + 1) = rgb
-              end do
+          end do
+          do in_ = inrange(1), inrange(2)
+            do plus = 0, pl % meshlines_width
+              data(:, outrange(1) + plus + 1, in_ + 1) = rgb
+              data(:, outrange(2) + plus + 1, in_ + 1) = rgb
+              data(:, outrange(1) - plus + 1, in_ + 1) = rgb
+              data(:, outrange(2) - plus + 1, in_ + 1) = rgb
             end do
+          end do
 
-          end if
-        end do
+        end if
       end do
-    end associate
+    end do
 
   end subroutine draw_mesh_lines
 

src/plot_header.F90

@@ -4,7 +4,6 @@ module plot_header
 
   use constants
   use dict_header, only: DictIntInt
-  use mesh_header, only: RegularMesh
 
   implicit none
 
@@ -31,7 +30,7 @@ module plot_header
     integer :: pixels(3)             ! pixel width/height of plot slice
     integer :: meshlines_width       ! pixel width of meshlines
     integer :: level                 ! universe depth to plot the cells of
-    type(RegularMesh), pointer :: meshlines_mesh => null() ! mesh to plot
+    integer :: index_meshlines_mesh = -1 ! index of  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

src/settings.cpp

@@ -1,5 +1,6 @@
 #include "openmc/settings.h"
 
+#include <cmath> // for ceil, pow
 #include <limits> // for numeric_limits
 #include <sstream>
 #include <string>
@@ -13,6 +14,8 @@
 #include "openmc/distribution_spatial.h"
 #include "openmc/error.h"
 #include "openmc/file_utils.h"
+#include "openmc/mesh.h"
+#include "openmc/output.h"
 #include "openmc/random_lcg.h"
 #include "openmc/source.h"
 #include "openmc/string_utils.h"
@@ -167,11 +170,7 @@ void get_run_parameters(pugi::xml_node node_base)
   }
 }
 
-extern "C" void title();
-extern "C" void read_settings_xml_f(pugi::xml_node_struct* root_ptr);
-
-extern "C" void
-read_settings_xml()
+void read_settings_xml()
 {
   using namespace settings;
   using namespace pugi;
@@ -484,7 +483,86 @@ read_settings_xml()
     //track_identifiers = reshape(temp_int_array, [3, n_tracks/3])
   }
 
-  // TODO: Read meshes
+  // Read meshes
+  read_meshes(&root);
+
+  // Shannon Entropy mesh
+  if (check_for_node(root, "entropy_mesh")) {
+    int temp = std::stoi(get_node_value(root, "entropy_mesh"));
+    if (mesh_map.find(temp) == mesh_map.end()) {
+      std::stringstream msg;
+      msg << "Mesh " << temp << " specified for Shannon entropy does not exist.";
+      fatal_error(msg);
+    }
+    index_entropy_mesh = mesh_map.at(temp);
+
+  } else if (check_for_node(root, "entropy")) {
+    warning("Specifying a Shannon entropy mesh via the <entropy> element "
+      "is deprecated. Please create a mesh using <mesh> and then reference "
+      "it by specifying its ID in an <entropy_mesh> element.");
+
+    // Read entropy mesh from <entropy>
+    auto node_entropy = root.child("entropy");
+    meshes.emplace_back(new RegularMesh{node_entropy});
+
+    // Set entropy mesh index
+    index_entropy_mesh = meshes.size() - 1;
+
+    // Assign ID and set mapping
+    meshes.back()->id_ = 10000;
+    mesh_map[10000] = index_entropy_mesh;
+  }
+
+  if (index_entropy_mesh >= 0) {
+    auto& m = *meshes[index_entropy_mesh];
+    if (m.shape_.dimension() == 0) {
+      // If the user did not specify how many mesh cells are to be used in
+      // each direction, we automatically determine an appropriate number of
+      // cells
+      int n = std::ceil(std::pow(settings::n_particles / 20.0, 1.0/3.0));
+      m.shape_ = {n, n, n};
+      m.n_dimension_ = 3;
+
+      // Calculate width
+      m.width_ = (m.upper_right_ - m.lower_left_) / m.shape_;
+    }
+
+    // Turn on Shannon entropy calculation
+    settings::entropy_on = true;
+  }
+
+  // Uniform fission source weighting mesh
+  if (check_for_node(root, "ufs_mesh")) {
+    auto temp = std::stoi(get_node_value(root, "ufs_mesh"));
+    if (mesh_map.find(temp) == mesh_map.end()) {
+      std::stringstream msg;
+      msg << "Mesh " << temp << " specified for uniform fission site method "
+        "does not exist.";
+      fatal_error(msg);
+    }
+    index_ufs_mesh = mesh_map.at(temp);
+
+  } else if (check_for_node(root, "uniform_fs")) {
+    warning("Specifying a UFS mesh via the <uniform_fs> element "
+      "is deprecated. Please create a mesh using <mesh> and then reference "
+      "it by specifying its ID in a <ufs_mesh> element.");
+
+    // Read entropy mesh from <entropy>
+    auto node_ufs = root.child("uniform_fs");
+    meshes.emplace_back(new RegularMesh{node_ufs});
+
+    // Set entropy mesh index
+    index_ufs_mesh = meshes.size() - 1;
+
+    // Assign ID and set mapping
+    meshes.back()->id_ = 10001;
+    mesh_map[10001] = index_entropy_mesh;
+  }
+
+  if (index_ufs_mesh >= 0) {
+    // Turn on uniform fission source weighting
+    settings::ufs_on = true;
+  }
 
   // TODO: Read <state_point>
 

src/simulation.F90

@@ -10,8 +10,7 @@ module simulation
   use cmfd_execute,    only: cmfd_init_batch, cmfd_tally_init, execute_cmfd
   use cmfd_header,     only: cmfd_on
   use constants,       only: ZERO
-  use eigenvalue,      only: count_source_for_ufs, calculate_average_keff, &
-                             calculate_generation_keff, shannon_entropy, &
+  use eigenvalue,      only: calculate_average_keff, calculate_generation_keff, &
                              synchronize_bank, keff_generation, k_sum
 #ifdef _OPENMP
   use eigenvalue,      only: join_bank_from_threads
@@ -234,12 +233,17 @@ contains
 
   subroutine initialize_generation()
 
+    interface
+      subroutine ufs_count_sites() bind(C)
+      end subroutine
+    end interface
+
     if (run_mode == MODE_EIGENVALUE) then
       ! Reset number of fission bank sites
       n_bank = 0
 
       ! Count source sites if using uniform fission source weighting
-      if (ufs) call count_source_for_ufs()
+      if (ufs) call ufs_count_sites()
 
       ! Store current value of tracklength k
       keff_generation = global_tallies(RESULT_VALUE, K_TRACKLENGTH)
@@ -256,6 +260,9 @@ contains
     interface
       subroutine fill_source_bank_fixedsource() bind(C)
       end subroutine
+
+      subroutine shannon_entropy() bind(C)
+      end subroutine
     end interface
 
     ! Update global tallies with the omp private accumulation variables
@@ -465,13 +472,17 @@ contains
 
     ! Allocate array for matching filter bins
     allocate(filter_matches(n_filters))
+    do i = 1, n_filters
+      allocate(filter_matches(i) % bins)
+      allocate(filter_matches(i) % weights)
+    end do
 !$omp end parallel
 
     ! Reset global variables -- this is done before loading state point (as that
     ! will potentially populate k_generation and entropy)
     current_batch = 0
     call k_generation % clear()
-    call entropy % clear()
+    call entropy_clear()
     need_depletion_rx = .false.
 
     ! If this is a restart run, load the state point data and binary source
@@ -550,6 +561,10 @@ contains
       deallocate(materials(i) % mat_nuclide_index)
     end do
 !$omp parallel
+    do i = 1, size(filter_matches)
+      deallocate(filter_matches(i) % bins)
+      deallocate(filter_matches(i) % weights)
+    end do
     deallocate(micro_xs, micro_photon_xs, filter_matches)
 !$omp end parallel
 

src/simulation_header.F90

@@ -47,13 +47,6 @@ module simulation_header
   real(8) :: k_col_tra = ZERO ! sum over batches of k_collision * k_tracklength
   real(8) :: k_abs_tra = ZERO ! sum over batches of k_absorption * k_tracklength
 
-  ! Shannon entropy
-  type(VectorReal)     :: entropy        ! shannon entropy at each generation
-  real(8), allocatable :: entropy_p(:,:) ! % of source sites in each cell
-
-  ! Uniform fission source weighting
-  real(8), allocatable :: source_frac(:,:)
-
   ! ============================================================================
   ! PARALLEL PROCESSING VARIABLES
 
@@ -71,6 +64,12 @@ module simulation_header
 
 !$omp threadprivate(trace, thread_id, current_work)
 
+  interface
+    subroutine entropy_clear() bind(C)
+    end subroutine
+  end interface
+
+
 contains
 
 !===============================================================================
@@ -87,14 +86,12 @@ contains
 !===============================================================================
 
   subroutine free_memory_simulation()
-    if (allocated(entropy_p)) deallocate(entropy_p)
-    if (allocated(source_frac)) deallocate(source_frac)
+
     if (allocated(work_index)) deallocate(work_index)
 
     call k_generation % clear()
     call k_generation % shrink_to_fit()
-    call entropy % clear()
-    call entropy % shrink_to_fit()
+    call entropy_clear()
   end subroutine free_memory_simulation
 
 end module simulation_header

src/state_point.F90

@@ -20,7 +20,6 @@ module state_point
   use endf,               only: reaction_name
   use error,              only: fatal_error, warning, write_message
   use hdf5_interface
-  use mesh_header,        only: RegularMesh, meshes, n_meshes
   use message_passing
   use mgxs_interface
   use nuclide_header,     only: nuclides
@@ -68,7 +67,7 @@ contains
     integer :: i_xs
     integer, allocatable :: id_array(:)
     integer(HID_T) :: file_id
-    integer(HID_T) :: cmfd_group, tallies_group, tally_group, meshes_group, &
+    integer(HID_T) :: cmfd_group, tallies_group, tally_group, &
                       filters_group, filter_group, derivs_group, &
                       deriv_group, runtime_group
     integer(C_INT) :: ignored_err
@@ -79,6 +78,17 @@ contains
     character(MAX_WORD_LEN, kind=C_CHAR) :: temp_name
     logical :: parallel
 
+    interface
+      subroutine meshes_to_hdf5(group) bind(C)
+        import HID_T
+        integer(HID_T), value :: group
+      end subroutine
+      subroutine entropy_to_hdf5(group) bind(C)
+        import HID_T
+        integer(HID_T), value :: group
+      end subroutine
+    end interface
+
     err = 0
 
     ! Set the filename
@@ -163,9 +173,7 @@ contains
         call write_dataset(file_id, "generations_per_batch", gen_per_batch)
         k = k_generation % size()
         call write_dataset(file_id, "k_generation", k_generation % data(1:k))
-        if (entropy_on) then
-          call write_dataset(file_id, "entropy", entropy % data(1:k))
-        end if
+        call entropy_to_hdf5(file_id)
         call write_dataset(file_id, "k_col_abs", k_col_abs)
         call write_dataset(file_id, "k_col_tra", k_col_tra)
         call write_dataset(file_id, "k_abs_tra", k_abs_tra)
@@ -192,26 +200,8 @@ contains
 
       tallies_group = create_group(file_id, "tallies")
 
-      ! Write number of meshes
-      meshes_group = create_group(tallies_group, "meshes")
-      call write_attribute(meshes_group, "n_meshes", n_meshes)
-
-      if (n_meshes > 0) then
-        ! Write IDs of meshes
-        allocate(id_array(n_meshes))
-        do i = 1, n_meshes
-          id_array(i) = meshes(i) % id
-        end do
-        call write_attribute(meshes_group, "ids", id_array)
-        deallocate(id_array)
-
-        ! Write information for meshes
-        MESH_LOOP: do i = 1, n_meshes
-          call meshes(i) % to_hdf5(meshes_group)
-        end do MESH_LOOP
-      end if
-
-      call close_group(meshes_group)
+      ! Write meshes
+      call meshes_to_hdf5(tallies_group)
 
       ! Write information for derivatives.
       if (size(tally_derivs) > 0) then
@@ -641,6 +631,11 @@ contains
     logical :: source_present
     character(MAX_WORD_LEN) :: word
 
+    interface
+      subroutine entropy_from_hdf5() bind(C)
+      end subroutine
+    end interface
+
     ! Write message
     call write_message("Loading state point " // trim(path_state_point) &
          // "...", 5)
@@ -722,10 +717,7 @@ contains
       call k_generation % resize(n)
       call read_dataset(k_generation % data(1:n), file_id, "k_generation")
 
-      if (entropy_on) then
-        call entropy % resize(n)
-        call read_dataset(entropy % data(1:n), file_id, "entropy")
-      end if
+      call entropy_from_hdf5()
       call read_dataset(k_col_abs, file_id, "k_col_abs")
       call read_dataset(k_col_tra, file_id, "k_col_tra")
       call read_dataset(k_abs_tra, file_id, "k_abs_tra")

src/stl_vector.F90

@@ -35,6 +35,8 @@ module stl_vector
   !
   ! size -- Returns the number of elements in the vector.
 
+  use, intrinsic :: ISO_C_BINDING
+
   implicit none
   private
 
@@ -521,4 +523,28 @@ contains
     size = this%size_
   end function size_char
 
+!===============================================================================
+! Procedures to be called from C++
+!===============================================================================
+
+  subroutine vector_int_push_back(ptr, val) bind(C)
+    type(C_PTR), value :: ptr
+    integer(C_INT), value :: val
+
+    type(VectorInt), pointer :: vec
+
+    call C_F_POINTER(ptr, vec)
+    call vec % push_back(val)
+  end subroutine
+
+  subroutine vector_real_push_back(ptr, val) bind(C)
+    type(C_PTR), value :: ptr
+    real(C_DOUBLE), value :: val
+
+    type(VectorReal), pointer :: vec
+
+    call C_F_POINTER(ptr, vec)
+    call vec % push_back(val)
+  end subroutine
+
 end module stl_vector

src/summary.F90

@@ -6,7 +6,6 @@ module summary
   use geometry_header
   use hdf5_interface
   use material_header, only: Material, n_materials, openmc_material_get_volume
-  use mesh_header,     only: RegularMesh
   use message_passing
   use mgxs_interface
   use nuclide_header

src/tallies/tally.F90

@@ -8,7 +8,6 @@ module tally
   use error,            only: fatal_error
   use geometry_header
   use math,             only: t_percentile
-  use mesh_header,      only: RegularMesh, meshes
   use message_passing
   use mgxs_interface
   use nuclide_header

src/tallies/tally_filter_header.F90

@@ -28,8 +28,8 @@ module tally_filter_header
   type, public :: TallyFilterMatch
     ! Index of the bin and weight being used in the current filter combination
     integer          :: i_bin
-    type(VectorInt)  :: bins
-    type(VectorReal) :: weights
+    type(VectorInt),  pointer :: bins
+    type(VectorReal), pointer :: weights
 
     ! Indicates whether all valid bins for this filter have been found
     logical          :: bins_present = .false.

src/tallies/tally_filter_mesh.F90

@@ -5,7 +5,7 @@ module tally_filter_mesh
   use constants
   use dict_header,         only: EMPTY
   use error
-  use mesh_header,         only: RegularMesh, meshes, n_meshes, mesh_dict
+  use mesh_header
   use hdf5_interface
   use particle_header,     only: Particle
   use string,              only: to_str
@@ -38,10 +38,11 @@ contains
     class(MeshFilter), intent(inout) :: this
     type(XMLNode), intent(in) :: node
 
-    integer :: i_mesh
+    integer :: i
     integer :: id
     integer :: n
-    integer :: val
+    integer(C_INT) :: err
+    type(RegularMesh) :: m
 
     n = node_word_count(node, "bins")
 
@@ -52,19 +53,18 @@ contains
     call get_node_value(node, "bins", id)
 
     ! Get pointer to mesh
-    val = mesh_dict % get(id)
-    if (val /= EMPTY) then
-      i_mesh = val
-    else
+    err = openmc_get_mesh_index(id, this % mesh)
+    if (err /= 0) then
       call fatal_error("Could not find mesh " // trim(to_str(id)) &
            // " specified on filter.")
     end if
 
     ! Determine number of bins
-    this % n_bins = product(meshes(i_mesh) % dimension)
-
-    ! Store the index of the mesh
-    this % mesh = i_mesh
+    m = meshes(this % mesh)
+    this % n_bins = 1
+    do i = 1, m % n_dimension()
+      this % n_bins = this % n_bins * m % dimension(i)
+    end do
   end subroutine from_xml
 
   subroutine get_all_bins_mesh(this, p, estimator, match)
@@ -73,184 +73,50 @@ contains
     integer,           intent(in)  :: estimator
     type(TallyFilterMatch), intent(inout) :: match
 
-    integer, parameter :: MAX_SEARCH_ITER = 100 ! Maximum number of times we can
-                                                !  can loop while trying to find
-                                                !  the first intersection.
-
-    integer :: j                    ! loop index for direction
-    integer :: n
-    integer :: ijk0(3)              ! indices of starting coordinates
-    integer :: ijk1(3)              ! indices of ending coordinates
-    integer :: search_iter          ! loop count for intersection search
     integer :: bin
-    real(8) :: uvw(3)               ! cosine of angle of particle
-    real(8) :: xyz0(3)              ! starting/intermediate coordinates
-    real(8) :: xyz1(3)              ! ending coordinates of particle
-    real(8) :: xyz_cross            ! coordinates of next boundary
-    real(8) :: d(3)                 ! distance to each bounding surface
-    real(8) :: total_distance       ! distance of entire particle track
-    real(8) :: distance             ! distance traveled in mesh cell
-    logical :: start_in_mesh        ! starting coordinates inside mesh?
-    logical :: end_in_mesh          ! ending coordinates inside mesh?
-    type(RegularMesh), pointer :: m
+    type(RegularMesh) :: m
+    type(C_PTR) :: ptr_bins, ptr_weights
+
+    interface
+      subroutine mesh_bins_crossed(m, p, bins, weights) bind(C)
+        import C_PTR, Particle
+        type(C_PTR), value :: m
+        type(Particle), intent(in) :: p
+        type(C_PTR), value :: bins
+        type(C_PTR), value :: weights
+      end subroutine
+    end interface
 
     ! Get a pointer to the mesh.
-    m => meshes(this % mesh)
-    n = m % n_dimension
+    m = meshes(this % mesh)
 
     if (estimator /= ESTIMATOR_TRACKLENGTH) then
       ! If this is an analog or collision tally, then there can only be one
       ! valid mesh bin.
       call m % get_bin(p % coord(1) % xyz, bin)
-      if (bin /= NO_BIN_FOUND) then
+      if (bin >= 0) then
         call match % bins % push_back(bin)
         call match % weights % push_back(ONE)
       end if
       return
+    else
+      ptr_bins = C_LOC(match % bins)
+      ptr_weights = C_LOC(match % weights)
+      call mesh_bins_crossed(m % ptr, p, ptr_bins, ptr_weights)
     end if
 
-    ! A track can span multiple mesh bins so we need to handle a lot of
-    ! intersection logic for tracklength tallies.
-
-    ! ========================================================================
-    ! Determine if the track intersects the tally mesh.
-
-    ! Copy the starting and ending coordinates of the particle.  Offset these
-    ! just a bit for the purposes of determining if there was an intersection
-    ! in case the mesh surfaces coincide with lattice/geometric surfaces which
-    ! might produce finite-precision errors.
-    xyz0 = p % last_xyz + TINY_BIT * p % coord(1) % uvw
-    xyz1 = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw
-
-    ! Determine indices for starting and ending location.
-    call m % get_indices(xyz0, ijk0(:n), start_in_mesh)
-    call m % get_indices(xyz1, ijk1(:n), end_in_mesh)
-
-    ! If this is the first iteration of the filter loop, check if the track
-    ! intersects any part of the mesh.
-    if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then
-      if (.not. m % intersects(xyz0, xyz1)) return
-    end if
-
-    ! ========================================================================
-    ! Figure out which mesh cell to tally.
-
-    ! Copy the un-modified coordinates the particle direction.
-    xyz0 = p % last_xyz
-    xyz1 = p % coord(1) % xyz
-    uvw = p % coord(1) % uvw
-
-    ! Compute the length of the entire track.
-    total_distance = sqrt(sum((xyz1 - xyz0)**2))
-
-    ! We are looking for the first valid mesh bin.  Check to see if the
-    ! particle starts inside the mesh.
-    if (any(ijk0(:n) < 1) .or. any(ijk0(:n) > m % dimension)) then
-      ! The particle does not start in the mesh.  Note that we nudged the
-      ! start and end coordinates by a TINY_BIT each so we will have
-      ! difficulty resolving tracks that are less than 2*TINY_BIT in length.
-      ! If the track is that short, it is also insignificant so we can
-      ! safely ignore it in the tallies.
-      if (total_distance < 2*TINY_BIT) return
-
-      ! The particle does not start in the mesh so keep iterating the ijk0
-      ! indices to cross the nearest mesh surface until we've found a valid
-      ! bin.  MAX_SEARCH_ITER prevents an infinite loop.
-      search_iter = 0
-      do while (any(ijk0(:n) < 1) .or. any(ijk0(:n) > m % dimension))
-        if (search_iter == MAX_SEARCH_ITER) then
-          call warning("Failed to find a mesh intersection on a tally mesh &
-               &filter.")
-          return
-        end if
-
-        do j = 1, n
-          if (abs(uvw(j)) < FP_PRECISION) then
-            d(j) = INFINITY
-          else if (uvw(j) > 0) then
-            xyz_cross = m % lower_left(j) + ijk0(j) * m % width(j)
-            d(j) = (xyz_cross - xyz0(j)) / uvw(j)
-          else
-            xyz_cross = m % lower_left(j) + (ijk0(j) - 1) * m % width(j)
-            d(j) = (xyz_cross - xyz0(j)) / uvw(j)
-          end if
-        end do
-        j = minloc(d(:n), 1)
-        if (uvw(j) > ZERO) then
-          ijk0(j) = ijk0(j) + 1
-        else
-          ijk0(j) = ijk0(j) - 1
-        end if
-
-        search_iter = search_iter + 1
-      end do
-      distance = d(j)
-      xyz0 = xyz0 + distance * uvw
-    end if
-
-    do
-      ! ========================================================================
-      ! Compute the length of the track segment in the appropiate mesh cell and
-      ! return.
-
-      if (all(ijk0(:n) == ijk1(:n))) then
-        ! The track ends in this cell.  Use the particle end location rather
-        ! than the mesh surface.
-        distance = sqrt(sum((xyz1 - xyz0)**2))
-      else
-        ! The track exits this cell.  Determine the distance to the closest mesh
-        ! surface.
-        do j = 1, n
-          if (abs(uvw(j)) < FP_PRECISION) then
-            d(j) = INFINITY
-          else if (uvw(j) > 0) then
-            xyz_cross = m % lower_left(j) + ijk0(j) * m % width(j)
-            d(j) = (xyz_cross - xyz0(j)) / uvw(j)
-          else
-            xyz_cross = m % lower_left(j) + (ijk0(j) - 1) * m % width(j)
-            d(j) = (xyz_cross - xyz0(j)) / uvw(j)
-          end if
-        end do
-        j = minloc(d(:n), 1)
-        distance = d(j)
-      end if
-
-      ! Assign the next tally bin and the score.
-      bin = m % get_bin_from_indices(ijk0(:n))
-      call match % bins % push_back(bin)
-      call match % weights % push_back(distance / total_distance)
-
-      ! Find the next mesh cell that the particle enters.
-
-      ! If the particle track ends in that bin, then we are done.
-      if (all(ijk0(:n) == ijk1(:n))) exit
-
-      ! Translate the starting coordintes by the distance to that face. This
-      ! should be the xyz that we computed the distance to in the last
-      ! iteration of the filter loop.
-      xyz0 = xyz0 + distance * uvw
-
-      ! Increment the indices into the next mesh cell.
-      if (uvw(j) > ZERO) then
-        ijk0(j) = ijk0(j) + 1
-      else
-        ijk0(j) = ijk0(j) - 1
-      end if
-
-      ! If the next indices are invalid, then the track has left the mesh and
-      ! we are done.
-      if (any(ijk0(:n) < 1) .or. any(ijk0(:n) > m % dimension)) exit
-    end do
-
   end subroutine get_all_bins_mesh
 
   subroutine to_statepoint_mesh(this, filter_group)
     class(MeshFilter), intent(in) :: this
     integer(HID_T),    intent(in) :: filter_group
 
+    type(RegularMesh) :: m
+
+    m = meshes(this % mesh)
     call write_dataset(filter_group, "type", "mesh")
     call write_dataset(filter_group, "n_bins", this % n_bins)
-    call write_dataset(filter_group, "bins", meshes(this % mesh) % id)
+    call write_dataset(filter_group, "bins", m % id())
   end subroutine to_statepoint_mesh
 
   function text_label_mesh(this, bin) result(label)
@@ -259,20 +125,20 @@ contains
     character(MAX_LINE_LEN)       :: label
 
     integer, allocatable       :: ijk(:)
+    type(RegularMesh) :: m
 
-    associate (m => meshes(this % mesh))
-      allocate(ijk(m % n_dimension))
-      call m % get_indices_from_bin(bin, ijk)
-      if (m % n_dimension == 1) then
-        label = "Mesh Index (" // trim(to_str(ijk(1))) // ")"
-      elseif (m % n_dimension == 2) then
-        label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // &
-             trim(to_str(ijk(2))) // ")"
-      elseif (m % n_dimension == 3) then
-        label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // &
-             trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")"
-      end if
-    end associate
+    m = meshes(this % mesh)
+    allocate(ijk(m % n_dimension()))
+    call m % get_indices_from_bin(bin, ijk)
+    if (m % n_dimension() == 1) then
+      label = "Mesh Index (" // trim(to_str(ijk(1))) // ")"
+    elseif (m % n_dimension() == 2) then
+      label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // &
+            trim(to_str(ijk(2))) // ")"
+    elseif (m % n_dimension() == 3) then
+      label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // &
+            trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")"
+    end if
   end function text_label_mesh
 
 !===============================================================================
@@ -304,18 +170,25 @@ contains
     integer(C_INT32_T), value, intent(in) :: index_mesh
     integer(C_INT) :: err
 
+    type(RegularMesh) :: m
+    integer :: i
+
     err = verify_filter(index)
     if (err == 0) then
       select type (f => filters(index) % obj)
       type is (MeshFilter)
-        if (index_mesh >= 1 .and. index_mesh <= n_meshes) then
+        if (index_mesh >= 0 .and. index_mesh < n_meshes()) then
           f % mesh = index_mesh
-          f % n_bins = product(meshes(index_mesh) % dimension)
+          f % n_bins = 1
+          m = meshes(index_mesh)
+          do i = 1, m % n_dimension()
+            f % n_bins = f % n_bins * m % dimension(i)
+          end do
         else
           err = E_OUT_OF_BOUNDS
           call set_errmsg("Index in 'meshes' array is out of bounds.")
         end if
-        class default
+      class default
         err = E_INVALID_TYPE
         call set_errmsg("Tried to set mesh on a non-mesh filter.")
       end select

src/tallies/tally_filter_meshsurface.F90

@@ -3,9 +3,8 @@ module tally_filter_meshsurface
   use, intrinsic :: ISO_C_BINDING
 
   use constants
-  use dict_header,         only: EMPTY
   use error
-  use mesh_header,         only: RegularMesh, meshes, n_meshes, mesh_dict
+  use mesh_header
   use hdf5_interface
   use particle_header,     only: Particle
   use string,              only: to_str
@@ -38,11 +37,11 @@ contains
     class(MeshSurfaceFilter), intent(inout) :: this
     type(XMLNode), intent(in) :: node
 
-    integer :: i_mesh
+    integer :: i
     integer :: id
     integer :: n
-    integer :: n_dim
-    integer :: val
+    integer(C_INT) :: err
+    type(RegularMesh) :: m
 
     n = node_word_count(node, "bins")
 
@@ -53,20 +52,18 @@ contains
     call get_node_value(node, "bins", id)
 
     ! Get pointer to mesh
-    val = mesh_dict % get(id)
-    if (val /= EMPTY) then
-      i_mesh = val
-    else
+    err = openmc_get_mesh_index(id, this % mesh)
+    if (err /= 0) then
       call fatal_error("Could not find mesh " // trim(to_str(id)) &
            // " specified on filter.")
     end if
 
     ! Determine number of bins
-    n_dim = meshes(i_mesh) % n_dimension
-    this % n_bins = 4*n_dim*product(meshes(i_mesh) % dimension)
-
-    ! Store the index of the mesh
-    this % mesh = i_mesh
+    m = meshes(this % mesh)
+    this % n_bins = 4 * m % n_dimension()
+    do i = 1, m % n_dimension()
+      this % n_bins = this % n_bins * m % dimension(i)
+    end do
   end subroutine from_xml
 
   subroutine get_all_bins(this, p, estimator, match)
@@ -75,150 +72,25 @@ contains
     integer,                  intent(in)  :: estimator
     type(TallyFilterMatch), intent(inout) :: match
 
-    integer :: j                    ! loop indices
-    integer :: n_dim                ! num dimensions of the mesh
-    integer :: d1                   ! dimension index
-    integer :: ijk0(3)              ! indices of starting coordinates
-    integer :: ijk1(3)              ! indices of ending coordinates
-    integer :: n_cross              ! number of surface crossings
-    integer :: i_mesh               ! flattened mesh bin index
-    integer :: i_surf               ! surface index (1--12)
-    integer :: i_bin                ! actual index for filter
-    real(8) :: uvw(3)               ! cosine of angle of particle
-    real(8) :: xyz0(3)              ! starting/intermediate coordinates
-    real(8) :: xyz1(3)              ! ending coordinates of particle
-    real(8) :: xyz_cross(3)         ! coordinates of bounding surfaces
-    real(8) :: d(3)                 ! distance to each bounding surface
-    real(8) :: distance             ! actual distance traveled
-    logical :: start_in_mesh        ! particle's starting xyz in mesh?
-    logical :: end_in_mesh          ! particle's ending xyz in mesh?
+    type(RegularMesh) :: m
+    type(C_PTR) :: ptr_bins, ptr_weights
 
-    ! Copy starting and ending location of particle
-    xyz0 = p % last_xyz_current
-    xyz1 = p % coord(1) % xyz
+    interface
+      subroutine mesh_surface_bins_crossed(m, p, bins, weights) bind(C)
+        import C_PTR, Particle
+        type(C_PTR), value :: m
+        type(Particle), intent(in) :: p
+        type(C_PTR), value :: bins
+        type(C_PTR), value :: weights
+      end subroutine
+    end interface
 
-    associate (m => meshes(this % mesh))
-      n_dim = m % n_dimension
+    ! Get a pointer to the mesh.
+    m = meshes(this % mesh)
 
-      ! Determine indices for starting and ending location
-      call m % get_indices(xyz0, ijk0, start_in_mesh)
-      call m % get_indices(xyz1, ijk1, end_in_mesh)
-
-      ! Check to see if start or end is in mesh -- if not, check if track still
-      ! intersects with mesh
-      if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then
-        if (.not. m % intersects(xyz0, xyz1)) return
-      end if
-
-      ! Calculate number of surface crossings
-      n_cross = sum(abs(ijk1(:n_dim) - ijk0(:n_dim)))
-      if (n_cross == 0) return
-
-      ! Copy particle's direction
-      uvw = p % coord(1) % uvw
-
-      ! Bounding coordinates
-      do d1 = 1, n_dim
-        if (uvw(d1) > 0) then
-          xyz_cross(d1) = m % lower_left(d1) + ijk0(d1) * m % width(d1)
-        else
-          xyz_cross(d1) = m % lower_left(d1) + (ijk0(d1) - 1) * m % width(d1)
-        end if
-      end do
-
-      do j = 1, n_cross
-        ! Set the distances to infinity
-        d = INFINITY
-
-        ! Calculate distance to each bounding surface. We need to treat
-        ! special case where the cosine of the angle is zero since this would
-        ! result in a divide-by-zero.
-        do d1 = 1, n_dim
-          if (uvw(d1) == 0) then
-            d(d1) = INFINITY
-          else
-            d(d1) = (xyz_cross(d1) - xyz0(d1))/uvw(d1)
-          end if
-        end do
-
-        ! Determine the closest bounding surface of the mesh cell by
-        ! calculating the minimum distance. Then use the minimum distance and
-        ! direction of the particle to determine which surface was crossed.
-        distance = minval(d)
-
-        ! Loop over the dimensions
-        do d1 = 1, n_dim
-
-          ! Check whether distance is the shortest distance
-          if (distance == d(d1)) then
-
-            ! Check whether particle is moving in positive d1 direction
-            if (uvw(d1) > 0) then
-
-              ! Outward current on d1 max surface
-              if (all(ijk0(:n_dim) >= 1) .and. &
-                   all(ijk0(:n_dim) <= m % dimension)) then
-                i_surf = d1 * 4 - 1
-                i_mesh = m % get_bin_from_indices(ijk0)
-                i_bin = 4*n_dim*(i_mesh - 1) + i_surf
-
-                call match % bins % push_back(i_bin)
-                call match % weights % push_back(ONE)
-              end if
-
-              ! Advance position
-              ijk0(d1) = ijk0(d1) + 1
-              xyz_cross(d1) = xyz_cross(d1) + m % width(d1)
-
-              ! If the particle crossed the surface, tally the inward current on
-              ! d1 min surface
-              if (all(ijk0(:n_dim) >= 1) .and. &
-                   all(ijk0(:n_dim) <= m % dimension)) then
-                i_surf = d1 * 4 - 2
-                i_mesh = m % get_bin_from_indices(ijk0)
-                i_bin = 4*n_dim*(i_mesh - 1) + i_surf
-
-                call match % bins % push_back(i_bin)
-                call match % weights % push_back(ONE)
-              end if
-
-            else
-              ! The particle is moving in the negative d1 direction
-
-              ! Outward current on d1 min surface
-              if (all(ijk0(:n_dim) >= 1) .and. &
-                   all(ijk0(:n_dim) <= m % dimension)) then
-                i_surf = d1 * 4 - 3
-                i_mesh = m % get_bin_from_indices(ijk0)
-                i_bin = 4*n_dim*(i_mesh - 1) + i_surf
-
-                call match % bins % push_back(i_bin)
-                call match % weights % push_back(ONE)
-              end if
-
-              ! Advance position
-              ijk0(d1) = ijk0(d1) - 1
-              xyz_cross(d1) = xyz_cross(d1) - m % width(d1)
-
-              ! If the particle crossed the surface, tally the inward current on
-              ! d1 max surface
-              if (all(ijk0(:n_dim) >= 1) .and. &
-                   all(ijk0(:n_dim) <= m % dimension)) then
-                i_surf = d1 * 4
-                i_mesh = m % get_bin_from_indices(ijk0)
-                i_bin = 4*n_dim*(i_mesh - 1) + i_surf
-
-                call match % bins % push_back(i_bin)
-                call match % weights % push_back(ONE)
-              end if
-            end if
-          end if
-        end do
-
-        ! Calculate new coordinates
-        xyz0 = xyz0 + distance * uvw
-      end do
-    end associate
+    ptr_bins = C_LOC(match % bins)
+    ptr_weights = C_LOC(match % weights)
+    call mesh_surface_bins_crossed(m % ptr, p, ptr_bins, ptr_weights)
 
   end subroutine get_all_bins
 
@@ -226,9 +98,12 @@ contains
     class(MeshSurfaceFilter), intent(in) :: this
     integer(HID_T),           intent(in) :: filter_group
 
+    type(RegularMesh) :: m
+
+    m = meshes(this % mesh)
     call write_dataset(filter_group, "type", "meshsurface")
     call write_dataset(filter_group, "n_bins", this % n_bins)
-    call write_dataset(filter_group, "bins", meshes(this % mesh) % id)
+    call write_dataset(filter_group, "bins", m % id())
   end subroutine to_statepoint
 
   function text_label(this, bin) result(label)
@@ -240,55 +115,55 @@ contains
     integer :: i_surf
     integer :: n_dim
     integer, allocatable :: ijk(:)
+    type(RegularMesh) :: m
 
-    associate (m => meshes(this % mesh))
-      n_dim = m % n_dimension
-      allocate(ijk(n_dim))
+    m = meshes(this % mesh)
+    n_dim = m % n_dimension()
+    allocate(ijk(n_dim))
 
-      ! Get flattend mesh index and surface index
-      i_mesh = (bin - 1) / (4*n_dim) + 1
-      i_surf = mod(bin - 1, 4*n_dim) + 1
+    ! Get flattend mesh index and surface index
+    i_mesh = (bin - 1) / (4*n_dim) + 1
+    i_surf = mod(bin - 1, 4*n_dim) + 1
 
-      ! Get mesh index part of label
-      call m % get_indices_from_bin(i_mesh, ijk)
-      if (m % n_dimension == 1) then
-        label = "Mesh Index (" // trim(to_str(ijk(1))) // ")"
-      elseif (m % n_dimension == 2) then
-        label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // &
-             trim(to_str(ijk(2))) // ")"
-      elseif (m % n_dimension == 3) then
-        label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // &
-             trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")"
-      end if
+    ! Get mesh index part of label
+    call m % get_indices_from_bin(i_mesh, ijk)
+    if (m % n_dimension() == 1) then
+      label = "Mesh Index (" // trim(to_str(ijk(1))) // ")"
+    elseif (m % n_dimension() == 2) then
+      label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // &
+            trim(to_str(ijk(2))) // ")"
+    elseif (m % n_dimension() == 3) then
+      label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // &
+            trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")"
+    end if
 
-      ! Get surface part of label
-      select case (i_surf)
-      case (OUT_LEFT)
-        label = trim(label) // " Outgoing, x-min"
-      case (IN_LEFT)
-        label = trim(label) // " Incoming, x-min"
-      case (OUT_RIGHT)
-        label = trim(label) // " Outgoing, x-max"
-      case (IN_RIGHT)
-        label = trim(label) // " Incoming, x-max"
-      case (OUT_BACK)
-        label = trim(label) // " Outgoing, y-min"
-      case (IN_BACK)
-        label = trim(label) // " Incoming, y-min"
-      case (OUT_FRONT)
-        label = trim(label) // " Outgoing, y-max"
-      case (IN_FRONT)
-        label = trim(label) // " Incoming, y-max"
-      case (OUT_BOTTOM)
-        label = trim(label) // " Outgoing, z-min"
-      case (IN_BOTTOM)
-        label = trim(label) // " Incoming, z-min"
-      case (OUT_TOP)
-        label = trim(label) // " Outgoing, z-max"
-      case (IN_TOP)
-        label = trim(label) // " Incoming, z-max"
-      end select
-    end associate
+    ! Get surface part of label
+    select case (i_surf)
+    case (OUT_LEFT)
+      label = trim(label) // " Outgoing, x-min"
+    case (IN_LEFT)
+      label = trim(label) // " Incoming, x-min"
+    case (OUT_RIGHT)
+      label = trim(label) // " Outgoing, x-max"
+    case (IN_RIGHT)
+      label = trim(label) // " Incoming, x-max"
+    case (OUT_BACK)
+      label = trim(label) // " Outgoing, y-min"
+    case (IN_BACK)
+      label = trim(label) // " Incoming, y-min"
+    case (OUT_FRONT)
+      label = trim(label) // " Outgoing, y-max"
+    case (IN_FRONT)
+      label = trim(label) // " Incoming, y-max"
+    case (OUT_BOTTOM)
+      label = trim(label) // " Outgoing, z-min"
+    case (IN_BOTTOM)
+      label = trim(label) // " Incoming, z-min"
+    case (OUT_TOP)
+      label = trim(label) // " Outgoing, z-max"
+    case (IN_TOP)
+      label = trim(label) // " Incoming, z-max"
+    end select
   end function text_label
 
 !===============================================================================
@@ -320,16 +195,22 @@ contains
     integer(C_INT32_T), value, intent(in) :: index_mesh
     integer(C_INT) :: err
 
+    integer :: i
     integer :: n_dim
+    type(RegularMesh) :: m
 
     err = verify_filter(index)
     if (err == 0) then
       select type (f => filters(index) % obj)
       type is (MeshSurfaceFilter)
-        if (index_mesh >= 1 .and. index_mesh <= n_meshes) then
+        if (index_mesh >= 0 .and. index_mesh < n_meshes()) then
           f % mesh = index_mesh
-          n_dim = meshes(index_mesh) % n_dimension
-          f % n_bins = 4*n_dim*product(meshes(index_mesh) % dimension)
+          m = meshes(index_mesh)
+          n_dim = m % n_dimension()
+          f % n_bins = 4*n_dim
+          do i = 1, n_dim
+            f % n_bins = f % n_bins * m % dimension(i)
+          end do
         else
           err = E_OUT_OF_BOUNDS
           call set_errmsg("Index in 'meshes' array is out of bounds.")

src/tallies/trigger.F90

@@ -257,14 +257,14 @@ contains
     logical :: print_ebin           ! should incoming energy bin be displayed?
     real(8) :: rel_err  = ZERO      ! temporary relative error of result
     real(8) :: std_dev  = ZERO      ! temporary standard deviration of result
-    type(RegularMesh), pointer :: m        ! surface current mesh
+    type(RegularMesh) :: m        ! surface current mesh
 
     ! Get pointer to mesh
     i_filter_mesh = t % filter(t % find_filter(FILTER_MESH))
     i_filter_surf = t % filter(t % find_filter(FILTER_SURFACE))
     select type(filt => filters(i_filter_mesh) % obj)
     type is (MeshFilter)
-      m => meshes(filt % mesh)
+      m = meshes(filt % mesh)
     end select
 
     ! initialize bins array
@@ -285,8 +285,11 @@ contains
     end if
 
     ! Get the dimensions and number of cells in the mesh
-    n_dim = m % n_dimension
-    n_cells = product(m % dimension)
+    n_dim = m % n_dimension()
+    n_cells = 1
+    do j = 1, n_dim
+      n_cells = n_cells * m % dimension(j)
+    end do
 
     ! Loop over all the mesh cells
     do i = 1, n_cells