diff --git a/docs/source/pythonapi/capi.rst b/docs/source/pythonapi/capi.rst index f35823399..38ce61fb3 100644 --- a/docs/source/pythonapi/capi.rst +++ b/docs/source/pythonapi/capi.rst @@ -44,5 +44,8 @@ Classes EnergyFilter MaterialFilter Material + Mesh + MeshFilter + MeshSurfaceFilter Nuclide Tally diff --git a/include/openmc.h b/include/openmc.h index 1c2106432..785aaeba5 100644 --- a/include/openmc.h +++ b/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); diff --git a/openmc/capi/__init__.py b/openmc/capi/__init__.py index bc173f994..217e782a8 100644 --- a/openmc/capi/__init__.py +++ b/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 diff --git a/openmc/capi/cell.py b/openmc/capi/cell.py index 0ab3f2583..8f0c0e1bb 100644 --- a/openmc/capi/cell.py +++ b/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 ---------- diff --git a/openmc/capi/filter.py b/openmc/capi/filter.py index 691ecc09c..a2b77fa92 100644 --- a/openmc/capi/filter.py +++ b/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' diff --git a/openmc/capi/mesh.py b/openmc/capi/mesh.py new file mode 100644 index 000000000..c7430f9f3 --- /dev/null +++ b/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() diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 86f4e2400..40b4abf57 100644 --- a/src/cmfd_input.F90 +++ b/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") diff --git a/src/mesh_header.F90 b/src/mesh_header.F90 index 8a2d61fd9..b1751a467 100644 --- a/src/mesh_header.F90 +++ b/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 diff --git a/src/tallies/tally_filter_mesh.F90 b/src/tallies/tally_filter_mesh.F90 index 3c0e6250c..d2acdc4c8 100644 --- a/src/tallies/tally_filter_mesh.F90 +++ b/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 diff --git a/src/tallies/tally_filter_meshsurface.F90 b/src/tallies/tally_filter_meshsurface.F90 index e1e8c4803..d90b63664 100644 --- a/src/tallies/tally_filter_meshsurface.F90 +++ b/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 diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index 8bc6c3d15..66a6f3bc1 100644 --- a/tests/unit_tests/test_capi.py +++ b/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