Add C API bindings to create/modify meshes (and filters)

docs/source/pythonapi/capi.rst

@@ -44,5 +44,8 @@ Classes
    EnergyFilter
    MaterialFilter
    Material
+   Mesh
+   MeshFilter
+   MeshSurfaceFilter
    Nuclide
    Tally

include/openmc.h

@@ -53,7 +53,15 @@ extern "C" {
   int openmc_material_set_id(int32_t index, int32_t id);
   int openmc_material_filter_get_bins(int32_t index, int32_t** bins, int32_t* n);
   int openmc_material_filter_set_bins(int32_t index, int32_t n, const int32_t* bins);
+  int openmc_mesh_filter_get_mesh(int32_t index, int32_t* index_mesh);
   int openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh);
+  int openmc_mesh_get_id(int32_t index, int32_t* id);
+  int openmc_mesh_get_dimension(int32_t index, int** id, int* n);
+  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_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 openmc_nuclide_name(int index, char** name);

openmc/capi/__init__.py

@@ -43,6 +43,7 @@ from .core import *
 from .nuclide import *
 from .material import *
 from .cell import *
+from .mesh import *
 from .filter import *
 from .tally import *
 from .settings import settings

openmc/capi/cell.py

@@ -44,7 +44,7 @@ class Cell(_FortranObjectWithID):
 
     This class exposes a cell that is stored internally in the OpenMC
     library. To obtain a view of a cell with a given ID, use the
-    :data:`openmc.capi.nuclides` mapping.
+    :data:`openmc.capi.cells` mapping.
 
     Parameters
     ----------

openmc/capi/filter.py

@@ -10,13 +10,14 @@ from . import _dll
 from .core import _FortranObjectWithID
 from .error import _error_handler, AllocationError, InvalidIDError
 from .material import Material
+from .mesh import Mesh
 
 
 __all__ = ['Filter', 'AzimuthalFilter', 'CellFilter',
            'CellbornFilter', 'CellfromFilter', 'DistribcellFilter',
            'DelayedGroupFilter', 'EnergyFilter', 'EnergyoutFilter',
            'EnergyFunctionFilter', 'MaterialFilter', 'MeshFilter',
-           'MuFilter', 'PolarFilter', 'SurfaceFilter',
+           'MeshSurfaceFilter', 'MuFilter', 'PolarFilter', 'SurfaceFilter',
            'UniverseFilter', 'filters']
 
 # Tally functions
@@ -52,9 +53,15 @@ _dll.openmc_material_filter_get_bins.errcheck = _error_handler
 _dll.openmc_material_filter_set_bins.argtypes = [c_int32, c_int32, POINTER(c_int32)]
 _dll.openmc_material_filter_set_bins.restype = c_int
 _dll.openmc_material_filter_set_bins.errcheck = _error_handler
+_dll.openmc_mesh_filter_get_mesh.argtypes = [c_int32, POINTER(c_int32)]
+_dll.openmc_mesh_filter_get_mesh.restype = c_int
+_dll.openmc_mesh_filter_get_mesh.errcheck = _error_handler
 _dll.openmc_mesh_filter_set_mesh.argtypes = [c_int32, c_int32]
 _dll.openmc_mesh_filter_set_mesh.restype = c_int
 _dll.openmc_mesh_filter_set_mesh.errcheck = _error_handler
+_dll.openmc_meshsurface_filter_get_mesh.argtypes = [c_int32, POINTER(c_int32)]
+_dll.openmc_meshsurface_filter_get_mesh.restype = c_int
+_dll.openmc_meshsurface_filter_get_mesh.errcheck = _error_handler
 _dll.openmc_meshsurface_filter_set_mesh.argtypes = [c_int32, c_int32]
 _dll.openmc_meshsurface_filter_set_mesh.restype = c_int
 _dll.openmc_meshsurface_filter_set_mesh.errcheck = _error_handler
@@ -190,10 +197,40 @@ class MaterialFilter(Filter):
 class MeshFilter(Filter):
     filter_type = 'mesh'
 
