Merge branch 'develop' into cpp-init-finalize

CMakeLists.txt

@@ -373,26 +373,15 @@ add_library(libopenmc SHARED
   src/tallies/tally_derivative_header.F90
   src/tallies/tally_filter.F90
   src/tallies/tally_filter_header.F90
-  src/tallies/tally_filter_azimuthal.F90
-  src/tallies/tally_filter_cell.F90
-  src/tallies/tally_filter_cellborn.F90
-  src/tallies/tally_filter_cellfrom.F90
+  src/tallies/tally_filter_cpp.F90
   src/tallies/tally_filter_delayedgroup.F90
   src/tallies/tally_filter_distribcell.F90
   src/tallies/tally_filter_energy.F90
-  src/tallies/tally_filter_energyfunc.F90
   src/tallies/tally_filter_legendre.F90
-  src/tallies/tally_filter_material.F90
   src/tallies/tally_filter_mesh.F90
   src/tallies/tally_filter_meshsurface.F90
-  src/tallies/tally_filter_mu.F90
   src/tallies/tally_filter_particle.F90
-  src/tallies/tally_filter_polar.F90
   src/tallies/tally_filter_sph_harm.F90
-  src/tallies/tally_filter_sptl_legendre.F90
-  src/tallies/tally_filter_surface.F90
-  src/tallies/tally_filter_universe.F90
-  src/tallies/tally_filter_zernike.F90
   src/tallies/tally_header.F90
   src/tallies/trigger.F90
   src/tallies/trigger_header.F90
@@ -442,6 +431,25 @@ add_library(libopenmc SHARED
   src/string_functions.cpp
   src/summary.cpp
   src/surface.cpp
+  src/tallies/filter.cpp
+  src/tallies/filter_azimuthal.cpp
+  src/tallies/filter_cellborn.cpp
+  src/tallies/filter_cellfrom.cpp
+  src/tallies/filter_cell.cpp
+  src/tallies/filter_distribcell.cpp
+  src/tallies/filter_energyfunc.cpp
+  src/tallies/filter_legendre.cpp
+  src/tallies/filter_material.cpp
+  src/tallies/filter_mesh.cpp
+  src/tallies/filter_meshsurface.cpp
+  src/tallies/filter_mu.cpp
+  src/tallies/filter_polar.cpp
+  src/tallies/filter_sph_harm.cpp
+  src/tallies/filter_sptl_legendre.cpp
+  src/tallies/filter_surface.cpp
+  src/tallies/filter_universe.cpp
+  src/tallies/filter_zernike.cpp
+  src/tallies/tally.cpp
   src/timer.cpp
   src/thermal.cpp
   src/xml_interface.cpp

docs/source/io_formats/data_wmp.rst

@@ -4,10 +4,15 @@
 Windowed Multipole Library Format
 =================================
 
-**/version** (*char[]*)
-  The format version of the file.  The current version is "v1.0"
+**/**
+
+:Attributes: - **filetype** (*char[]*) -- String indicating the type of file
+             - **version** (*int[2]*) -- Major and minor version of the data
+
+**/<nuclide name>/**
+
+:Datasets:
 
-**/nuclide/**
     - **broaden_poly** (*int[]*)
         If 1, Doppler broaden curve fit for window with corresponding index.
         If 0, do not.

examples/jupyter/nuclear-data.ipynb

@@ -1466,7 +1466,7 @@
    },
    "outputs": [],
    "source": [
-    "url = 'https://github.com/mit-crpg/WMP_Library/releases/download/v1.0/092238.h5'\n",
+    "url = 'https://github.com/mit-crpg/WMP_Library/releases/download/v1.1/092238.h5'\n",
     "filename, headers = urllib.request.urlretrieve(url, '092238.h5')"
    ]
   },

include/openmc/capi.h

@@ -136,7 +136,6 @@ extern "C" {
   // Global variables
   extern char openmc_err_msg[256];
   extern int32_t n_cells;
-  extern int32_t n_filters;
   extern int32_t n_lattices;
   extern int32_t n_materials;
   extern int n_nuclides;

include/openmc/constants.h

@@ -62,6 +62,7 @@ constexpr int MAX_SAMPLE {100000};
 
 // Maximum number of words in a single line, length of line, and length of
 // single word
+constexpr int MAX_LINE_LEN {250};
 constexpr int MAX_WORD_LEN {150};
 
 // Maximum number of external source spatial resamples to encounter before an

include/openmc/geometry_aux.h

@@ -5,6 +5,7 @@
 #define OPENMC_GEOMETRY_AUX_H
 
 #include <cstdint>
+#include <string>
 
 
 namespace openmc {
@@ -26,7 +27,7 @@ extern "C" void assign_temperatures();
 //!
 //! This function looks for a universe that is not listed in a Cell::fill or in
 //! a Lattice.
-//! @return The index of the root universe.
+//! \return The index of the root universe.
 //==============================================================================
 
 extern "C" int32_t find_root_universe();
@@ -41,14 +42,14 @@ extern "C" void neighbor_lists();
 //! Populate all data structures needed for distribcells.
 //==============================================================================
 
-extern "C" void prepare_distribcell(int32_t* filter_cell_list, int n);
+extern "C" void prepare_distribcell();
 
 //==============================================================================
 //! Recursively search through the geometry and count cell instances.
 //!
 //! This function will update the Cell::n_instances value for each cell in the
 //! geometry.
-//! @param univ_indx The index of the universe to begin searching from (probably
+//! \param univ_indx The index of the universe to begin searching from (probably
 //!   the root universe).
 //==============================================================================
 
@@ -56,51 +57,33 @@ extern "C" void count_cell_instances(int32_t univ_indx);
 
 //==============================================================================
 //! Recursively search through universes and count universe instances.
-//! @param search_univ The index of the universe to begin searching from.
-//! @param target_univ_id The ID of the universe to be counted.
-//! @return The number of instances of target_univ_id in the geometry tree under
+//! \param search_univ The index of the universe to begin searching from.
+//! \param target_univ_id The ID of the universe to be counted.
+//! \return The number of instances of target_univ_id in the geometry tree under
 //!   search_univ.
 //==============================================================================
 
 extern "C" int
 count_universe_instances(int32_t search_univ, int32_t target_univ_id);
 
-//==============================================================================
-//! Find the length necessary for a string to contain a distribcell path.
-//! @param target_cell The index of the Cell in the global Cell array.
-//! @param map The index of the distribcell mapping corresponding to the target
-//!   cell.
-//! @param target_offset An instance number for a distributed cell.
-//! @param root_univ The index of the root Universe in the global Universe
-//!   array.
-//! @return The size of a character array needed to fit the distribcell path.
-//==============================================================================
-
-extern "C" int
-distribcell_path_len(int32_t target_cell, int32_t map, int32_t target_offset,
-                     int32_t root_univ);
-
 //==============================================================================
 //! Build a character array representing the path to a distribcell instance.
-//! @param target_cell The index of the Cell in the global Cell array.
-//! @param map The index of the distribcell mapping corresponding to the target
+//! \param target_cell The index of the Cell in the global Cell array.
+//! \param map The index of the distribcell mapping corresponding to the target
 //!   cell.
-//! @param target_offset An instance number for a distributed cell.
-//! @param root_univ The index of the root Universe in the global Universe
-//!   array.
-//! @param[out] path The unique traversal through the geometry tree that leads
-//!   to the desired instance of the target cell.
+//! \param target_offset An instance number for a distributed cell.
+//! \return The unique traversal through the geometry tree that leads to the
+//!   desired instance of the target cell.
 //==============================================================================
 
-extern "C" void
-distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset,
-                 int32_t root_univ, char *path);
+std::string
+distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset);
 
 //==============================================================================
 //! Determine the maximum number of nested coordinate levels in the geometry.
-//! @param univ The index of the universe to begin seraching from (probably the
+//! \param univ The index of the universe to begin seraching from (probably the
 //!   root universe).
-//! @return The number of coordinate levels.
+//! \return The number of coordinate levels.
 //==============================================================================
 
 extern "C" int maximum_levels(int32_t univ);

include/openmc/hdf5_interface.h

@@ -380,5 +380,11 @@ write_dataset(hid_t obj_id, const char* name, Position r)
   write_dataset(obj_id, name, buffer);
 }
 
+inline void
+write_dataset(hid_t obj_id, const char* name, std::string buffer)
+{
+  write_string(obj_id, name, buffer.c_str(), false);
+}
+
 } // namespace openmc
 #endif // OPENMC_HDF5_INTERFACE_H

include/openmc/material.h

@@ -24,7 +24,7 @@ extern std::unordered_map<int32_t, int32_t> material_map;
 class Material
 {
 public:
-  int32_t id; //!< Unique ID
+  int32_t id_; //!< Unique ID
   double volume_ {-1.0}; //!< Volume in [cm^3]
   bool fissionable {false}; //!< Does this material contain fissionable nuclides
 

include/openmc/math_functions.h

@@ -90,7 +90,7 @@ extern "C" void calc_rn_c(int n, const double uvw[3], double rn[]);
 //!   evaluated at rho and phi.
 //==============================================================================
 
-extern "C" void calc_zn_c(int n, double rho, double phi, double zn[]);
+extern "C" void calc_zn(int n, double rho, double phi, double zn[]);
 
 //==============================================================================
 //! Calculate only the even radial components of n-th order modified Zernike
@@ -113,7 +113,7 @@ extern "C" void calc_zn_c(int n, double rho, double phi, double zn[]);
 //!   evaluated at rho and phi when m = 0.
 //==============================================================================
 
-extern "C" void calc_zn_rad_c(int n, double rho, double zn_rad[]);
+extern "C" void calc_zn_rad(int n, double rho, double zn_rad[]);
 
 //==============================================================================
 //! Rotate the direction cosines through a polar angle whose cosine is mu and

include/openmc/tallies/filter.h

@@ -0,0 +1,91 @@
+#ifndef OPENMC_TALLIES_FILTER_H
+#define OPENMC_TALLIES_FILTER_H
+
+#include <cstdint>
+#include <string>
+#include <vector>
+
+#include "openmc/hdf5_interface.h"
+#include "openmc/particle.h"
+#include "pugixml.hpp"
+
+
+namespace openmc {
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+extern "C" int32_t n_filters;
+
+class FilterMatch;
+extern std::vector<FilterMatch> filter_matches;
+#pragma omp threadprivate(filter_matches)
+
+class Filter;
+extern std::vector<Filter*> tally_filters;
+
+//==============================================================================
+//! Stores bins and weights for filtered tally events.
+//==============================================================================
+
+class FilterMatch
+{
+public:
+  std::vector<int> bins_;
+  std::vector<double> weights_;
+};
+
+//==============================================================================
+//! Modifies tally score events.
+//==============================================================================
+
+class Filter
+{
+public:
+  virtual ~Filter() = default;
+
+  virtual std::string type() const = 0;
+
+  //! Uses an XML input to fill the filter's data fields.
+  virtual void from_xml(pugi::xml_node node) = 0;
+
+  //! Matches a tally event to a set of filter bins and weights.
+  //!
+  //! \param[out] match will contain the matching bins and corresponding
+  //!   weights; note that there may be zero matching bins
+  virtual void
+  get_all_bins(const Particle* p, int estimator, FilterMatch& match) const = 0;
+
+  //! Writes data describing this filter to an HDF5 statepoint group.
+  virtual void
+  to_statepoint(hid_t filter_group) const
+  {
+    write_dataset(filter_group, "type", type());
+    write_dataset(filter_group, "n_bins", n_bins_);
+  }
+
+  //! Return a string describing a filter bin for the tallies.out file.
+  //
+  //! For example, an `EnergyFilter` might return the string
+  //! "Incoming Energy [0.625E-6, 20.0)".
+  virtual std::string text_label(int bin) const = 0;
+
+  virtual void initialize() {}
+
+  int n_bins_;
+};
+
+//==============================================================================
+
+extern "C" void free_memory_tally_c();
+
+//==============================================================================
+
+// Filter-related Fortran functions that will be called from C++
+extern "C" int verify_filter(int32_t index);
+extern "C" Filter* filter_from_f(int32_t index);
+extern "C" void filter_update_n_bins(int32_t index);
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_H

include/openmc/tallies/filter_azimuthal.h

@@ -0,0 +1,35 @@
+#ifndef OPENMC_TALLIES_FILTER_AZIMUTHAL_H
+#define OPENMC_TALLIES_FILTER_AZIMUTHAL_H
+
+#include <string>
+#include <vector>
+
+#include "openmc/tallies/filter.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Bins the incident neutron azimuthal angle (relative to the global xy-plane).
+//==============================================================================
+
+class AzimuthalFilter : public Filter
+{
+public:
+  ~AzimuthalFilter() = default;
+
+  std::string type() const override {return "azimuthal";}
+
+  void from_xml(pugi::xml_node node) override;
+
+  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  std::vector<double> bins_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_AZIMUTHAL_H

include/openmc/tallies/filter_cell.h

@@ -0,0 +1,42 @@
+#ifndef OPENMC_TALLIES_FILTER_CELL_H
+#define OPENMC_TALLIES_FILTER_CELL_H
+
+#include <cstdint>
+#include <unordered_map>
+#include <vector>
+
+#include "openmc/tallies/filter.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Specifies which geometric cells tally events reside in.
+//==============================================================================
+
+class CellFilter : public Filter
+{
+public:
+  ~CellFilter() = default;
+
+  std::string type() const override {return "cell";}
+
+  void from_xml(pugi::xml_node node) override;
+
+  void initialize() override;
+
+  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  //! The indices of the cells binned by this filter.
+  std::vector<int32_t> cells_;
+
+  //! A map from cell indices to filter bin indices.
+  std::unordered_map<int32_t, int> map_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_CELL_H

include/openmc/tallies/filter_cellborn.h

@@ -0,0 +1,26 @@
+#ifndef OPENMC_TALLIES_FILTER_CELLBORN_H
+#define OPENMC_TALLIES_FILTER_CELLBORN_H
+
+#include <string>
+
+#include "openmc/tallies/filter_cell.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Specifies which cell the particle was born in.
+//==============================================================================
+
+class CellbornFilter : public CellFilter
+{
+public:
+  std::string type() const override {return "cellborn";}
+
+  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  const override;
+
+  std::string text_label(int bin) const override;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_CELLBORN_H

include/openmc/tallies/filter_cellfrom.h

@@ -0,0 +1,26 @@
+#ifndef OPENMC_TALLIES_FILTER_CELLFROM_H
+#define OPENMC_TALLIES_FILTER_CELLFROM_H
+
+#include <string>
+
+#include "openmc/tallies/filter_cell.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Specifies which geometric cells particles exit when crossing a surface.
+//==============================================================================
+
+class CellFromFilter : public CellFilter
+{
+public:
+  std::string type() const override {return "cellfrom";}
+
+  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  const override;
+
+  std::string text_label(int bin) const override;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_CELLFROM_H

include/openmc/tallies/filter_distribcell.h

@@ -0,0 +1,36 @@
+#ifndef OPENMC_TALLIES_FILTER_DISTRIBCELL_H
+#define OPENMC_TALLIES_FILTER_DISTRIBCELL_H
+
+#include <string>
+
+#include "openmc/tallies/filter.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Specifies which distributed geometric cells tally events reside in.
+//==============================================================================
+
+class DistribcellFilter : public Filter
+{
+public:
+  ~DistribcellFilter() = default;
+
+  std::string type() const override {return "distribcell";}
+
+  void from_xml(pugi::xml_node node) override;
+
+  void initialize() override;
+
+  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  int32_t cell_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_DISTRIBCELL_H

include/openmc/tallies/filter_energyfunc.h

@@ -0,0 +1,45 @@
+#ifndef OPENMC_TALLIES_FILTER_ENERGYFUNC_H
+#define OPENMC_TALLIES_FILTER_ENERGYFUNC_H
+
+#include <vector>
+
+#include "openmc/tallies/filter.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Multiplies tally scores by an arbitrary function of incident energy
+//! described by a piecewise linear-linear interpolation.
+//==============================================================================
+
+class EnergyFunctionFilter : public Filter
+{
+public:
+  EnergyFunctionFilter()
+    : Filter {}
+  {
+    n_bins_ = 1;
+  }
+
+  ~EnergyFunctionFilter() = default;
+
+  std::string type() const override {return "energyfunction";}
+
+  void from_xml(pugi::xml_node node) override;
+
+  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  //! Incident neutron energy interpolation grid.
+  std::vector<double> energy_;
+
+  //! Interpolant values.
+  std::vector<double> y_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_ENERGYFUNC_H

include/openmc/tallies/filter_legendre.h

@@ -0,0 +1,34 @@
+#ifndef OPENMC_TALLIES_FILTER_LEGENDRE_H
+#define OPENMC_TALLIES_FILTER_LEGENDRE_H
+
+#include <string>
+
+#include "openmc/tallies/filter.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Gives Legendre moments of the change in scattering angle
+//==============================================================================
+
+class LegendreFilter : public Filter
+{
+public:
+  ~LegendreFilter() = default;
+
+  std::string type() const override {return "legendre";}
+
+  void from_xml(pugi::xml_node node) override;
+
+  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  int order_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_LEGENDRE_H

include/openmc/tallies/filter_material.h

@@ -0,0 +1,42 @@
+#ifndef OPENMC_TALLIES_FILTER_MATERIAL_H
+#define OPENMC_TALLIES_FILTER_MATERIAL_H
+
+#include <cstdint>
+#include <unordered_map>
+#include <vector>
+
+#include "openmc/tallies/filter.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Specifies which material tally events reside in.
+//==============================================================================
+
+class MaterialFilter : public Filter
+{
+public:
+  ~MaterialFilter() = default;
+
+  std::string type() const override {return "material";}
+
+  void from_xml(pugi::xml_node node) override;
+
+  void initialize() override;
+
+  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  //! The indices of the materials binned by this filter.
+  std::vector<int32_t> materials_;
+
+  //! A map from material indices to filter bin indices.
+  std::unordered_map<int32_t, int> map_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_MATERIAL_H

include/openmc/tallies/filter_mesh.h

@@ -0,0 +1,41 @@
+#ifndef OPENMC_TALLIES_FILTER_MESH_H
+#define OPENMC_TALLIES_FILTER_MESH_H
+
+#include <cstdint>
+
+#include "openmc/tallies/filter.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Indexes the location of particle events to a regular mesh.  For tracklength
+//! tallies, it will produce multiple valid bins and the bin weight will
+//! correspond to the fraction of the track length that lies in that bin.
+//==============================================================================
+
+class MeshFilter : public Filter
+{
+public:
+  ~MeshFilter() = default;
+
+  std::string type() const override {return "mesh";}
+
+  void from_xml(pugi::xml_node node) override;
+
+  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  virtual int32_t mesh() const {return mesh_;}
+
+  virtual void set_mesh(int32_t mesh);
+
+protected:
+  int32_t mesh_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_MESH_H

include/openmc/tallies/filter_meshsurface.h

@@ -0,0 +1,22 @@
+#ifndef OPENMC_TALLIES_FILTER_MESHSURFACE_H
+#define OPENMC_TALLIES_FILTER_MESHSURFACE_H
+
+#include "openmc/tallies/filter_mesh.h"
+
+namespace openmc {
+
+class MeshSurfaceFilter : public MeshFilter
+{
+public:
+  std::string type() const override {return "meshsurface";}
+
+  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  const override;
+
+  std::string text_label(int bin) const override;
+
+  void set_mesh(int32_t mesh) override;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_MESHSURFACE_H

include/openmc/tallies/filter_mu.h

@@ -0,0 +1,35 @@
+#ifndef OPENMC_TALLIES_FILTER_MU_H
+#define OPENMC_TALLIES_FILTER_MU_H
+
+#include <vector>
+
+#include "openmc/tallies/filter.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Bins the incoming-outgoing direction cosine.  This is only used for scatter
+//! reactions.
+//==============================================================================
+
+class MuFilter : public Filter
+{
+public:
+  ~MuFilter() = default;
+
+  std::string type() const override {return "mu";}
+
+  void from_xml(pugi::xml_node node) override;
+
+  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  std::vector<double> bins_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_MU_H

include/openmc/tallies/filter_polar.h

@@ -0,0 +1,35 @@
+#ifndef OPENMC_TALLIES_FILTER_POLAR_H
+#define OPENMC_TALLIES_FILTER_POLAR_H
+
+#include <cmath>
+#include <vector>
+
+#include "openmc/tallies/filter.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Bins the incident neutron polar angle (relative to the global z-axis).
+//==============================================================================
+
+class PolarFilter : public Filter
+{
+public:
+  ~PolarFilter() = default;
+
+  std::string type() const override {return "polar";}
+
+  void from_xml(pugi::xml_node node) override;
+
+  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  std::vector<double> bins_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_POLAR_H

include/openmc/tallies/filter_sph_harm.h

@@ -0,0 +1,42 @@
+#ifndef OPENMC_TALLIES_FILTER_SPH_HARM_H
+#define OPENMC_TALLIES_FILTER_SPH_HARM_H
+
+#include <string>
+
+#include "openmc/tallies/filter.h"
+
+namespace openmc {
+
+//TODO: those integer values are not needed when Fortran interop is removed
+enum class SphericalHarmonicsCosine {
+  scatter = 1, particle = 2
+};
+
+//==============================================================================
+//! Gives spherical harmonics expansion moments of a tally score
+//==============================================================================
+
+class SphericalHarmonicsFilter : public Filter
+{
+public:
+  ~SphericalHarmonicsFilter() = default;
+
+  std::string type() const override {return "sphericalharmonics";}
+
+  void from_xml(pugi::xml_node node) override;
+
+  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  int order_;
+
+  //! The type of angle that this filter measures when binning events.
+  SphericalHarmonicsCosine cosine_ {SphericalHarmonicsCosine::particle};
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_SPH_HARM_H

include/openmc/tallies/filter_sptl_legendre.h

@@ -0,0 +1,48 @@
+#ifndef OPENMC_TALLIES_FILTER_SPTL_LEGENDRE_H
+#define OPENMC_TALLIES_FILTER_SPTL_LEGENDRE_H
+
+#include <string>
+
+#include "openmc/tallies/filter.h"
+
+namespace openmc {
+
+//TODO: those integer values are not needed when Fortran interop is removed
+enum class LegendreAxis {
+  x = 1, y = 2, z = 3
+};
+
+//==============================================================================
+//! Gives Legendre moments of the particle's normalized position along an axis
+//==============================================================================
+
+class SpatialLegendreFilter : public Filter
+{
+public:
+  ~SpatialLegendreFilter() = default;
+
+  std::string type() const override {return "spatiallegendre";}
+
+  void from_xml(pugi::xml_node node) override;
+
+  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  int order_;
+
+  //! The Cartesian coordinate axis that the Legendre expansion is applied to.
+  LegendreAxis axis_;
+
+  //! The minimum coordinate along the reference axis that the expansion covers.
+  double min_;
+
+  //! The maximum coordinate along the reference axis that the expansion covers.
+  double max_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_SPTL_LEGENDRE_H

include/openmc/tallies/filter_surface.h

@@ -0,0 +1,42 @@
+#ifndef OPENMC_TALLIES_FILTER_SURFACE_H
+#define OPENMC_TALLIES_FILTER_SURFACE_H
+
+#include <cstdint>
+#include <unordered_map>
+#include <vector>
+
+#include "openmc/tallies/filter.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Specifies which surface particles are crossing
+//==============================================================================
+
+class SurfaceFilter : public Filter
+{
+public:
+  ~SurfaceFilter() = default;
+
+  std::string type() const override {return "surface";}
+
+  void from_xml(pugi::xml_node node) override;
+
+  void initialize() override;
+
+  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  //! The indices of the surfaces binned by this filter.
+  std::vector<int32_t> surfaces_;
+
+  //! A map from surface indices to filter bin indices.
+  std::unordered_map<int32_t, int> map_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_SURFACE_H

include/openmc/tallies/filter_universe.h

@@ -0,0 +1,42 @@
+#ifndef OPENMC_TALLIES_FILTER_UNIVERSE_H
+#define OPENMC_TALLIES_FILTER_UNIVERSE_H
+
+#include <cstdint>
+#include <unordered_map>
+#include <vector>
+
+#include "openmc/tallies/filter.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Specifies which geometric universes tally events reside in.
+//==============================================================================
+
+class UniverseFilter : public Filter
+{
+public:
+  ~UniverseFilter() = default;
+
+  std::string type() const override {return "universe";}
+
+  void from_xml(pugi::xml_node node) override;
+
+  void initialize() override;
+
+  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  //! The indices of the universes binned by this filter.
+  std::vector<int32_t> universes_;
+
+  //! A map from universe indices to filter bin indices.
+  std::unordered_map<int32_t, int> map_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_UNIVERSE_H

include/openmc/tallies/filter_zernike.h

@@ -0,0 +1,65 @@
+#ifndef OPENMC_TALLIES_FILTER_ZERNIKE_H
+#define OPENMC_TALLIES_FILTER_ZERNIKE_H
+
+#include <string>
+
+#include "openmc/tallies/filter.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Gives Zernike polynomial moments of a particle's position
+//==============================================================================
+
+class ZernikeFilter : public Filter
+{
+public:
+  std::string type() const override {return "zernike";}
+
+  ~ZernikeFilter() = default;
+
+  void from_xml(pugi::xml_node node) override;
+
+  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  const override;
+
+  void to_statepoint(hid_t filter_group) const override;
+
+  std::string text_label(int bin) const override;
+
+  int order() const {return order_;}
+
+  virtual void set_order(int order);
+
+  //! Cartesian x coordinate for the origin of this expansion.
+  double x_;
+
+  //! Cartesian y coordinate for the origin of this expansion.
+  double y_;
+
+  //! Maximum radius from the origin covered by this expansion.
+  double r_;
+
+protected:
+  int order_;
+};
+
+//==============================================================================
+//! Gives even order radial Zernike polynomial moments of a particle's position
+//==============================================================================
+
+class ZernikeRadialFilter : public ZernikeFilter
+{
+public:
+  std::string type() const override {return "zernikeradial";}
+
+  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  const override;
+
+  std::string text_label(int bin) const override;
+
+  void set_order(int order) override;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_ZERNIKE_H

openmc/capi/filter.py

@@ -238,13 +238,15 @@ class MaterialFilter(Filter):
         materials = POINTER(c_int32)()
         n = c_int32()
         _dll.openmc_material_filter_get_bins(self._index, materials, n)
-        return [Material(index=materials[i]) for i in range(n.value)]
+        #TODO: fix this off-by-one when materials become 0-indexed
+        return [Material(index=materials[i]+1) for i in range(n.value)]
 
     @bins.setter
     def bins(self, materials):
         # Get material indices as int32_t[]
         n = len(materials)
-        bins = (c_int32*n)(*(m._index for m in materials))
+        #TODO: fix this off-by-one when materials become 0-indexed
+        bins = (c_int32*n)(*(m._index-1 for m in materials))
 
         _dll.openmc_material_filter_set_bins(self._index, n, bins)
 

openmc/capi/math.py

@@ -18,11 +18,11 @@ _dll.evaluate_legendre.argtypes = [c_int, POINTER(c_double), c_double]
 _dll.calc_rn_c.restype = None
 _dll.calc_rn_c.argtypes = [c_int, ndpointer(c_double), ndpointer(c_double)]
 
-_dll.calc_zn_c.restype = None
-_dll.calc_zn_c.argtypes = [c_int, c_double, c_double, ndpointer(c_double)]
+_dll.calc_zn.restype = None
+_dll.calc_zn.argtypes = [c_int, c_double, c_double, ndpointer(c_double)]
 
-_dll.calc_zn_rad_c.restype = None
-_dll.calc_zn_rad_c.argtypes = [c_int, c_double, ndpointer(c_double)]
+_dll.calc_zn_rad.restype = None
+_dll.calc_zn_rad.argtypes = [c_int, c_double, ndpointer(c_double)]
 
 _dll.rotate_angle_c.restype = None
 _dll.rotate_angle_c.argtypes = [ndpointer(c_double), c_double,
@@ -153,7 +153,7 @@ def calc_zn(n, rho, phi):
 
     num_bins = ((n + 1) * (n + 2)) // 2
     zn = np.zeros(num_bins, dtype=np.float64)
-    _dll.calc_zn_c(n, rho, phi, zn)
+    _dll.calc_zn(n, rho, phi, zn)
     return zn
 
 
@@ -179,7 +179,7 @@ def calc_zn_rad(n, rho):
 
     num_bins = n // 2 + 1
     zn_rad = np.zeros(num_bins, dtype=np.float64)
-    _dll.calc_zn_rad_c(n, rho, zn_rad)
+    _dll.calc_zn_rad(n, rho, zn_rad)
     return zn_rad
 
 

openmc/data/__init__.py

@@ -4,7 +4,9 @@ HDF5_VERSION_MINOR = 0
 HDF5_VERSION = (HDF5_VERSION_MAJOR, HDF5_VERSION_MINOR)
 
 # Version of WMP nuclear data format
-WMP_VERSION = 'v1.0'
+WMP_VERSION_MAJOR = 1
+WMP_VERSION_MINOR = 1
+WMP_VERSION = (WMP_VERSION_MAJOR, WMP_VERSION_MINOR)
 
 
 from .data import *

openmc/data/multipole.py

@@ -4,7 +4,7 @@ from math import exp, erf, pi, sqrt
 import h5py
 import numpy as np
 
-from . import WMP_VERSION
+from . import WMP_VERSION, WMP_VERSION_MAJOR
 from .data import K_BOLTZMANN
 import openmc.checkvalue as cv
 from openmc.mixin import EqualityMixin
@@ -133,6 +133,8 @@ class WindowedMultipole(EqualityMixin):
 
     Parameters
     ----------
+    name : str
+        Name of the nuclide using the GND naming convention
 
     Attributes
     ----------
@@ -171,7 +173,8 @@ class WindowedMultipole(EqualityMixin):
         a/E + b/sqrt(E) + c + d sqrt(E) + ...
 
     """
-    def __init__(self):
+    def __init__(self, name):
+        self.name = name
         self.spacing = None
         self.sqrtAWR = None
         self.E_min = None
@@ -181,6 +184,10 @@ class WindowedMultipole(EqualityMixin):
         self.broaden_poly = None
         self.curvefit = None
 
+    @property
+    def name(self):
+        return self._name
+
     @property
     def fit_order(self):
         return self.curvefit.shape[1] - 1
@@ -221,6 +228,11 @@ class WindowedMultipole(EqualityMixin):
     def curvefit(self):
         return self._curvefit
 
+    @name.setter
+    def name(self, name):
+        cv.check_type('name', name, str)
+        self._name = name
+
     @spacing.setter
     def spacing(self, spacing):
         if spacing is not None:
@@ -322,17 +334,25 @@ class WindowedMultipole(EqualityMixin):
             group = group_or_filename
         else:
             h5file = h5py.File(group_or_filename, 'r')
-            try:
-                version = h5file['version'].value.decode()
-            except AttributeError:
-                version = h5file['version'].value[0].decode()
-            if version != WMP_VERSION:
-                raise ValueError('The given WMP data uses version '
-                    + version + ' whereas your installation of the OpenMC '
-                    'Python API expects version ' + WMP_VERSION)
-            group = h5file['nuclide']
 
-        out = cls()
+            # Make sure version matches
+            if 'version' in h5file.attrs:
+                major, minor = h5file.attrs['version']
+                if major != WMP_VERSION_MAJOR:
+                    raise IOError(
+                        'WMP data format uses version {}. {} whereas your '
+                        'installation of the OpenMC Python API expects version '
+                        '{}.x.'.format(major, minor, WMP_VERSION_MAJOR))
+            else:
+                raise IOError(
+                    'WMP data does not indicate a version. Your installation of '
+                    'the OpenMC Python API expects version {}.x data.'
+                    .format(WMP_VERSION_MAJOR))
+
+            group = list(h5file.values())[0]
+
+        name = group.name[1:]
+        out = cls(name)
 
         # Read scalars.
 
@@ -478,13 +498,16 @@ class WindowedMultipole(EqualityMixin):
         fun = np.vectorize(lambda x: self._evaluate(x, T))
         return fun(E)
 
-    def export_to_hdf5(self, path, libver='earliest'):
+    def export_to_hdf5(self, path, mode='a', libver='earliest'):
         """Export windowed multipole data to an HDF5 file.
 
         Parameters
         ----------
         path : str
             Path to write HDF5 file to
+        mode : {'r', r+', 'w', 'x', 'a'}
+            Mode that is used to open the HDF5 file. This is the second argument
+            to the :class:`h5py.File` constructor.
         libver : {'earliest', 'latest'}
             Compatibility mode for the HDF5 file. 'latest' will produce files
             that are less backwards compatible but have performance benefits.
@@ -492,12 +515,11 @@ class WindowedMultipole(EqualityMixin):
         """
 
         # Open file and write version.
-        with h5py.File(path, 'w', libver=libver) as f:
-            f.create_dataset('version', (1, ), dtype='S10')
-            f['version'][:] = WMP_VERSION.encode('ASCII')
+        with h5py.File(path, mode, libver=libver) as f:
+            f.attrs['filetype'] = np.string_('data_wmp')
+            f.attrs['version'] = np.array(WMP_VERSION)
 
-            # Make a nuclide group.
-            g = f.create_group('nuclide')
+            g = f.create_group(self.name)
 
             # Write scalars.
             g.create_dataset('spacing', data=np.array(self.spacing))

openmc/polynomial.py

@@ -1,6 +1,5 @@
 import numpy as np
 import openmc
-import openmc.capi as capi
 from collections.abc import Iterable
 
 
@@ -74,6 +73,7 @@ class ZernikeRadial(Polynomial):
         return self._order
 
     def __call__(self, r):
+        import openmc.capi as capi
         if isinstance(r, Iterable):
             return [np.sum(self._norm_coef * capi.calc_zn_rad(self.order, r_i))
                     for r_i in r]

scripts/openmc-get-multipole-data

@@ -29,9 +29,9 @@ parser.add_argument('-b', '--batch', action='store_true',
 args = parser.parse_args()
 
 
-baseUrl = 'https://github.com/mit-crpg/WMP_Library/releases/download/v1.0/'
-files = ['WMP_Library_v1.0.tar.gz']
-checksums = ['22cb675734cfccb278dffd40dcfbf26a']
+baseUrl = 'https://github.com/mit-crpg/WMP_Library/releases/download/v1.1/'
+files = ['WMP_Library_v1.1.tar.gz']
+checksums = ['8523895928dd6ba63fba803e3a45d4f3']
 block_size = 16384
 
 # ==============================================================================

scripts/openmc-plot-mesh-tally

@@ -11,8 +11,10 @@ import tkinter.font as font
 import tkinter.messagebox as messagebox
 import tkinter.ttk as ttk
 
+import matplotlib
+matplotlib.use("TkAgg")
 from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
-from matplotlib.backends.backend_tkagg import NavigationToolbar2TkAgg
+from matplotlib.backends.backend_tkagg import NavigationToolbar2Tk
 from matplotlib.figure import Figure
 import matplotlib.pyplot as plt
 import numpy as np
@@ -56,7 +58,7 @@ class MeshPlotter(tk.Frame):
         self.canvas.get_tk_widget().pack(side=tk.TOP, fill=tk.BOTH, expand=1)
 
         # Create the navigation toolbar, tied to the canvas
-        self.mpl_toolbar = NavigationToolbar2TkAgg(self.canvas, figureFrame)
+        self.mpl_toolbar = NavigationToolbar2Tk(self.canvas, figureFrame)
         self.mpl_toolbar.update()
         self.canvas._tkcanvas.pack(side=tk.TOP, fill=tk.BOTH, expand=1)
 

src/cell.cpp

@@ -476,7 +476,7 @@ CSGCell::to_hdf5(hid_t cell_group) const
     std::vector<int32_t> mat_ids;
     for (auto i_mat : material_) {
       if (i_mat != MATERIAL_VOID) {
-        mat_ids.push_back(materials[i_mat]->id);
+        mat_ids.push_back(materials[i_mat]->id_);
       } else {
         mat_ids.push_back(MATERIAL_VOID);
       }
@@ -649,6 +649,19 @@ read_cells(pugi::xml_node* node)
     cells.push_back(new CSGCell(cell_node));
   }
 
+  // Fill the cell map.
+  for (int i = 0; i < cells.size(); i++) {
+    int32_t id = cells[i]->id_;
+    auto search = cell_map.find(id);
+    if (search == cell_map.end()) {
+      cell_map[id] = i;
+    } else {
+      std::stringstream err_msg;
+      err_msg << "Two or more cells use the same unique ID: " << id;
+      fatal_error(err_msg);
+    }
+  }
+
   // Populate the Universe vector and map.
   for (int i = 0; i < cells.size(); i++) {
     int32_t uid = cells[i]->universe_;
@@ -766,7 +779,19 @@ extern "C" {
 
   int32_t cell_id(Cell* c) {return c->id_;}
 
-  void cell_set_id(Cell* c, int32_t id) {c->id_ = id;}
+  void
+  cell_set_id(Cell* c, int32_t id)
+  {
+    c->id_ = id;
+
+    // Find the index of this cell and update the cell map.
+    for (int i = 0; i < cells.size(); i++) {
+      if (cells[i] == c) {
+        cell_map[id] = i;
+        break;
+      }
+    }
+  }
 
   int cell_type(Cell* c) {return c->type_;}
 
@@ -789,10 +814,6 @@ extern "C" {
 
   int32_t cell_fill(Cell* c) {return c->fill_;}
 
-  int32_t cell_n_instances(Cell* c) {return c->n_instances_;}
-
-  int cell_distribcell_index(Cell* c) {return c->distribcell_index_;}
-
   int cell_material_size(Cell* c) {return c->material_.size();}
 
   //TODO: off-by-one
@@ -807,8 +828,6 @@ extern "C" {
 
   double cell_sqrtkT(Cell* c, int i) {return c->sqrtkT_[i];}
 
-  int32_t cell_offset(Cell* c, int map) {return c->offset_[map];}
-
   void extend_cells_c(int32_t n)
   {
     cells.reserve(cells.size() + n);

src/cmfd_data.F90

@@ -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
@@ -154,27 +154,27 @@ contains
 
                 ! Reset all bins to 1
                 do l = 1, size(t % filter)
-                  call filter_matches(t % filter(l)) % bins % clear()
-                  call filter_matches(t % filter(l)) % bins % push_back(1)
+                  call filter_matches(t % filter(l)) % bins_clear()
+                  call filter_matches(t % filter(l)) % bins_push_back(1)
                 end do
 
                 ! Set ijk as mesh indices
                 ijk = (/ i, j, k /)
 
                 ! Get bin number for mesh indices
-                filter_matches(i_filter_mesh) % bins % data(1) = &
-                     m % get_bin_from_indices(ijk)
+                call filter_matches(i_filter_mesh) &
+                     % bins_set_data(1, m % get_bin_from_indices(ijk))
 
                 ! Apply energy in filter
                 if (energy_filters) then
-                  filter_matches(i_filter_ein) % bins % data(1) = ng - h + 1
+                  call filter_matches(i_filter_ein) % bins_set_data(1, ng - h + 1)
                 end if
 
                 ! Calculate score index from bins
                 score_index = 1
                 do l = 1, size(t % filter)
                   score_index = score_index + (filter_matches(t % filter(l)) &
-                       % bins % data(1) - 1) * t % stride(l)
+                       % bins_data(1) - 1) * t % stride(l)
                 end do
 
                 ! Get flux
@@ -198,30 +198,30 @@ contains
 
                   ! Reset all bins to 1
                   do l = 1, size(t % filter)
-                    call filter_matches(t % filter(l)) % bins % clear()
-                    call filter_matches(t % filter(l)) % bins % push_back(1)
+                    call filter_matches(t % filter(l)) % bins_clear()
+                    call filter_matches(t % filter(l)) % bins_push_back(1)
                   end do
 
                   ! Set ijk as mesh indices
                   ijk = (/ i, j, k /)
 
                   ! Get bin number for mesh indices
-                  filter_matches(i_filter_mesh) % bins % data(1) = &
-                       m % get_bin_from_indices(ijk)
+                  call filter_matches(i_filter_mesh) &
+                       % bins_set_data(1, m % get_bin_from_indices(ijk))
 
                   if (energy_filters) then
                     ! Apply energy in filter
-                    filter_matches(i_filter_ein) % bins % data(1) = ng - h + 1
+                    call filter_matches(i_filter_ein) % bins_set_data(1, ng - h + 1)
 
                     ! Set energy out bin
-                    filter_matches(i_filter_eout) % bins % data(1) = ng - g + 1
+                    call filter_matches(i_filter_eout) % bins_set_data(1, ng - g + 1)
                   end if
 
                   ! Calculate score index from bins
                   score_index = 1
                   do l = 1, size(t % filter)
                     score_index = score_index + (filter_matches(t % filter(l)) &
-                         % bins % data(1) - 1) * t % stride(l)
+                         % bins_data(1) - 1) * t % stride(l)
                   end do
 
                   ! Get scattering
@@ -244,23 +244,23 @@ contains
 
                 ! Initialize and filter for energy
                 do l = 1, size(t % filter)
-                  call filter_matches(t % filter(l)) % bins % clear()
-                  call filter_matches(t % filter(l)) % bins % push_back(1)
+                  call filter_matches(t % filter(l)) % bins_clear()
+                  call filter_matches(t % filter(l)) % bins_push_back(1)
                 end do
 
                 ! Set the bin for this mesh cell
                 i_mesh = m % get_bin_from_indices([ i, j, k ])
-                filter_matches(i_filter_mesh) % bins % data(1) = 12*(i_mesh - 1) + 1
+                call filter_matches(i_filter_mesh) % bins_set_data(1, 12*(i_mesh - 1) + 1)
 
                 ! Set the energy bin if needed
                 if (energy_filters) then
-                  filter_matches(i_filter_ein) % bins % data(1) = ng - h + 1
+                  call filter_matches(i_filter_ein) % bins_set_data(1, ng - h + 1)
                 end if
 
                 score_index = 0
                 do l = 1, size(t % filter)
                   score_index = score_index + (filter_matches(t % filter(l)) &
-                       % bins % data(1) - 1) * t % stride(l)
+                       % bins_data(1) - 1) * t % stride(l)
                 end do
 
                 ! Left surface
@@ -303,30 +303,30 @@ contains
 
                 ! Reset all bins to 1
                 do l = 1, size(t % filter)
-                  call filter_matches(t % filter(l)) % bins % clear()
-                  call filter_matches(t % filter(l)) % bins % push_back(1)
+                  call filter_matches(t % filter(l)) % bins_clear()
+                  call filter_matches(t % filter(l)) % bins_push_back(1)
                 end do
 
                 ! Set ijk as mesh indices
                 ijk = (/ i, j, k /)
 
                 ! Get bin number for mesh indices
-                filter_matches(i_filter_mesh) % bins % data(1) = &
-                     m % get_bin_from_indices(ijk)
+                call filter_matches(i_filter_mesh) &
+                     % bins_set_data(1, m % get_bin_from_indices(ijk))
 
                 ! Apply energy in filter
                 if (energy_filters) then
-                  filter_matches(i_filter_ein) % bins % data(1) = ng - h + 1
+                  call filter_matches(i_filter_ein) % bins_set_data(1, ng - h + 1)
                 end if
 
                 ! Apply Legendre filter
-                filter_matches(i_filter_legendre) % bins % data(1) = 2
+                call filter_matches(i_filter_legendre) % bins_set_data(1, 2)
 
                 ! Calculate score index from bins
                 score_index = 1
                 do l = 1, size(t % filter)
                   score_index = score_index + (filter_matches(t % filter(l)) &
-                       % bins % data(1) - 1) * t % stride(l)
+                       % bins_data(1) - 1) * t % stride(l)
                 end do
 
                 ! Get p1 scatter rr and convert to p1 scatter xs

src/constants.F90

@@ -17,6 +17,9 @@ module constants
   ! HDF5 data format
   integer, parameter :: HDF5_VERSION(2) = [1, 0]
 
+  ! WMP data format
+  integer, parameter :: WMP_VERSION(2) = [1, 1]
+
   ! Version numbers for binary files
   integer, parameter :: VERSION_STATEPOINT(2)       = [17, 0]
   integer, parameter :: VERSION_PARTICLE_RESTART(2) = [2, 0]
@@ -25,7 +28,6 @@ module constants
   integer, parameter :: VERSION_VOLUME(2)           = [1, 0]
   integer, parameter :: VERSION_VOXEL(2)            = [1, 0]
   integer, parameter :: VERSION_MGXS_LIBRARY(2)     = [1, 0]
-  character(10), parameter :: VERSION_MULTIPOLE     = "v1.0"
 
   ! ============================================================================
   ! ADJUSTABLE PARAMETERS

src/dagmc.cpp

@@ -54,7 +54,7 @@ void load_dagmc_geometry()
     c->fill_ = C_NONE; // no fill, single universe
 
     cells.push_back(c);
-    cell_map[c->id_] = c->id_;
+    cell_map[c->id_] = i;
 
     // Populate the Universe vector and dict
     auto it = universe_map.find(dagmc_univ_id);

src/geometry_aux.cpp

@@ -12,6 +12,8 @@
 #include "openmc/material.h"
 #include "openmc/settings.h"
 #include "openmc/surface.h"
+#include "openmc/tallies/filter.h"
+#include "openmc/tallies/filter_distribcell.h"
 
 
 namespace openmc {
@@ -181,12 +183,15 @@ neighbor_lists()
 //==============================================================================
 
 void
-prepare_distribcell(int32_t* filter_cell_list, int n)
+prepare_distribcell()
 {
-  // Read the list of cells contained in distribcell filters from Fortran.
+  // Find all cells listed in a DistribcellFilter.
   std::unordered_set<int32_t> distribcells;
-  for (int i = 0; i < n; i++) {
-    distribcells.insert(filter_cell_list[i]);
+  for (auto* filt : tally_filters) {
+    if (filt->type() == "distribcell") {
+      auto* distrib_filt = static_cast<DistribcellFilter*>(filt);
+      distribcells.insert(distrib_filt->cell_);
+    }
   }
 
   // Find all cells with distributed materials or temperatures.  Make sure that
@@ -393,29 +398,11 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset,
   }
 }
 
-//==============================================================================
-
-int
-distribcell_path_len(int32_t target_cell, int32_t map, int32_t target_offset,
-                     int32_t root_univ)
+std::string
+distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset)
 {
-  Universe& root = *universes[root_univ];
-  std::string path_ {distribcell_path_inner(target_cell, map, target_offset,
-                                            root, 0)};
-  return path_.size() + 1;
-}
-
-//==============================================================================
-
-void
-distribcell_path(int32_t target_cell, int32_t map, int32_t target_offset,
-                 int32_t root_univ, char* path)
-{
-  Universe& root = *universes[root_univ];
-  std::string path_ {distribcell_path_inner(target_cell, map, target_offset,
-                                            root, 0)};
-  path_.copy(path, path_.size());
-  path[path_.size()] = '\0';
+  auto& root_univ = *universes[openmc_root_universe];
+  return distribcell_path_inner(target_cell, map, target_offset, root_univ, 0);
 }
 
 //==============================================================================

src/geometry_header.F90

@@ -59,20 +59,6 @@ module geometry_header
       integer(C_INT32_T)             :: fill
     end function cell_fill_c
 
-    function cell_n_instances_c(cell_ptr) bind(C, name='cell_n_instances') &
-         result(n_instances)
-      import C_PTR, C_INT32_T
-      type(C_PTR), intent(in), value :: cell_ptr
-      integer(C_INT32_T)             :: n_instances
-    end function cell_n_instances_c
-
-    function cell_distribcell_index_c(cell_ptr) &
-         bind(C, name='cell_distribcell_index') result(distribcell_index)
-      import C_PTR, C_INT
-      type(C_PTR), intent(in), value :: cell_ptr
-      integer(C_INT)                 :: distribcell_index
-    end function cell_distribcell_index_c
-
     function cell_material_size_c(cell_ptr) bind(C, name='cell_material_size') &
          result(n)
       import C_PTR, C_INT
@@ -103,14 +89,6 @@ module geometry_header
       real(C_DOUBLE)                    :: sqrtkT
     end function cell_sqrtkT_c
 
-    function cell_offset_c(cell_ptr, map) bind(C, name="cell_offset") &
-         result(offset)
-      import C_PTR, C_INT, C_INT32_T
-      type(C_PTR),    intent(in), value :: cell_ptr
-      integer(C_INT), intent(in), value :: map
-      integer(C_INT32_T)                :: offset
-    end function cell_offset_c
-
     function lattice_pointer(lat_ind) bind(C) result(ptr)
       import C_PTR, C_INT32_T
       integer(C_INT32_T), intent(in), value :: lat_ind
@@ -123,23 +101,6 @@ module geometry_header
       integer(C_INT32_T)             :: id
     end function lattice_id_c
 
-    function lattice_are_valid_indices_c(lat_ptr, i_xyz) &
-         bind(C, name='lattice_are_valid_indices') result (is_valid)
-      import C_PTR, C_INT, C_BOOL
-      type(C_PTR),    intent(in), value :: lat_ptr
-      integer(C_INT), intent(in)        :: i_xyz(3)
-      logical(C_BOOL)                   :: is_valid
-    end function lattice_are_valid_indices_c
-
-    function lattice_offset_c(lat_ptr, map, i_xyz) &
-         bind(C, name='lattice_offset') result(offset)
-      import C_PTR, C_INT, C_INT32_T
-      type(C_PTR),    intent(in), value :: lat_ptr
-      integer(C_INT), intent(in), value :: map
-      integer(C_INT), intent(in)        :: i_xyz(3)
-      integer(C_INT32_T)                :: offset
-    end function lattice_offset_c
-
     subroutine extend_cells_c(n) bind(C)
       import C_INT32_t
       integer(C_INT32_T), intent(in), value :: n
@@ -154,8 +115,6 @@ module geometry_header
     type(C_PTR) :: ptr
   contains
     procedure :: id => lattice_id
-    procedure :: are_valid_indices => lattice_are_valid_indices
-    procedure :: offset => lattice_offset
   end type Lattice
 
 !===============================================================================
@@ -176,13 +135,10 @@ module geometry_header
     procedure :: type => cell_type
     procedure :: universe => cell_universe
     procedure :: fill => cell_fill
-    procedure :: n_instances => cell_n_instances
-    procedure :: distribcell_index => cell_distribcell_index
     procedure :: material_size => cell_material_size
     procedure :: material => cell_material
     procedure :: sqrtkT_size => cell_sqrtkT_size
     procedure :: sqrtkT => cell_sqrtkT
-    procedure :: offset => cell_offset
 
   end type Cell
 
@@ -197,7 +153,6 @@ module geometry_header
 
   ! Dictionaries which map user IDs to indices in the global arrays
   type(DictIntInt) :: cell_dict
-  type(DictIntInt) :: universe_dict
   type(DictIntInt) :: lattice_dict
 
 contains
@@ -208,21 +163,6 @@ contains
     id = lattice_id_c(this % ptr)
   end function lattice_id
 
-  function lattice_are_valid_indices(this, i_xyz) result (is_valid)
-    class(Lattice), intent(in) :: this
-    integer(C_INT), intent(in) :: i_xyz(3)
-    logical(C_BOOL)            :: is_valid
-    is_valid = lattice_are_valid_indices_c(this % ptr, i_xyz)
-  end function lattice_are_valid_indices
-
-  function lattice_offset(this, map, i_xyz) result(offset)
-    class(Lattice), intent(in) :: this
-    integer(C_INT), intent(in) :: map
-    integer(C_INT), intent(in) :: i_xyz(3)
-    integer(C_INT32_T)         :: offset
-    offset = lattice_offset_c(this % ptr, map, i_xyz)
-  end function lattice_offset
-
 !===============================================================================
 
   function cell_id(this) result(id)
@@ -255,18 +195,6 @@ contains
     fill = cell_fill_c(this % ptr)
   end function cell_fill
 
-  function cell_n_instances(this) result(n_instances)
-    class(Cell), intent(in) :: this
-    integer(C_INT32_T)      :: n_instances
-    n_instances = cell_n_instances_c(this % ptr)
-  end function cell_n_instances
-
-  function cell_distribcell_index(this) result(distribcell_index)
-    class(Cell), intent(in) :: this
-    integer(C_INT)          :: distribcell_index
-    distribcell_index = cell_distribcell_index_c(this % ptr)
-  end function cell_distribcell_index
-
   function cell_material_size(this) result(n)
     class(Cell), intent(in) :: this
     integer(C_INT)          :: n
@@ -293,13 +221,6 @@ contains
     sqrtkT = cell_sqrtkT_c(this % ptr, i)
   end function cell_sqrtkT
 
-  function cell_offset(this, map) result(offset)
-    class(Cell), intent(in)    :: this
-    integer(C_INT), intent(in) :: map
-    integer(C_INT32_T)         :: offset
-    offset = cell_offset_c(this % ptr, map)
-  end function cell_offset
-
 !===============================================================================
 ! GET_TEMPERATURES returns a list of temperatures that each nuclide/S(a,b) table
 ! appears at in the model. Later, this list is used to determine the actual
@@ -383,7 +304,6 @@ contains
     if (allocated(lattices)) deallocate(lattices)
 
     call cell_dict % clear()
-    call universe_dict % clear()
     call lattice_dict % clear()
 
   end subroutine free_memory_geometry

src/input_xml.F90

@@ -22,6 +22,7 @@ module input_xml
   use mgxs_data,        only: create_macro_xs, read_mgxs
   use mgxs_interface
   use nuclide_header
+  use multipole_header
   use output,           only: title, header, print_plot
   use photon_header
   use plot_header
@@ -59,12 +60,8 @@ module input_xml
       integer(C_INT32_T), intent(in), value :: univ_indx
     end subroutine count_cell_instances
 
-    subroutine prepare_distribcell_c(cell_list, n) &
-         bind(C, name="prepare_distribcell")
-      import C_INT32_T, C_INT
-      integer(C_INT),     intent(in), value :: n
-      integer(C_INT32_T), intent(in)        :: cell_list(n)
-    end subroutine prepare_distribcell_c
+    subroutine prepare_distribcell() bind(C)
+    end subroutine prepare_distribcell
 
     subroutine read_surfaces(node_ptr) bind(C)
       import C_PTR
@@ -253,7 +250,7 @@ contains
 
   subroutine read_geometry_dagmc()
 
-    integer :: i, j
+    integer :: i
     integer :: univ_id
     integer :: n_cells_in_univ
     logical :: file_exists
@@ -285,7 +282,6 @@ contains
       if (.not. cells_in_univ_dict % has(univ_id)) then
         n_universes = n_universes + 1
         n_cells_in_univ = 1
-        call universe_dict % set(univ_id, n_universes)
         call univ_ids % push_back(univ_id)
       else
         n_cells_in_univ = 1 + cells_in_univ_dict % get(univ_id)
@@ -363,9 +359,6 @@ contains
       surfaces(i) % ptr = surface_pointer(i - 1);
 
       if (surfaces(i) % bc() /= BC_TRANSMIT) boundary_exists = .true.
-
-      ! Add surface to dictionary
-      call surface_dict % set(surfaces(i) % id(), i)
     end do
 
     ! Check to make sure a boundary condition was applied to at least one
@@ -455,7 +448,6 @@ contains
       if (.not. cells_in_univ_dict % has(univ_id)) then
         n_universes = n_universes + 1
         n_cells_in_univ = 1
-        call universe_dict % set(univ_id, n_universes - 1)
         call univ_ids % push_back(univ_id)
       else
         n_cells_in_univ = 1 + cells_in_univ_dict % get(univ_id)
@@ -524,8 +516,6 @@ contains
 
     do i = 1, n_surfaces
       surfaces(i) % ptr = surface_pointer(i - 1);
-      ! Add surface to dictionary
-      call surface_dict % set(surfaces(i) % id(), i)
     end do
 
   end subroutine allocate_surfaces
@@ -3014,7 +3004,9 @@ contains
 
     logical :: file_exists                 ! Does multipole library exist?
     character(7) :: readable               ! Is multipole library readable?
-    character(MAX_FILE_LEN) :: filename  ! Path to multipole xs library
+    character(MAX_FILE_LEN) :: filename    ! Path to multipole xs library
+    integer(HID_T) :: file_id
+    integer(HID_T) :: group_id
 
     ! For the time being, and I know this is a bit hacky, we just assume
     ! that the file will be ZZZAAAmM.h5.
@@ -3040,45 +3032,24 @@ contains
       end if
 
       ! Display message
-      call write_message("Loading Multipole XS table: " // filename, 6)
+      call write_message("Loading Windowed Multipole XS from " // filename, 6)
 
+      ! Open file and make sure version is sufficient
+      file_id = file_open(filename, 'r')
+      call check_wmp_version(file_id)
+
+      ! Read nuclide data from HDF5
+      group_id = open_group(file_id, nuc % name)
       allocate(nuc % multipole)
-
-      ! Call the read routine
-      call nuc % multipole % from_hdf5(filename)
+      call nuc % multipole % from_hdf5(group_id)
       nuc % mp_present = .true.
 
+      ! Close the group and file.
+      call close_group(group_id)
+      call file_close(file_id)
+
     end associate
 
   end subroutine read_multipole_data
 
-!===============================================================================
-! PREPARE_DISTRIBCELL initializes any distribcell filters present and sets the
-! offsets for distribcells
-!===============================================================================
-
-  subroutine prepare_distribcell()
-
-    integer :: i, j
-    type(SetInt)  :: cell_list  ! distribcells to track
-    integer(C_INT32_T), allocatable :: cell_list_c(:)
-
-    ! Find all cells listed in a distribcell filter.
-    do i = 1, n_tallies
-      do j = 1, size(tallies(i) % obj % filter)
-        select type(filt => filters(tallies(i) % obj % filter(j)) % obj)
-        type is (DistribcellFilter)
-          call cell_list % add(filt % cell)
-        end select
-      end do
-    end do
-
-    allocate(cell_list_c(cell_list % size()))
-    do i = 1, cell_list % size()
-      cell_list_c(i) = cell_list % get_item(i) - 1
-    end do
-    call prepare_distribcell_c(cell_list_c, cell_list % size())
-
-  end subroutine prepare_distribcell
-
 end module input_xml

src/lattice.cpp

@@ -893,12 +893,6 @@ extern "C" {
   Lattice* lattice_pointer(int lat_ind) {return lattices[lat_ind];}
 
   int32_t lattice_id(Lattice *lat) {return lat->id_;}
-
-  bool lattice_are_valid_indices(Lattice *lat, const int i_xyz[3])
-  {return lat->are_valid_indices(i_xyz);}
-
-  int32_t lattice_offset(Lattice *lat, int map, const int i_xyz[3])
-  {return lat->offset(map, i_xyz);}
 }
 
 } // namespace openmc

src/material.cpp

@@ -23,7 +23,7 @@ std::unordered_map<int32_t, int32_t> material_map;
 Material::Material(pugi::xml_node node)
 {
   if (check_for_node(node, "id")) {
-    id = std::stoi(get_node_value(node, "id"));
+    id_ = std::stoi(get_node_value(node, "id"));
   } else {
     fatal_error("Must specify id of material in materials XML file.");
   }
@@ -52,7 +52,7 @@ read_materials(pugi::xml_node* node)
 
   // Populate the material map.
   for (int i = 0; i < materials.size(); i++) {
-    int32_t mid = materials[i]->id;
+    int32_t mid = materials[i]->id_;
     auto search = material_map.find(mid);
     if (search == material_map.end()) {
       material_map[mid] = i;
@@ -78,7 +78,7 @@ openmc_material_get_volume(int32_t index, double* volume)
       return 0;
     } else {
       std::stringstream msg;
-      msg << "Volume for material with ID=" << m->id << " not set.";
+      msg << "Volume for material with ID=" << m->id_ << " not set.";
       set_errmsg(msg);
       return OPENMC_E_UNASSIGNED;
     }
@@ -113,11 +113,11 @@ openmc_material_set_volume(int32_t index, double volume)
 extern "C" {
   Material* material_pointer(int32_t indx) {return materials[indx];}
 
-  int32_t material_id(Material* mat) {return mat->id;}
+  int32_t material_id(Material* mat) {return mat->id_;}
 
   void material_set_id(Material* mat, int32_t id, int32_t index)
   {
-    mat->id = id;
+    mat->id_ = id;
     //TODO: off-by-one
     material_map[id] = index - 1;
   }

src/math.F90

@@ -10,8 +10,6 @@ module math
   public :: t_percentile
   public :: calc_pn
   public :: calc_rn
-  public :: calc_zn
-  public :: calc_zn_rad
   public :: rotate_angle
   public :: faddeeva
   public :: w_derivative
@@ -47,23 +45,6 @@ module math
       real(C_DOUBLE), intent(out) :: rn(2 * n + 1)
     end subroutine calc_rn
 
-    pure subroutine calc_zn(n, rho, phi, zn) bind(C, name='calc_zn_c')
-      use ISO_C_BINDING
-      implicit none
-      integer(C_INT), value, intent(in) :: n
-      real(C_DOUBLE), value, intent(in) :: rho
-      real(C_DOUBLE), value, intent(in) :: phi
-      real(C_DOUBLE), intent(out) :: zn(((n + 1) * (n + 2)) / 2)
-    end subroutine calc_zn
-
-    pure subroutine calc_zn_rad(n, rho, zn_rad) bind(C, name='calc_zn_rad_c')
-      use ISO_C_BINDING
-      implicit none
-      integer(C_INT), value, intent(in) :: n
-      real(C_DOUBLE), value, intent(in) :: rho
-      real(C_DOUBLE), intent(out) :: zn_rad((n / 2) + 1)
-    end subroutine calc_zn_rad
-
     subroutine rotate_angle_c_intfc(uvw, mu, phi) bind(C, name='rotate_angle_c')
       use ISO_C_BINDING
       implicit none

src/math_functions.cpp

@@ -510,7 +510,7 @@ void calc_rn_c(int n, const double uvw[3], double rn[]){
 }
 
 
-void calc_zn_c(int n, double rho, double phi, double zn[]) {
+void calc_zn(int n, double rho, double phi, double zn[]) {
   // ===========================================================================
   // Determine vector of sin(n*phi) and cos(n*phi). This takes advantage of the
   // following recurrence relations so that only a single sin/cos have to be
@@ -587,7 +587,7 @@ void calc_zn_c(int n, double rho, double phi, double zn[]) {
 
 }
 
-void calc_zn_rad_c(int n, double rho, double zn_rad[]) {
+void calc_zn_rad(int n, double rho, double zn_rad[]) {
   // Calculate R_p0(rho) as Zn_p0(rho)
   // Set up the array of the coefficients
 

src/mesh.cpp

@@ -20,6 +20,7 @@
 #include "openmc/hdf5_interface.h"
 #include "openmc/message_passing.h"
 #include "openmc/search.h"
+#include "openmc/tallies/filter.h"
 #include "openmc/xml_interface.h"
 
 namespace openmc {
@@ -899,10 +900,6 @@ void meshes_to_hdf5(hid_t group)
 //==============================================================================
 
 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(); }
@@ -941,35 +938,6 @@ extern "C" {
     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();

src/multipole_header.F90

@@ -1,8 +1,9 @@
 module multipole_header
 
   use constants
-  use error,            only: fatal_error
   use hdf5_interface
+  use error,            only: fatal_error
+  use string,           only: to_str
 
   implicit none
 
@@ -30,17 +31,15 @@ module multipole_header
 
   type MultipoleArray
 
-    !=========================================================================
+    character(20) :: name                   ! name of nuclide, e.g. U235
+
     ! Isotope Properties
+    logical                 :: fissionable  ! Is this isotope fissionable?
+    complex(8), allocatable :: data(:,:)    ! Poles and residues
+    real(8)                 :: sqrtAWR      ! Square root of the atomic
+                                            ! weight ratio
 
-    logical                 :: fissionable     ! Is this isotope fissionable?
-    complex(8), allocatable :: data(:,:)       ! Poles and residues
-    real(8)                 :: sqrtAWR         ! Square root of the atomic
-                                               ! weight ratio
-
-    !=========================================================================
     ! Windows
-
     integer :: fit_order                    ! Order of the fit. 1 linear,
                                             !   2 quadratic, etc.
     real(8) :: E_min                        ! Start energy for the windows
@@ -53,11 +52,9 @@ module multipole_header
     real(8), allocatable :: curvefit(:,:,:) ! Contains the fitting function.
                                             !   (reaction type, coeff index,
                                             !   window index)
-
     integer, allocatable :: broaden_poly(:) ! if 1, broaden, if 0, don't.
 
   contains
-
     procedure :: from_hdf5 => multipole_from_hdf5
 
   end type MultipoleArray
@@ -65,29 +62,22 @@ module multipole_header
 contains
 
 !===============================================================================
-! FROM_HDF5 loads multipole data from an HDF5 file.
+! FROM_HDF5 loads multipole data from a HDF5 file.
 !===============================================================================
 
-  subroutine multipole_from_hdf5(this, filename)
+  subroutine multipole_from_hdf5(this, group_id)
     class(MultipoleArray), intent(inout) :: this
-    character(len=*),      intent(in)    :: filename
+    integer(HID_T),        intent(in)    :: group_id
 
-    character(len=10) :: version
     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
 
-    ! Open file for reading and move into the /isotope group
-    file_id = file_open(filename, 'r', parallel=.true.)
-    group_id = open_group(file_id, "/nuclide")
+    ! Get name of nuclide from group
+    this % name = get_name(group_id)
 
-    ! Check the file version number.
-    call read_dataset(version, file_id, "version")
-    if (version /= VERSION_MULTIPOLE) call fatal_error("The current multipole&
-         & format version is " // trim(VERSION_MULTIPOLE) // " but the file "&
-         // trim(filename) // " uses version " // trim(version) // ".")
+    ! Get rid of leading '/'
+    this % name = trim(this % name(2:))
 
     ! Read scalar values.
     call read_dataset(this % spacing, group_id, "spacing")
@@ -121,8 +111,8 @@ contains
     dset = open_dataset(group_id, "broaden_poly")
     call get_shape(dset, dims_1d)
     if (dims_1d(1) /= n_windows) call fatal_error("broaden_poly array shape is&
-         &not consistent with the windows array shape in multipole library"&
-         // trim(filename) // ".")
+         &not consistent with the windows array shape in WMP library for"&
+         // trim(this % name) // ".")
     allocate(this % broaden_poly(n_windows))
     call read_dataset(this % broaden_poly, dset)
     call close_dataset(dset)
@@ -131,15 +121,38 @@ contains
     dset = open_dataset(group_id, "curvefit")
     call get_shape(dset, dims_3d)
     if (dims_3d(3) /= n_windows) call fatal_error("curvefit array shape is not&
-         &consistent with the windows array shape in multipole library"&
-         // trim(filename) // ".")
+         &consistent with the windows array shape in WMP library for"&
+         // trim(this % name) // ".")
     allocate(this % curvefit(dims_3d(1), dims_3d(2), dims_3d(3)))
     call read_dataset(this % curvefit, dset)
     call close_dataset(dset)
     this % fit_order = int(dims_3d(2), 4) - 1
 
-    ! Close the group and file.
-    call close_group(group_id)
-    call file_close(file_id)
   end subroutine multipole_from_hdf5
+
+!===============================================================================
+! CHECK_WMP_VERSION checks for the right version of WMP data within HDF5
+! files
+!===============================================================================
+
+  subroutine check_wmp_version(file_id)
+    integer(HID_T), intent(in) :: file_id
+
+    integer, allocatable :: version(:)
+
+    if (attribute_exists(file_id, 'version')) then
+      call read_attribute(version, file_id, 'version')
+      if (version(1) /= WMP_VERSION(1)) then
+        call fatal_error("WMP data format uses version " // trim(to_str(&
+             version(1))) // "." // trim(to_str(version(2))) // " whereas &
+             &your installation of OpenMC expects version " // trim(to_str(&
+             WMP_VERSION(1))) // ".x data.")
+      end if
+    else
+      call fatal_error("WMP data does not indicate a version. Your &
+           &installation of OpenMC expects version " // trim(to_str(&
+           WMP_VERSION(1))) // ".x data.")
+    end if
+  end subroutine check_wmp_version
+
 end module multipole_header

src/nuclide_header.F90

@@ -53,7 +53,7 @@ module nuclide_header
 
   type :: Nuclide
     ! Nuclide meta-data
-    character(20) :: name    ! name of nuclide, e.g. U235.71c
+    character(20) :: name    ! name of nuclide, e.g. U235
     integer       :: Z       ! atomic number
     integer       :: A       ! mass number
     integer       :: metastable ! metastable state

src/output.F90

@@ -653,17 +653,13 @@ contains
     character(36)           :: score_names(N_SCORE_TYPES)  ! names of scoring function
     character(36)           :: score_name                  ! names of scoring function
                                                            ! to be applied at write-time
-    type(TallyFilterMatch), allocatable :: matches(:)
+    integer, allocatable :: filter_bins(:)
     character(MAX_WORD_LEN) :: temp_name
 
     ! Skip if there are no tallies
     if (n_tallies == 0) return
 
-    allocate(matches(n_filters))
-    do i = 1, n_filters
-      allocate(matches(i) % bins)
-      allocate(matches(i) % weights)
-    end do
+    allocate(filter_bins(n_filters))
 
     ! Initialize names for scores
     score_names(abs(SCORE_FLUX))               = "Flux"
@@ -747,8 +743,7 @@ contains
 
       ! Initialize bins, filter level, and indentation
       do h = 1, size(t % filter)
-        call matches(t % filter(h)) % bins % clear()
-        call matches(t % filter(h)) % bins % push_back(0)
+        filter_bins(t % filter(h)) = 0
       end do
       j = 1
       indent = 0
@@ -759,18 +754,17 @@ contains
           if (size(t % filter) == 0) exit find_bin
 
           ! Increment bin combination
-          matches(t % filter(j)) % bins % data(1) = &
-               matches(t % filter(j)) % bins % data(1) + 1
+          filter_bins(t % filter(j)) = filter_bins(t % filter(j)) + 1
 
           ! =================================================================
           ! REACHED END OF BINS FOR THIS FILTER, MOVE TO NEXT FILTER
 
-          if (matches(t % filter(j)) % bins % data(1) > &
+          if (filter_bins(t % filter(j)) > &
                filters(t % filter(j)) % obj % n_bins) then
             ! If this is the first filter, then exit
             if (j == 1) exit print_bin
 
-            matches(t % filter(j)) % bins % data(1) = 0
+            filter_bins(t % filter(j)) = 0
             j = j - 1
             indent = indent - 2
 
@@ -784,7 +778,7 @@ contains
             ! Print current filter information
             write(UNIT=unit_tally, FMT='(1X,2A)') repeat(" ", indent), &
                  trim(filters(t % filter(j)) % obj % &
-                 text_label(matches(t % filter(j)) % bins % data(1)))
+                 text_label(filter_bins(t % filter(j))))
             indent = indent + 2
             j = j + 1
           end if
@@ -795,7 +789,7 @@ contains
         if (size(t % filter) > 0) then
           write(UNIT=unit_tally, FMT='(1X,2A)') repeat(" ", indent), &
                trim(filters(t % filter(j)) % obj % &
-               text_label(matches(t % filter(j)) % bins % data(1)))
+               text_label(filter_bins(t % filter(j))))
         end if
 
         ! Determine scoring index for this bin combination -- note that unlike
@@ -804,8 +798,8 @@ contains
 
         filter_index = 1
         do h = 1, size(t % filter)
-          filter_index = filter_index + (max(matches(t % filter(h)) &
-               % bins % data(1),1) - 1) * t % stride(h)
+          filter_index = filter_index &
+               + (max(filter_bins(t % filter(h)) ,1) - 1) * t % stride(h)
         end do
 
         ! Write results for this filter bin combination
@@ -859,11 +853,6 @@ 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/simulation.F90

@@ -16,7 +16,7 @@ module simulation
   use settings
   use simulation_header
   use tally_header
-  use tally_filter_header, only: filter_matches, n_filters
+  use tally_filter_header, only: filter_matches, n_filters, filter_match_pointer
 
   implicit none
   private
@@ -44,8 +44,7 @@ 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)
+      filter_matches(i) % ptr = filter_match_pointer(i - 1)
     end do
 !$omp end parallel
 
@@ -65,10 +64,6 @@ 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.cpp

@@ -12,6 +12,7 @@
 #include "openmc/source.h"
 #include "openmc/state_point.h"
 #include "openmc/timer.h"
+#include "openmc/tallies/filter.h"
 #include "openmc/tallies/tally.h"
 
 #include <algorithm>
@@ -53,7 +54,8 @@ extern "C" void write_tallies();
 // over the batches, generations, and histories in a fixed source or k-eigenvalue
 // calculation.
 
-int openmc_run() {
+int openmc_run()
+{
   openmc_simulation_init();
 
   int err = 0;
@@ -76,6 +78,12 @@ int openmc_simulation_init()
   // Determine how much work each process should do
   calculate_work();
 
+  // Allocate array for matching filter bins
+#pragma omp parallel
+  {
+    filter_matches.resize(n_filters);
+  }
+
   // Set up tally procedure pointers
   init_tally_routines();
 
@@ -134,6 +142,11 @@ int openmc_simulation_finalize()
   time_active.stop();
   time_finalize.start();
 
+#pragma omp parallel
+  {
+    filter_matches.clear();
+  }
+
   // Deallocate Fortran variables, set tallies to inactive
   simulation_finalize_f();
 

src/surface_header.F90

@@ -2,7 +2,6 @@ module surface_header
 
   use, intrinsic :: ISO_C_BINDING
 
-  use dict_header, only: DictIntInt
   use hdf5_interface
 
   implicit none
@@ -83,9 +82,6 @@ module surface_header
 
   type(Surface), allocatable, target :: surfaces(:)
 
-  ! Dictionary that maps user IDs to indices in 'surfaces'
-  type(DictIntInt) :: surface_dict
-
 contains
 
   pure function surface_id(this) result(id)
@@ -128,7 +124,6 @@ contains
 
   subroutine free_memory_surfaces()
     if (allocated(surfaces)) deallocate(surfaces)
-    call surface_dict % clear()
     call free_memory_surfaces_c()
   end subroutine free_memory_surfaces
 

src/tallies/filter.cpp

@@ -0,0 +1,170 @@
+#include "openmc/tallies/filter.h"
+
+#include <string>
+
+#include "openmc/capi.h"
+#include "openmc/constants.h"  // for MAX_LINE_LEN;
+#include "openmc/error.h"
+#include "openmc/tallies/filter_azimuthal.h"
+#include "openmc/tallies/filter_cell.h"
+#include "openmc/tallies/filter_cellborn.h"
+#include "openmc/tallies/filter_cellfrom.h"
+#include "openmc/tallies/filter_distribcell.h"
+#include "openmc/tallies/filter_energyfunc.h"
+#include "openmc/tallies/filter_legendre.h"
+#include "openmc/tallies/filter_material.h"
+#include "openmc/tallies/filter_mesh.h"
+#include "openmc/tallies/filter_meshsurface.h"
+#include "openmc/tallies/filter_mu.h"
+#include "openmc/tallies/filter_polar.h"
+#include "openmc/tallies/filter_sph_harm.h"
+#include "openmc/tallies/filter_sptl_legendre.h"
+#include "openmc/tallies/filter_surface.h"
+#include "openmc/tallies/filter_universe.h"
+#include "openmc/tallies/filter_zernike.h"
+
+
+namespace openmc {
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+std::vector<FilterMatch> filter_matches;
+std::vector<Filter*> tally_filters;
+
+//==============================================================================
+// Non-member functions
+//==============================================================================
+
+void
+free_memory_tally_c()
+{
+  #pragma omp parallel
+  {
+    filter_matches.clear();
+  }
+
+  for (Filter* filt : tally_filters) {delete filt;}
+  tally_filters.clear();
+}
+
+//==============================================================================
+// Fortran compatibility functions
+//==============================================================================
+
+extern "C" {
+  FilterMatch* filter_match_pointer(int indx)
+  {return &filter_matches[indx];}
+
+  void
+  filter_match_bins_push_back(FilterMatch* match, int val)
+  {match->bins_.push_back(val);}
+
+  void
+  filter_match_weights_push_back(FilterMatch* match, double val)
+  {match->weights_.push_back(val);}
+
+  void
+  filter_match_bins_clear(FilterMatch* match)
+  {match->bins_.clear();}
+
+  void
+  filter_match_weights_clear(FilterMatch* match)
+  {match->weights_.clear();}
+
+  int
+  filter_match_bins_size(FilterMatch* match)
+  {return match->bins_.size();}
+
+  int
+  filter_match_bins_data(FilterMatch* match, int indx)
+  {return match->bins_[indx-1];}
+
+  double
+  filter_match_weights_data(FilterMatch* match, int indx)
+  {return match->weights_[indx-1];}
+
+  void
+  filter_match_bins_set_data(FilterMatch* match, int indx, int val)
+  {match->bins_[indx-1] = val;}
+
+  Filter*
+  allocate_filter(const char* type)
+  {
+    std::string type_ {type};
+    if (type_ == "azimuthal") {
+      tally_filters.push_back(new AzimuthalFilter());
+    } else if (type_ == "cell") {
+      tally_filters.push_back(new CellFilter());
+    } else if (type_ == "cellborn") {
+      tally_filters.push_back(new CellbornFilter());
+    } else if (type_ == "cellfrom") {
+      tally_filters.push_back(new CellFromFilter());
+    } else if (type_ == "distribcell") {
+      tally_filters.push_back(new DistribcellFilter());
+    } else if (type_ == "energyfunction") {
+      tally_filters.push_back(new EnergyFunctionFilter());
+    } else if (type_ == "legendre") {
+      tally_filters.push_back(new LegendreFilter());
+    } else if (type_ == "material") {
+      tally_filters.push_back(new MaterialFilter());
+    } else if (type_ == "mesh") {
+      tally_filters.push_back(new MeshFilter());
+    } else if (type_ == "meshsurface") {
+      tally_filters.push_back(new MeshSurfaceFilter());
+    } else if (type_ == "mu") {
+      tally_filters.push_back(new MuFilter());
+    } else if (type_ == "polar") {
+      tally_filters.push_back(new PolarFilter());
+    } else if (type_ == "surface") {
+      tally_filters.push_back(new SurfaceFilter());
+    } else if (type_ == "spatiallegendre") {
+      tally_filters.push_back(new SpatialLegendreFilter());
+    } else if (type_ == "sphericalharmonics") {
+      tally_filters.push_back(new SphericalHarmonicsFilter());
+    } else if (type_ == "universe") {
+      tally_filters.push_back(new UniverseFilter());
+    } else if (type_ == "zernike") {
+      tally_filters.push_back(new ZernikeFilter());
+    } else if (type_ == "zernikeradial") {
+      tally_filters.push_back(new ZernikeRadialFilter());
+    } else {
+      return nullptr;
+    }
+    return tally_filters.back();
+  }
+
+  void filter_from_xml(Filter* filt, pugi::xml_node* node)
+  {filt->from_xml(*node);}
+
+  void
+  filter_get_all_bins(Filter* filt, Particle* p, int estimator,
+                      FilterMatch* match)
+  {
+    filt->get_all_bins(p, estimator, *match);
+  }
+
+  void filter_to_statepoint(Filter* filt, hid_t group)
+  {filt->to_statepoint(group);}
+
+  void filter_text_label(Filter* filt, int bin, char* label)
+  {
+    std::string label_str = filt->text_label(bin);
+    int i = 0;
+    for (; i < label_str.size() && i < MAX_LINE_LEN; i++)
+      label[i] = label_str[i];
+    label[i] = '\0';
+  }
+
+  void filter_initialize(Filter* filt) {filt->initialize();}
+
+  int filter_n_bins(Filter* filt) {return filt->n_bins_;}
+
+  int mesh_filter_get_mesh(MeshFilter* filt) {return filt->mesh();}
+
+  int sphharm_filter_get_cosine(SphericalHarmonicsFilter* filt)
+  {return static_cast<int>(filt->cosine_);}
+}
+
+} // namespace openmc

src/tallies/filter_azimuthal.cpp

@@ -0,0 +1,73 @@
+#include "openmc/tallies/filter_azimuthal.h"
+
+#include <cmath>
+#include <sstream>
+
+#include "openmc/constants.h"
+#include "openmc/error.h"
+#include "openmc/search.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+void
+AzimuthalFilter::from_xml(pugi::xml_node node)
+{
+  auto bins = get_node_array<double>(node, "bins");
+
+  if (bins.size() > 1) {
+    bins_ = bins;
+
+  } else {
+    // Allow a user to input a lone number which will mean that you subdivide
+    // [-pi,pi) evenly with the input being the number of bins
+
+    int n_angle = bins[0];
+
+    if (n_angle <= 1) fatal_error("Number of bins for azimuthal filter must "
+                                  "be greater than 1.");
+
+    double d_angle = 2.0 * PI / n_angle;
+    bins_.resize(n_angle + 1);
+    for (int i = 0; i < n_angle; i++) bins_[i] = -PI + i * d_angle;
+    bins_[n_angle] = PI;
+  }
+
+  n_bins_ = bins_.size() - 1;
+}
+
+void
+AzimuthalFilter::get_all_bins(const Particle* p, int estimator,
+                              FilterMatch& match) const
+{
+  double phi;
+  if (estimator == ESTIMATOR_TRACKLENGTH) {
+    phi = std::atan2(p->coord[0].uvw[1], p->coord[0].uvw[0]);
+  } else {
+    phi = std::atan2(p->last_uvw[1], p->last_uvw[0]);
+  }
+
+  if (phi >= bins_.front() && phi <= bins_.back()) {
+    //TODO: off-by-one
+    auto bin = lower_bound_index(bins_.begin(), bins_.end(), phi) + 1;
+    match.bins_.push_back(bin);
+    match.weights_.push_back(1.0);
+  }
+}
+
+void
+AzimuthalFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  write_dataset(filter_group, "bins", bins_);
+}
+
+std::string
+AzimuthalFilter::text_label(int bin) const
+{
+  std::stringstream out;
+  out << "Azimuthal Angle [" << bins_[bin-1] << ", " << bins_[bin] << ")";
+  return out.str();
+}
+
+} // namespace openmc

src/tallies/filter_cell.cpp

@@ -0,0 +1,91 @@
+#include "openmc/tallies/filter_cell.h"
+
+#include <sstream>
+
+#include "openmc/capi.h"
+#include "openmc/cell.h"
+#include "openmc/error.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+void
+CellFilter::from_xml(pugi::xml_node node)
+{
+  cells_ = get_node_array<int32_t>(node, "bins");
+  n_bins_ = cells_.size();
+}
+
+void
+CellFilter::initialize()
+{
+  // Convert cell IDs to indices of the global array.
+  for (auto& c : cells_) {
+    auto search = cell_map.find(c);
+    if (search != cell_map.end()) {
+      c = search->second;
+    } else {
+      std::stringstream err_msg;
+      err_msg << "Could not find cell " << c << " specified on tally filter.";
+      fatal_error(err_msg);
+    }
+  }
+
+  // Populate the index->bin map.
+  for (int i = 0; i < cells_.size(); i++) {
+    map_[cells_[i]] = i;
+  }
+}
+
+void
+CellFilter::get_all_bins(const Particle* p, int estimator,
+                         FilterMatch& match) const
+{
+  for (int i = 0; i < p->n_coord; i++) {
+    auto search = map_.find(p->coord[i].cell);
+    if (search != map_.end()) {
+      //TODO: off-by-one
+      match.bins_.push_back(search->second + 1);
+      match.weights_.push_back(1.0);
+    }
+  }
+}
+
+void
+CellFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  std::vector<int32_t> cell_ids;
+  for (auto c : cells_) cell_ids.push_back(cells[c]->id_);
+  write_dataset(filter_group, "bins", cell_ids);
+}
+
+std::string
+CellFilter::text_label(int bin) const
+{
+  return "Cell " + std::to_string(cells[cells_[bin-1]]->id_);
+}
+
+//==============================================================================
+// C-API functions
+//==============================================================================
+
+extern "C" int
+openmc_cell_filter_get_bins(int32_t index, int32_t** cells, int32_t* n)
+{
+  int err = verify_filter(index);
+  if (err) return err;
+
+  auto filt = filter_from_f(index);
+  if (filt->type() != "cell") {
+    set_errmsg("Tried to get cells from a non-cell filter.");
+    return OPENMC_E_INVALID_TYPE;
+  }
+
+  auto cell_filt = static_cast<CellFilter*>(filt);
+  *cells = cell_filt->cells_.data();
+  *n = cell_filt->cells_.size();
+  return 0;
+}
+
+} // namespace openmc

src/tallies/filter_cellborn.cpp

@@ -0,0 +1,25 @@
+#include "openmc/tallies/filter_cellborn.h"
+
+#include "openmc/cell.h"
+
+namespace openmc {
+
+void
+CellbornFilter::get_all_bins(const Particle* p, int estimator,
+                             FilterMatch& match) const
+{
+  auto search = map_.find(p->cell_born);
+  if (search != map_.end()) {
+    //TODO: off-by-one
+    match.bins_.push_back(search->second + 1);
+    match.weights_.push_back(1.0);
+  }
+}
+
+std::string
+CellbornFilter::text_label(int bin) const
+{
+  return "Birth Cell " + std::to_string(cells[cells_[bin-1]]->id_);
+}
+
+} // namespace openmc

src/tallies/filter_cellfrom.cpp

@@ -0,0 +1,27 @@
+#include "openmc/tallies/filter_cellfrom.h"
+
+#include "openmc/cell.h"
+
+namespace openmc {
+
+void
+CellFromFilter::get_all_bins(const Particle* p, int estimator,
+                             FilterMatch& match) const
+{
+  for (int i = 0; i < p->last_n_coord; i++) {
+    auto search = map_.find(p->last_cell[i]);
+    if (search != map_.end()) {
+      //TODO: off-by-one
+      match.bins_.push_back(search->second + 1);
+      match.weights_.push_back(1.0);
+    }
+  }
+}
+
+std::string
+CellFromFilter::text_label(int bin) const
+{
+  return "Cell from " + std::to_string(cells[cells_[bin-1]]->id_);
+}
+
+} // namespace openmc

src/tallies/filter_distribcell.cpp

@@ -0,0 +1,80 @@
+#include "openmc/tallies/filter_distribcell.h"
+
+#include "openmc/cell.h"
+#include "openmc/error.h"
+#include "openmc/geometry_aux.h" // For distribcell_path
+#include "openmc/lattice.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+void
+DistribcellFilter::from_xml(pugi::xml_node node)
+{
+  auto cells = get_node_array<int32_t>(node, "bins");
+  if (cells.size() != 1) {
+    fatal_error("Only one cell can be specified per distribcell filter.");
+  }
+  cell_ = cells[0];
+}
+
+void
+DistribcellFilter::initialize()
+{
+  // Convert the cell ID to an index of the global array.
+  auto search = cell_map.find(cell_);
+  if (search != cell_map.end()) {
+    cell_ = search->second;
+    n_bins_ = cells[cell_]->n_instances_;
+  } else {
+    std::stringstream err_msg;
+    err_msg << "Could not find cell " << cell_
+            << " specified on tally filter.";
+    fatal_error(err_msg);
+  }
+}
+
+void
+DistribcellFilter::get_all_bins(const Particle* p, int estimator,
+                                FilterMatch& match) const
+{
+  int offset = 0;
+  auto distribcell_index = cells[cell_]->distribcell_index_;
+  for (int i = 0; i < p->n_coord; i++) {
+    auto& c {*cells[p->coord[i].cell]};
+    if (c.type_ == FILL_UNIVERSE) {
+      offset += c.offset_[distribcell_index];
+    } else if (c.type_ == FILL_LATTICE) {
+      auto& lat {*lattices[p->coord[i+1].lattice-1]};
+      int i_xyz[3] {p->coord[i+1].lattice_x,
+                    p->coord[i+1].lattice_y,
+                    p->coord[i+1].lattice_z};
+      if (lat.are_valid_indices(i_xyz)) {
+        offset += lat.offset(distribcell_index, i_xyz);
+      }
+    }
+    if (cell_ == p->coord[i].cell) {
+      //TODO: off-by-one
+      match.bins_.push_back(offset + 1);
+      match.weights_.push_back(1.0);
+      return;
+    }
+  }
+}
+
+void
+DistribcellFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  write_dataset(filter_group, "bins", cells[cell_]->id_);
+}
+
+std::string
+DistribcellFilter::text_label(int bin) const
+{
+  auto map = cells[cell_]->distribcell_index_;
+  auto path = distribcell_path(cell_, map, bin-1);
+  return "Distributed Cell " + path;
+}
+
+} // namespace openmc

src/tallies/filter_energyfunc.cpp

@@ -0,0 +1,69 @@
+#include "openmc/tallies/filter_energyfunc.h"
+
+#include <iomanip>  // for setprecision
+#include <ios>  // for scientific
+#include <sstream>
+
+#include "openmc/error.h"
+#include "openmc/search.h"
+#include "openmc/settings.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+void
+EnergyFunctionFilter::from_xml(pugi::xml_node node)
+{
+  if (!settings::run_CE)
+    fatal_error("EnergyFunction filters are only supported for "
+                "continuous-energy transport calculations");
+
+  if (!check_for_node(node, "energy"))
+    fatal_error("Energy grid not specified for EnergyFunction filter.");
+
+  energy_ = get_node_array<double>(node, "energy");
+
+  if (!check_for_node(node, "y"))
+    fatal_error("y values not specified for EnergyFunction filter.");
+
+  y_ = get_node_array<double>(node, "y");
+}
+
+void
+EnergyFunctionFilter::get_all_bins(const Particle* p, int estimator,
+                                   FilterMatch& match) const
+{
+  if (p->last_E >= energy_.front() && p->last_E <= energy_.back()) {
+    // Search for the incoming energy bin.
+    auto i = lower_bound_index(energy_.begin(), energy_.end(), p->last_E);
+
+    // Compute the interpolation factor between the nearest bins.
+    double f = (p->last_E - energy_[i]) / (energy_[i+1] - energy_[i]);
+
+    // Interpolate on the lin-lin grid.
+    //TODO: off-by-one
+    match.bins_.push_back(1);
+    match.weights_.push_back((1-f) * y_[i] + f * y_[i+1]);
+  }
+}
+
+void
+EnergyFunctionFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  write_dataset(filter_group, "energy", energy_);
+  write_dataset(filter_group, "y", y_);
+}
+
+std::string
+EnergyFunctionFilter::text_label(int bin) const
+{
+  std::stringstream out;
+  out << std::scientific << std::setprecision(1)
+      << "Energy Function f"
+      << "([ " << energy_.front() << ", ..., " << energy_.back() << "]) = "
+      << "[" << y_.front() << ", ..., " << y_.back() << "]";
+  return out.str();
+}
+
+} // namespace openmc

src/tallies/filter_legendre.cpp

@@ -0,0 +1,93 @@
+#include "openmc/tallies/filter_legendre.h"
+
+#include "openmc/capi.h"
+#include "openmc/error.h"
+#include "openmc/math_functions.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+void
+LegendreFilter::from_xml(pugi::xml_node node)
+{
+  order_ = std::stoi(get_node_value(node, "order"));
+  n_bins_ = order_ + 1;
+}
+
+void
+LegendreFilter::get_all_bins(const Particle* p, int estimator,
+                             FilterMatch& match) const
+{
+  double wgt[n_bins_];
+  calc_pn_c(order_, p->mu, wgt);
+  for (int i = 0; i < n_bins_; i++) {
+    //TODO: off-by-one
+    match.bins_.push_back(i + 1);
+    match.weights_.push_back(wgt[i]);
+  }
+}
+
+void
+LegendreFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  write_dataset(filter_group, "order", order_);
+}
+
+std::string
+LegendreFilter::text_label(int bin) const
+{
+  return "Legendre expansion, P" + std::to_string(bin - 1);
+}
+
+//==============================================================================
+// C-API functions
+//==============================================================================
+
+extern "C" int
+openmc_legendre_filter_get_order(int32_t index, int* order)
+{
+  // Make sure this is a valid index to an allocated filter.
+  int err = verify_filter(index);
+  if (err) return err;
+
+  // Get a pointer to the filter and downcast.
+  auto* filt_base = filter_from_f(index);
+  auto* filt = dynamic_cast<LegendreFilter*>(filt_base);
+
+  // Check the filter type.
+  if (!filt) {
+    set_errmsg("Not a legendre filter.");
+    return OPENMC_E_INVALID_TYPE;
+  }
+
+  // Output the order.
+  *order = filt->order_;
+  return 0;
+}
+
+extern "C" int
+openmc_legendre_filter_set_order(int32_t index, int order)
+{
+  // Make sure this is a valid index to an allocated filter.
+  int err = verify_filter(index);
+  if (err) return err;
+
+  // Get a pointer to the filter and downcast.
+  auto* filt_base = filter_from_f(index);
+  auto* filt = dynamic_cast<LegendreFilter*>(filt_base);
+
+  // Check the filter type.
+  if (!filt) {
+    set_errmsg("Not a legendre filter.");
+    return OPENMC_E_INVALID_TYPE;
+  }
+
+  // Update the filter.
+  filt->order_ = order;
+  filt->n_bins_ = order + 1;
+  filter_update_n_bins(index);
+  return 0;
+}
+
+} // namespace openmc

src/tallies/filter_material.cpp

@@ -0,0 +1,123 @@
+#include "openmc/tallies/filter_material.h"
+
+#include <sstream>
+
+#include "openmc/capi.h"
+#include "openmc/error.h"
+#include "openmc/material.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+void
+MaterialFilter::from_xml(pugi::xml_node node)
+{
+  materials_ = get_node_array<int32_t>(node, "bins");
+  n_bins_ = materials_.size();
+}
+
+void
+MaterialFilter::initialize()
+{
+  // Convert material IDs to indices of the global array.
+  for (auto& m : materials_) {
+    auto search = material_map.find(m);
+    if (search != material_map.end()) {
+      m = search->second;
+    } else {
+      std::stringstream err_msg;
+      err_msg << "Could not find material " << m
+              << " specified on tally filter.";
+      fatal_error(err_msg);
+    }
+  }
+
+  // Populate the index->bin map.
+  for (int i = 0; i < materials_.size(); i++) {
+    map_[materials_[i]] = i;
+  }
+}
+
+void
+MaterialFilter::get_all_bins(const Particle* p, int estimator,
+                             FilterMatch& match) const
+{
+  auto search = map_.find(p->material - 1);
+  if (search != map_.end()) {
+    //TODO: off-by-one
+    match.bins_.push_back(search->second + 1);
+    match.weights_.push_back(1.0);
+  }
+}
+
+void
+MaterialFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  std::vector<int32_t> material_ids;
+  for (auto c : materials_) material_ids.push_back(materials[c]->id_);
+  write_dataset(filter_group, "bins", material_ids);
+}
+
+std::string
+MaterialFilter::text_label(int bin) const
+{
+  return "Material " + std::to_string(materials[materials_[bin-1]]->id_);
+}
+
+//==============================================================================
+// C-API functions
+//==============================================================================
+
+extern "C" int
+openmc_material_filter_get_bins(int32_t index, int32_t** bins, int32_t* n)
+{
+  // Make sure this is a valid index to an allocated filter.
+  int err = verify_filter(index);
+  if (err) return err;
+
+  // Get a pointer to the filter and downcast.
+  auto* filt_base = filter_from_f(index);
+  auto* filt = dynamic_cast<MaterialFilter*>(filt_base);
+
+  // Check the filter type.
+  if (!filt) {
+    set_errmsg("Tried to get material filter bins on a non-material filter.");
+    return OPENMC_E_INVALID_TYPE;
+  }
+
+  // Output the bins.
+  *bins = filt->materials_.data();
+  *n = filt->materials_.size();
+  return 0;
+}
+
+extern "C" int
+openmc_material_filter_set_bins(int32_t index, int32_t n, const int32_t* bins)
+{
+  // Make sure this is a valid index to an allocated filter.
+  int err = verify_filter(index);
+  if (err) return err;
+
+  // Get a pointer to the filter and downcast.
+  auto* filt_base = filter_from_f(index);
+  auto* filt = dynamic_cast<MaterialFilter*>(filt_base);
+
+  // Check the filter type.
+  if (!filt) {
+    set_errmsg("Tried to set material filter bins on a non-material filter.");
+    return OPENMC_E_INVALID_TYPE;
+  }
+
+  // Update the filter.
+  filt->materials_.clear();
+  filt->materials_.resize(n);
+  for (int i = 0; i < n; i++) filt->materials_[i] = bins[i];
+  filt->n_bins_ = filt->materials_.size();
+  filt->map_.clear();
+  for (int i = 0; i < n; i++) filt->map_[filt->materials_[i]] = i;
+  filter_update_n_bins(index);
+  return 0;
+}
+
+} // namespace openmc

src/tallies/filter_mesh.cpp

@@ -0,0 +1,137 @@
+#include "openmc/tallies/filter_mesh.h"
+
+#include <sstream>
+
+#include "openmc/capi.h"
+#include "openmc/constants.h"
+#include "openmc/error.h"
+#include "openmc/mesh.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+void
+MeshFilter::from_xml(pugi::xml_node node)
+{
+  auto bins_ = get_node_array<int32_t>(node, "bins");
+  if (bins_.size() != 1) {
+    fatal_error("Only one mesh can be specified per " + type()
+                + " mesh filter.");
+  }
+
+  auto id = bins_[0];
+  auto search = mesh_map.find(id);
+  if (search != mesh_map.end()) {
+    set_mesh(search->second);
+  } else{
+    std::stringstream err_msg;
+    err_msg << "Could not find cell " << id << " specified on tally filter.";
+    fatal_error(err_msg);
+  }
+}
+
+void
+MeshFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+const
+{
+  if (estimator != ESTIMATOR_TRACKLENGTH) {
+    auto bin = meshes[mesh_]->get_bin(p->coord[0].xyz);
+    if (bin >= 0) {
+      //TODO: off-by-one
+      match.bins_.push_back(bin);
+      match.weights_.push_back(1.0);
+    }
+  } else {
+    meshes[mesh_]->bins_crossed(p, match.bins_, match.weights_);
+  }
+}
+
+void
+MeshFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  write_dataset(filter_group, "bins", meshes[mesh_]->id_);
+}
+
+std::string
+MeshFilter::text_label(int bin) const
+{
+  auto& mesh = *meshes[mesh_];
+  int n_dim = mesh.n_dimension_;
+
+  int ijk[n_dim];
+  mesh.get_indices_from_bin(bin, ijk);
+
+  std::stringstream out;
+  out << "Mesh Index (" << ijk[0];
+  if (n_dim > 1) out << ", " << ijk[1];
+  if (n_dim > 2) out << ", " << ijk[2];
+  out << ")";
+
+  return out.str();
+}
+
+void
+MeshFilter::set_mesh(int32_t mesh)
+{
+  mesh_ = mesh;
+  n_bins_ = 1;
+  for (auto dim : meshes[mesh_]->shape_) n_bins_ *= dim;
+}
+
+//==============================================================================
+// C-API functions
+//==============================================================================
+
+extern "C" int
+openmc_mesh_filter_get_mesh(int32_t index, int32_t* index_mesh)
+{
+  // Make sure this is a valid index to an allocated filter.
+  int err = verify_filter(index);
+  if (err) return err;
+
+  // Get a pointer to the filter and downcast.
+  auto* filt_base = filter_from_f(index);
+  auto* filt = dynamic_cast<MeshFilter*>(filt_base);
+
+  // Check the filter type.
+  if (!filt) {
+    set_errmsg("Tried to get mesh on a non-mesh filter.");
+    return OPENMC_E_INVALID_TYPE;
+  }
+
+  // Output the mesh.
+  *index_mesh = filt->mesh();
+  return 0;
+}
+
+extern "C" int
+openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh)
+{
+  // Make sure this is a valid index to an allocated filter.
+  int err = verify_filter(index);
+  if (err) return err;
+
+  // Get a pointer to the filter and downcast.
+  auto* filt_base = filter_from_f(index);
+  auto* filt = dynamic_cast<MeshFilter*>(filt_base);
+
+  // Check the filter type.
+  if (!filt) {
+    set_errmsg("Tried to set mesh on a non-mesh filter.");
+    return OPENMC_E_INVALID_TYPE;
+  }
+
+  // Check the mesh index.
+  if (index_mesh < 0 || index_mesh >= meshes.size()) {
+    set_errmsg("Index in 'meshes' array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+
+  // Update the filter.
+  filt->set_mesh(index_mesh);
+  filter_update_n_bins(index);
+  return 0;
+}
+
+} // namespace openmc

src/tallies/filter_meshsurface.cpp

@@ -0,0 +1,94 @@
+#include "openmc/tallies/filter_meshsurface.h"
+
+#include "openmc/capi.h"
+#include "openmc/constants.h"
+#include "openmc/error.h"
+#include "openmc/mesh.h"
+
+namespace openmc {
+
+void
+MeshSurfaceFilter::get_all_bins(const Particle* p, int estimator,
+                                FilterMatch& match) const
+{
+  meshes[mesh_]->surface_bins_crossed(p, match.bins_);
+  for (auto b : match.bins_) match.weights_.push_back(1.0);
+}
+
+std::string
+MeshSurfaceFilter::text_label(int bin) const
+{
+  auto& mesh = *meshes[mesh_];
+  int n_dim = mesh.n_dimension_;
+
+  // Get flattend mesh index and surface index.
+  int i_mesh = (bin - 1) / (4 * n_dim) + 1;
+  int i_surf = ((bin - 1) % (4 * n_dim)) + 1;
+
+  // Get mesh index part of label.
+  std::string out = MeshFilter::text_label(i_mesh);
+
+  // Get surface part of label.
+  switch (i_surf) {
+    case OUT_LEFT:
+      out += " Outgoing, x-min";
+      break;
+    case IN_LEFT:
+      out += " Incoming, x-min";
+      break;
+    case OUT_RIGHT:
+      out += " Outgoing, x-max";
+      break;
+    case IN_RIGHT:
+      out += " Incoming, x-max";
+      break;
+    case OUT_BACK:
+      out += " Outgoing, y-min";
+      break;
+    case IN_BACK:
+      out += " Incoming, y-min";
+      break;
+    case OUT_FRONT:
+      out += " Outgoing, y-max";
+      break;
+    case IN_FRONT:
+      out += " Incoming, y-max";
+      break;
+    case OUT_BOTTOM:
+      out += " Outgoing, z-min";
+      break;
+    case IN_BOTTOM:
+      out += " Incoming, z-min";
+      break;
+    case OUT_TOP:
+      out += " Outgoing, z-max";
+      break;
+    case IN_TOP:
+      out += " Incoming, z-max";
+      break;
+  }
+
+  return out;
+}
+
+void
+MeshSurfaceFilter::set_mesh(int32_t mesh)
+{
+  mesh_ = mesh;
+  n_bins_ = 4 * meshes[mesh_]->n_dimension_;
+  for (auto dim : meshes[mesh_]->shape_) n_bins_ *= dim;
+}
+
+//==============================================================================
+// C-API functions
+//==============================================================================
+
+extern"C" int
+openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh)
+{return openmc_mesh_filter_get_mesh(index, index_mesh);}
+
+extern"C" int
+openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh)
+{return openmc_mesh_filter_set_mesh(index, index_mesh);}
+
+} // namespace openmc

src/tallies/filter_mu.cpp

@@ -0,0 +1,64 @@
+#include "openmc/tallies/filter_mu.h"
+
+#include <sstream>
+
+#include "openmc/error.h"
+#include "openmc/search.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+void
+MuFilter::from_xml(pugi::xml_node node)
+{
+  auto bins = get_node_array<double>(node, "bins");
+
+  if (bins.size() > 1) {
+    bins_ = bins;
+
+  } else {
+    // Allow a user to input a lone number which will mean that you subdivide
+    // [-1,1) evenly with the input being the number of bins
+
+    int n_angle = bins[0];
+
+    if (n_angle <= 1) fatal_error("Number of bins for mu filter must "
+                                  "be greater than 1.");
+
+    double d_angle = 2.0 / n_angle;
+    bins_.resize(n_angle + 1);
+    for (int i = 0; i < n_angle; i++) bins_[i] = -1 + i * d_angle;
+    bins_[n_angle] = 1;
+  }
+
+  n_bins_ = bins_.size() - 1;
+}
+
+void
+MuFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+const
+{
+  if (p->mu >= bins_.front() && p->mu <= bins_.back()) {
+    //TODO: off-by-one
+    auto bin = lower_bound_index(bins_.begin(), bins_.end(), p->mu) + 1;
+    match.bins_.push_back(bin);
+    match.weights_.push_back(1.0);
+  }
+}
+
+void
+MuFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  write_dataset(filter_group, "bins", bins_);
+}
+
+std::string
+MuFilter::text_label(int bin) const
+{
+  std::stringstream out;
+  out << "Change-in-Angle [" << bins_[bin-1] << ", " << bins_[bin] << ")";
+  return out.str();
+}
+
+} // namespace openmc

src/tallies/filter_polar.cpp

@@ -0,0 +1,72 @@
+#include "openmc/tallies/filter_polar.h"
+
+#include <sstream>
+
+#include "openmc/constants.h"
+#include "openmc/error.h"
+#include "openmc/search.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+void
+PolarFilter::from_xml(pugi::xml_node node)
+{
+  auto bins = get_node_array<double>(node, "bins");
+
+  if (bins.size() > 1) {
+    bins_ = bins;
+
+  } else {
+    // Allow a user to input a lone number which will mean that you subdivide
+    // [0,pi] evenly with the input being the number of bins
+
+    int n_angle = bins[0];
+
+    if (n_angle <= 1) fatal_error("Number of bins for polar filter must "
+                                  "be greater than 1.");
+
+    double d_angle = PI / n_angle;
+    bins_.resize(n_angle + 1);
+    for (int i = 0; i < n_angle; i++) bins_[i] = i * d_angle;
+    bins_[n_angle] = PI;
+  }
+
+  n_bins_ = bins_.size() - 1;
+}
+
+void
+PolarFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+const
+{
+  double theta;
+  if (estimator == ESTIMATOR_TRACKLENGTH) {
+    theta = std::acos(p->coord[0].uvw[2]);
+  } else {
+    theta = std::acos(p->last_uvw[2]);
+  }
+
+  if (theta >= bins_.front() && theta <= bins_.back()) {
+    //TODO: off-by-one
+    auto bin = lower_bound_index(bins_.begin(), bins_.end(), theta) + 1;
+    match.bins_.push_back(bin);
+    match.weights_.push_back(1.0);
+  }
+}
+
+void
+PolarFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  write_dataset(filter_group, "bins", bins_);
+}
+
+std::string
+PolarFilter::text_label(int bin) const
+{
+  std::stringstream out;
+  out << "Polar Angle [" << bins_[bin-1] << ", " << bins_[bin] << ")";
+  return out.str();
+}
+
+} // namespace openmc

src/tallies/filter_sph_harm.cpp

@@ -0,0 +1,182 @@
+#include "openmc/tallies/filter_sph_harm.h"
+
+#include <utility>  // For pair
+
+#include "openmc/capi.h"
+#include "openmc/error.h"
+#include "openmc/math_functions.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+void
+SphericalHarmonicsFilter::from_xml(pugi::xml_node node)
+{
+  order_ = std::stoi(get_node_value(node, "order"));
+  n_bins_ = (order_ + 1) * (order_ + 1);
+
+  if (check_for_node(node, "cosine")) {
+    auto cos = get_node_value(node, "cosine", true);
+    if (cos == "scatter") {
+      cosine_ = SphericalHarmonicsCosine::scatter;
+    } else if (cos == "particle") {
+      cosine_ = SphericalHarmonicsCosine::particle;
+    } else {
+      std::stringstream err_msg;
+      err_msg << "Unrecognized cosine type, \"" << cos
+              << "\" in spherical harmonics filter";
+      fatal_error(err_msg);
+    }
+  }
+}
+
+void
+SphericalHarmonicsFilter::get_all_bins(const Particle* p, int estimator,
+                                       FilterMatch& match) const
+{
+  // Determine cosine term for scatter expansion if necessary
+  double wgt[order_ + 1];
+  if (cosine_ == SphericalHarmonicsCosine::scatter) {
+    calc_pn_c(order_, p->mu, wgt);
+  } else {
+    for (int i = 0; i < order_ + 1; i++) wgt[i] = 1;
+  }
+
+  // Find the Rn,m values
+  double rn[n_bins_];
+  calc_rn_c(order_, p->last_uvw, rn);
+
+  int j = 0;
+  for (int n = 0; n < order_ + 1; n++) {
+    // Calculate n-th order spherical harmonics for (u,v,w)
+    int num_nm = 2*n + 1;
+
+    // Append the matching (bin,weight) for each moment
+    for (int i = 0; i < num_nm; i++) {
+      match.weights_.push_back(wgt[n] * rn[j]);
+      //TODO: off-by-one
+      match.bins_.push_back(++j);
+    }
+  }
+}
+
+void
+SphericalHarmonicsFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  write_dataset(filter_group, "order", order_);
+  if (cosine_ == SphericalHarmonicsCosine::scatter) {
+    write_dataset(filter_group, "cosine", "scatter");
+  } else {
+    write_dataset(filter_group, "cosine", "particle");
+  }
+}
+
+std::string
+SphericalHarmonicsFilter::text_label(int bin) const
+{
+  std::stringstream out;
+  for (int n = 0; n < order_ + 1; n++) {
+    if (bin <= (n + 1) * (n + 1)) {
+      int m = (bin - n*n - 1) - n;
+      out << "Spherical harmonic expansion, Y" << n << "," << m;
+      return out.str();
+    }
+  }
+}
+
+//==============================================================================
+// C-API functions
+//==============================================================================
+
+std::pair<int, SphericalHarmonicsFilter*>
+check_sphharm_filter(int32_t index)
+{
+  // Make sure this is a valid index to an allocated filter.
+  int err = verify_filter(index);
+  if (err) {
+    return {err, nullptr};
+  }
+
+  // Get a pointer to the filter and downcast.
+  auto* filt_base = filter_from_f(index);
+  auto* filt = dynamic_cast<SphericalHarmonicsFilter*>(filt_base);
+
+  // Check the filter type.
+  if (!filt) {
+    set_errmsg("Not a spherical harmonics filter.");
+    err = OPENMC_E_INVALID_TYPE;
+  }
+  return {err, filt};
+}
+
+extern "C" int
+openmc_sphharm_filter_get_order(int32_t index, int* order)
+{
+  // Check the filter.
+  auto check_result = check_sphharm_filter(index);
+  auto err = check_result.first;
+  auto filt = check_result.second;
+  if (err) return err;
+
+  // Output the order.
+  *order = filt->order_;
+  return 0;
+}
+
+extern "C" int
+openmc_sphharm_filter_get_cosine(int32_t index, char cosine[])
+{
+  // Check the filter.
+  auto check_result = check_sphharm_filter(index);
+  auto err = check_result.first;
+  auto filt = check_result.second;
+  if (err) return err;
+
+  // Output the cosine.
+  if (filt->cosine_ == SphericalHarmonicsCosine::scatter) {
+    strcpy(cosine, "scatter");
+  } else {
+    strcpy(cosine, "particle");
+  }
+  return 0;
+}
+
+extern "C" int
+openmc_sphharm_filter_set_order(int32_t index, int order)
+{
+  // Check the filter.
+  auto check_result = check_sphharm_filter(index);
+  auto err = check_result.first;
+  auto filt = check_result.second;
+  if (err) return err;
+
+  // Update the filter.
+  filt->order_ = order;
+  filt->n_bins_ = (order + 1) * (order + 1);
+  filter_update_n_bins(index);
+  return 0;
+}
+
+extern "C" int
+openmc_sphharm_filter_set_cosine(int32_t index, const char cosine[])
+{
+  // Check the filter.
+  auto check_result = check_sphharm_filter(index);
+  auto err = check_result.first;
+  auto filt = check_result.second;
+  if (err) return err;
+
+  // Update the filter.
+  if (strcmp(cosine, "scatter") == 0) {
+    filt->cosine_ = SphericalHarmonicsCosine::scatter;
+  } else if (strcmp(cosine, "particle") == 0) {
+    filt->cosine_ = SphericalHarmonicsCosine::particle;
+  } else {
+    set_errmsg("Invalid spherical harmonics cosine.");
+    return OPENMC_E_INVALID_ARGUMENT;
+  }
+  return 0;
+}
+
+} // namespace openmc

src/tallies/filter_sptl_legendre.cpp

@@ -0,0 +1,184 @@
+#include "openmc/tallies/filter_sptl_legendre.h"
+
+#include <utility>  // For pair
+
+#include "openmc/capi.h"
+#include "openmc/error.h"
+#include "openmc/math_functions.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+void
+SpatialLegendreFilter::from_xml(pugi::xml_node node)
+{
+  order_ = std::stoi(get_node_value(node, "order"));
+
+  auto axis = get_node_value(node, "axis");
+  if (axis == "x") {
+    axis_ = LegendreAxis::x;
+  } else if (axis == "y") {
+    axis_ = LegendreAxis::y;
+  } else if (axis == "z") {
+    axis_ = LegendreAxis::z;
+  } else {
+    fatal_error("Unrecognized axis on SpatialLegendreFilter");
+  }
+
+  min_ = std::stod(get_node_value(node, "min"));
+  max_ = std::stod(get_node_value(node, "max"));
+
+  n_bins_ = order_ + 1;
+}
+
+void
+SpatialLegendreFilter::get_all_bins(const Particle* p, int estimator,
+                                    FilterMatch& match) const
+{
+  // Get the coordinate along the axis of interest.
+  double x;
+  if (axis_ == LegendreAxis::x) {
+    x = p->coord[0].xyz[0];
+  } else if (axis_ == LegendreAxis::y) {
+    x = p->coord[0].xyz[1];
+  } else {
+    x = p->coord[0].xyz[2];
+  }
+
+  if (x >= min_ && x <= max_) {
+    // Compute the normalized coordinate value.
+    double x_norm = 2.0*(x - min_) / (max_ - min_) - 1.0;
+
+    // Compute and return the Legendre weights.
+    double wgt[order_ + 1];
+    calc_pn_c(order_, x_norm, wgt);
+    for (int i = 0; i < order_ + 1; i++) {
+      //TODO: off-by-one
+      match.bins_.push_back(i + 1);
+      match.weights_.push_back(wgt[i]);
+    }
+  }
+}
+
+void
+SpatialLegendreFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  write_dataset(filter_group, "order", order_);
+  if (axis_ == LegendreAxis::x) {
+    write_dataset(filter_group, "axis", "x");
+  } else if (axis_ == LegendreAxis::y) {
+    write_dataset(filter_group, "axis", "y");
+  } else {
+    write_dataset(filter_group, "axis", "z");
+  }
+  write_dataset(filter_group, "min", min_);
+  write_dataset(filter_group, "max", max_);
+}
+
+std::string
+SpatialLegendreFilter::text_label(int bin) const
+{
+  std::stringstream out;
+  out << "Legendre expansion, ";
+  if (axis_ == LegendreAxis::x) {
+    out << "x";
+  } else if (axis_ == LegendreAxis::y) {
+    out << "y";
+  } else {
+    out << "z";
+  }
+  out << " axis, P" << std::to_string(bin - 1);
+  return out.str();
+}
+
+//==============================================================================
+// C-API functions
+//==============================================================================
+
+std::pair<int, SpatialLegendreFilter*>
+check_sptl_legendre_filter(int32_t index)
+{
+  // Make sure this is a valid index to an allocated filter.
+  int err = verify_filter(index);
+  if (err) {
+    return {err, nullptr};
+  }
+
+  // Get a pointer to the filter and downcast.
+  auto* filt_base = filter_from_f(index);
+  auto* filt = dynamic_cast<SpatialLegendreFilter*>(filt_base);
+
+  // Check the filter type.
+  if (!filt) {
+    set_errmsg("Not a spatial Legendre filter.");
+    err = OPENMC_E_INVALID_TYPE;
+  }
+  return {err, filt};
+}
+
+extern "C" int
+openmc_spatial_legendre_filter_get_order(int32_t index, int* order)
+{
+  // Check the filter.
+  auto check_result = check_sptl_legendre_filter(index);
+  auto err = check_result.first;
+  auto filt = check_result.second;
+  if (err) return err;
+
+  // Output the order.
+  *order = filt->order_;
+  return 0;
+}
+
+extern "C" int
+openmc_spatial_legendre_filter_get_params(int32_t index, int* axis,
+                                          double* min, double* max)
+{
+  // Check the filter.
+  auto check_result = check_sptl_legendre_filter(index);
+  auto err = check_result.first;
+  auto filt = check_result.second;
+  if (err) return err;
+
+  // Output the params.
+  *axis = static_cast<int>(filt->axis_);
+  *min = filt->min_;
+  *max = filt->max_;
+  return 0;
+}
+
+extern "C" int
+openmc_spatial_legendre_filter_set_order(int32_t index, int order)
+{
+  // Check the filter.
+  auto check_result = check_sptl_legendre_filter(index);
+  auto err = check_result.first;
+  auto filt = check_result.second;
+  if (err) return err;
+
+  // Update the filter.
+  filt->order_ = order;
+  filt->n_bins_ = order + 1;
+  filter_update_n_bins(index);
+  return 0;
+}
+
+extern "C" int
+openmc_spatial_legendre_filter_set_params(int32_t index, const int* axis,
+  const double* min, const double* max)
+{
+  // Check the filter.
+  auto check_result = check_sptl_legendre_filter(index);
+  auto err = check_result.first;
+  auto filt = check_result.second;
+  if (err) return err;
+
+  // Update the filter.
+  if (axis) filt->axis_ = static_cast<LegendreAxis>(*axis);
+  if (min) filt->min_ = *min;
+  if (max) filt->max_ = *max;
+  return 0;
+}
+
+} // namespace openmc

src/tallies/filter_surface.cpp

@@ -0,0 +1,71 @@
+#include "openmc/tallies/filter_surface.h"
+
+#include <sstream>
+
+#include "openmc/error.h"
+#include "openmc/surface.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+void
+SurfaceFilter::from_xml(pugi::xml_node node)
+{
+  surfaces_ = get_node_array<int32_t>(node, "bins");
+  n_bins_ = surfaces_.size();
+}
+
+void
+SurfaceFilter::initialize()
+{
+  // Convert surface IDs to indices of the global array.
+  for (auto& s : surfaces_) {
+    auto search = surface_map.find(s);
+    if (search != surface_map.end()) {
+      s = search->second;
+    } else {
+      std::stringstream err_msg;
+      err_msg << "Could not find surface " << s
+              << " specified on tally filter.";
+      fatal_error(err_msg);
+    }
+  }
+
+  // Populate the index->bin map.
+  for (int i = 0; i < surfaces_.size(); i++) {
+    map_[surfaces_[i]] = i;
+  }
+}
+
+void
+SurfaceFilter::get_all_bins(const Particle* p, int estimator,
+                            FilterMatch& match) const
+{
+  auto search = map_.find(std::abs(p->surface)-1);
+  if (search != map_.end()) {
+    //TODO: off-by-one
+    match.bins_.push_back(search->second + 1);
+    if (p->surface < 0) {
+      match.weights_.push_back(-1.0);
+    } else {
+      match.weights_.push_back(1.0);
+    }
+  }
+}
+
+void
+SurfaceFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  std::vector<int32_t> surface_ids;
+  for (auto c : surfaces_) surface_ids.push_back(surfaces[c]->id_);
+  write_dataset(filter_group, "bins", surface_ids);
+}
+
+std::string
+SurfaceFilter::text_label(int bin) const
+{
+  return "Surface " + std::to_string(surfaces[surfaces_[bin-1]]->id_);
+}
+
+} // namespace openmc

src/tallies/filter_universe.cpp

@@ -0,0 +1,69 @@
+#include "openmc/tallies/filter_universe.h"
+
+#include <sstream>
+
+#include "openmc/cell.h"
+#include "openmc/error.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+void
+UniverseFilter::from_xml(pugi::xml_node node)
+{
+  universes_ = get_node_array<int32_t>(node, "bins");
+  n_bins_ = universes_.size();
+}
+
+void
+UniverseFilter::initialize()
+{
+  // Convert universe IDs to indices of the global array.
+  for (auto& u : universes_) {
+    auto search = universe_map.find(u);
+    if (search != universe_map.end()) {
+      u = search->second;
+    } else {
+      std::stringstream err_msg;
+      err_msg << "Could not find universe " << u
+              << " specified on tally filter.";
+      fatal_error(err_msg);
+    }
+  }
+
+  // Populate the index->bin map.
+  for (int i = 0; i < universes_.size(); i++) {
+    map_[universes_[i]] = i;
+  }
+}
+
+void
+UniverseFilter::get_all_bins(const Particle* p, int estimator,
+                             FilterMatch& match) const
+{
+  for (int i = 0; i < p->n_coord; i++) {
+    auto search = map_.find(p->coord[i].universe);
+    if (search != map_.end()) {
+      //TODO: off-by-one
+      match.bins_.push_back(search->second + 1);
+      match.weights_.push_back(1.0);
+    }
+  }
+}
+
+void
+UniverseFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  std::vector<int32_t> universe_ids;
+  for (auto u : universes_) universe_ids.push_back(universes[u]->id_);
+  write_dataset(filter_group, "bins", universe_ids);
+}
+
+std::string
+UniverseFilter::text_label(int bin) const
+{
+  return "Universe " + std::to_string(universes[universes_[bin-1]]->id_);
+}
+
+} // namespace openmc

src/tallies/filter_zernike.cpp

@@ -0,0 +1,207 @@
+#include "openmc/tallies/filter_zernike.h"
+
+#include <cmath>
+#include <sstream>
+#include <utility>  // For pair
+
+#include "openmc/capi.h"
+#include "openmc/error.h"
+#include "openmc/math_functions.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+//==============================================================================
+// ZernikeFilter implementation
+//==============================================================================
+
+void
+ZernikeFilter::from_xml(pugi::xml_node node)
+{
+  set_order(std::stoi(get_node_value(node, "order")));
+  x_ = std::stod(get_node_value(node, "x"));
+  y_ = std::stod(get_node_value(node, "y"));
+  r_ = std::stod(get_node_value(node, "r"));
+}
+
+void
+ZernikeFilter::get_all_bins(const Particle* p, int estimator,
+                            FilterMatch& match) const
+{
+  // Determine the normalized (r,theta) coordinates.
+  double x = p->coord[0].xyz[0] - x_;
+  double y = p->coord[0].xyz[1] - y_;
+  double r = std::sqrt(x*x + y*y) / r_;
+  double theta = std::atan2(y, x);
+
+  if (r <= 1.0) {
+    // Compute and return the Zernike weights.
+    double zn[n_bins_];
+    calc_zn(order_, r, theta, zn);
+    for (int i = 0; i < n_bins_; i++) {
+      //TODO: off-by-one
+      match.bins_.push_back(i+1);
+      match.weights_.push_back(zn[i]);
+    }
+  }
+}
+
+void
+ZernikeFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  write_dataset(filter_group, "order", order_);
+  write_dataset(filter_group, "x", x_);
+  write_dataset(filter_group, "y", y_);
+  write_dataset(filter_group, "r", r_);
+}
+
+std::string
+ZernikeFilter::text_label(int bin) const
+{
+  std::stringstream out;
+  for (int n = 0; n < order_+1; n++) {
+    int last = (n + 1) * (n + 2) / 2;
+    if (bin <= last) {
+      int first = last - n;
+      int m = -n + (bin - first) * 2;
+      out << "Zernike expansion, Z" << n << "," << m;
+      return out.str();
+    }
+  }
+}
+
+void
+ZernikeFilter::set_order(int order)
+{
+  order_ = order;
+  n_bins_ = ((order+1) * (order+2)) / 2;
+}
+
+//==============================================================================
+// ZernikeRadialFilter implementation
+//==============================================================================
+
+void
+ZernikeRadialFilter::get_all_bins(const Particle* p, int estimator,
+                                  FilterMatch& match) const
+{
+  // Determine the normalized radius coordinate.
+  double x = p->coord[0].xyz[0] - x_;
+  double y = p->coord[0].xyz[1] - y_;
+  double r = std::sqrt(x*x + y*y) / r_;
+
+  if (r <= 1.0) {
+    // Compute and return the Zernike weights.
+    double zn[n_bins_];
+    calc_zn_rad(order_, r, zn);
+    for (int i = 0; i < n_bins_; i++) {
+      //TODO: off-by-one
+      match.bins_.push_back(i+1);
+      match.weights_.push_back(zn[i]);
+    }
+  }
+}
+
+std::string
+ZernikeRadialFilter::text_label(int bin) const
+{
+  return "Zernike expansion, Z" + std::to_string(2*(bin-1)) + ",0";
+}
+
+void
+ZernikeRadialFilter::set_order(int order)
+{
+  order_ = order;
+  n_bins_ = order / 2 + 1;
+}
+
+//==============================================================================
+// C-API functions
+//==============================================================================
+
+std::pair<int, ZernikeFilter*>
+check_zernike_filter(int32_t index)
+{
+  // Make sure this is a valid index to an allocated filter.
+  int err = verify_filter(index);
+  if (err) {
+    return {err, nullptr};
+  }
+
+  // Get a pointer to the filter and downcast.
+  auto* filt_base = filter_from_f(index);
+  auto* filt = dynamic_cast<ZernikeFilter*>(filt_base);
+
+  // Check the filter type.
+  if (!filt) {
+    set_errmsg("Not a Zernike filter.");
+    err = OPENMC_E_INVALID_TYPE;
+  }
+  return {err, filt};
+}
+
+extern "C" int
+openmc_zernike_filter_get_order(int32_t index, int* order)
+{
+  // Check the filter.
+  auto check_result = check_zernike_filter(index);
+  auto err = check_result.first;
+  auto filt = check_result.second;
+  if (err) return err;
+
+  // Output the order.
+  *order = filt->order();
+  return 0;
+}
+
+extern "C" int
+openmc_zernike_filter_get_params(int32_t index, double* x, double* y,
+                                 double* r)
+{
+  // Check the filter.
+  auto check_result = check_zernike_filter(index);
+  auto err = check_result.first;
+  auto filt = check_result.second;
+  if (err) return err;
+
+  // Output the params.
+  *x = filt->x_;
+  *y = filt->y_;
+  *r = filt->r_;
+  return 0;
+}
+
+extern "C" int
+openmc_zernike_filter_set_order(int32_t index, int order)
+{
+  // Check the filter.
+  auto check_result = check_zernike_filter(index);
+  auto err = check_result.first;
+  auto filt = check_result.second;
+  if (err) return err;
+
+  // Update the filter.
+  filt->set_order(order);
+  filter_update_n_bins(index);
+  return 0;
+}
+
+extern "C" int
+openmc_zernike_filter_set_params(int32_t index, const double* x,
+                                 const double* y, const double* r)
+{
+  // Check the filter.
+  auto check_result = check_zernike_filter(index);
+  auto err = check_result.first;
+  auto filt = check_result.second;
+  if (err) return err;
+
+  // Update the filter.
+  if (x) filt->x_ = *x;
+  if (y) filt->y_ = *y;
+  if (r) filt->r_ = *r;
+  return 0;
+}
+
+} // namespace openmc

src/tallies/tally.F90

@@ -2079,15 +2079,15 @@ contains
       do j = 1, size(t % filter)
         i_filt = t % filter(j)
         if (.not. filter_matches(i_filt) % bins_present) then
-          call filter_matches(i_filt) % bins % clear()
-          call filter_matches(i_filt) % weights % clear()
+          call filter_matches(i_filt) % bins_clear()
+          call filter_matches(i_filt) % weights_clear()
           call filters(i_filt) % obj % get_all_bins(p, t % estimator, &
                filter_matches(i_filt))
           filter_matches(i_filt) % bins_present = .true.
         end if
         ! If there are no valid bins for this filter, then there is nothing to
         ! score and we can move on to the next tally.
-        if (filter_matches(i_filt) % bins % size() == 0) cycle TALLY_LOOP
+        if (filter_matches(i_filt) % bins_size() == 0) cycle TALLY_LOOP
 
         ! Set the index of the bin used in the first filter combination
         filter_matches(i_filt) % i_bin = 1
@@ -2106,10 +2106,10 @@ contains
         do j = 1, size(t % filter)
           i_filt = t % filter(j)
           i_bin = filter_matches(i_filt) % i_bin
-          filter_index = filter_index + (filter_matches(i_filt) % bins % &
-               data(i_bin) - 1) * t % stride(j)
-          filter_weight = filter_weight * filter_matches(i_filt) % weights % &
-               data(i_bin)
+          filter_index = filter_index + (filter_matches(i_filt) &
+               % bins_data(i_bin) - 1) * t % stride(j)
+          filter_weight = filter_weight * filter_matches(i_filt) &
+               % weights_data(i_bin)
         end do
 
         ! ======================================================================
@@ -2166,7 +2166,7 @@ contains
         do j = size(t % filter), 1, -1
           i_filt = t % filter(j)
           if (filter_matches(i_filt) % i_bin < filter_matches(i_filt) % &
-               bins % size()) then
+               bins_size()) then
             filter_matches(i_filt) % i_bin = filter_matches(i_filt) % i_bin + 1
             finished = .false.
             exit
@@ -2225,15 +2225,15 @@ contains
       do j = 1, size(t % filter)
         i_filt = t % filter(j)
         if (.not. filter_matches(i_filt) % bins_present) then
-          call filter_matches(i_filt) % bins % clear()
-          call filter_matches(i_filt) % weights % clear()
+          call filter_matches(i_filt) % bins_clear()
+          call filter_matches(i_filt) % weights_clear()
           call filters(i_filt) % obj % get_all_bins(p, t % estimator, &
                filter_matches(i_filt))
           filter_matches(i_filt) % bins_present = .true.
         end if
         ! If there are no valid bins for this filter, then there is nothing to
         ! score and we can move on to the next tally.
-        if (filter_matches(i_filt) % bins % size() == 0) cycle TALLY_LOOP
+        if (filter_matches(i_filt) % bins_size() == 0) cycle TALLY_LOOP
 
         ! Set the index of the bin used in the first filter combination
         filter_matches(i_filt) % i_bin = 1
@@ -2252,10 +2252,10 @@ contains
         do j = 1, size(t % filter)
           i_filt = t % filter(j)
           i_bin = filter_matches(i_filt) % i_bin
-          filter_index = filter_index + (filter_matches(i_filt) % bins % &
-               data(i_bin) - 1) * t % stride(j)
-          filter_weight = filter_weight * filter_matches(i_filt) % weights % &
-               data(i_bin)
+          filter_index = filter_index + (filter_matches(i_filt) &
+               % bins_data(i_bin) - 1) * t % stride(j)
+          filter_weight = filter_weight * filter_matches(i_filt) &
+               % weights_data(i_bin)
         end do
 
         ! ======================================================================
@@ -2295,7 +2295,7 @@ contains
         do j = size(t % filter), 1, -1
           i_filt = t % filter(j)
           if (filter_matches(i_filt) % i_bin < filter_matches(i_filt) % &
-               bins % size()) then
+               bins_size()) then
             filter_matches(i_filt) % i_bin = filter_matches(i_filt) % i_bin + 1
             finished = .false.
             exit
@@ -2360,7 +2360,7 @@ contains
     ! save original outgoing energy bin and score index
     i = t % filter(t % find_filter(FILTER_ENERGYOUT))
     i_bin = filter_matches(i) % i_bin
-    bin_energyout = filter_matches(i) % bins % data(i_bin)
+    bin_energyout = filter_matches(i) % bins_data(i_bin)
 
     ! declare the energyout filter type
     select type(eo_filt => filters(i) % obj)
@@ -2400,7 +2400,7 @@ contains
           g_out = size(eo_filt % bins) - g_out
 
           ! change outgoing energy bin
-          filter_matches(i) % bins % data(i_bin) = g_out
+          call filter_matches(i) % bins_set_data(i_bin, g_out)
 
         else
 
@@ -2416,8 +2416,8 @@ contains
           if (E_out < eo_filt % bins(1) .or. E_out > eo_filt % bins(n)) cycle
 
           ! change outgoing energy bin
-          filter_matches(i) % bins % data(i_bin) = &
-               binary_search(eo_filt % bins, n, E_out)
+          call filter_matches(i) % bins_set_data(i_bin, &
+               binary_search(eo_filt % bins, n, E_out))
 
         end if
 
@@ -2428,8 +2428,8 @@ contains
           ! determine scoring index and weight for this filter combination
           i_filter = 1
           do l = 1, size(t % filter)
-            i_filter = i_filter + (filter_matches(t % filter(l)) % bins % &
-                 data(filter_matches(t % filter(l)) % i_bin) - 1) * &
+            i_filter = i_filter + (filter_matches(t % filter(l)) &
+                 % bins_data(filter_matches(t % filter(l)) % i_bin) - 1) * &
                  t % stride(l)
           end do
 
@@ -2470,10 +2470,10 @@ contains
                   do l = 1, size(t % filter)
                     f = t % filter(l)
                     b = filter_matches(f) % i_bin
-                    i_filter = i_filter + (filter_matches(f) % bins % &
-                         data(b) - 1) * t % stride(l)
-                    filter_weight = filter_weight * filter_matches(f) % &
-                         weights % data(b)
+                    i_filter = i_filter + (filter_matches(f) &
+                         % bins_data(b) - 1) * t % stride(l)
+                    filter_weight = filter_weight * filter_matches(f) &
+                         % weights_data(b)
                   end do
 
                   call score_fission_delayed_dg(t, d_bin, &
@@ -2493,10 +2493,10 @@ contains
             do l = 1, size(t % filter)
               f = t % filter(l)
               b = filter_matches(f) % i_bin
-              i_filter = i_filter + (filter_matches(f) % bins % data(b) - 1) &
+              i_filter = i_filter + (filter_matches(f) % bins_data(b) - 1) &
                    * t % stride(l)
-              filter_weight = filter_weight * filter_matches(f) % weights % &
-                   data(b)
+              filter_weight = filter_weight * filter_matches(f) &
+                   % weights_data(b)
             end do
 
             ! Add score to tally
@@ -2509,7 +2509,7 @@ contains
     end select
 
     ! reset outgoing energy bin and score index
-    filter_matches(i) % bins % data(i_bin) = bin_energyout
+    call filter_matches(i) % bins_set_data(i_bin, bin_energyout)
 
   end subroutine score_fission_eout
 
@@ -2534,14 +2534,14 @@ contains
     ! save original delayed group bin
     i_filt = t % filter(t % find_filter(FILTER_DELAYEDGROUP))
     i_bin = filter_matches(i_filt) % i_bin
-    bin_original = filter_matches(i_filt) % bins % data(i_bin)
-    filter_matches(i_filt) % bins % data(i_bin) = d_bin
+    bin_original = filter_matches(i_filt) % bins_data(i_bin)
+    call filter_matches(i_filt) % bins_set_data(i_bin, d_bin)
 
     ! determine scoring index and weight on the modified matching bins
     filter_index = 1
     do i = 1, size(t % filter)
-      filter_index = filter_index + (filter_matches(t % filter(i)) % bins % &
-           data(filter_matches(t % filter(i)) % i_bin) - 1) * t % stride(i)
+      filter_index = filter_index + (filter_matches(t % filter(i)) % &
+           bins_data(filter_matches(t % filter(i)) % i_bin) - 1) * t % stride(i)
     end do
 
 !$omp atomic
@@ -2549,7 +2549,7 @@ contains
          t % results(RESULT_VALUE, score_index, filter_index) + score
 
     ! reset original delayed group bin
-    filter_matches(i_filt) % bins % data(i_bin) = bin_original
+    call filter_matches(i_filt) % bins_set_data(i_bin, bin_original)
 
   end subroutine score_fission_delayed_dg
 
@@ -2595,15 +2595,15 @@ contains
       do j = 1, size(t % filter)
         i_filt = t % filter(j)
         if (.not. filter_matches(i_filt) % bins_present) then
-          call filter_matches(i_filt) % bins % clear()
-          call filter_matches(i_filt) % weights % clear()
+          call filter_matches(i_filt) % bins_clear()
+          call filter_matches(i_filt) % weights_clear()
           call filters(i_filt) % obj % get_all_bins(p, t % estimator, &
                filter_matches(i_filt))
           filter_matches(i_filt) % bins_present = .true.
         end if
         ! If there are no valid bins for this filter, then there is nothing to
         ! score and we can move on to the next tally.
-        if (filter_matches(i_filt) % bins % size() == 0) cycle TALLY_LOOP
+        if (filter_matches(i_filt) % bins_size() == 0) cycle TALLY_LOOP
 
         ! Set the index of the bin used in the first filter combination
         filter_matches(i_filt) % i_bin = 1
@@ -2622,10 +2622,10 @@ contains
         do j = 1, size(t % filter)
           i_filt = t % filter(j)
           i_bin = filter_matches(i_filt) % i_bin
-          filter_index = filter_index + (filter_matches(i_filt) % bins % &
-               data(i_bin) - 1) * t % stride(j)
-          filter_weight = filter_weight * filter_matches(i_filt) % weights % &
-               data(i_bin)
+          filter_index = filter_index + (filter_matches(i_filt) &
+               % bins_data(i_bin) - 1) * t % stride(j)
+          filter_weight = filter_weight * filter_matches(i_filt) &
+               % weights_data(i_bin)
         end do
 
         ! ======================================================================
@@ -2674,7 +2674,7 @@ contains
         do j = size(t % filter), 1, -1
           i_filt = t % filter(j)
           if (filter_matches(i_filt) % i_bin < filter_matches(i_filt) % &
-               bins % size()) then
+               bins_size()) then
             filter_matches(i_filt) % i_bin = filter_matches(i_filt) % i_bin + 1
             finished = .false.
             exit
@@ -2752,15 +2752,15 @@ contains
       do j = 1, size(t % filter)
         i_filt = t % filter(j)
         if (.not. filter_matches(i_filt) % bins_present) then
-          call filter_matches(i_filt) % bins % clear()
-          call filter_matches(i_filt) % weights % clear()
+          call filter_matches(i_filt) % bins_clear()
+          call filter_matches(i_filt) % weights_clear()
           call filters(i_filt) % obj % get_all_bins(p, t % estimator, &
                filter_matches(i_filt))
           filter_matches(i_filt) % bins_present = .true.
         end if
         ! If there are no valid bins for this filter, then there is nothing to
         ! score and we can move on to the next tally.
-        if (filter_matches(i_filt) % bins % size() == 0) cycle TALLY_LOOP
+        if (filter_matches(i_filt) % bins_size() == 0) cycle TALLY_LOOP
 
         ! Set the index of the bin used in the first filter combination
         filter_matches(i_filt) % i_bin = 1
@@ -2779,10 +2779,10 @@ contains
         do j = 1, size(t % filter)
           i_filt = t % filter(j)
           i_bin = filter_matches(i_filt) % i_bin
-          filter_index = filter_index + (filter_matches(i_filt) % bins % &
-               data(i_bin) - 1) * t % stride(j)
-          filter_weight = filter_weight * filter_matches(i_filt) % weights % &
-               data(i_bin)
+          filter_index = filter_index + (filter_matches(i_filt) &
+               % bins_data(i_bin) - 1) * t % stride(j)
+          filter_weight = filter_weight * filter_matches(i_filt) &
+               % weights_data(i_bin)
         end do
 
         ! ======================================================================
@@ -2831,7 +2831,7 @@ contains
         do j = size(t % filter), 1, -1
           i_filt = t % filter(j)
           if (filter_matches(i_filt) % i_bin < filter_matches(i_filt) % &
-               bins % size()) then
+               bins_size()) then
             filter_matches(i_filt) % i_bin = filter_matches(i_filt) % i_bin + 1
             finished = .false.
             exit
@@ -2898,15 +2898,15 @@ contains
       do j = 1, size(t % filter)
         i_filt = t % filter(j)
         if (.not. filter_matches(i_filt) % bins_present) then
-          call filter_matches(i_filt) % bins % clear()
-          call filter_matches(i_filt) % weights % clear()
+          call filter_matches(i_filt) % bins_clear()
+          call filter_matches(i_filt) % weights_clear()
           call filters(i_filt) % obj % get_all_bins(p, t % estimator, &
                filter_matches(i_filt))
           filter_matches(i_filt) % bins_present = .true.
         end if
         ! If there are no valid bins for this filter, then there is nothing to
         ! score and we can move on to the next tally.
-        if (filter_matches(i_filt) % bins % size() == 0) cycle TALLY_LOOP
+        if (filter_matches(i_filt) % bins_size() == 0) cycle TALLY_LOOP
 
         ! Set the index of the bin used in the first filter combination
         filter_matches(i_filt) % i_bin = 1
@@ -2925,10 +2925,10 @@ contains
         do j = 1, size(t % filter)
           i_filt = t % filter(j)
           i_bin = filter_matches(i_filt) % i_bin
-          filter_index = filter_index + (filter_matches(i_filt) % bins % &
-               data(i_bin) - 1) * t % stride(j)
-          filter_weight = filter_weight * filter_matches(i_filt) % weights % &
-               data(i_bin)
+          filter_index = filter_index + (filter_matches(i_filt) &
+               % bins_data(i_bin) - 1) * t % stride(j)
+          filter_weight = filter_weight * filter_matches(i_filt) &
+               % weights_data(i_bin)
         end do
 
         ! Determine score
@@ -2961,7 +2961,7 @@ contains
         do j = size(t % filter), 1, -1
           i_filt = t % filter(j)
           if (filter_matches(i_filt) % i_bin < filter_matches(i_filt) % &
-               bins % size()) then
+               bins_size()) then
             filter_matches(i_filt) % i_bin = filter_matches(i_filt) % i_bin + 1
             finished = .false.
             exit

src/tallies/tally_filter.F90

@@ -7,26 +7,15 @@ module tally_filter
   use tally_filter_header
 
   ! Inherit other filters
-  use tally_filter_azimuthal
-  use tally_filter_cell
-  use tally_filter_cellborn
-  use tally_filter_cellfrom
+  use tally_filter_cpp
   use tally_filter_delayedgroup
   use tally_filter_distribcell
   use tally_filter_energy
-  use tally_filter_energyfunc
   use tally_filter_legendre
-  use tally_filter_material
   use tally_filter_mesh
   use tally_filter_meshsurface
-  use tally_filter_mu
   use tally_filter_particle
-  use tally_filter_polar
   use tally_filter_sph_harm
-  use tally_filter_sptl_legendre
-  use tally_filter_surface
-  use tally_filter_universe
-  use tally_filter_zernike
 
   implicit none
 
@@ -114,6 +103,14 @@ contains
 
     character(:), allocatable :: type_
 
+    interface
+      function allocate_filter(type) result(ptr) bind(C)
+        import C_CHAR, C_PTR
+        character(kind=C_CHAR), intent(in)        :: type(*)
+        type(C_PTR)                               :: ptr
+      end function allocate_filter
+    end interface
+
     ! Convert C string to Fortran string
     type_ = to_f_string(type)
 
@@ -131,7 +128,7 @@ contains
         case ('cellborn')
           allocate(CellbornFilter :: filters(index) % obj)
         case ('cellfrom')
-          allocate(CellfromFilter :: filters(index) % obj)
+          allocate(CellFromFilter :: filters(index) % obj)
         case ('delayedgroup')
           allocate(DelayedGroupFilter :: filters(index) % obj)
         case ('distribcell')
@@ -172,6 +169,16 @@ contains
           err = E_UNASSIGNED
           call set_errmsg("Unknown filter type: " // trim(type_))
         end select
+
+        select type(filt => filters(index) % obj)
+        class is (CppTallyFilter)
+          filt % ptr = allocate_filter(type)
+          if (.not. c_associated(filt % ptr)) then
+            err = E_UNASSIGNED
+            call set_errmsg("Could not allocate C++ tally filter")
+          end if
+        end select
+
       end if
     else
       err = E_OUT_OF_BOUNDS

src/tallies/tally_filter_azimuthal.F90

@@ -1,120 +0,0 @@
-module tally_filter_azimuthal
-
-  use, intrinsic :: ISO_C_BINDING
-
-  use algorithm,          only: binary_search
-  use constants
-  use error,              only: fatal_error
-  use hdf5_interface
-  use particle_header,    only: Particle
-  use string,             only: to_str
-  use tally_filter_header
-  use xml_interface
-
-  implicit none
-  private
-
-!===============================================================================
-! AZIMUTHALFILTER bins the incident neutron azimuthal angle (relative to the
-! global xy-plane).
-!===============================================================================
-
-  type, public, extends(TallyFilter) :: AzimuthalFilter
-    real(8), allocatable :: bins(:)
-  contains
-    procedure :: from_xml
-    procedure :: get_all_bins => get_all_bins_azimuthal
-    procedure :: to_statepoint => to_statepoint_azimuthal
-    procedure :: text_label => text_label_azimuthal
-  end type AzimuthalFilter
-
-contains
-
-  subroutine from_xml(this, node)
-    class(AzimuthalFilter), intent(inout) :: this
-    type(XMLNode), intent(in) :: node
-
-    integer :: i
-    integer :: n
-    integer :: n_angle
-    real(8) :: d_angle
-
-    n = node_word_count(node, "bins")
-
-    ! Allocate and store bins
-    this % n_bins = n - 1
-    allocate(this % bins(n))
-    call get_node_array(node, "bins", this % bins)
-
-    ! Allow a user to input a lone number which will mean that you
-    ! subdivide [-pi,pi) evenly with the input being the number of
-    ! bins
-    if (n == 1) then
-      n_angle = int(this % bins(1))
-      if (n_angle > 1) then
-        this % n_bins = n_angle
-        d_angle = TWO * PI / n_angle
-        deallocate(this % bins)
-        allocate(this % bins(n_angle + 1))
-        do i = 1, n_angle
-          this % bins(i) = -PI + (i - 1) * d_angle
-        end do
-        this % bins(n_angle + 1) = PI
-      else
-        call fatal_error("Number of bins for azimuthal filter must be&
-             & greater than 1.")
-      end if
-    end if
-  end subroutine from_xml
-
-  subroutine get_all_bins_azimuthal(this, p, estimator, match)
-    class(AzimuthalFilter), intent(in)  :: this
-    type(Particle),         intent(in)  :: p
-    integer,                intent(in)  :: estimator
-    type(TallyFilterMatch),      intent(inout) :: match
-
-    integer :: n
-    integer :: bin
-    real(8) :: phi
-
-    n = this % n_bins
-
-    ! Make sure the correct direction vector is used.
-    if (estimator == ESTIMATOR_TRACKLENGTH) then
-      phi = atan2(p % coord(1) % uvw(2), p % coord(1) % uvw(1))
-    else
-      phi = atan2(p % last_uvw(2), p % last_uvw(1))
-    end if
-
-    ! Search to find azimuthal angle bin.
-    bin = binary_search(this % bins, n + 1, phi)
-    if (bin /= NO_BIN_FOUND) then
-      call match % bins % push_back(bin)
-      call match % weights % push_back(ONE)
-    end if
-
-  end subroutine get_all_bins_azimuthal
-
-  subroutine to_statepoint_azimuthal(this, filter_group)
-    class(AzimuthalFilter), intent(in) :: this
-    integer(HID_T),         intent(in) :: filter_group
-
-    call write_dataset(filter_group, "type", "azimuthal")
-    call write_dataset(filter_group, "n_bins", this % n_bins)
-    call write_dataset(filter_group, "bins", this % bins)
-  end subroutine to_statepoint_azimuthal
-
-  function text_label_azimuthal(this, bin) result(label)
-    class(AzimuthalFilter), intent(in) :: this
-    integer,                intent(in) :: bin
-    character(MAX_LINE_LEN)            :: label
-
-    real(8) :: E0, E1
-
-    E0 = this % bins(bin)
-    E1 = this % bins(bin + 1)
-    label = "Azimuthal Angle [" // trim(to_str(E0)) // ", " &
-         // trim(to_str(E1)) // ")"
-  end function text_label_azimuthal
-
-end module tally_filter_azimuthal

src/tallies/tally_filter_cell.F90

@@ -1,144 +0,0 @@
-module tally_filter_cell
-
-  use, intrinsic :: ISO_C_BINDING
-
-  use constants,          only: ONE, MAX_LINE_LEN
-  use dict_header,        only: EMPTY
-  use error,              only: fatal_error
-  use hdf5_interface
-  use geometry_header
-  use particle_header,    only: Particle
-  use string,             only: to_str
-  use tally_filter_header
-  use xml_interface
-
-  implicit none
-  private
-  public :: openmc_cell_filter_get_bins
-
-!===============================================================================
-! CELLFILTER specifies which geometric cells tally events reside in.
-!===============================================================================
-
-  type, public, extends(TallyFilter) :: CellFilter
-    integer(C_INT32_T), allocatable :: cells(:)
-    type(DictIntInt)     :: map
-  contains
-    procedure :: from_xml
-    procedure :: get_all_bins => get_all_bins_cell
-    procedure :: to_statepoint => to_statepoint_cell
-    procedure :: text_label => text_label_cell
-    procedure :: initialize => initialize_cell
-  end type CellFilter
-
-contains
-
-  subroutine from_xml(this, node)
-    class(CellFilter), intent(inout) :: this
-    type(XMLNode), intent(in) :: node
-
-    integer :: n
-
-    ! Determine how many bins were given
-    n = node_word_count(node, "bins")
-
-    ! Allocate and store bins
-    this % n_bins = n
-    allocate(this % cells(n))
-    call get_node_array(node, "bins", this % cells)
-  end subroutine from_xml
-
-  subroutine get_all_bins_cell(this, p, estimator, match)
-    class(CellFilter), intent(in)  :: this
-    type(Particle),    intent(in)  :: p
-    integer,           intent(in)  :: estimator
-    type(TallyFilterMatch), intent(inout) :: match
-
-    integer :: i
-    integer :: val
-
-    ! Iterate over coordinate levels to see with cells match
-    do i = 1, p % n_coord
-      val = this % map % get(p % coord(i) % cell + 1)
-      if (val /= EMPTY) then
-        call match % bins % push_back(val)
-        call match % weights % push_back(ONE)
-      end if
-    end do
-
-  end subroutine get_all_bins_cell
-
-  subroutine to_statepoint_cell(this, filter_group)
-    class(CellFilter), intent(in) :: this
-    integer(HID_T),    intent(in) :: filter_group
-
-    integer :: i
-    integer, allocatable :: cell_ids(:)
-
-    call write_dataset(filter_group, "type", "cell")
-    call write_dataset(filter_group, "n_bins", this % n_bins)
-
-    allocate(cell_ids(size(this % cells)))
-    do i = 1, size(this % cells)
-      cell_ids(i) = cells(this % cells(i)) % id()
-    end do
-    call write_dataset(filter_group, "bins", cell_ids)
-  end subroutine to_statepoint_cell
-
-  subroutine initialize_cell(this)
-    class(CellFilter), intent(inout) :: this
-
-    integer :: i, id
-    integer :: val
-
-    ! Convert ids to indices.
-    do i = 1, this % n_bins
-      id = this % cells(i)
-      val = cell_dict % get(id)
-      if (val /= EMPTY) then
-        this % cells(i) = val
-      else
-        call fatal_error("Could not find cell " // trim(to_str(id)) &
-             &// " specified on tally filter.")
-      end if
-    end do
-
-    ! Generate mapping from cell indices to filter bins.
-    do i = 1, this % n_bins
-      call this % map % set(this % cells(i), i)
-    end do
-  end subroutine initialize_cell
-
-  function text_label_cell(this, bin) result(label)
-    class(CellFilter), intent(in) :: this
-    integer,           intent(in) :: bin
-    character(MAX_LINE_LEN)       :: label
-
-    label = "Cell " // to_str(cells(this % cells(bin)) % id())
-  end function text_label_cell
-
-!===============================================================================
-!                               C API FUNCTIONS
-!===============================================================================
-
-  function openmc_cell_filter_get_bins(index, cells, n) result(err) bind(C)
-    ! Return the cells associated with a cell filter
-    integer(C_INT32_T), value :: index
-    type(C_PTR), intent(out) :: cells
-    integer(C_INT32_T), intent(out) :: n
-    integer(C_INT) :: err
-
-    err = verify_filter(index)
-    if (err == 0) then
-       select type (f => filters(index) % obj)
-       type is (CellFilter)
-          cells = C_LOC(f % cells)
-          n = size(f % cells)
-       class default
-          err = E_INVALID_TYPE
-          call set_errmsg("Tried to get cells from a non-cell filter.")
-       end select
-    end if
-  end function openmc_cell_filter_get_bins
-
-end module tally_filter_cell

src/tallies/tally_filter_cellborn.F90

@@ -1,114 +0,0 @@
-module tally_filter_cellborn
-
-  use, intrinsic :: ISO_C_BINDING
-
-  use constants,          only: ONE, MAX_LINE_LEN
-  use dict_header,        only: EMPTY
-  use error,              only: fatal_error
-  use hdf5_interface
-  use geometry_header
-  use particle_header,    only: Particle
-  use string,             only: to_str
-  use tally_filter_header
-  use xml_interface
-
-  implicit none
-  private
-
-!===============================================================================
-! CELLBORNFILTER specifies which cell the particle was born in.
-!===============================================================================
-
-  type, public, extends(TallyFilter) :: CellbornFilter
-    integer, allocatable :: cells(:)
-    type(DictIntInt)     :: map
-  contains
-    procedure :: from_xml
-    procedure :: get_all_bins => get_all_bins_cellborn
-    procedure :: to_statepoint => to_statepoint_cellborn
-    procedure :: text_label => text_label_cellborn
-    procedure :: initialize => initialize_cellborn
-  end type CellbornFilter
-
-contains
-
-  subroutine from_xml(this, node)
-    class(CellbornFilter), intent(inout) :: this
-    type(XMLNode), intent(in) :: node
-
-    integer :: n
-
-    ! Determine number of bins
-    n = node_word_count(node, "bins")
-
-    ! Allocate and store bins
-    this % n_bins = n
-    allocate(this % cells(n))
-    call get_node_array(node, "bins", this % cells)
-  end subroutine from_xml
-
-  subroutine get_all_bins_cellborn(this, p, estimator, match)
-    class(CellbornFilter), intent(in)  :: this
-    type(Particle),        intent(in)  :: p
-    integer,               intent(in)  :: estimator
-    type(TallyFilterMatch),     intent(inout) :: match
-
-    integer :: val
-
-      val = this % map % get(p % cell_born + 1)
-      if (val /= EMPTY) then
-        call match % bins % push_back(val)
-        call match % weights % push_back(ONE)
-      end if
-
-  end subroutine get_all_bins_cellborn
-
-  subroutine to_statepoint_cellborn(this, filter_group)
-    class(CellbornFilter), intent(in) :: this
-    integer(HID_T),        intent(in) :: filter_group
-
-    integer :: i
-    integer, allocatable :: cell_ids(:)
-
-    call write_dataset(filter_group, "type", "cellborn")
-    call write_dataset(filter_group, "n_bins", this % n_bins)
-    allocate(cell_ids(size(this % cells)))
-    do i = 1, size(this % cells)
-      cell_ids(i) = cells(this % cells(i)) % id()
-    end do
-    call write_dataset(filter_group, "bins", cell_ids)
-  end subroutine to_statepoint_cellborn
-
-  subroutine initialize_cellborn(this)
-    class(CellbornFilter), intent(inout) :: this
-
-    integer :: i, id
-    integer :: val
-
-    ! Convert ids to indices.
-    do i = 1, this % n_bins
-      id = this % cells(i)
-      val = cell_dict % get(id)
-      if (val /= EMPTY) then
-        this % cells(i) = val
-      else
-        call fatal_error("Could not find cell " // trim(to_str(id)) &
-             &// " specified on tally filter.")
-      end if
-    end do
-
-    ! Generate mapping from cell indices to filter bins.
-    do i = 1, this % n_bins
-      call this % map % set(this % cells(i), i)
-    end do
-  end subroutine initialize_cellborn
-
-  function text_label_cellborn(this, bin) result(label)
-    class(CellbornFilter), intent(in) :: this
-    integer,               intent(in) :: bin
-    character(MAX_LINE_LEN)           :: label
-
-    label = "Birth Cell " // to_str(cells(this % cells(bin)) % id())
-  end function text_label_cellborn
-
-end module tally_filter_cellborn

src/tallies/tally_filter_cellfrom.F90

@@ -1,80 +0,0 @@
-module tally_filter_cellfrom
-
-  use, intrinsic :: ISO_C_BINDING
-
-  use constants,          only: ONE, MAX_LINE_LEN
-  use dict_header,        only: EMPTY
-  use error,              only: fatal_error
-  use hdf5_interface
-  use geometry_header
-  use particle_header,    only: Particle
-  use string,             only: to_str
-  use tally_filter_header
-  use tally_filter_cell
-  use xml_interface
-
-  implicit none
-  private
-
-!===============================================================================
-! CELLFROMFILTER specifies which geometric cells particles exit when crossing a
-! surface.
-!===============================================================================
-
-  type, public, extends(CellFilter) :: CellFromFilter
-  contains
-    ! Inherit from_xml from CellFilter
-    procedure :: get_all_bins => get_all_bins_cell_from
-    procedure :: to_statepoint => to_statepoint_cell_from
-    procedure :: text_label => text_label_cell_from
-  end type CellFromFilter
-
-contains
-
-  subroutine get_all_bins_cell_from(this, p, estimator, match)
-    class(CellFromFilter), intent(in)  :: this
-    type(Particle),    intent(in)  :: p
-    integer,           intent(in)  :: estimator
-    type(TallyFilterMatch), intent(inout) :: match
-
-    integer :: i
-    integer :: val
-
-    ! Starting one coordinate level deeper, find the next bin.
-    do i = 1, p % last_n_coord
-      val = this % map % get(p % last_cell(i) + 1)
-      if (val /= EMPTY) then
-        call match % bins % push_back(val)
-        call match % weights % push_back(ONE)
-        exit
-      end if
-    end do
-
-  end subroutine get_all_bins_cell_from
-
-  subroutine to_statepoint_cell_from(this, filter_group)
-    class(CellFromFilter), intent(in) :: this
-    integer(HID_T),        intent(in) :: filter_group
-
-    integer :: i
-    integer, allocatable :: cell_ids(:)
-
-    call write_dataset(filter_group, "type", "cellfrom")
-    call write_dataset(filter_group, "n_bins", this % n_bins)
-
-    allocate(cell_ids(size(this % cells)))
-    do i = 1, size(this % cells)
-      cell_ids(i) = cells(this % cells(i)) % id()
-    end do
-    call write_dataset(filter_group, "bins", cell_ids)
-  end subroutine to_statepoint_cell_from
-
-  function text_label_cell_from(this, bin) result(label)
-    class(CellFromFilter), intent(in) :: this
-    integer,               intent(in) :: bin
-    character(MAX_LINE_LEN)           :: label
-
-    label = "Cell from " // to_str(cells(this % cells(bin)) % id())
-  end function text_label_cell_from
-
-end module tally_filter_cellfrom

src/tallies/tally_filter_cpp.F90

@@ -0,0 +1,235 @@
+module tally_filter_cpp
+  use, intrinsic :: ISO_C_BINDING
+
+  use constants,       only: MAX_LINE_LEN
+  use hdf5_interface,  only: HID_T
+  use particle_header, only: Particle
+  use tally_filter_header
+  use xml_interface,   only: XMLNode
+
+  implicit none
+
+!===============================================================================
+! CPPTALLYFILTER is a TallyFilter that is at least partially implemented in C++.
+!
+! Moving a tally to C++ is easier if done piece-by-piece, so this filter has
+! *_cpp_inner procedures that allows the Fortran code to override e.g.
+! get_all_bins and still call the C++ implementation for side-by-side comparison
+!===============================================================================
+
+  type, abstract, extends(TallyFilter) :: CppTallyFilter
+    type(C_PTR) :: ptr
+  contains
+    procedure :: n_bins_cpp
+    procedure :: from_xml_cpp_inner
+    procedure :: get_all_bins_cpp_inner
+    procedure :: to_statepoint_cpp_inner
+    procedure :: text_label_cpp_inner
+    procedure :: initialize_cpp_inner
+    procedure :: from_xml => from_xml_cpp_default
+    procedure :: get_all_bins => get_all_bins_cpp_default
+    procedure :: to_statepoint => to_statepoint_cpp_default
+    procedure :: text_label => text_label_cpp_default
+    procedure :: initialize => initialize_cpp_default
+  end type CppTallyFilter
+
+!===============================================================================
+! Pure C++ filters
+!===============================================================================
+
+  type, extends(CppTallyFilter) :: AzimuthalFilter
+  end type
+
+  type, extends(CppTallyFilter) :: CellFilter
+  end type
+
+  type, extends(CppTallyFilter) :: CellbornFilter
+  end type
+
+  type, extends(CellFilter) :: CellFromFilter
+  end type
+
+  type, extends(CppTallyFilter) :: EnergyFunctionFilter
+  end type
+
+  type, extends(CppTallyFilter) :: MaterialFilter
+  end type
+
+  type, extends(CppTallyFilter) :: MuFilter
+  end type
+
+  type, extends(CppTallyFilter) :: PolarFilter
+  end type
+
+  type, extends(CppTallyFilter) :: SpatialLegendreFilter
+  end type
+
+  type, extends(CppTallyFilter) :: SurfaceFilter
+    ! True if this filter is used for surface currents
+    logical :: current = .false.
+  end type
+
+  type, extends(CppTallyFilter) :: UniverseFilter
+  end type
+
+  type, extends(CppTallyFilter) :: ZernikeFilter
+  end type
+
+  type, extends(ZernikeFilter) :: ZernikeRadialFilter
+  end type
+
+contains
+
+!===============================================================================
+! CppTallyFilter implementation
+!===============================================================================
+
+  function n_bins_cpp(this) result(n_bins)
+    class(CppTallyFilter), intent(inout) :: this
+    integer                              :: n_bins
+    interface
+      function filter_n_bins(filt) result(n_bins) bind(C)
+        import C_PTR, C_INT
+        type(C_PTR), value :: filt
+        integer(C_INT)     :: n_bins
+      end function filter_n_bins
+    end interface
+    n_bins = filter_n_bins(this % ptr)
+  end function n_bins_cpp
+
+  subroutine from_xml_cpp_inner(this, node)
+    class(CppTallyFilter), intent(inout) :: this
+    class(XMLNode),        intent(in)    :: node
+    interface
+      subroutine filter_from_xml(filt, node) bind(C)
+        import C_PTR
+        type(C_PTR), value :: filt
+        type(C_PTR) :: node
+      end subroutine filter_from_xml
+    end interface
+    call filter_from_xml(this % ptr, node % ptr)
+  end subroutine from_xml_cpp_inner
+
+  subroutine get_all_bins_cpp_inner(this, p, estimator, match)
+    class(CppTallyFilter), intent(in) :: this
+    type(Particle),     intent(in)    :: p
+    integer,            intent(in)    :: estimator
+    type(TallyFilterMatch), intent(inout) :: match
+    interface
+      subroutine filter_get_all_bins(filt, p, estimator, match) bind(C)
+        import C_PTR, Particle, C_INT
+        type(C_PTR),                value :: filt
+        type(Particle), intent(in)        :: p
+        integer(C_INT), intent(in), value :: estimator
+        type(C_PTR),                value :: match
+      end subroutine filter_get_all_bins
+    end interface
+    call filter_get_all_bins(this % ptr, p, estimator, match % ptr)
+  end subroutine get_all_bins_cpp_inner
+
+  subroutine to_statepoint_cpp_inner(this, filter_group)
+    class(CppTallyFilter), intent(in) :: this
+    integer(HID_T),        intent(in) :: filter_group
+    interface
+      subroutine filter_to_statepoint(filt, filter_group) bind(C)
+        import C_PTR, HID_T
+        type(C_PTR),                value :: filt
+        integer(HID_T), intent(in), value :: filter_group
+      end subroutine filter_to_statepoint
+    end interface
+    call filter_to_statepoint(this % ptr, filter_group)
+  end subroutine to_statepoint_cpp_inner
+
+  function text_label_cpp_inner(this, bin) result(label)
+    class(CppTallyFilter), intent(in) :: this
+    integer,               intent(in) :: bin
+    character(MAX_LINE_LEN)           :: label
+    character(kind=C_CHAR)            :: label_(MAX_LINE_LEN+1)
+    integer :: i
+    interface
+      subroutine filter_text_label(filt, bin, label) bind(C)
+        import C_PTR, C_INT, C_CHAR
+        type(C_PTR), value     :: filt
+        integer(C_INT), value  :: bin
+        character(kind=C_CHAR) :: label(*)
+      end subroutine filter_text_label
+    end interface
+    call filter_text_label(this % ptr, bin, label_)
+    label = " "
+    do i = 1, MAX_LINE_LEN
+      if (label_(i) == C_NULL_CHAR) exit
+      label(i:i) = label_(i)
+    end do
+  end function text_label_cpp_inner
+
+  subroutine initialize_cpp_inner(this)
+    class(CppTallyFilter), intent(inout) :: this
+    interface
+      subroutine filter_initialize(filt) bind(C)
+        import C_PTR
+        type(C_PTR), value :: filt
+      end subroutine filter_initialize
+    end interface
+    call filter_initialize(this % ptr)
+  end subroutine initialize_cpp_inner
+
+  subroutine from_xml_cpp_default(this, node)
+    class(CppTallyFilter), intent(inout) :: this
+    type(XMLNode),         intent(in)    :: node
+    call this % from_xml_cpp_inner(node)
+    this % n_bins = this % n_bins_cpp()
+  end subroutine from_xml_cpp_default
+
+  subroutine get_all_bins_cpp_default(this, p, estimator, match)
+    class(CppTallyFilter), intent(in) :: this
+    type(Particle),     intent(in)    :: p
+    integer,            intent(in)    :: estimator
+    type(TallyFilterMatch), intent(inout) :: match
+    call this % get_all_bins_cpp_inner(p, estimator, match)
+  end subroutine get_all_bins_cpp_default
+
+  subroutine to_statepoint_cpp_default(this, filter_group)
+    class(CppTallyFilter), intent(in) :: this
+    integer(HID_T),        intent(in) :: filter_group
+    call this % to_statepoint_cpp_inner(filter_group)
+  end subroutine to_statepoint_cpp_default
+
+  function text_label_cpp_default(this, bin) result(label)
+    class(CppTallyFilter), intent(in) :: this
+    integer,               intent(in) :: bin
+    character(MAX_LINE_LEN)           :: label
+    label = this % text_label_cpp_inner(bin)
+  end function text_label_cpp_default
+
+  subroutine initialize_cpp_default(this)
+    class(CppTallyFilter), intent(inout) :: this
+    call this % initialize_cpp_inner()
+  end subroutine initialize_cpp_default
+
+!===============================================================================
+! FILTER_FROM_F given a Fortran index, return a pointer to a C++ filter.
+!===============================================================================
+
+  function filter_from_f(index) result(filt) bind(C)
+    integer(C_INT32_T), intent(in), value :: index
+    type(C_PTR) :: filt
+
+    filt = C_NULL_PTR
+    select type(f => filters(index) % obj)
+    class is (CppTallyFilter)
+      filt = f % ptr
+    end select
+  end function
+
+!===============================================================================
+! FILTER_UPDATE_N_BINS given a Fortran index, updates filt % n_bins using C++.
+!===============================================================================
+
+  subroutine filter_update_n_bins(index) bind(C)
+    integer(C_INT32_T), intent(in), value :: index
+    select type(f => filters(index) % obj)
+    class is (CppTallyFilter)
+      f % n_bins = f % n_bins_cpp()
+    end select
+  end subroutine
+end module tally_filter_cpp

src/tallies/tally_filter_delayedgroup.F90

@@ -62,8 +62,8 @@ contains
     integer,                   intent(in)  :: estimator
     type(TallyFilterMatch),         intent(inout) :: match
 
-    call match % bins % push_back(1)
-    call match % weights % push_back(ONE)
+    call match % bins_push_back(1)
+    call match % weights_push_back(ONE)
   end subroutine get_all_bins_dg
 
   subroutine to_statepoint_dg(this, filter_group)

src/tallies/tally_filter_distribcell.F90

@@ -1,169 +1,21 @@
 module tally_filter_distribcell
 
-  use, intrinsic :: ISO_C_BINDING
-
-  use constants
-  use dict_header,     only: EMPTY
-  use error
-  use geometry_header
-  use hdf5_interface
-  use particle_header, only: Particle
-  use string,          only: to_str
-  use tally_filter_header
-  use xml_interface
+  use tally_filter_cpp
 
   implicit none
   private
 
-!===============================================================================
-! DISTRIBCELLFILTER specifies which distributed geometric cells tally events
-! reside in.
-!===============================================================================
-
-  type, public, extends(TallyFilter) :: DistribcellFilter
-    integer :: cell
+  type, public, extends(CppTallyFilter) :: DistribcellFilter
   contains
-    procedure :: from_xml
-    procedure :: get_all_bins => get_all_bins_distribcell
-    procedure :: to_statepoint => to_statepoint_distribcell
-    procedure :: text_label => text_label_distribcell
     procedure :: initialize => initialize_distribcell
   end type DistribcellFilter
 
 contains
 
-  subroutine from_xml(this, node)
-    class(DistribcellFilter), intent(inout) :: this
-    type(XMLNode), intent(in) :: node
-
-    integer :: n
-
-    n = node_word_count(node, "bins")
-    if (n /= 1) call fatal_error("Only one cell can be &
-         &specified per distribcell filter.")
-
-    ! Store bins
-    call get_node_value(node, "bins", this % cell)
-  end subroutine from_xml
-
-  subroutine get_all_bins_distribcell(this, p, estimator, match)
-    class(DistribcellFilter), intent(in)  :: this
-    type(Particle),           intent(in)  :: p
-    integer,                  intent(in)  :: estimator
-    type(TallyFilterMatch),   intent(inout) :: match
-
-    integer :: distribcell_index, offset, i
-
-      distribcell_index = cells(this % cell) % distribcell_index()
-      offset = 0
-      do i = 1, p % n_coord
-        if (cells(p % coord(i) % cell + 1) % type() == FILL_UNIVERSE) then
-          offset = offset + cells(p % coord(i) % cell + 1) &
-               % offset(distribcell_index)
-        elseif (cells(p % coord(i) % cell + 1) % type() == FILL_LATTICE) then
-          if (lattices(p % coord(i + 1) % lattice) % are_valid_indices([&
-               p % coord(i + 1) % lattice_x, &
-               p % coord(i + 1) % lattice_y, &
-               p % coord(i + 1) % lattice_z])) then
-            offset = offset + lattices(p % coord(i + 1) % lattice) &
-                 % offset(distribcell_index, &
-                 [p % coord(i + 1) % lattice_x, &
-                 p % coord(i + 1) % lattice_y, &
-                 p % coord(i + 1) % lattice_z])
-          end if
-        end if
-        if (this % cell == p % coord(i) % cell + 1) then
-          call match % bins % push_back(offset + 1)
-          call match % weights % push_back(ONE)
-          return
-        end if
-      end do
-  end subroutine get_all_bins_distribcell
-
-  subroutine to_statepoint_distribcell(this, filter_group)
-    class(DistribcellFilter), intent(in) :: this
-    integer(HID_T),           intent(in) :: filter_group
-
-    call write_dataset(filter_group, "type", "distribcell")
-    call write_dataset(filter_group, "n_bins", this % n_bins)
-    call write_dataset(filter_group, "bins", cells(this % cell) % id())
-  end subroutine to_statepoint_distribcell
-
   subroutine initialize_distribcell(this)
     class(DistribcellFilter), intent(inout) :: this
-
-    integer :: id
-    integer :: val
-
-    ! Convert id to index.
-    id = this % cell
-    val = cell_dict % get(id)
-    if (val /= EMPTY) then
-      this % cell = val
-      this % n_bins = cells(this % cell) % n_instances()
-    else
-      call fatal_error("Could not find cell " // trim(to_str(id)) &
-           &// " specified on tally filter.")
-    end if
+    call this % initialize_cpp_inner()
+    this % n_bins = this % n_bins_cpp()
   end subroutine initialize_distribcell
 
-  function text_label_distribcell(this, bin) result(label)
-    class(DistribcellFilter), intent(in) :: this
-    integer,                  intent(in) :: bin
-    character(MAX_LINE_LEN)              :: label
-
-    label = ''
-    call find_offset(this % cell, bin-1, label)
-    label = "Distributed Cell " // label
-  end function text_label_distribcell
-
-!===============================================================================
-! FIND_OFFSET (for distribcell) uses a given map number, a target cell ID, and
-! a target offset to build a string which is the path from the base universe to
-! the target cell with the given offset
-!===============================================================================
-
-  subroutine find_offset(i_cell, target_offset, path)
-    integer, intent(in) :: i_cell         ! The target cell index
-    integer, intent(in) :: target_offset  ! Target offset
-    character(*), intent(inout) :: path   ! Path to offset
-
-    integer :: map                  ! Index in maps vector
-    integer :: i                    ! Index over cells
-
-    integer(C_INT) :: path_len
-    character(kind=C_CHAR), allocatable, target :: path_c(:)
-
-    interface
-      function distribcell_path_len(target_cell, map, target_offset, root_univ)&
-           bind(C) result(len)
-        import C_INT32_T, C_INT
-        integer(C_INT32_T), intent(in), value :: target_cell, map, &
-                                                 target_offset, root_univ
-        integer(C_INT) :: len
-      end function distribcell_path_len
-
-      subroutine distribcell_path(target_cell, map, target_offset, root_univ, &
-                                  path) bind(C)
-        import C_INT32_T, C_CHAR
-        integer(C_INT32_T), intent(in), value :: target_cell, map, &
-                                                 target_offset, root_univ
-        character(kind=C_CHAR), intent(out)   :: path(*)
-      end subroutine distribcell_path
-    end interface
-
-    ! Get the distribcell index for this cell
-    map = cells(i_cell) % distribcell_index()
-
-    path_len = distribcell_path_len(i_cell-1, map, target_offset, &
-                                    root_universe)
-    allocate(path_c(path_len))
-    call distribcell_path(i_cell-1, map, target_offset, root_universe, &
-                          path_c)
-    do i = 1, min(path_len, MAX_LINE_LEN)
-      if (path_c(i) == C_NULL_CHAR) exit
-      path(i:i) = path_c(i)
-    end do
-  end subroutine find_offset
-
 end module tally_filter_distribcell

src/tallies/tally_filter_energy.F90

@@ -95,11 +95,11 @@ contains
 
     if (p % g /= NONE .and. this % matches_transport_groups) then
       if (estimator == ESTIMATOR_TRACKLENGTH) then
-        call match % bins % push_back(num_energy_groups - p % g + 1)
-        call match % weights % push_back(ONE)
+        call match % bins_push_back(num_energy_groups - p % g + 1)
+        call match % weights_push_back(ONE)
       else
-        call match % bins % push_back(num_energy_groups - p % last_g + 1)
-        call match % weights % push_back(ONE)
+        call match % bins_push_back(num_energy_groups - p % last_g + 1)
+        call match % weights_push_back(ONE)
       end if
 
     else
@@ -109,8 +109,8 @@ contains
       ! Search to find incoming energy bin.
       bin = binary_search(this % bins, n + 1, E)
       if (bin /= NO_BIN_FOUND) then
-        call match % bins % push_back(bin)
-        call match % weights % push_back(ONE)
+        call match % bins_push_back(bin)
+        call match % weights_push_back(ONE)
       end if
     end if
   end subroutine get_all_bins_energy
@@ -155,16 +155,16 @@ contains
     if (p % g /= NONE .and. this % matches_transport_groups) then
       ! Tallies are ordered in increasing groups, group indices
       ! however are the opposite, so switch
-      call match % bins % push_back(num_energy_groups - p % g + 1)
-      call match % weights % push_back(ONE)
+      call match % bins_push_back(num_energy_groups - p % g + 1)
+      call match % weights_push_back(ONE)
 
     else
 
       ! Search to find incoming energy bin.
       bin = binary_search(this % bins, n + 1, p % E)
       if (bin /= NO_BIN_FOUND) then
-        call match % bins % push_back(bin)
-        call match % weights % push_back(ONE)
+        call match % bins_push_back(bin)
+        call match % weights_push_back(ONE)
       end if
     end if
   end subroutine get_all_bins_energyout

src/tallies/tally_filter_energyfunc.F90

@@ -1,124 +0,0 @@
-module tally_filter_energyfunc
-
-  use, intrinsic :: ISO_C_BINDING
-
-  use algorithm,          only: binary_search
-  use constants
-  use error,              only: fatal_error
-  use hdf5_interface
-  use particle_header,    only: Particle
-  use settings,           only: run_CE
-  use string,             only: to_str
-  use tally_filter_header
-  use xml_interface
-
-  implicit none
-  private
-
-!===============================================================================
-! EnergyFunctionFilter multiplies tally scores by an arbitrary function of
-! incident energy described by a piecewise linear-linear interpolation.
-!===============================================================================
-
-  type, public, extends(TallyFilter) :: EnergyFunctionFilter
-    real(8), allocatable :: energy(:)
-    real(8), allocatable :: y(:)
-
-  contains
-    procedure :: from_xml
-    procedure :: get_all_bins => get_all_bins_energyfunction
-    procedure :: to_statepoint => to_statepoint_energyfunction
-    procedure :: text_label => text_label_energyfunction
-  end type EnergyFunctionFilter
-
-contains
-
-  subroutine from_xml(this, node)
-    class(EnergyFunctionFilter), intent(inout) :: this
-    type(XMLNode), intent(in) :: node
-
-    integer :: n
-
-    this % n_bins = 1
-    ! Make sure this is continuous-energy mode.
-    if (.not. run_CE) then
-      call fatal_error("EnergyFunction filters are only supported for &
-           &continuous-energy transport calculations")
-    end if
-
-    ! Allocate and store energy grid.
-    if (.not. check_for_node(node, "energy")) then
-      call fatal_error("Energy grid not specified for EnergyFunction &
-           &filter.")
-    end if
-    n = node_word_count(node, "energy")
-    allocate(this % energy(n))
-    call get_node_array(node, "energy", this % energy)
-
-    ! Allocate and store interpolant values.
-    if (.not. check_for_node(node, "y")) then
-      call fatal_error("y values not specified for EnergyFunction &
-           &filter.")
-    end if
-    n = node_word_count(node, "y")
-    allocate(this % y(n))
-    call get_node_array(node, "y", this % y)
-  end subroutine from_xml
-
-  subroutine get_all_bins_energyfunction(this, p, estimator, match)
-    class(EnergyFunctionFilter), intent(in)  :: this
-    type(Particle),              intent(in)  :: p
-    integer,                     intent(in)  :: estimator
-    type(TallyFilterMatch),      intent(inout) :: match
-
-    integer :: n, indx
-    real(8) :: E, f, weight
-
-    select type(this)
-    type is (EnergyFunctionFilter)
-      n = size(this % energy)
-
-      ! Get pre-collision energy of particle
-      E = p % last_E
-
-      ! Search to find incoming energy bin.
-      indx = binary_search(this % energy, n, E)
-
-      ! Compute an interpolation factor between nearest bins.
-      f = (E - this % energy(indx)) &
-           / (this % energy(indx+1) - this % energy(indx))
-
-      ! Interpolate on the lin-lin grid.
-      call match % bins % push_back(1)
-      weight = (ONE - f) * this % y(indx) + f * this % y(indx+1)
-      call match % weights % push_back(weight)
-    end select
-  end subroutine get_all_bins_energyfunction
-
-  subroutine to_statepoint_energyfunction(this, filter_group)
-    class(EnergyFunctionFilter), intent(in) :: this
-    integer(HID_T),              intent(in) :: filter_group
-
-    select type(this)
-    type is (EnergyFunctionFilter)
-      call write_dataset(filter_group, "type", "energyfunction")
-      call write_dataset(filter_group, "energy", this % energy)
-      call write_dataset(filter_group, "y", this % y)
-    end select
-  end subroutine to_statepoint_energyfunction
-
-  function text_label_energyfunction(this, bin) result(label)
-    class(EnergyFunctionFilter), intent(in) :: this
-    integer,                     intent(in) :: bin
-    character(MAX_LINE_LEN)                 :: label
-
-    select type(this)
-    type is (EnergyFunctionFilter)
-      write(label, FMT="(A, ES8.1, A, ES8.1, A, ES8.1, A, ES8.1, A)") &
-           "Energy Function f([", this % energy(1), ", ..., ", &
-           this % energy(size(this % energy)), "]) = [", this % y(1), &
-           ", ..., ", this % y(size(this % y)), "]"
-    end select
-  end function text_label_energyfunction
-
-end module tally_filter_energyfunc

src/tallies/tally_filter_header.F90

@@ -7,7 +7,6 @@ module tally_filter_header
   use error
   use hdf5_interface,  only: HID_T
   use particle_header, only: Particle
-  use stl_vector,      only: VectorInt, VectorReal
   use string,          only: to_str
   use xml_interface,   only: XMLNode
 
@@ -20,19 +19,38 @@ module tally_filter_header
   public :: openmc_filter_set_id
   public :: openmc_get_filter_index
   public :: openmc_get_filter_next_id
+  public :: filter_match_pointer
+
+  interface
+    function filter_match_pointer(indx) bind(C) result(ptr)
+      import C_PTR, C_INT
+      integer(C_INT), intent(in), value :: indx
+      type(C_PTR)                       :: ptr
+    end function filter_match_pointer
+  end interface
 
 !===============================================================================
 ! TALLYFILTERMATCH stores every valid bin and weight for a filter
 !===============================================================================
 
   type, public :: TallyFilterMatch
+    type(C_PTR) :: ptr
+
     ! Index of the bin and weight being used in the current filter combination
     integer          :: i_bin
-    type(VectorInt),  pointer :: bins
-    type(VectorReal), pointer :: weights
 
     ! Indicates whether all valid bins for this filter have been found
     logical          :: bins_present = .false.
+
+  contains
+    procedure :: bins_push_back
+    procedure :: weights_push_back
+    procedure :: bins_clear
+    procedure :: weights_clear
+    procedure :: bins_size
+    procedure :: bins_data
+    procedure :: weights_data
+    procedure :: bins_set_data
   end type TallyFilterMatch
 
 !===============================================================================
@@ -129,6 +147,116 @@ module tally_filter_header
 
 contains
 
+!===============================================================================
+! TallyFilterMatch implementation
+!===============================================================================
+
+  subroutine bins_push_back(this, val)
+    class(TallyFilterMatch), intent(inout) :: this
+    integer,                 intent(in)    :: val
+    interface
+      subroutine filter_match_bins_push_back(ptr, val) bind(C)
+        import C_PTR, C_INT
+        type(C_PTR),                value :: ptr
+        integer(C_INT), intent(in), value :: val
+      end subroutine
+    end interface
+    call filter_match_bins_push_back(this % ptr, val)
+  end subroutine bins_push_back
+
+  subroutine weights_push_back(this, val)
+    class(TallyFilterMatch), intent(inout) :: this
+    real(8),                 intent(in)    :: val
+    interface
+      subroutine filter_match_weights_push_back(ptr, val) bind(C)
+        import C_PTR, C_DOUBLE
+        type(C_PTR),                value :: ptr
+        real(C_DOUBLE), intent(in), value :: val
+      end subroutine
+    end interface
+    call filter_match_weights_push_back(this % ptr, val)
+  end subroutine weights_push_back
+
+  subroutine bins_clear(this)
+    class(TallyFilterMatch), intent(inout) :: this
+    interface
+      subroutine filter_match_bins_clear(ptr) bind(C)
+        import C_PTR
+        type(C_PTR), value :: ptr
+      end subroutine
+    end interface
+    call filter_match_bins_clear(this % ptr)
+  end subroutine bins_clear
+
+  subroutine weights_clear(this)
+    class(TallyFilterMatch), intent(inout) :: this
+    interface
+      subroutine filter_match_weights_clear(ptr) bind(C)
+        import C_PTR
+        type(C_PTR), value :: ptr
+      end subroutine
+    end interface
+    call filter_match_weights_clear(this % ptr)
+  end subroutine weights_clear
+
+  function bins_size(this) result(len)
+    class(TallyFilterMatch), intent(inout) :: this
+    integer                                :: len
+    interface
+      function filter_match_bins_size(ptr) bind(C) result(len)
+        import C_PTR, C_INT
+        type(C_PTR), value :: ptr
+        integer(C_INT)     :: len
+      end function
+    end interface
+    len = filter_match_bins_size(this % ptr)
+  end function bins_size
+
+  function bins_data(this, indx) result(val)
+    class(TallyFilterMatch), intent(inout) :: this
+    integer,                 intent(in)    :: indx
+    integer                                :: val
+    interface
+      function filter_match_bins_data(ptr, indx) bind(C) result(val)
+        import C_PTR, C_INT
+        type(C_PTR),                value :: ptr
+        integer(C_INT), intent(in), value :: indx
+        integer(C_INT)                    :: val
+      end function
+    end interface
+    val = filter_match_bins_data(this % ptr, indx)
+  end function bins_data
+
+  function weights_data(this, indx) result(val)
+    class(TallyFilterMatch), intent(inout) :: this
+    integer,                 intent(in)    :: indx
+    real(8)                                :: val
+    interface
+      function filter_match_weights_data(ptr, indx) bind(C) result(val)
+        import C_PTR, C_INT, C_DOUBLE
+        type(C_PTR),                value :: ptr
+        integer(C_INT), intent(in), value :: indx
+        real(C_DOUBLE)                    :: val
+      end function
+    end interface
+    val = filter_match_weights_data(this % ptr, indx)
+  end function weights_data
+
+  subroutine bins_set_data(this, indx, val)
+    class(TallyFilterMatch), intent(inout) :: this
+    integer,                 intent(in)    :: indx
+    integer,                 intent(in)    :: val
+    interface
+      subroutine filter_match_bins_set_data(ptr, indx, val) bind(C)
+        import C_PTR, C_INT
+        type(C_PTR),    value             :: ptr
+        integer(C_INT), value, intent(in) :: indx
+        integer(C_INT), value, intent(in) :: val
+      end subroutine
+    end interface
+    call filter_match_bins_set_data(this % ptr, indx, val)
+  end subroutine bins_set_data
+
 !===============================================================================
 ! INITIALIZE sets up any internal data, as necessary.  If this procedure is not
 ! overriden by the derived class, then it will do nothing by default.
@@ -153,8 +281,8 @@ contains
 ! array is sufficient and a filter object has already been allocated.
 !===============================================================================
 
-  function verify_filter(index) result(err)
-    integer(C_INT32_T), intent(in) :: index
+  function verify_filter(index) result(err) bind(C)
+    integer(C_INT32_T), intent(in), value :: index
     integer(C_INT) :: err
 
     err = 0

src/tallies/tally_filter_legendre.F90

@@ -2,121 +2,20 @@ module tally_filter_legendre
 
   use, intrinsic :: ISO_C_BINDING
 
-  use constants
-  use error
-  use hdf5_interface
-  use math,                only: calc_pn
-  use particle_header,     only: Particle
-  use string,              only: to_str
-  use tally_filter_header
-  use xml_interface
+  use tally_filter_cpp
 
   implicit none
-  private
-  public :: openmc_legendre_filter_get_order
-  public :: openmc_legendre_filter_set_order
 
-!===============================================================================
-! LEGENDREFILTER gives Legendre moments of the change in scattering angle
-!===============================================================================
+  interface
+    function openmc_legendre_filter_set_order(index, order) result(err) bind(C)
+      import C_INT32_T, C_INT
+      integer(C_INT32_T), value :: index
+      integer(C_INT),     value :: order
+      integer(C_INT) :: err
+    end function
+  end interface
 
-  type, public, extends(TallyFilter) :: LegendreFilter
-    integer(C_INT) :: order
-  contains
-    procedure :: from_xml
-    procedure :: get_all_bins
-    procedure :: to_statepoint
-    procedure :: text_label
+  type, extends(CppTallyFilter) :: LegendreFilter
   end type LegendreFilter
 
-contains
-
-!===============================================================================
-! LegendreFilter methods
-!===============================================================================
-
-  subroutine from_xml(this, node)
-    class(LegendreFilter), intent(inout) :: this
-    type(XMLNode), intent(in) :: node
-
-    ! Get specified order
-    call get_node_value(node, "order", this % order)
-    this % n_bins = this % order + 1
-  end subroutine from_xml
-
-  subroutine get_all_bins(this, p, estimator, match)
-    class(LegendreFilter), intent(in)  :: this
-    type(Particle),      intent(in)  :: p
-    integer,             intent(in)  :: estimator
-    type(TallyFilterMatch),   intent(inout) :: match
-
-    integer :: i
-    real(C_DOUBLE) :: wgt(this % n_bins)
-
-    call calc_pn(this % order, p % mu, wgt)
-    do i = 1, this % n_bins
-      call match % bins % push_back(i)
-      call match % weights % push_back(wgt(i))
-    end do
-  end subroutine get_all_bins
-
-  subroutine to_statepoint(this, filter_group)
-    class(LegendreFilter), intent(in) :: this
-    integer(HID_T),      intent(in) :: filter_group
-
-    call write_dataset(filter_group, "type", "legendre")
-    call write_dataset(filter_group, "n_bins", this % n_bins)
-    call write_dataset(filter_group, "order", this % order)
-  end subroutine to_statepoint
-
-  function text_label(this, bin) result(label)
-    class(LegendreFilter), intent(in) :: this
-    integer,             intent(in) :: bin
-    character(MAX_LINE_LEN)         :: label
-
-    label = "Legendre expansion, P" // trim(to_str(bin - 1))
-  end function text_label
-
-!===============================================================================
-!                               C API FUNCTIONS
-!===============================================================================
-
-  function openmc_legendre_filter_get_order(index, order) result(err) bind(C)
-    ! Get the order of an expansion filter
-    integer(C_INT32_T), value       :: index
-    integer(C_INT),     intent(out) :: order
-    integer(C_INT) :: err
-
-    err = verify_filter(index)
-    if (err == 0) then
-      select type (f => filters(index) % obj)
-      type is (LegendreFilter)
-        order = f % order
-      class default
-        err = E_INVALID_TYPE
-        call set_errmsg("Tried to get order on a non-expansion filter.")
-      end select
-    end if
-  end function openmc_legendre_filter_get_order
-
-
-  function openmc_legendre_filter_set_order(index, order) result(err) bind(C)
-    ! Set the order of an expansion filter
-    integer(C_INT32_T), value :: index
-    integer(C_INT),     value :: order
-    integer(C_INT) :: err
-
-    err = verify_filter(index)
-    if (err == 0) then
-      select type (f => filters(index) % obj)
-      type is (LegendreFilter)
-        f % order = order
-        f % n_bins = order + 1
-      class default
-        err = E_INVALID_TYPE
-        call set_errmsg("Tried to set order on a non-expansion filter.")
-      end select
-    end if
-  end function openmc_legendre_filter_set_order
-
 end module tally_filter_legendre

src/tallies/tally_filter_material.F90

@@ -1,175 +0,0 @@
-module tally_filter_material
-
-  use, intrinsic :: ISO_C_BINDING
-
-  use constants
-  use dict_header,     only: DictIntInt, EMPTY
-  use error
-  use hdf5_interface
-  use material_header, only: materials, material_dict
-  use particle_header, only: Particle
-  use string,          only: to_str
-  use tally_filter_header
-  use xml_interface
-
-  implicit none
-  private
-  public :: openmc_material_filter_get_bins
-  public :: openmc_material_filter_set_bins
-
-!===============================================================================
-! MATERIAL specifies which material tally events reside in.
-!===============================================================================
-
-  type, public, extends(TallyFilter) :: MaterialFilter
-    integer, allocatable :: materials(:)
-    type(DictIntInt)     :: map
-  contains
-    procedure :: from_xml
-    procedure :: get_all_bins => get_all_bins_material
-    procedure :: to_statepoint => to_statepoint_material
-    procedure :: text_label => text_label_material
-    procedure :: initialize => initialize_material
-  end type MaterialFilter
-
-contains
-
-  subroutine from_xml(this, node)
-    class(MaterialFilter), intent(inout) :: this
-    type(XMLNode), intent(in) :: node
-
-    integer :: n
-
-    n = node_word_count(node, "bins")
-
-    ! Allocate and store bins
-    this % n_bins = n
-    allocate(this % materials(n))
-    call get_node_array(node, "bins", this % materials)
-  end subroutine from_xml
-
-  subroutine get_all_bins_material(this, p, estimator, match)
-    class(MaterialFilter), intent(in)  :: this
-    type(Particle),        intent(in)  :: p
-    integer,               intent(in)  :: estimator
-    type(TallyFilterMatch),     intent(inout) :: match
-
-    integer :: val
-
-    val = this % map % get(p % material)
-    if (val /= EMPTY) then
-      call match % bins % push_back(val)
-      call match % weights % push_back(ONE)
-    end if
-
-  end subroutine get_all_bins_material
-
-  subroutine to_statepoint_material(this, filter_group)
-    class(MaterialFilter), intent(in) :: this
-    integer(HID_T),        intent(in) :: filter_group
-
-    integer :: i
-    integer, allocatable :: material_ids(:)
-
-    call write_dataset(filter_group, "type", "material")
-    call write_dataset(filter_group, "n_bins", this % n_bins)
-
-    allocate(material_ids(size(this % materials)))
-    do i = 1, size(this % materials)
-      material_ids(i) = materials(this % materials(i)) % id()
-    end do
-    call write_dataset(filter_group, "bins", material_ids)
-  end subroutine to_statepoint_material
-
-  subroutine initialize_material(this)
-    class(MaterialFilter), intent(inout) :: this
-
-    integer :: i, id
-    integer :: val
-
-    ! Convert ids to indices.
-    do i = 1, this % n_bins
-      id = this % materials(i)
-      val = material_dict % get(id)
-      if (val /= EMPTY) then
-        this % materials(i) = val
-      else
-        call fatal_error("Could not find material " // trim(to_str(id)) &
-             &// " specified on a tally filter.")
-      end if
-    end do
-
-    ! Generate mapping from material indices to filter bins.
-    do i = 1, this % n_bins
-      call this % map % set(this % materials(i), i)
-    end do
-  end subroutine initialize_material
-
-  function text_label_material(this, bin) result(label)
-    class(MaterialFilter), intent(in) :: this
-    integer,               intent(in) :: bin
-    character(MAX_LINE_LEN)           :: label
-
-    label = "Material " // to_str(materials(this % materials(bin)) % id())
-  end function text_label_material
-
-!===============================================================================
-!                               C API FUNCTIONS
-!===============================================================================
-
-  function openmc_material_filter_get_bins(index, bins, n) result(err) bind(C)
-    ! Return the bins for a material filter
-    integer(C_INT32_T), value :: index
-    type(C_PTR), intent(out) :: bins
-    integer(C_INT32_T), intent(out) :: n
-    integer(C_INT) :: err
-
-    err = verify_filter(index)
-    if (err == 0) then
-      select type (f => filters(index) % obj)
-      type is (MaterialFilter)
-        bins = C_LOC(f % materials)
-        n = size(f % materials)
-        err = 0
-        class default
-        err = E_INVALID_TYPE
-        call set_errmsg("Tried to get material filter bins on a &
-             &non-material filter.")
-      end select
-    end if
-  end function openmc_material_filter_get_bins
-
-
-  function openmc_material_filter_set_bins(index, n, bins) result(err) bind(C)
-    ! Set the materials for the filter
-    integer(C_INT32_T), value, intent(in) :: index
-    integer(C_INT32_T), value, intent(in) :: n
-    integer(C_INT32_T), intent(in) :: bins(n)
-    integer(C_INT) :: err
-
-    integer :: i
-
-    err = verify_filter(index)
-    if (err == 0) then
-      select type (f => filters(index) % obj)
-      type is (MaterialFilter)
-        f % n_bins = n
-        if (allocated(f % materials)) deallocate(f % materials)
-        allocate(f % materials(n))
-        f % materials(:) = bins
-
-        ! Generate mapping from material indices to filter bins.
-        call f % map % clear()
-        do i = 1, n
-          call f % map % set(f % materials(i), i)
-        end do
-
-        class default
-        err = E_INVALID_TYPE
-        call set_errmsg("Tried to set material filter bins on a &
-             &non-material filter.")
-      end select
-    end if
-  end function openmc_material_filter_set_bins
-
-end module tally_filter_material

src/tallies/tally_filter_mesh.F90

@@ -2,197 +2,37 @@ module tally_filter_mesh
 
   use, intrinsic :: ISO_C_BINDING
 
-  use constants
-  use dict_header,         only: EMPTY
-  use error
-  use mesh_header
-  use hdf5_interface
-  use particle_header,     only: Particle
-  use string,              only: to_str
-  use tally_filter_header
-  use xml_interface
+  use tally_filter_cpp
 
   implicit none
-  private
-  public :: openmc_mesh_filter_get_mesh
-  public :: openmc_mesh_filter_set_mesh
 
-!===============================================================================
-! MESHFILTER indexes the location of particle events to a regular mesh.  For
-! tracklength tallies, it will produce multiple valid bins and the bin weight
-! will correspond to the fraction of the track length that lies in that bin.
-!===============================================================================
+  interface
+    function openmc_mesh_filter_set_mesh(index, index_mesh) 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) :: index_mesh
+      integer(C_INT) :: err
+    end function
+  end interface
 
-  type, public, extends(TallyFilter) :: MeshFilter
-    integer :: mesh
+  type, public, extends(CppTallyFilter) :: MeshFilter
   contains
-    procedure :: from_xml
-    procedure :: get_all_bins => get_all_bins_mesh
-    procedure :: to_statepoint => to_statepoint_mesh
-    procedure :: text_label => text_label_mesh
+    procedure :: mesh => get_mesh
   end type MeshFilter
 
 contains
 
-  subroutine from_xml(this, node)
-    class(MeshFilter), intent(inout) :: this
-    type(XMLNode), intent(in) :: node
-
-    integer :: i
-    integer :: id
-    integer :: n
-    integer(C_INT) :: err
-    type(RegularMesh) :: m
-
-    n = node_word_count(node, "bins")
-
-    if (n /= 1) call fatal_error("Only one mesh can be &
-         &specified per mesh filter.")
-
-    ! Determine id of mesh
-    call get_node_value(node, "bins", id)
-
-    ! Get pointer to mesh
-    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
-    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)
-    class(MeshFilter), intent(in)  :: this
-    type(Particle),    intent(in)  :: p
-    integer,           intent(in)  :: estimator
-    type(TallyFilterMatch), intent(inout) :: match
-
-    integer :: bin
-    type(RegularMesh) :: m
-    type(C_PTR) :: ptr_bins, ptr_weights
-
+  function get_mesh(this) result(mesh)
+    class(MeshFilter) :: this
+    integer :: mesh
     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
+      function mesh_filter_get_mesh(filt) result(index_mesh) bind(C)
+        import C_PTR, C_INT
+        type(C_PTR), value :: filt
+        integer(C_INT)     :: index_mesh
+      end function mesh_filter_get_mesh
     end interface
-
-    ! Get a pointer to the mesh.
-    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 >= 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
-
-  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", m % id())
-  end subroutine to_statepoint_mesh
-
-  function text_label_mesh(this, bin) result(label)
-    class(MeshFilter), intent(in) :: this
-    integer,           intent(in) :: bin
-    character(MAX_LINE_LEN)       :: label
-
-    integer, allocatable       :: ijk(:)
-    type(RegularMesh) :: m
-
-    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
-
-!===============================================================================
-!                               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
-
-    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 >= 0 .and. index_mesh < n_meshes()) then
-          f % mesh = index_mesh
-          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
-        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
+    mesh = mesh_filter_get_mesh(this % ptr)
+  end function get_mesh
 
 end module tally_filter_mesh

src/tallies/tally_filter_meshsurface.F90

@@ -2,224 +2,21 @@ module tally_filter_meshsurface
 
   use, intrinsic :: ISO_C_BINDING
 
-  use constants
-  use error
-  use mesh_header
-  use hdf5_interface
-  use particle_header,     only: Particle
-  use string,              only: to_str
-  use tally_filter_header
-  use xml_interface
+  use tally_filter_cpp
 
   implicit none
-  private
-  public :: openmc_meshsurface_filter_get_mesh
-  public :: openmc_meshsurface_filter_set_mesh
 
-!===============================================================================
-! MESHFILTER indexes the location of particle events to a regular mesh.  For
-! tracklength tallies, it will produce multiple valid bins and the bin weight
-! will correspond to the fraction of the track length that lies in that bin.
-!===============================================================================
+  interface
+    function openmc_meshsurface_filter_set_mesh(index, index_mesh) 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) :: index_mesh
+      integer(C_INT) :: err
+    end function
+  end interface
 
-  type, public, extends(TallyFilter) :: MeshSurfaceFilter
-    integer :: mesh
-  contains
-    procedure :: from_xml
-    procedure :: get_all_bins
-    procedure :: to_statepoint
-    procedure :: text_label
+  type, extends(CppTallyFilter) :: MeshSurfaceFilter
   end type MeshSurfaceFilter
 
-contains
-
-  subroutine from_xml(this, node)
-    class(MeshSurfaceFilter), intent(inout) :: this
-    type(XMLNode), intent(in) :: node
-
-    integer :: i
-    integer :: id
-    integer :: n
-    integer(C_INT) :: err
-    type(RegularMesh) :: m
-
-    n = node_word_count(node, "bins")
-
-    if (n /= 1) call fatal_error("Only one mesh can be &
-         &specified per meshsurface filter.")
-
-    ! Determine id of mesh
-    call get_node_value(node, "bins", id)
-
-    ! Get pointer to mesh
-    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
-    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)
-    class(MeshSurfaceFilter), intent(in)  :: this
-    type(Particle),           intent(in)  :: p
-    integer,                  intent(in)  :: estimator
-    type(TallyFilterMatch), intent(inout) :: match
-
-    type(RegularMesh) :: m
-    type(C_PTR) :: ptr_bins, ptr_weights
-
-    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
-
-    ! Get a pointer to the mesh.
-    m = meshes(this % mesh)
-
-    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
-
-  subroutine to_statepoint(this, filter_group)
-    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", m % id())
-  end subroutine to_statepoint
-
-  function text_label(this, bin) result(label)
-    class(MeshSurfaceFilter), intent(in) :: this
-    integer,                  intent(in) :: bin
-    character(MAX_LINE_LEN)              :: label
-
-    integer :: i_mesh
-    integer :: i_surf
-    integer :: n_dim
-    integer, allocatable :: ijk(:)
-    type(RegularMesh) :: m
-
-    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 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 function text_label
-
-!===============================================================================
-!                               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
-    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 >= 0 .and. index_mesh < n_meshes()) then
-          f % mesh = index_mesh
-          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.")
-        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
-
 end module tally_filter_meshsurface

src/tallies/tally_filter_mu.F90

@@ -1,110 +0,0 @@
-module tally_filter_mu
-
-  use, intrinsic :: ISO_C_BINDING
-
-  use algorithm,          only: binary_search
-  use constants,          only: ONE, TWO, MAX_LINE_LEN, NO_BIN_FOUND
-  use error,              only: fatal_error
-  use hdf5_interface
-  use particle_header,    only: Particle
-  use string,             only: to_str
-  use tally_filter_header
-  use xml_interface
-
-  implicit none
-  private
-
-!===============================================================================
-! MUFILTER bins the incoming-outgoing direction cosine.  This is only used for
-! scatter reactions.
-!===============================================================================
-
-  type, public, extends(TallyFilter) :: MuFilter
-    real(8), allocatable :: bins(:)
-  contains
-    procedure :: from_xml
-    procedure :: get_all_bins => get_all_bins_mu
-    procedure :: to_statepoint => to_statepoint_mu
-    procedure :: text_label => text_label_mu
-  end type MuFilter
-
-contains
-
-  subroutine from_xml(this, node)
-    class(MuFilter), intent(inout) :: this
-    type(XMLNode), intent(in) :: node
-
-    integer :: i
-    integer :: n_angle
-    integer :: n
-    real(8) :: d_angle
-
-    n = node_word_count(node, "bins")
-
-    ! Allocate and store bins
-    this % n_bins = n - 1
-    allocate(this % bins(n))
-    call get_node_array(node, "bins", this % bins)
-
-    ! Allow a user to input a lone number which will mean that you
-    ! subdivide [-1,1] evenly with the input being the number of bins
-    if (n == 1) then
-      n_angle = int(this % bins(1))
-      if (n_angle > 1) then
-        this % n_bins = n_angle
-        d_angle = TWO / n_angle
-        deallocate(this % bins)
-        allocate(this % bins(n_angle + 1))
-        do i = 1, n_angle
-          this % bins(i) = -ONE + (i - 1) * d_angle
-        end do
-        this % bins(n_angle + 1) = ONE
-      else
-        call fatal_error("Number of bins for mu filter must be&
-             & greater than 1.")
-      end if
-    end if
-  end subroutine from_xml
-
-  subroutine get_all_bins_mu(this, p, estimator, match)
-    class(MuFilter), intent(in)  :: this
-    type(Particle),  intent(in)  :: p
-    integer,         intent(in)  :: estimator
-    type(TallyFilterMatch), intent(inout) :: match
-
-    integer :: n
-    integer :: bin
-
-    n = this % n_bins
-
-    ! Search to find incoming energy bin.
-    bin = binary_search(this % bins, n + 1, p % mu)
-    if (bin /= NO_BIN_FOUND) then
-      call match % bins % push_back(bin)
-      call match % weights % push_back(ONE)
-    end if
-  end subroutine get_all_bins_mu
-
-  subroutine to_statepoint_mu(this, filter_group)
-    class(MuFilter), intent(in) :: this
-    integer(HID_T),  intent(in) :: filter_group
-
-    call write_dataset(filter_group, "type", "mu")
-    call write_dataset(filter_group, "n_bins", this % n_bins)
-    call write_dataset(filter_group, "bins", this % bins)
-  end subroutine to_statepoint_mu
-
-  function text_label_mu(this, bin) result(label)
-    class(MuFilter), intent(in) :: this
-    integer,         intent(in) :: bin
-    character(MAX_LINE_LEN)     :: label
-
-    real(8) :: E0, E1
-
-    E0 = this % bins(bin)
-    E1 = this % bins(bin + 1)
-    label = "Change-in-Angle [" // trim(to_str(E0)) // ", " &
-         // trim(to_str(E1)) // ")"
-  end function text_label_mu
-
-end module tally_filter_mu

src/tallies/tally_filter_particle.F90

@@ -52,8 +52,8 @@ contains
 
     do i = 1, this % n_bins
       if (this % particles(i) == p % type) then
-        call match % bins % push_back(i)
-        call match % weights % push_back(ONE)
+        call match % bins_push_back(i)
+        call match % weights_push_back(ONE)
       end if
     end do
   end subroutine get_all_bins

src/tallies/tally_filter_polar.F90

@@ -1,118 +0,0 @@
-module tally_filter_polar
-
-  use, intrinsic :: ISO_C_BINDING
-
-  use algorithm,          only: binary_search
-  use constants
-  use error,              only: fatal_error
-  use hdf5_interface
-  use particle_header,    only: Particle
-  use string,             only: to_str
-  use tally_filter_header
-  use xml_interface
-
-  implicit none
-  private
-
-!===============================================================================
-! POLARFILTER bins the incident neutron polar angle (relative to the global
-! z-axis).
-!===============================================================================
-
-  type, public, extends(TallyFilter) :: PolarFilter
-    real(8), allocatable :: bins(:)
-  contains
-    procedure :: from_xml
-    procedure :: get_all_bins => get_all_bins_polar
-    procedure :: to_statepoint => to_statepoint_polar
-    procedure :: text_label => text_label_polar
-  end type PolarFilter
-
-contains
-
-  subroutine from_xml(this, node)
-    class(PolarFilter), intent(inout) :: this
-    type(XMLNode), intent(in) :: node
-
-    integer :: i
-    integer :: n_angle
-    integer :: n
-    real(8) :: d_angle
-
-    n = node_word_count(node, "bins")
-
-    ! Allocate and store bins
-    this % n_bins = n - 1
-    allocate(this % bins(n))
-    call get_node_array(node, "bins", this % bins)
-
-    ! Allow a user to input a lone number which will mean that you
-    ! subdivide [0,pi] evenly with the input being the number of bins
-    if (n == 1) then
-      n_angle = int(this % bins(1))
-      if (n_angle > 1) then
-        this % n_bins = n_angle
-        d_angle = PI / real(n_angle,8)
-        deallocate(this % bins)
-        allocate(this % bins(n_angle + 1))
-        do i = 1, n_angle
-          this % bins(i) = (i - 1) * d_angle
-        end do
-        this % bins(n_angle + 1) = PI
-      else
-        call fatal_error("Number of bins for polar filter must be&
-             & greater than 1.")
-      end if
-    end if
-  end subroutine from_xml
-
-  subroutine get_all_bins_polar(this, p, estimator, match)
-    class(PolarFilter), intent(in)  :: this
-    type(Particle),     intent(in)  :: p
-    integer,            intent(in)  :: estimator
-    type(TallyFilterMatch),  intent(inout) :: match
-
-    integer :: n
-    integer :: bin
-    real(8) :: theta
-
-    n = this % n_bins
-
-    ! Make sure the correct direction vector is used.
-    if (estimator == ESTIMATOR_TRACKLENGTH) then
-      theta = acos(p % coord(1) % uvw(3))
-    else
-      theta = acos(p % last_uvw(3))
-    end if
-
-    ! Search to find polar angle bin.
-    bin = binary_search(this % bins, n + 1, theta)
-    if (bin /= NO_BIN_FOUND) then
-      call match % bins % push_back(bin)
-      call match % weights % push_back(ONE)
-    end if
-  end subroutine get_all_bins_polar
-
-  subroutine to_statepoint_polar(this, filter_group)
-    class(PolarFilter), intent(in) :: this
-    integer(HID_T),     intent(in) :: filter_group
-
-    call write_dataset(filter_group, "type", "polar")
-    call write_dataset(filter_group, "n_bins", this % n_bins)
-    call write_dataset(filter_group, "bins", this % bins)
-  end subroutine to_statepoint_polar
-
-  function text_label_polar(this, bin) result(label)
-    class(PolarFilter), intent(in) :: this
-    integer,            intent(in) :: bin
-    character(MAX_LINE_LEN)        :: label
-
-    real(8) :: E0, E1
-
-    E0 = this % bins(bin)
-    E1 = this % bins(bin + 1)
-    label = "Polar Angle [" // trim(to_str(E0)) // ", " // trim(to_str(E1)) &
-         // ")"
-  end function text_label_polar
-
-end module tally_filter_polar

src/tallies/tally_filter_sph_harm.F90

@@ -2,241 +2,32 @@ module tally_filter_sph_harm
 
   use, intrinsic :: ISO_C_BINDING
 
-  use constants
-  use error
-  use hdf5_interface
-  use math,                only: calc_pn, calc_rn
-  use particle_header,     only: Particle
-  use string,              only: to_str, to_lower, to_f_string
-  use tally_filter_header
-  use xml_interface
+  use tally_filter_cpp
 
   implicit none
   private
-  public :: openmc_sphharm_filter_get_order
-  public :: openmc_sphharm_filter_get_cosine
-  public :: openmc_sphharm_filter_set_order
-  public :: openmc_sphharm_filter_set_cosine
 
   integer, public, parameter :: COSINE_SCATTER = 1
   integer, public, parameter :: COSINE_PARTICLE = 2
 
-!===============================================================================
-! SPHERICALHARMONICSFILTER gives spherical harmonics expansion moments of a
-! tally score
-!===============================================================================
-
-  type, public, extends(TallyFilter) :: SphericalHarmonicsFilter
-    integer :: order
-    integer :: cosine
+  type, public, extends(CppTallyFilter) :: SphericalHarmonicsFilter
   contains
-    procedure :: from_xml
-    procedure :: get_all_bins
-    procedure :: to_statepoint
-    procedure :: text_label
+    procedure :: cosine
   end type SphericalHarmonicsFilter
 
 contains
 
-!===============================================================================
-! SphericalHarmonicsFilter methods
-!===============================================================================
-
-  subroutine from_xml(this, node)
-    class(SphericalHarmonicsFilter), intent(inout) :: this
-    type(XMLNode), intent(in) :: node
-
-    character(MAX_WORD_LEN) :: temp_str
-
-    ! Get specified order
-    call get_node_value(node, "order", this % order)
-    this % n_bins = (this % order + 1)**2
-
-    ! Determine how cosine term is to be treated
-    if (check_for_node(node, "cosine")) then
-      call get_node_value(node, "cosine", temp_str)
-      select case (to_lower(temp_str))
-      case ('scatter')
-        this % cosine = COSINE_SCATTER
-      case ('particle')
-        this % cosine = COSINE_PARTICLE
-      end select
-    else
-      this % cosine = COSINE_PARTICLE
-    end if
-  end subroutine from_xml
-
-  subroutine get_all_bins(this, p, estimator, match)
-    class(SphericalHarmonicsFilter), intent(in)  :: this
-    type(Particle),      intent(in)  :: p
-    integer,             intent(in)  :: estimator
-    type(TallyFilterMatch),   intent(inout) :: match
-
-    integer :: i, j, n
-    integer :: num_nm
-    real(C_DOUBLE) :: wgt(this % order + 1)
-    real(C_DOUBLE) :: rn(this % n_bins)
-
-    ! Determine cosine term for scatter expansion if necessary
-    if (this % cosine == COSINE_SCATTER) then
-      call calc_pn(this % order, p % mu, wgt)
-    else
-      wgt = ONE
-    end if
-
-    ! Find the Rn,m values
-    call calc_rn(this % order, p % last_uvw, rn)
-
-    j = 0
-    do n = 0, this % order
-      ! Calculate n-th order spherical harmonics for (u,v,w)
-      num_nm = 2*n + 1
-
-      ! Append matching (bin,weight) for each moment
-      do i = 1, num_nm
-        j = j + 1
-        call match % bins % push_back(j)
-        call match % weights % push_back(wgt(n + 1) * rn(j))
-      end do
-    end do
-
-  end subroutine get_all_bins
-
-  subroutine to_statepoint(this, filter_group)
-    class(SphericalHarmonicsFilter), intent(in) :: this
-    integer(HID_T),      intent(in) :: filter_group
-
-    call write_dataset(filter_group, "type", "sphericalharmonics")
-    call write_dataset(filter_group, "n_bins", this % n_bins)
-    call write_dataset(filter_group, "order", this % order)
-    if (this % cosine == COSINE_SCATTER) then
-      call write_dataset(filter_group, "cosine", "scatter")
-    else
-      call write_dataset(filter_group, "cosine", "particle")
-    end if
-  end subroutine to_statepoint
-
-  function text_label(this, bin) result(label)
-    class(SphericalHarmonicsFilter), intent(in) :: this
-    integer,             intent(in) :: bin
-    character(MAX_LINE_LEN)         :: label
-
-    integer :: n, m
-
-    do n = 0, this % order
-      if (bin <= (n + 1)**2) then
-        m = (bin - n**2 - 1) - n
-        label = "Spherical harmonic expansion, Y" // trim(to_str(n)) // &
-             "," // trim(to_str(m))
-        exit
-      end if
-    end do
-  end function text_label
-
-!===============================================================================
-!                               C API FUNCTIONS
-!===============================================================================
-
-  function openmc_sphharm_filter_get_order(index, order) result(err) bind(C)
-    ! Get the order of an expansion filter
-    integer(C_INT32_T), value       :: index
-    integer(C_INT),     intent(out) :: order
-    integer(C_INT) :: err
-
-    err = verify_filter(index)
-    if (err == 0) then
-      select type (f => filters(index) % obj)
-      type is (SphericalHarmonicsFilter)
-        order = f % order
-      class default
-        err = E_INVALID_TYPE
-        call set_errmsg("Not a spherical harmonics filter.")
-      end select
-    end if
-  end function openmc_sphharm_filter_get_order
-
-
-  function openmc_sphharm_filter_get_cosine(index, cosine) result(err) bind(C)
-    ! Get the order of an expansion filter
-    integer(C_INT32_T), value :: index
-    character(kind=C_CHAR), intent(out) :: cosine(*)
-    integer(C_INT) :: err
-
-    integer :: i
-    character(10) :: cosine_
-
-    err = verify_filter(index)
-    if (err == 0) then
-      select type (f => filters(index) % obj)
-      type is (SphericalHarmonicsFilter)
-        select case (f % cosine)
-        case (COSINE_SCATTER)
-          cosine_ = 'scatter'
-        case (COSINE_PARTICLE)
-          cosine_ = 'particle'
-        end select
-
-        ! Convert to C string
-        do i = 1, len_trim(cosine_)
-          cosine(i) = cosine_(i:i)
-        end do
-        cosine(len_trim(cosine_) + 1) = C_NULL_CHAR
-
-      class default
-        err = E_INVALID_TYPE
-        call set_errmsg("Not a spherical harmonics filter.")
-      end select
-    end if
-  end function openmc_sphharm_filter_get_cosine
-
-
-  function openmc_sphharm_filter_set_order(index, order) result(err) bind(C)
-    ! Set the order of an expansion filter
-    integer(C_INT32_T), value :: index
-    integer(C_INT),     value :: order
-    integer(C_INT) :: err
-
-    err = verify_filter(index)
-    if (err == 0) then
-      select type (f => filters(index) % obj)
-      type is (SphericalHarmonicsFilter)
-        f % order = order
-        f % n_bins = (order + 1)**2
-      class default
-        err = E_INVALID_TYPE
-        call set_errmsg("Not a spherical harmonics filter.")
-      end select
-    end if
-  end function openmc_sphharm_filter_set_order
-
-
-  function openmc_sphharm_filter_set_cosine(index, cosine) result(err) bind(C)
-    ! Set the cosine parameter
-    integer(C_INT32_T), value :: index
-    character(kind=C_CHAR), intent(in) :: cosine(*)
-    integer(C_INT) :: err
-
-    character(:), allocatable :: cosine_
-
-    ! Convert C string to Fortran string
-    cosine_ = to_f_string(cosine)
-
-    err = verify_filter(index)
-    if (err == 0) then
-      select type (f => filters(index) % obj)
-      type is (SphericalHarmonicsFilter)
-        select case (cosine_)
-        case ('scatter')
-          f % cosine = COSINE_SCATTER
-        case ('particle')
-          f % cosine = COSINE_PARTICLE
-        end select
-
-      class default
-        err = E_INVALID_TYPE
-        call set_errmsg("Not a spherical harmonics filter.")
-      end select
-    end if
-  end function openmc_sphharm_filter_set_cosine
+  function cosine(this) result(val)
+    class(SphericalHarmonicsFilter) :: this
+    integer                         :: val
+    interface
+      function sphharm_filter_get_cosine(filt) result(val) bind(C)
+        import C_PTR, C_INT
+        type(C_PTR), value :: filt
+        integer(C_INT) :: val
+      end function
+    end interface
+    val = sphharm_filter_get_cosine(this % ptr)
+  end function cosine
 
 end module tally_filter_sph_harm

src/tallies/tally_filter_sptl_legendre.F90

@@ -1,225 +0,0 @@
-module tally_filter_sptl_legendre
-
-  use, intrinsic :: ISO_C_BINDING
-
-  use constants
-  use error
-  use hdf5_interface
-  use math,                only: calc_pn
-  use particle_header,     only: Particle
-  use string,              only: to_str, to_lower
-  use tally_filter_header
-  use xml_interface
-
-  implicit none
-  private
-  public :: openmc_spatial_legendre_filter_get_order
-  public :: openmc_spatial_legendre_filter_get_params
-  public :: openmc_spatial_legendre_filter_set_order
-  public :: openmc_spatial_legendre_filter_set_params
-
-  integer, parameter :: AXIS_X = 1
-  integer, parameter :: AXIS_Y = 2
-  integer, parameter :: AXIS_Z = 3
-
-!===============================================================================
-! SPATIALLEGENDREFILTER gives Legendre moments of the particle's normalized
-! position along an axis
-!===============================================================================
-
-  type, public, extends(TallyFilter) :: SpatialLegendreFilter
-    integer(C_INT) :: order
-    integer(C_INT) :: axis
-    real(C_DOUBLE) :: min
-    real(C_DOUBLE) :: max
-  contains
-    procedure :: from_xml
-    procedure :: get_all_bins
-    procedure :: to_statepoint
-    procedure :: text_label
-  end type SpatialLegendreFilter
-
-contains
-
-!===============================================================================
-! SpatialLegendreFilter methods
-!===============================================================================
-
-  subroutine from_xml(this, node)
-    class(SpatialLegendreFilter), intent(inout) :: this
-    type(XMLNode), intent(in) :: node
-
-    character(MAX_WORD_LEN) :: axis
-
-    ! Get attributes from XML
-    call get_node_value(node, "order", this % order)
-    call get_node_value(node, "axis", axis)
-    select case (to_lower(axis))
-    case ('x')
-      this % axis = AXIS_X
-    case ('y')
-      this % axis = AXIS_Y
-    case ('z')
-      this % axis = AXIS_Z
-    end select
-    call get_node_value(node, "min", this % min)
-    call get_node_value(node, "max", this % max)
-
-    this % n_bins = this % order + 1
-  end subroutine from_xml
-
-  subroutine get_all_bins(this, p, estimator, match)
-    class(SpatialLegendreFilter), intent(in)  :: this
-    type(Particle),      intent(in)  :: p
-    integer,             intent(in)  :: estimator
-    type(TallyFilterMatch),   intent(inout) :: match
-
-    integer :: i
-    real(C_DOUBLE) :: wgt(this % n_bins)
-    real(C_DOUBLE) :: x       ! Position on specified axis
-    real(C_DOUBLE) :: x_norm  ! Normalized position
-
-    x = p % coord(1) % xyz(this % axis)
-    if (this % min <= x .and. x <= this % max) then
-      ! Calculate normalized position between min and max
-      x_norm = TWO*(x - this % min)/(this % max - this % min) - ONE
-
-      call calc_pn(this % order, x_norm, wgt)
-      do i = 1, this % n_bins
-        call match % bins % push_back(i)
-        call match % weights % push_back(wgt(i))
-      end do
-    end if
-  end subroutine get_all_bins
-
-  subroutine to_statepoint(this, filter_group)
-    class(SpatialLegendreFilter), intent(in) :: this
-    integer(HID_T),      intent(in) :: filter_group
-    character(kind=C_CHAR) :: axis
-
-    call write_dataset(filter_group, "type", "spatiallegendre")
-    call write_dataset(filter_group, "n_bins", this % n_bins)
-    call write_dataset(filter_group, "order", this % order)
-    select case (this % axis)
-    case (AXIS_X)
-      axis = 'x'
-    case (AXIS_Y)
-      axis = 'y'
-    case (AXIS_Z)
-      axis = 'z'
-    end select
-    call write_dataset(filter_group, 'axis', axis)
-    call write_dataset(filter_group, 'min', this % min)
-    call write_dataset(filter_group, 'max', this % max)
-  end subroutine to_statepoint
-
-  function text_label(this, bin) result(label)
-    class(SpatialLegendreFilter), intent(in) :: this
-    integer,             intent(in) :: bin
-    character(MAX_LINE_LEN)         :: label
-
-    character(1) :: axis
-
-    select case (this % axis)
-    case (AXIS_X)
-      axis = 'x'
-    case (AXIS_Y)
-      axis = 'y'
-    case (AXIS_Z)
-      axis = 'z'
-    end select
-    label = "Legendre expansion, " // axis // " axis, P" // &
-         trim(to_str(bin - 1))
-  end function text_label
-
-!===============================================================================
-!                               C API FUNCTIONS
-!===============================================================================
-
-  function openmc_spatial_legendre_filter_get_order(index, order) result(err) bind(C)
-    ! Get the order of an expansion filter
-    integer(C_INT32_T), value       :: index
-    integer(C_INT),     intent(out) :: order
-    integer(C_INT) :: err
-
-    err = verify_filter(index)
-    if (err == 0) then
-      select type (f => filters(index) % obj)
-      type is (SpatialLegendreFilter)
-        order = f % order
-      class default
-        err = E_INVALID_TYPE
-        call set_errmsg("Not a spatial Legendre filter.")
-      end select
-    end if
-  end function openmc_spatial_legendre_filter_get_order
-
-
-  function openmc_spatial_legendre_filter_set_order(index, order) result(err) bind(C)
-    ! Set the order of an expansion filter
-    integer(C_INT32_T), value :: index
-    integer(C_INT),     value :: order
-    integer(C_INT) :: err
-
-    err = verify_filter(index)
-    if (err == 0) then
-      select type (f => filters(index) % obj)
-      type is (SpatialLegendreFilter)
-        f % order = order
-        f % n_bins = order + 1
-      class default
-        err = E_INVALID_TYPE
-        call set_errmsg("Not a spatial Legendre filter.")
-      end select
-    end if
-  end function openmc_spatial_legendre_filter_set_order
-
-
-  function openmc_spatial_legendre_filter_get_params(index, axis, min, max) &
-       result(err) bind(C)
-    ! Get the parameters for a spatial Legendre filter
-    integer(C_INT32_T), value :: index
-    integer(C_INT), intent(out) :: axis
-    real(C_DOUBLE), intent(out) :: min
-    real(C_DOUBLE), intent(out) :: max
-    integer(C_INT) :: err
-
-    err = verify_filter(index)
-    if (err == 0) then
-      select type(f => filters(index) % obj)
-      type is (SpatialLegendreFilter)
-        axis = f % axis
-        min = f % min
-        max = f % max
-      class default
-        err = E_INVALID_TYPE
-        call set_errmsg("Not a spatial Legendre filter.")
-      end select
-    end if
-  end function openmc_spatial_legendre_filter_get_params
-
-
-  function openmc_spatial_legendre_filter_set_params(index, axis, min, max) &
-       result(err) bind(C)
-    ! Set the parameters for a spatial Legendre filter
-    integer(C_INT32_T), value :: index
-    integer(C_INT), intent(in), optional :: axis
-    real(C_DOUBLE), intent(in), optional :: min
-    real(C_DOUBLE), intent(in), optional :: max
-    integer(C_INT) :: err
-
-    err = verify_filter(index)
-    if (err == 0) then
-      select type(f => filters(index) % obj)
-      type is (SpatialLegendreFilter)
-        if (present(axis)) f % axis = axis
-        if (present(min)) f % min = min
-        if (present(max)) f % max = max
-      class default
-        err = E_INVALID_TYPE
-        call set_errmsg("Not a spatial Legendre filter.")
-      end select
-    end if
-  end function openmc_spatial_legendre_filter_set_params
-
-end module tally_filter_sptl_legendre

src/tallies/tally_filter_surface.F90

@@ -1,109 +0,0 @@
-module tally_filter_surface
-
-  use, intrinsic :: ISO_C_BINDING
-
-  use constants,          only: ONE, MAX_LINE_LEN
-  use dict_header,        only: EMPTY
-  use error,              only: fatal_error
-  use hdf5_interface
-  use surface_header
-  use particle_header,    only: Particle
-  use string,             only: to_str
-  use tally_filter_header
-  use xml_interface
-
-  implicit none
-  private
-
-!===============================================================================
-! SURFACEFILTER specifies which surface particles are crossing
-!===============================================================================
-
-  type, public, extends(TallyFilter) :: SurfaceFilter
-    integer, allocatable :: surfaces(:)
-
-    ! True if this filter is used for surface currents
-    logical              :: current = .false.
-  contains
-    procedure :: from_xml
-    procedure :: get_all_bins => get_all_bins_surface
-    procedure :: to_statepoint => to_statepoint_surface
-    procedure :: text_label => text_label_surface
-    procedure :: initialize => initialize_surface
-  end type SurfaceFilter
-
-contains
-
-  subroutine from_xml(this, node)
-    class(SurfaceFilter), intent(inout) :: this
-    type(XMLNode), intent(in) :: node
-
-    integer :: n
-
-    n = node_word_count(node, "bins")
-
-    ! Allocate and store bins
-    this % n_bins = n
-    allocate(this % surfaces(n))
-    call get_node_array(node, "bins", this % surfaces)
-  end subroutine from_xml
-
-  subroutine get_all_bins_surface(this, p, estimator, match)
-    class(SurfaceFilter), intent(in)  :: this
-    type(Particle),       intent(in)  :: p
-    integer,              intent(in)  :: estimator
-    type(TallyFilterMatch),    intent(inout) :: match
-
-    integer :: i
-
-      do i = 1, this % n_bins
-        if (abs(p % surface) == this % surfaces(i)) then
-          call match % bins % push_back(i)
-          if (p % surface < 0) then
-            call match % weights % push_back(-ONE)
-          else
-            call match % weights % push_back(ONE)
-          end if
-          exit
-        end if
-      end do
-
-  end subroutine get_all_bins_surface
-
-  subroutine to_statepoint_surface(this, filter_group)
-    class(SurfaceFilter), intent(in) :: this
-    integer(HID_T),       intent(in) :: filter_group
-
-    call write_dataset(filter_group, "type", "surface")
-    call write_dataset(filter_group, "n_bins", this % n_bins)
-    call write_dataset(filter_group, "bins", this % surfaces)
-  end subroutine to_statepoint_surface
-
-  subroutine initialize_surface(this)
-    class(SurfaceFilter), intent(inout) :: this
-
-    integer :: i, id
-    integer :: val
-
-    ! Convert ids to indices.
-    do i = 1, this % n_bins
-      id = this % surfaces(i)
-      val = surface_dict % get(id)
-      if (val /= EMPTY) then
-        this % surfaces(i) = val
-      else
-        call fatal_error("Could not find surface " // trim(to_str(id)) &
-             &// " specified on tally filter.")
-      end if
-    end do
-  end subroutine initialize_surface
-
-  function text_label_surface(this, bin) result(label)
-    class(SurfaceFilter), intent(in) :: this
-    integer,              intent(in) :: bin
-    character(MAX_LINE_LEN)          :: label
-
-    label = "Surface " // to_str(surfaces(this % surfaces(bin)) % id())
-  end function text_label_surface
-
-end module tally_filter_surface

src/tallies/tally_filter_universe.F90

@@ -1,118 +0,0 @@
-module tally_filter_universe
-
-  use, intrinsic :: ISO_C_BINDING
-
-  use constants,          only: ONE, MAX_LINE_LEN
-  use dict_header,        only: EMPTY
-  use error,              only: fatal_error
-  use hdf5_interface
-  use geometry_header
-  use particle_header,    only: Particle
-  use string,             only: to_str
-  use tally_filter_header
-  use xml_interface
-
-  implicit none
-  private
-
-!===============================================================================
-! UNIVERSEFILTER specifies which geometric universes tally events reside in.
-!===============================================================================
-
-  type, public, extends(TallyFilter) :: UniverseFilter
-    integer, allocatable :: universes(:)
-    type(DictIntInt)     :: map
-  contains
-    procedure :: from_xml
-    procedure :: get_all_bins => get_all_bins_universe
-    procedure :: to_statepoint => to_statepoint_universe
-    procedure :: text_label => text_label_universe
-    procedure :: initialize => initialize_universe
-  end type UniverseFilter
-
-contains
-
-  subroutine from_xml(this, node)
-    class(UniverseFilter), intent(inout) :: this
-    type(XMLNode), intent(in) :: node
-
-    integer :: n
-
-    n = node_word_count(node, "bins")
-
-    ! Allocate and store bins
-    this % n_bins = n
-    allocate(this % universes(n))
-    call get_node_array(node, "bins", this % universes)
-  end subroutine from_xml
-
-  subroutine get_all_bins_universe(this, p, estimator, match)
-    class(UniverseFilter), intent(in)  :: this
-    type(Particle),        intent(in)  :: p
-    integer,               intent(in)  :: estimator
-    type(TallyFilterMatch),     intent(inout) :: match
-
-    integer :: i
-    integer :: val
-
-    ! Iterate over coordinate levels to see which universes match
-    do i = 1, p % n_coord
-      val = this % map % get(p % coord(i) % universe)
-      if (val /= EMPTY) then
-        call match % bins % push_back(val)
-        call match % weights % push_back(ONE)
-      end if
-    end do
-
-  end subroutine get_all_bins_universe
-
-  subroutine to_statepoint_universe(this, filter_group)
-    class(UniverseFilter), intent(in) :: this
-    integer(HID_T),        intent(in) :: filter_group
-
-    integer :: i
-    integer, allocatable :: universe_ids(:)
-
-    call write_dataset(filter_group, "type", "universe")
-    call write_dataset(filter_group, "n_bins", this % n_bins)
-
-    allocate(universe_ids(size(this % universes)))
-    do i = 1, size(this % universes)
-      universe_ids(i) = universe_id(this % universes(i))
-    end do
-    call write_dataset(filter_group, "bins", universe_ids)
-  end subroutine to_statepoint_universe
-
-  subroutine initialize_universe(this)
-    class(UniverseFilter), intent(inout) :: this
-
-    integer :: i, id
-    integer :: val
-
-    ! Convert ids to indices.
-    do i = 1, this % n_bins
-      id = this % universes(i)
-      val = universe_dict % get(id)
-      if (val /= EMPTY) then
-        this % universes(i) = val
-      else
-        call fatal_error("Could not find universe " // trim(to_str(id)) &
-             &// " specified on a tally filter.")
-      end if
-    end do
-
-    ! Generate mapping from universe indices to filter bins.
-    do i = 1, this % n_bins
-      call this % map % set(this % universes(i), i)
-    end do
-  end subroutine initialize_universe
-
-  function text_label_universe(this, bin) result(label)
-    class(UniverseFilter), intent(in) :: this
-    integer,               intent(in) :: bin
-    character(MAX_LINE_LEN)           :: label
-
-    label = "Universe " // to_str(universe_id(this % universes(bin)))
-  end function text_label_universe
-
-end module tally_filter_universe

src/tallies/tally_filter_zernike.F90

@@ -1,295 +0,0 @@
-module tally_filter_zernike
-
-  use, intrinsic :: ISO_C_BINDING
-
-  use constants
-  use error
-  use hdf5_interface
-  use math,                only: calc_zn, calc_zn_rad
-  use particle_header,     only: Particle
-  use string,              only: to_str
-  use tally_filter_header
-  use xml_interface
-
-  implicit none
-  private
-  public :: openmc_zernike_filter_get_order
-  public :: openmc_zernike_filter_get_params
-  public :: openmc_zernike_filter_set_order
-  public :: openmc_zernike_filter_set_params
-
-!===============================================================================
-! ZERNIKEFILTER gives Zernike polynomial moments of a particle's position
-!===============================================================================
-
-  type, public, extends(TallyFilter) :: ZernikeFilter
-    integer :: order
-    real(8) :: x
-    real(8) :: y
-    real(8) :: r
-  contains
-    procedure :: calc_n_bins => calc_n_bins_zn
-    procedure :: from_xml => from_xml_zn
-    procedure :: get_all_bins => get_all_bins_zn
-    procedure :: to_statepoint => to_statepoint_zn
-    procedure :: text_label => text_label_zn
-  end type ZernikeFilter
-
-!===============================================================================
-! ZERNIKERADIALFILTER gives even order radial Zernike polynomial moments of a
-! particle's position
-!===============================================================================
-
-  type, public, extends(ZernikeFilter) :: ZernikeRadialFilter
-  contains
-    procedure :: calc_n_bins => calc_n_bins_zn_rad
-    ! Inherit from_xml from ZernikeFilter
-    procedure :: get_all_bins => get_all_bins_zn_rad
-    procedure :: to_statepoint => to_statepoint_zn_rad
-    procedure :: text_label => text_label_zn_rad
-  end type ZernikeRadialFilter
-
-contains
-
-!===============================================================================
-! ZernikeFilter methods
-!===============================================================================
-
-  function calc_n_bins_zn(this) result(n_bins)
-    class(ZernikeFilter), intent(in) :: this
-    integer :: n
-    integer :: n_bins
-
-    n = this % order
-    n_bins = ((n+1) * (n+2))/2
-  end function calc_n_bins_zn
-
-  subroutine from_xml_zn(this, node)
-    class(ZernikeFilter), intent(inout) :: this
-    type(XMLNode), intent(in) :: node
-
-    integer :: n
-
-    ! Get center of cylinder and radius
-    call get_node_value(node, "x", this % x)
-    call get_node_value(node, "y", this % y)
-    call get_node_value(node, "r", this % r)
-
-    ! Get specified order
-    call get_node_value(node, "order", n)
-    this % order = n
-    this % n_bins = this % calc_n_bins()
-  end subroutine from_xml_zn
-
-  subroutine get_all_bins_zn(this, p, estimator, match)
-    class(ZernikeFilter), intent(in)  :: this
-    type(Particle),      intent(in)  :: p
-    integer,             intent(in)  :: estimator
-    type(TallyFilterMatch),   intent(inout) :: match
-
-    integer :: i
-    real(8) :: x, y, r, theta
-    real(C_DOUBLE) :: zn(this % n_bins)
-
-    ! Determine normalized (r,theta) positions
-    x = p % coord(1) % xyz(1) - this % x
-    y = p % coord(1) % xyz(2) - this % y
-    r = sqrt(x*x + y*y)/this % r
-    if (r <= ONE) then
-      theta = atan2(y, x)
-      
-      ! Get moments for Zernike polynomial orders 0..n
-      call calc_zn(this % order, r, theta, zn)
-
-      do i = 1, this % n_bins
-        call match % bins % push_back(i)
-        call match % weights % push_back(zn(i))
-      end do
-    endif
-  end subroutine get_all_bins_zn
-
-  subroutine to_statepoint_zn(this, filter_group)
-    class(ZernikeFilter), intent(in) :: this
-    integer(HID_T),      intent(in) :: filter_group
-
-    call write_dataset(filter_group, "type", "zernike")
-    call write_dataset(filter_group, "n_bins", this % n_bins)
-    call write_dataset(filter_group, "order", this % order)
-    call write_dataset(filter_group, "x", this % x)
-    call write_dataset(filter_group, "y", this % y)
-    call write_dataset(filter_group, "r", this % r)
-  end subroutine to_statepoint_zn
-
-  function text_label_zn(this, bin) result(label)
-    class(ZernikeFilter), intent(in) :: this
-    integer,             intent(in) :: bin
-    character(MAX_LINE_LEN)         :: label
-
-    integer :: n, m
-    integer :: first, last
-
-    do n = 0, this % order
-      last = (n + 1)*(n + 2)/2
-      if (bin <= last) then
-        first = last - n
-        m = -n + (bin - first)*2
-        label = "Zernike expansion, Z" // trim(to_str(n)) // "," &
-             // trim(to_str(m))
-        exit
-      end if
-    end do
-  end function text_label_zn
-
-!===============================================================================
-! ZernikeRadialFilter methods
-!===============================================================================
-
-  function calc_n_bins_zn_rad(this) result(num_n_bins)
-    class(ZernikeRadialFilter), intent(in) :: this
-    integer :: n
-    integer :: num_n_bins
-
-    n = this % order
-    num_n_bins = n/2 + 1
-  end function calc_n_bins_zn_rad
-
-  subroutine get_all_bins_zn_rad(this, p, estimator, match)
-    class(ZernikeRadialFilter), intent(in)  :: this
-    type(Particle),      intent(in)  :: p
-    integer,             intent(in)  :: estimator
-    type(TallyFilterMatch),   intent(inout) :: match
-
-    integer :: i
-    real(8) :: x, y, r
-    real(C_DOUBLE) :: zn_rad(this % n_bins)
-
-    ! Determine normalized (r,theta) positions
-    x = p % coord(1) % xyz(1) - this % x
-    y = p % coord(1) % xyz(2) - this % y
-    r = sqrt(x*x + y*y)/this % r
-    if (r <= ONE) then
-
-      ! Get moments for even order Zernike polynomial orders 0..n
-      call calc_zn_rad(this % order, r, zn_rad)
-
-      do i = 1, this % n_bins
-        call match % bins % push_back(i)
-        call match % weights % push_back(zn_rad(i))
-      end do
-    endif
-  end subroutine get_all_bins_zn_rad
-
-  subroutine to_statepoint_zn_rad(this, filter_group)
-    class(ZernikeRadialFilter), intent(in) :: this
-    integer(HID_T),      intent(in) :: filter_group
-
-    call write_dataset(filter_group, "type", "zernikeradial")
-    call write_dataset(filter_group, "n_bins", this % n_bins)
-    call write_dataset(filter_group, "order", this % order)
-    call write_dataset(filter_group, "x", this % x)
-    call write_dataset(filter_group, "y", this % y)
-    call write_dataset(filter_group, "r", this % r)
-  end subroutine to_statepoint_zn_rad
-
-  function text_label_zn_rad(this, bin) result(label)
-    class(ZernikeRadialFilter), intent(in) :: this
-    integer,             intent(in) :: bin
-    character(MAX_LINE_LEN)         :: label
-
-    label = "Zernike expansion, Z" // trim(to_str(2*(bin-1))) // ",0"
-  end function text_label_zn_rad
-
-!===============================================================================
-!                               C API FUNCTIONS
-!===============================================================================
-
-  function openmc_zernike_filter_get_order(index, order) result(err) bind(C)
-    ! Get the order of an expansion filter
-    integer(C_INT32_T), value       :: index
-    integer(C_INT),     intent(out) :: order
-    integer(C_INT) :: err
-
-    err = verify_filter(index)
-    if (err == 0) then
-      select type (f => filters(index) % obj)
-      type is (ZernikeFilter)
-        order = f % order
-      type is (ZernikeRadialFilter)
-        order = f % order
-      class default
-        err = E_INVALID_TYPE
-        call set_errmsg("Not a Zernike filter.")
-      end select
-    end if
-  end function openmc_zernike_filter_get_order
-
-
-  function openmc_zernike_filter_get_params(index, x, y, r) result(err) bind(C)
-    ! Get the Zernike filter parameters
-    integer(C_INT32_T), value :: index
-    real(C_DOUBLE), intent(out) :: x
-    real(C_DOUBLE), intent(out) :: y
-    real(C_DOUBLE), intent(out) :: r
-    integer(C_INT) :: err
-
-    err = verify_filter(index)
-    if (err == 0) then
-      select type (f => filters(index) % obj)
-      class is (ZernikeFilter)
-        x = f % x
-        y = f % y
-        r = f % r
-      class default
-        err = E_INVALID_TYPE
-        call set_errmsg("Not a Zernike filter.")
-      end select
-    end if
-  end function openmc_zernike_filter_get_params
-
-
-  function openmc_zernike_filter_set_order(index, order) result(err) bind(C)
-    ! Set the order of an expansion filter
-    integer(C_INT32_T), value :: index
-    integer(C_INT),     value :: order
-    integer(C_INT) :: err
-
-    err = verify_filter(index)
-    if (err == 0) then
-      select type (f => filters(index) % obj)
-      type is (ZernikeFilter)
-        f % order = order
-        f % n_bins = ((order + 1)*(order + 2))/2
-      type is (ZernikeRadialFilter)
-        f % order = order
-        f % n_bins = order/2 + 1
-      class default
-        err = E_INVALID_TYPE
-        call set_errmsg("Not a Zernike filter.")
-      end select
-    end if
-  end function openmc_zernike_filter_set_order
-
-
-  function openmc_zernike_filter_set_params(index, x, y, r) result(err) bind(C)
-    ! Set the Zernike filter parameters
-    integer(C_INT32_T), value :: index
-    real(C_DOUBLE), intent(in), optional :: x
-    real(C_DOUBLE), intent(in), optional :: y
-    real(C_DOUBLE), intent(in), optional :: r
-    integer(C_INT) :: err
-
-    err = verify_filter(index)
-    if (err == 0) then
-      select type (f => filters(index) % obj)
-      class is (ZernikeFilter)
-        if (present(x)) f % x = x
-        if (present(y)) f % y = y
-        if (present(r)) f % r = r
-      class default
-        err = E_INVALID_TYPE
-        call set_errmsg("Not a Zernike filter.")
-      end select
-    end if
-  end function openmc_zernike_filter_set_params
-
-end module tally_filter_zernike

src/tallies/tally_header.F90

@@ -344,7 +344,7 @@ contains
         this % estimator = ESTIMATOR_ANALOG
       type is (SphericalHarmonicsFilter)
         j = FILTER_SPH_HARMONICS
-        if (filt % cosine == COSINE_SCATTER) then
+        if (filt % cosine() == COSINE_SCATTER) then
           this % estimator = ESTIMATOR_ANALOG
         end if
       type is (SpatialLegendreFilter)
@@ -413,6 +413,13 @@ contains
 !===============================================================================
 
   subroutine free_memory_tally()
+    interface
+      subroutine free_memory_tally_c() bind(C)
+      end subroutine free_memory_tally_c
+    end interface
+
+    call free_memory_tally_c()
+
     n_tallies = 0
     if (allocated(tallies)) deallocate(tallies)
     call tally_dict % clear()

src/tallies/trigger.F90

@@ -173,8 +173,8 @@ contains
 
             ! Initialize bins, filter level
             do j = 1, size(t % filter)
-              call filter_matches(t % filter(j)) % bins % clear()
-              call filter_matches(t % filter(j)) % bins % push_back(0)
+              call filter_matches(t % filter(j)) % bins_clear()
+              call filter_matches(t % filter(j)) % bins_push_back(0)
             end do
 
             FILTER_LOOP: do filter_index = 1, t % n_filter_bins
@@ -264,13 +264,13 @@ contains
     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
     do j = 1, size(t % filter)
-      call filter_matches(t % filter(j)) % bins % clear()
-      call filter_matches(t % filter(j)) % bins % push_back(1)
+      call filter_matches(t % filter(j)) % bins_clear()
+      call filter_matches(t % filter(j)) % bins_push_back(1)
     end do
 
     ! determine how many energyin bins there are
@@ -296,20 +296,20 @@ contains
 
       ! Get the indices for this cell
       call m % get_indices_from_bin(i, ijk)
-      filter_matches(i_filter_mesh) % bins % data(1) = i
+      call filter_matches(i_filter_mesh) % bins_set_data(1, i)
 
       do l = 1, n
 
         if (print_ebin) then
-          filter_matches(i_filter_ein) % bins % data(1) = l
+          call filter_matches(i_filter_ein) % bins_set_data(1, l)
         end if
 
         ! Left Surface
-        filter_matches(i_filter_surf) % bins % data(1) = OUT_LEFT
+        call filter_matches(i_filter_surf) % bins_set_data(1, OUT_LEFT)
         filter_index = 1
         do j = 1, size(t % filter)
           filter_index = filter_index + (filter_matches(t % filter(j)) % &
-               bins % data(1) - 1) * t % stride(j)
+               bins_data(1) - 1) * t % stride(j)
         end do
         call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
         if (trigger % std_dev < std_dev) then
@@ -321,11 +321,11 @@ contains
         trigger % variance = std_dev**2
 
         ! Right Surface
-        filter_matches(i_filter_surf) % bins % data(1) = OUT_RIGHT
+        call filter_matches(i_filter_surf) % bins_set_data(1, OUT_RIGHT)
         filter_index = 1
         do j = 1, size(t % filter)
           filter_index = filter_index + (filter_matches(t % filter(j)) % &
-               bins % data(1) - 1) * t % stride(j)
+               bins_data(1) - 1) * t % stride(j)
         end do
         call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
         if (trigger % std_dev < std_dev) then
@@ -337,11 +337,11 @@ contains
         trigger % variance = trigger % std_dev**2
 
         ! Back Surface
-        filter_matches(i_filter_surf) % bins % data(1) = OUT_BACK
+        call filter_matches(i_filter_surf) % bins_set_data(1, OUT_BACK)
         filter_index = 1
         do j = 1, size(t % filter)
           filter_index = filter_index + (filter_matches(t % filter(j)) % &
-               bins % data(1) - 1) * t % stride(j)
+               bins_data(1) - 1) * t % stride(j)
         end do
         call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
         if (trigger % std_dev < std_dev) then
@@ -353,11 +353,11 @@ contains
         trigger % variance = trigger % std_dev**2
 
         ! Front Surface
-        filter_matches(i_filter_surf) % bins % data(1) = OUT_FRONT
+        call filter_matches(i_filter_surf) % bins_set_data(1, OUT_FRONT)
         filter_index = 1
         do j = 1, size(t % filter)
           filter_index = filter_index + (filter_matches(t % filter(j)) % &
-               bins % data(1) - 1) * t % stride(j)
+               bins_data(1) - 1) * t % stride(j)
         end do
         call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
         if (trigger % std_dev < std_dev) then
@@ -369,11 +369,11 @@ contains
         trigger % variance = trigger % std_dev**2
 
         ! Bottom Surface
-        filter_matches(i_filter_surf) % bins % data(1) = OUT_BOTTOM
+        call filter_matches(i_filter_surf) % bins_set_data(1, OUT_BOTTOM)
         filter_index = 1
         do j = 1, size(t % filter)
           filter_index = filter_index + (filter_matches(t % filter(j)) % &
-               bins % data(1) - 1) * t % stride(j)
+               bins_data(1) - 1) * t % stride(j)
         end do
         call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
         if (trigger % std_dev < std_dev) then
@@ -385,11 +385,11 @@ contains
         trigger % variance = trigger % std_dev**2
 
         ! Top Surface
-        filter_matches(i_filter_surf) % bins % data(1) = OUT_TOP
+        call filter_matches(i_filter_surf) % bins_set_data(1, OUT_TOP)
         filter_index = 1
         do j = 1, size(t % filter)
           filter_index = filter_index + (filter_matches(t % filter(j)) % &
-               bins % data(1) - 1) * t % stride(j)
+               bins_data(1) - 1) * t % stride(j)
         end do
         call get_trigger_uncertainty(std_dev, rel_err, 1, filter_index, t)
         if (trigger % std_dev < std_dev) then

tests/regression_tests/dagmc/test.py

@@ -1,7 +1,7 @@
 from tests.testing_harness import TestHarness
 import os
 import pytest
-import openmc
+import openmc.capi
 
 pytestmark = pytest.mark.skipif(
     not openmc.capi.dagmc_enabled,
@@ -9,4 +9,4 @@ pytestmark = pytest.mark.skipif(
 
 def test_dagmc():
     harness = TestHarness('statepoint.5.h5')
-    harness.main()    
+    harness.main()

tests/unit_tests/test_capi.py

@@ -168,6 +168,8 @@ def test_tally(capi_init):
         openmc.capi.MaterialFilter(uid=1)
     mats = openmc.capi.materials
     f = openmc.capi.MaterialFilter([mats[2], mats[1]])
+    assert f.bins[0] == mats[2]
+    assert f.bins[1] == mats[1]
     t.filters = [f]
     assert t.filters == [f]
 

tools/ci/download-xs.sh

@@ -14,6 +14,6 @@ fi
 
 # Download multipole library
 if [[ ! -e $HOME/WMP_Library/092235.h5 ]]; then
-    wget -q https://github.com/mit-crpg/WMP_Library/releases/download/v1.0/WMP_Library_v1.0.tar.gz
-    tar -C $HOME -xzf WMP_Library_v1.0.tar.gz
+    wget -q https://github.com/mit-crpg/WMP_Library/releases/download/v1.1/WMP_Library_v1.1.tar.gz
+    tar -C $HOME -xzf WMP_Library_v1.1.tar.gz
 fi