Add openmc_find function, move openmc_tally_results

openmc/libopenmc.py

@@ -7,26 +7,30 @@ from numpy.ctypeslib import ndpointer, as_array
 
 import pkg_resources
 
+_double3 = c_double*3
+
 
 class _OpenMCLibrary(object):
     def __init__(self, filename):
         self._dll = CDLL(filename)
 
         # Set argument/return types
+        self._dll.openmc_calculate_volumes.restype = None
+        self._dll.openmc_cell_set_temperature.argtypes = [
+            c_int, c_double]
+        self._dll.openmc_cell_set_temperature.restype = c_int
+        self._dll.openmc_finalize.restype = None
+        self._dll.openmc_find.argtypes = [POINTER(_double3), c_int]
+        self._dll.openmc_find.restype = c_int
         self._dll.openmc_init.argtypes = [POINTER(c_int)]
         self._dll.openmc_init.restype = None
-        self._dll.openmc_run.restype = None
         self._dll.openmc_plot_geometry.restype = None
-        self._dll.openmc_calculate_volumes.restype = None
-        self._dll.openmc_finalize.restype = None
+        self._dll.openmc_run.restype = None
         self._dll.openmc_reset.restype = None
         self._dll.openmc_tally_results.argtypes = [
             c_int, POINTER(POINTER(c_double)), ndpointer(
                 np.intc, shape=(3,))]
         self._dll.openmc_tally_results.restype = None
-        self._dll.openmc_cell_set_temperature.argtypes = [
-            c_int, c_double]
-        self._dll.openmc_cell_set_temperature.restype = c_int
 
     def init(self, intracomm=None):
         """Initialize OpenMC
@@ -94,6 +98,31 @@ class _OpenMCLibrary(object):
         """Set the temperature of a cell"""
         return self._dll.openmc_cell_set_temperature(cell_id, temperature)
 
+    def find(self, xyz, rtype='cell'):
+        """Find the cell or material at a given point
+
+        Parameters
+        ----------
+        xyz : iterable of float
+            Cartesian coordinates of position
+        rtype : {'cell', 'material'}
+            Whether to return the cell or material ID
+
+        Returns
+        -------
+        int or None
+            ID of the cell or material. If 'material' is requested and no
+            material exists at the given coordinate, None is returned.
+
+        """
+        if rtype == 'cell':
+            return self._dll.openmc_find(_double3(*xyz), 1)
+        elif rtype == 'material':
+            uid = self._dll.openmc_find(_double3(*xyz), 2)
+            return uid if uid != 0 else None
+        else:
+            raise ValueError('Unknown return type: {}'.format(rtype))
+
     def __getattr__(self, key):
         # Fall-back for other functions that may be available from library
         try:

src/capi.F90

@@ -2,13 +2,27 @@ module openmc_capi
 
   use, intrinsic :: ISO_C_BINDING
 
-  use constants, only: K_BOLTZMANN
-  use eigenvalue, only: k_sum
+  use constants,       only: K_BOLTZMANN
+  use eigenvalue,      only: k_sum
+  use finalize,        only: openmc_finalize
+  use geometry,        only: find_cell
   use global
+  use initialize,      only: openmc_init
+  use particle_header, only: Particle
+  use plot,            only: openmc_plot_geometry
+  use simulation,      only: openmc_run
+  use volume_calc,     only: openmc_calculate_volumes
 
   private
+  public :: openmc_calculate_volumes
   public :: openmc_cell_set_temperature
+  public :: openmc_finalize
+  public :: openmc_find
+  public :: openmc_init
+  public :: openmc_plot_geometry
   public :: openmc_reset
+  public :: openmc_run
+  public :: openmc_tally_results
 
 contains
 
@@ -37,6 +51,37 @@ contains
     end if
   end function openmc_cell_set_temperature
 
+!===============================================================================
+! OPENMC_FIND determines the ID or a cell or material at a given point in space
+!===============================================================================
+
+  function openmc_find(xyz, rtype) result(id) bind(C)
+    real(C_DOUBLE), intent(in)        :: xyz(3) ! Cartesian point
+    integer(C_INT), intent(in), value :: rtype  ! 1 for cell, 2 for material
+    integer(C_INT) :: id
+
+    logical :: found
+    type(Particle) :: p
+
+    call p % initialize()
+    p % coord(1) % xyz(:) = xyz
+    p % coord(1) % uvw(:) = [ZERO, ZERO, ONE]
+    call find_cell(p, found)
+
+    id = -1
+    if (found) then
+      if (rtype == 1) then
+        id = cells(p % coord(p % n_coord) % cell) % id
+      elseif (rtype == 2) then
+        if (p % material == MATERIAL_VOID) then
+          id = 0
+        else
+          id = materials(p % material) % id
+        end if
+      end if
+    end if
+  end function openmc_find
+
 !===============================================================================
 ! OPENMC_RESET resets all tallies
 !===============================================================================
@@ -62,4 +107,26 @@ contains
 
   end subroutine openmc_reset
 
+!===============================================================================
+! OPENMC_TALLY_RESULTS returns a pointer to a tally results array along with its
+! shape. This allows a user to obtain in-memory tally results from Python
+! directly.
+!===============================================================================
+
+  subroutine openmc_tally_results(i, ptr, shape_) bind(C)
+    integer(C_INT), intent(in), value :: i
+    type(C_PTR),    intent(out) :: ptr
+    integer(C_INT), intent(out) :: shape_(3)
+
+    ptr = C_NULL_PTR
+    if (allocated(tallies)) then
+      if (i >= 1 .and. i <= size(tallies)) then
+        if (allocated(tallies(i) % results)) then
+          ptr = C_LOC(tallies(i) % results(1,1,1))
+          shape_(:) = shape(tallies(i) % results)
+        end if
+      end if
+    end if
+  end subroutine openmc_tally_results
+
 end module openmc_capi

src/tally.F90

@@ -4430,20 +4430,4 @@ contains
 
   end subroutine setup_active_cmfdtallies
 
-  subroutine openmc_tally_results(i, ptr, shape_) bind(C)
-    integer(C_INT), intent(in), value :: i
-    type(C_PTR),    intent(out) :: ptr
-    integer(C_INT), intent(out) :: shape_(3)
-
-    ptr = C_NULL_PTR
-    if (allocated(tallies)) then
-      if (i >= 1 .and. i <= size(tallies)) then
-        if (allocated(tallies(i) % results)) then
-          ptr = C_LOC(tallies(i) % results(1,1,1))
-          shape_(:) = shape(tallies(i) % results)
-        end if
-      end if
-    end if
-  end subroutine openmc_tally_results
-
 end module tally