+    def __init__(self, mesh=None, uid=None, new=True, index=None):
+        super().__init__(uid, new, index)
+        if mesh is not None:
+            self.mesh = mesh
+
+    @property
+    def mesh(self):
+        index_mesh = c_int32()
+        _dll.openmc_mesh_filter_get_mesh(self._index, index_mesh)
+        return Mesh(index=index_mesh.value)
+
+    @mesh.setter
+    def mesh(self, mesh):
+        _dll.openmc_mesh_filter_set_mesh(self._index, mesh._index)
+
 
 class MeshSurfaceFilter(Filter):
     filter_type = 'meshsurface'
 
+    def __init__(self, mesh=None, uid=None, new=True, index=None):
+        super().__init__(uid, new, index)
+        if mesh is not None:
+            self.mesh = mesh
+
+    @property
+    def mesh(self):
+        index_mesh = c_int32()
+        _dll.openmc_meshsurface_filter_get_mesh(self._index, index_mesh)
+        return Mesh(index=index_mesh.value)
+
+    @mesh.setter
+    def mesh(self, mesh):
+        _dll.openmc_meshsurface_filter_set_mesh(self._index, mesh._index)
+
 
 class MuFilter(Filter):
     filter_type = 'mu'

openmc/capi/mesh.py

@@ -0,0 +1,172 @@
+from collections.abc import Mapping, Iterable
+from ctypes import c_int, c_int32, c_double, POINTER
+from weakref import WeakValueDictionary
+
+import numpy as np
+from numpy.ctypeslib import as_array
+
+from . import _dll
+from .core import _FortranObjectWithID
+from .error import _error_handler, AllocationError, InvalidIDError
+from .material import Material
+
+__all__ = ['Mesh', 'meshes']
+
+# Mesh functions
+_dll.openmc_extend_meshes.argtypes = [c_int32, POINTER(c_int32), POINTER(c_int32)]
+_dll.openmc_extend_meshes.restype = c_int
+_dll.openmc_extend_meshes.errcheck = _error_handler
+_dll.openmc_mesh_get_id.argtypes = [c_int32, POINTER(c_int32)]
+_dll.openmc_mesh_get_id.restype = c_int
+_dll.openmc_mesh_get_id.errcheck = _error_handler
+_dll.openmc_mesh_get_dimension.argtypes = [c_int32, POINTER(POINTER(c_int)), POINTER(c_int)]
+_dll.openmc_mesh_get_dimension.restype = c_int
+_dll.openmc_mesh_get_dimension.errcheck = _error_handler
+_dll.openmc_mesh_get_params.argtypes = [
+    c_int32, POINTER(POINTER(c_double)), POINTER(POINTER(c_double)),
+    POINTER(POINTER(c_double)), POINTER(c_int)]
+_dll.openmc_mesh_get_params.restype = c_int
+_dll.openmc_mesh_get_params.errcheck = _error_handler
+_dll.openmc_mesh_set_id.argtypes = [c_int32, c_int32]
+_dll.openmc_mesh_set_id.restype = c_int
+_dll.openmc_mesh_set_id.errcheck = _error_handler
+_dll.openmc_mesh_set_dimension.argtypes = [c_int32, c_int, POINTER(c_int)]
+_dll.openmc_mesh_set_dimension.restype = c_int
+_dll.openmc_mesh_set_dimension.errcheck = _error_handler
+_dll.openmc_mesh_set_params.argtypes = [
+    c_int32, c_int, POINTER(c_double), POINTER(c_double), POINTER(c_double)]
+_dll.openmc_mesh_set_params.restype = c_int
+_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
+
+
+class Mesh(_FortranObjectWithID):
+    """Mesh stored internally.
+
+    This class exposes a mesh that is stored internally in the OpenMC
+    library. To obtain a view of a mesh with a given ID, use the
+    :data:`openmc.capi.meshes` mapping.
+
+    Parameters
+    ----------
+    index : int
+         Index in the `meshes` array.
+
+    Attributes
+    ----------
+    id : int
+        ID of the mesh
+
+    """
+    __instances = WeakValueDictionary()
+
+    def __new__(cls, uid=None, new=True, index=None):
+        mapping = meshes
+        if index is None:
+            if new:
+                # Determine ID to assign
+                if uid is None:
+                    uid = max(mapping, default=0) + 1
+                else:
+                    if uid in mapping:
+                        raise AllocationError('A mesh with ID={} has already '
+                                              'been allocated.'.format(uid))
+
+                index = c_int32()
+                _dll.openmc_extend_meshes(1, index, None)
+                index = index.value
+            else:
+                index = mapping[uid]._index
+
+        if index not in cls.__instances:
+            instance = super().__new__(cls)
+            instance._index = index
+            if uid is not None:
+                instance.id = uid
+            cls.__instances[index] = instance
+
+        return cls.__instances[index]
+
+    @property
+    def id(self):
+        mesh_id = c_int32()
+        _dll.openmc_mesh_get_id(self._index, mesh_id)
+        return mesh_id.value
+
+    @id.setter
+    def id(self, mesh_id):
+        _dll.openmc_mesh_set_id(self._index, mesh_id)
+
+    @property
+    def dimension(self):
+        dims = POINTER(c_int)()
+        n = c_int()
+        _dll.openmc_mesh_get_dimension(self._index, dims, n)
+        return tuple(as_array(dims, (n.value,)))
+
+    @dimension.setter
+    def dimension(self, dimension):
+        n = len(dimension)
+        dimension = (c_int*n)(*dimension)
+        _dll.openmc_mesh_set_dimension(self._index, n, dimension)
+
+    @property
+    def lower_left(self):
+        return self._get_parameters()[0]
+
+    @property
+    def upper_right(self):
+        return self._get_parameters()[1]
+
+    @property
+    def width(self):
+        return self._get_parameters()[2]
+
+    def _get_parameters(self):
+        ll = POINTER(c_double)()
+        ur = POINTER(c_double)()
+        w = POINTER(c_double)()
+        n = c_int()
+        _dll.openmc_mesh_get_params(self._index, ll, ur, w, n)
+        return (
+            as_array(ll, (n.value,)),
+            as_array(ur, (n.value,)),
+            as_array(w, (n.value,))
+        )
+
+    def set_parameters(self, lower_left=None, upper_right=None, width=None):
+        if lower_left is not None:
+            n = len(lower_left)
+            lower_left = (c_double*n)(*lower_left)
+        if upper_right is not None:
+            n = len(upper_right)
+            upper_right = (c_double*n)(*upper_right)
+        if width is not None:
+            n = len(width)
+            width = (c_double*n)(*width)
+        _dll.openmc_mesh_set_params(self._index, n, lower_left, upper_right, width)
+
+
+class _MeshMapping(Mapping):
+    def __getitem__(self, key):
+        index = c_int32()
+        try:
+            _dll.openmc_get_mesh_index(key, index)
+        except (AllocationError, InvalidIDError) as e:
+            # __contains__ expects a KeyError to work correctly
+            raise KeyError(str(e))
+        return Mesh(index=index.value)
+
+    def __iter__(self):
+        for i in range(len(self)):
+            yield Mesh(index=i + 1).id
+
+    def __len__(self):
+        return c_int32.in_dll(_dll, 'n_meshes').value
+
+    def __repr__(self):
+        return repr(dict(self))
+
+meshes = _MeshMapping()

src/cmfd_input.F90

@@ -278,7 +278,7 @@ contains
     m % id = i_start
 
     ! Set mesh type to rectangular
-    m % type = LATTICE_RECT
+    m % type = MESH_REGULAR
 
     ! Get pointer to mesh XML node
     node_mesh = root % child("mesh")

src/mesh_header.F90

@@ -6,7 +6,7 @@ module mesh_header
 
   use constants
   use dict_header, only: DictIntInt
-  use error, only: warning, fatal_error
+  use error
   use hdf5_interface
   use string, only: to_str, to_lower
   use xml_interface
@@ -15,6 +15,13 @@ module mesh_header
   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
@@ -23,13 +30,13 @@ module mesh_header
 
   type, public :: RegularMesh
     integer :: id = -1                     ! user-specified id
-    integer :: type                        ! rectangular, hexagonal
-    integer :: n_dimension                 ! rank of mesh
+    integer :: type = MESH_REGULAR         ! rectangular, hexagonal
+    integer(C_INT) :: n_dimension          ! rank of mesh
     real(8) :: volume_frac                 ! volume fraction of each cell
-    integer, allocatable :: dimension(:)   ! number of cells in each direction
-    real(8), allocatable :: lower_left(:)  ! lower-left corner of mesh
-    real(8), allocatable :: upper_right(:) ! upper-right corner of mesh
-    real(8), allocatable :: width(:)       ! width of each mesh 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
   contains
     procedure :: from_xml => regular_from_xml
     procedure :: get_bin => regular_get_bin
@@ -592,4 +599,180 @@ contains
     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)
+          ur = C_LOC(m % upper_right)
+          width = C_LOC(m % width)
+          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/tallies/tally_filter_mesh.F90

@@ -16,6 +16,7 @@ module tally_filter_mesh
 
   implicit none
   private
+  public :: openmc_mesh_filter_get_mesh
   public :: openmc_mesh_filter_set_mesh
 
 !===============================================================================
@@ -280,35 +281,46 @@ contains
 !                               C API FUNCTIONS
 !===============================================================================
 
+  function openmc_mesh_filter_get_mesh(index, index_mesh) result(err) bind(C)
+    ! Get the mesh for a mesh filter
+    integer(C_INT32_T), value, intent(in) :: index
+    integer(C_INT32_T), intent(out)       :: index_mesh
+    integer(C_INT) :: err
+
+    err = verify_filter(index)
+    if (err == 0) then
+      select type (f => filters(index) % obj)
+      type is (MeshFilter)
+        index_mesh = f % mesh
+      class default
+        err = E_INVALID_TYPE
+        call set_errmsg("Tried to set mesh on a non-mesh filter.")
+      end select
+    end if
+  end function openmc_mesh_filter_get_mesh
+
+
   function openmc_mesh_filter_set_mesh(index, index_mesh) result(err) bind(C)
     ! Set the mesh for a mesh filter
     integer(C_INT32_T), value, intent(in) :: index
     integer(C_INT32_T), value, intent(in) :: index_mesh
     integer(C_INT) :: err
 
-    err = 0
-    if (index >= 1 .and. index <= n_filters) then
-      if (allocated(filters(index) % obj)) then
-        select type (f => filters(index) % obj)
-        type is (MeshFilter)
-          if (index_mesh >= 1 .and. index_mesh <= n_meshes) then
-            f % mesh = index_mesh
-            f % n_bins = product(meshes(index_mesh) % dimension)
-          else
-            err = E_OUT_OF_BOUNDS
-            call set_errmsg("Index in 'meshes' array is out of bounds.")
-          end if
+    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
+          f % mesh = index_mesh
+          f % n_bins = product(meshes(index_mesh) % dimension)
+        else
+          err = E_OUT_OF_BOUNDS
+          call set_errmsg("Index in 'meshes' array is out of bounds.")
+        end if
         class default
-          err = E_INVALID_TYPE
-          call set_errmsg("Tried to set mesh on a non-mesh filter.")
-        end select
-      else
-        err = E_ALLOCATE
-        call set_errmsg("Filter type has not been set yet.")
-      end if
-    else
-      err = E_OUT_OF_BOUNDS
-      call set_errmsg("Index in filters array out of bounds.")
+        err = E_INVALID_TYPE
+        call set_errmsg("Tried to set mesh on a non-mesh filter.")
+      end select
     end if
   end function openmc_mesh_filter_set_mesh
 

src/tallies/tally_filter_meshsurface.F90

@@ -16,6 +16,7 @@ module tally_filter_meshsurface
 
   implicit none
   private
+  public :: openmc_meshsurface_filter_get_mesh
   public :: openmc_meshsurface_filter_set_mesh
 
 !===============================================================================
@@ -296,6 +297,25 @@ contains
 !                               C API FUNCTIONS
 !===============================================================================
 
+  function openmc_meshsurface_filter_get_mesh(index, index_mesh) result(err) bind(C)
+    ! Get the mesh for a mesh surface filter
+    integer(C_INT32_T), value, intent(in) :: index
+    integer(C_INT32_T), intent(out)       :: index_mesh
+    integer(C_INT) :: err
+
+    err = verify_filter(index)
+    if (err == 0) then
+      select type (f => filters(index) % obj)
+      type is (MeshSurfaceFilter)
+        index_mesh = f % mesh
+      class default
+        err = E_INVALID_TYPE
+        call set_errmsg("Tried to set mesh on a non-mesh filter.")
+      end select
+    end if
+  end function openmc_meshsurface_filter_get_mesh
+
+
   function openmc_meshsurface_filter_set_mesh(index, index_mesh) result(err) bind(C)
     ! Set the mesh for a mesh surface filter
     integer(C_INT32_T), value, intent(in) :: index
@@ -304,30 +324,22 @@ contains
 
     integer :: n_dim
 
-    err = 0
-    if (index >= 1 .and. index <= n_filters) then
-      if (allocated(filters(index) % obj)) then
-        select type (f => filters(index) % obj)
-        type is (MeshSurfaceFilter)
-          if (index_mesh >= 1 .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)
-          else
-            err = E_OUT_OF_BOUNDS
-            call set_errmsg("Index in 'meshes' array is out of bounds.")
-          end if
-        class default
-          err = E_INVALID_TYPE
-          call set_errmsg("Tried to set mesh on a non-mesh filter.")
-        end select
-      else
-        err = E_ALLOCATE
-        call set_errmsg("Filter type has not been set yet.")
-      end if
-    else
-      err = E_OUT_OF_BOUNDS
-      call set_errmsg("Index in filters array out of bounds.")
+    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
+          f % mesh = index_mesh
+          n_dim = meshes(index_mesh) % n_dimension
+          f % n_bins = 4*n_dim*product(meshes(index_mesh) % dimension)
+        else
+          err = E_OUT_OF_BOUNDS
+          call set_errmsg("Index in 'meshes' array is out of bounds.")
+        end if
+      class default
+        err = E_INVALID_TYPE
+        call set_errmsg("Tried to set mesh on a non-mesh filter.")
+      end select
     end if
   end function openmc_meshsurface_filter_set_mesh
 

tests/unit_tests/test_capi.py

@@ -250,3 +250,38 @@ def test_find_material(capi_init):
     assert mat is openmc.capi.materials[1]
     mat = openmc.capi.find_material((0.4, 0., 0.))
     assert mat is openmc.capi.materials[2]
+
+
+def test_mesh(capi_init):
+    mesh = openmc.capi.Mesh()
+    mesh.dimension = (2, 3, 4)
+    assert mesh.dimension == (2, 3, 4)
+    with pytest.raises(openmc.capi.AllocationError):
+        mesh2 = openmc.capi.Mesh(mesh.id)
+
+    # Make sure each combination of parameters works
+    ll = (0., 0., 0.)
+    ur = (10., 10., 10.)
+    width = (1., 1., 1.)
+    mesh.set_parameters(lower_left=ll, upper_right=ur)
+    assert mesh.lower_left == pytest.approx(ll)
+    assert mesh.upper_right == pytest.approx(ur)
+    mesh.set_parameters(lower_left=ll, width=width)
+    assert mesh.lower_left == pytest.approx(ll)
+    assert mesh.width == pytest.approx(width)
+    mesh.set_parameters(upper_right=ur, width=width)
+    assert mesh.upper_right == pytest.approx(ur)
+    assert mesh.width == pytest.approx(width)
+
+    meshes = openmc.capi.meshes
+    assert isinstance(meshes, Mapping)
+    assert len(meshes) == 1
+    for mesh_id, mesh in meshes.items():
+        assert isinstance(mesh, openmc.capi.Mesh)
+        assert mesh_id == mesh.id
+
+    mf = openmc.capi.MeshFilter(mesh)
+    assert mf.mesh == mesh
+
+    msf = openmc.capi.MeshSurfaceFilter(mesh)
+    assert msf.mesh == mesh