Merge pull request #1110 from smharper/cpp_tallies

Move remaining tally filters to C++

CMakeLists.txt

@@ -373,7 +373,6 @@ 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_cpp.F90
   src/tallies/tally_filter_delayedgroup.F90
   src/tallies/tally_filter_distribcell.F90
   src/tallies/tally_filter_energy.F90
@@ -436,13 +435,16 @@ add_library(libopenmc SHARED
   src/tallies/filter_cellborn.cpp
   src/tallies/filter_cellfrom.cpp
   src/tallies/filter_cell.cpp
+  src/tallies/filter_delayedgroup.cpp
   src/tallies/filter_distribcell.cpp
   src/tallies/filter_energyfunc.cpp
+  src/tallies/filter_energy.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_particle.cpp
   src/tallies/filter_polar.cpp
   src/tallies/filter_sph_harm.cpp
   src/tallies/filter_sptl_legendre.cpp

include/openmc/constants.h

@@ -420,6 +420,7 @@ constexpr int DIFF_NUCLIDE_DENSITY {2};
 constexpr int DIFF_TEMPERATURE {3};
 
 constexpr int C_NONE {-1};
+constexpr int F90_NONE {0}; //TODO: replace usage of this with C_NONE
 
 // Interpolation rules
 enum class Interpolation {

include/openmc/tallies/filter.h

@@ -2,6 +2,7 @@
 #define OPENMC_TALLIES_FILTER_H
 
 #include <cstdint>
+#include <memory>
 #include <string>
 #include <vector>
 
@@ -23,7 +24,7 @@ extern std::vector<FilterMatch> filter_matches;
 #pragma omp threadprivate(filter_matches)
 
 class Filter;
-extern std::vector<Filter*> tally_filters;
+extern std::vector<std::unique_ptr<Filter>> tally_filters;
 
 //==============================================================================
 //! Stores bins and weights for filtered tally events.

include/openmc/tallies/filter_delayedgroup.h

@@ -0,0 +1,37 @@
+#ifndef OPENMC_TALLIES_FILTER_DELAYEDGROUP_H
+#define OPENMC_TALLIES_FILTER_DELAYEDGROUP_H
+
+#include <vector>
+
+#include "openmc/tallies/filter.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Bins outgoing fission neutrons in their delayed groups.
+//!
+//! The get_all_bins functionality is not actually used.  The bins are manually
+//! iterated over in the scoring subroutines.
+//==============================================================================
+
+class DelayedGroupFilter : public Filter
+{
+public:
+  ~DelayedGroupFilter() = default;
+
+  std::string type() const override {return "delayedgroup";}
+
+  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<int> groups_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_DELAYEDGROUP_H

include/openmc/tallies/filter_energy.h

@@ -0,0 +1,55 @@
+#ifndef OPENMC_TALLIES_FILTER_ENERGY_H
+#define OPENMC_TALLIES_FILTER_ENERGY_H
+
+#include <vector>
+
+#include "openmc/tallies/filter.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Bins the incident neutron energy.
+//==============================================================================
+
+class EnergyFilter : public Filter
+{
+public:
+  ~EnergyFilter() = default;
+
+  std::string type() const override {return "energy";}
+
+  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_;
+
+  //! True if transport group number can be used directly to get bin number
+  bool matches_transport_groups_ {false};
+};
+
+//==============================================================================
+//! Bins the outgoing neutron energy.
+//!
+//! Only scattering events use the get_all_bins functionality.  Nu-fission
+//! tallies manually iterate over the filter bins.
+//==============================================================================
+
+class EnergyoutFilter : public EnergyFilter
+{
+public:
+  std::string type() const override {return "energyout";}
+
+  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_ENERGY_H

include/openmc/tallies/filter_particle.h

@@ -0,0 +1,34 @@
+#ifndef OPENMC_TALLIES_FILTER_PARTICLE_H
+#define OPENMC_TALLIES_FILTER_PARTICLE_H
+
+#include <vector>
+
+#include "openmc/tallies/filter.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Bins by type of particle (e.g. neutron, photon).
+//==============================================================================
+
+class ParticleFilter : public Filter
+{
+public:
+  ~ParticleFilter() = default;
+
+  std::string type() const override {return "particle";}
+
+  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<int> particles_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLIES_FILTER_PARTICLE_H

src/geometry.cpp

@@ -16,8 +16,6 @@ namespace openmc {
 
 std::vector<int64_t> overlap_check_count;
 
-constexpr int F90_NONE {0}; //TODO: replace usage of this with C_NONE
-
 //==============================================================================
 
 extern "C" bool

src/geometry_aux.cpp

@@ -187,9 +187,9 @@ prepare_distribcell()
 {
   // Find all cells listed in a DistribcellFilter.
   std::unordered_set<int32_t> distribcells;
-  for (auto* filt : tally_filters) {
-    if (filt->type() == "distribcell") {
-      auto* distrib_filt = static_cast<DistribcellFilter*>(filt);
+  for (auto& filt : tally_filters) {
+    auto* distrib_filt = dynamic_cast<DistribcellFilter*>(filt.get());
+    if (distrib_filt) {
       distribcells.insert(distrib_filt->cell_);
     }
   }

src/mgxs_interface.cpp

@@ -99,11 +99,11 @@ void read_mg_cross_sections_header_c(hid_t file_id)
   read_attribute(file_id, "energy_groups", num_energy_groups);
 
   ensure_exists(file_id, "group structure", true);
-  read_attribute(file_id, "group structure", energy_bins);
+  read_attribute(file_id, "group structure", rev_energy_bins);
 
-  // Create reverse energy bins
-  std::copy(energy_bins.crbegin(), energy_bins.crend(),
-    std::back_inserter(rev_energy_bins));
+  // Reverse energy bins
+  std::copy(rev_energy_bins.crbegin(), rev_energy_bins.crend(),
+    std::back_inserter(energy_bins));
 
   // Create average energies
   for (int i = 0; i < energy_bins.size() - 1; ++i) {

src/source.cpp

@@ -292,8 +292,6 @@ void initialize_source()
   }
 }
 
-extern "C" double* rev_energy_bins_ptr();
-
 Bank sample_external_source()
 {
   // Set the random number generator to the source stream.
@@ -320,11 +318,8 @@ Bank sample_external_source()
 
   // If running in MG, convert site % E to group
   if (!settings::run_CE) {
-    // Get pointer to rev_energy_bins array on Fortran side
-    double* rev_energy_bins = rev_energy_bins_ptr();
-
-    int n = num_energy_groups + 1;
-    site.E = lower_bound_index(rev_energy_bins, rev_energy_bins + n, site.E);
+    site.E = lower_bound_index(rev_energy_bins.begin(), rev_energy_bins.end(),
+                               site.E);
     site.E = num_energy_groups - site.E;
   }
 

src/tallies/filter.cpp

@@ -9,13 +9,16 @@
 #include "openmc/tallies/filter_cell.h"
 #include "openmc/tallies/filter_cellborn.h"
 #include "openmc/tallies/filter_cellfrom.h"
+#include "openmc/tallies/filter_delayedgroup.h"
 #include "openmc/tallies/filter_distribcell.h"
 #include "openmc/tallies/filter_energyfunc.h"
+#include "openmc/tallies/filter_energy.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_particle.h"
 #include "openmc/tallies/filter_polar.h"
 #include "openmc/tallies/filter_sph_harm.h"
 #include "openmc/tallies/filter_sptl_legendre.h"
@@ -31,7 +34,7 @@ namespace openmc {
 //==============================================================================
 
 std::vector<FilterMatch> filter_matches;
-std::vector<Filter*> tally_filters;
+std::vector<std::unique_ptr<Filter>> tally_filters;
 
 //==============================================================================
 // Non-member functions
@@ -45,7 +48,6 @@ free_memory_tally_c()
     filter_matches.clear();
   }
 
-  for (Filter* filt : tally_filters) {delete filt;}
   tally_filters.clear();
 }
 
@@ -94,45 +96,53 @@ extern "C" {
   {
     std::string type_ {type};
     if (type_ == "azimuthal") {
-      tally_filters.push_back(new AzimuthalFilter());
+      tally_filters.push_back(std::make_unique<AzimuthalFilter>());
     } else if (type_ == "cell") {
-      tally_filters.push_back(new CellFilter());
+      tally_filters.push_back(std::make_unique<CellFilter>());
     } else if (type_ == "cellborn") {
-      tally_filters.push_back(new CellbornFilter());
+      tally_filters.push_back(std::make_unique<CellbornFilter>());
     } else if (type_ == "cellfrom") {
-      tally_filters.push_back(new CellFromFilter());
+      tally_filters.push_back(std::make_unique<CellFromFilter>());
     } else if (type_ == "distribcell") {
-      tally_filters.push_back(new DistribcellFilter());
+      tally_filters.push_back(std::make_unique<DistribcellFilter>());
+    } else if (type_ == "delayedgroup") {
+      tally_filters.push_back(std::make_unique<DelayedGroupFilter>());
     } else if (type_ == "energyfunction") {
-      tally_filters.push_back(new EnergyFunctionFilter());
+      tally_filters.push_back(std::make_unique<EnergyFunctionFilter>());
+    } else if (type_ == "energy") {
+      tally_filters.push_back(std::make_unique<EnergyFilter>());
+    } else if (type_ == "energyout") {
+      tally_filters.push_back(std::make_unique<EnergyoutFilter>());
     } else if (type_ == "legendre") {
-      tally_filters.push_back(new LegendreFilter());
+      tally_filters.push_back(std::make_unique<LegendreFilter>());
     } else if (type_ == "material") {
-      tally_filters.push_back(new MaterialFilter());
+      tally_filters.push_back(std::make_unique<MaterialFilter>());
     } else if (type_ == "mesh") {
-      tally_filters.push_back(new MeshFilter());
+      tally_filters.push_back(std::make_unique<MeshFilter>());
     } else if (type_ == "meshsurface") {
-      tally_filters.push_back(new MeshSurfaceFilter());
+      tally_filters.push_back(std::make_unique<MeshSurfaceFilter>());
     } else if (type_ == "mu") {
-      tally_filters.push_back(new MuFilter());
+      tally_filters.push_back(std::make_unique<MuFilter>());
+    } else if (type_ == "particle") {
+      tally_filters.push_back(std::make_unique<ParticleFilter>());
     } else if (type_ == "polar") {
-      tally_filters.push_back(new PolarFilter());
+      tally_filters.push_back(std::make_unique<PolarFilter>());
     } else if (type_ == "surface") {
-      tally_filters.push_back(new SurfaceFilter());
+      tally_filters.push_back(std::make_unique<SurfaceFilter>());
     } else if (type_ == "spatiallegendre") {
-      tally_filters.push_back(new SpatialLegendreFilter());
+      tally_filters.push_back(std::make_unique<SpatialLegendreFilter>());
     } else if (type_ == "sphericalharmonics") {
-      tally_filters.push_back(new SphericalHarmonicsFilter());
+      tally_filters.push_back(std::make_unique<SphericalHarmonicsFilter>());
     } else if (type_ == "universe") {
-      tally_filters.push_back(new UniverseFilter());
+      tally_filters.push_back(std::make_unique<UniverseFilter>());
     } else if (type_ == "zernike") {
-      tally_filters.push_back(new ZernikeFilter());
+      tally_filters.push_back(std::make_unique<ZernikeFilter>());
     } else if (type_ == "zernikeradial") {
-      tally_filters.push_back(new ZernikeRadialFilter());
+      tally_filters.push_back(std::make_unique<ZernikeRadialFilter>());
     } else {
       return nullptr;
     }
-    return tally_filters.back();
+    return tally_filters.back().get();
   }
 
   void filter_from_xml(Filter* filt, pugi::xml_node* node)

src/tallies/filter_azimuthal.cpp

@@ -66,6 +66,7 @@ std::string
 AzimuthalFilter::text_label(int bin) const
 {
   std::stringstream out;
+  //TODO: off-by-one
   out << "Azimuthal Angle [" << bins_[bin-1] << ", " << bins_[bin] << ")";
   return out.str();
 }

src/tallies/filter_cell.cpp

@@ -63,6 +63,7 @@ CellFilter::to_statepoint(hid_t filter_group) const
 std::string
 CellFilter::text_label(int bin) const
 {
+  //TODO: off-by-one
   return "Cell " + std::to_string(cells[cells_[bin-1]]->id_);
 }
 

src/tallies/filter_cellborn.cpp

@@ -19,6 +19,7 @@ CellbornFilter::get_all_bins(const Particle* p, int estimator,
 std::string
 CellbornFilter::text_label(int bin) const
 {
+  //TODO: off-by-one
   return "Birth Cell " + std::to_string(cells[cells_[bin-1]]->id_);
 }
 

src/tallies/filter_cellfrom.cpp

@@ -21,6 +21,7 @@ CellFromFilter::get_all_bins(const Particle* p, int estimator,
 std::string
 CellFromFilter::text_label(int bin) const
 {
+  //TODO: off-by-one
   return "Cell from " + std::to_string(cells[cells_[bin-1]]->id_);
 }
 

src/tallies/filter_delayedgroup.cpp

@@ -0,0 +1,58 @@
+#include "openmc/tallies/filter_delayedgroup.h"
+
+#include "openmc/error.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+void
+DelayedGroupFilter::from_xml(pugi::xml_node node)
+{
+  groups_ = get_node_array<int>(node, "bins");
+  n_bins_ = groups_.size();
+
+  // Make sure all the group index values are valid.
+  // TODO: do these need to be decremented for zero-based indexing?
+  for (auto group : groups_) {
+    if (group < 1) {
+      fatal_error("Encountered delayedgroup bin with index "
+        + std::to_string(group) + " which is less than 1");
+    } else if (group > MAX_DELAYED_GROUPS) {
+      fatal_error("Encountered delayedgroup bin with index "
+        + std::to_string(group) + " which is greater than MAX_DELATED_GROUPS ("
+        + std::to_string(MAX_DELAYED_GROUPS) + ")");
+    }
+  }
+}
+
+void
+DelayedGroupFilter::get_all_bins(const Particle* p, int estimator,
+                                 FilterMatch& match) const
+{
+  //TODO: off-by-one
+  match.bins_.push_back(1);
+  match.weights_.push_back(1.0);
+}
+
+void
+DelayedGroupFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  write_dataset(filter_group, "bins", groups_);
+}
+
+std::string
+DelayedGroupFilter::text_label(int bin) const
+{
+  //TODO: off-by-one
+  return "Delayed Group " + std::to_string(groups_[bin-1]);
+}
+
+//==============================================================================
+// Fortran interoperability
+//==============================================================================
+
+extern "C" int delayedgroup_filter_groups(DelayedGroupFilter* filt, int i)
+{return filt->groups_[i-1];}
+
+} // namespace openmc

src/tallies/filter_distribcell.cpp

@@ -73,6 +73,7 @@ std::string
 DistribcellFilter::text_label(int bin) const
 {
   auto map = cells[cell_]->distribcell_index_;
+  //TODO: off-by-one
   auto path = distribcell_path(cell_, map, bin-1);
   return "Distributed Cell " + path;
 }

src/tallies/filter_energy.cpp

@@ -0,0 +1,185 @@
+#include "openmc/tallies/filter_energy.h"
+
+#include "openmc/capi.h"
+#include "openmc/constants.h"  // For F90_NONE
+#include "openmc/mgxs_interface.h"
+#include "openmc/search.h"
+#include "openmc/settings.h"
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+//==============================================================================
+// EnergyFilter implementation
+//==============================================================================
+
+void
+EnergyFilter::from_xml(pugi::xml_node node)
+{
+  bins_ = get_node_array<double>(node, "bins");
+  n_bins_ = bins_.size() - 1;
+
+  // In MG mode, check if the filter bins match the transport bins.
+  // We can save tallying time if we know that the tally bins match the energy
+  // group structure.  In that case, the matching bin index is simply the group
+  // (after flipping for the different ordering of the library and tallying
+  // systems).
+  if (!settings::run_CE) {
+    if (n_bins_ == num_energy_groups) {
+      matches_transport_groups_ = true;
+      for (auto i = 0; i < n_bins_ + 1; i++) {
+        if (rev_energy_bins[i] != bins_[i]) {
+          matches_transport_groups_ = false;
+          break;
+        }
+      }
+    }
+  }
+}
+
+void
+EnergyFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+const
+{
+  if (p->g != F90_NONE && matches_transport_groups_) {
+    if (estimator == ESTIMATOR_TRACKLENGTH) {
+      //TODO: off-by-one
+      match.bins_.push_back(num_energy_groups - p->g + 1);
+    } else {
+      //TODO: off-by-one
+      match.bins_.push_back(num_energy_groups - p->last_g + 1);
+    }
+    match.weights_.push_back(1.0);
+
+  } else {
+    // Get the pre-collision energy of the particle.
+    auto E = p->last_E;
+
+    // Bin the energy.
+    if (E >= bins_.front() && E <= bins_.back()) {
+      //TODO: off-by-one
+      auto bin = lower_bound_index(bins_.begin(), bins_.end(), E) + 1;
+      match.bins_.push_back(bin);
+      match.weights_.push_back(1.0);
+    }
+  }
+}
+
+void
+EnergyFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  write_dataset(filter_group, "bins", bins_);
+}
+
+std::string
+EnergyFilter::text_label(int bin) const
+{
+  std::stringstream out;
+  //TODO: off-by-one
+  out << "Incoming Energy [" << bins_[bin-1] << ", " << bins_[bin] << ")";
+  return out.str();
+}
+
+//==============================================================================
+// EnergyoutFilter implementation
+//==============================================================================
+
+void
+EnergyoutFilter::get_all_bins(const Particle* p, int estimator,
+                              FilterMatch& match) const
+{
+  if (p->g != F90_NONE && matches_transport_groups_) {
+    match.bins_.push_back(num_energy_groups - p->g + 1);
+    match.weights_.push_back(1.0);
+
+  } else {
+    if (p->E >= bins_.front() && p->E <= bins_.back()) {
+      //TODO: off-by-one
+      auto bin = lower_bound_index(bins_.begin(), bins_.end(), p->E) + 1;
+      match.bins_.push_back(bin);
+      match.weights_.push_back(1.0);
+    }
+  }
+}
+
+std::string
+EnergyoutFilter::text_label(int bin) const
+{
+  std::stringstream out;
+  //TODO: off-by-one
+  out << "Outgoing Energy [" << bins_[bin-1] << ", " << bins_[bin] << ")";
+  return out.str();
+}
+
+//==============================================================================
+// Fortran interoperability
+//==============================================================================
+
+extern "C" bool energy_filter_matches_transport_groups(EnergyFilter* filt)
+{return filt->matches_transport_groups_;}
+
+extern "C" int
+energy_filter_search(EnergyFilter* filt, double val)
+{
+  if (val < filt->bins_.front() || val > filt->bins_.back()) {
+    return -1;
+  } else {
+    //TODO: off-by-one
+    return lower_bound_index(filt->bins_.begin(), filt->bins_.end(), val) + 1;
+  }
+}
+
+//==============================================================================
+// C-API functions
+//==============================================================================
+
+extern"C" int
+openmc_energy_filter_get_bins(int32_t index, double** energies, 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<EnergyFilter*>(filt_base);
+
+  // Check the filter type.
+  if (!filt) {
+    set_errmsg("Tried to get energy bins on a non-energy filter.");
+    return OPENMC_E_INVALID_TYPE;
+  }
+
+  // Output the bins.
+  *energies = filt->bins_.data();
+  *n = filt->bins_.size();
+  return 0;
+}
+
+extern "C" int
+openmc_energy_filter_set_bins(int32_t index, int32_t n, const double* energies)
+{
+  // 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<EnergyFilter*>(filt_base);
+
+  // Check the filter type.
+  if (!filt) {
+    set_errmsg("Tried to set energy bins on a non-energy filter.");
+    return OPENMC_E_INVALID_TYPE;
+  }
+
+  // Update the filter.
+  filt->bins_.clear();
+  filt->bins_.resize(n);
+  for (int i = 0; i < n; i++) filt->bins_[i] = energies[i];
+  filt->n_bins_ = n - 1;
+  filter_update_n_bins(index);
+}
+
+}// namespace openmc

src/tallies/filter_legendre.cpp

@@ -37,6 +37,7 @@ LegendreFilter::to_statepoint(hid_t filter_group) const
 std::string
 LegendreFilter::text_label(int bin) const
 {
+  //TODO: off-by-one
   return "Legendre expansion, P" + std::to_string(bin - 1);
 }
 

src/tallies/filter_material.cpp

@@ -62,6 +62,7 @@ MaterialFilter::to_statepoint(hid_t filter_group) const
 std::string
 MaterialFilter::text_label(int bin) const
 {
+  //TODO: off-by-one
   return "Material " + std::to_string(materials[materials_[bin-1]]->id_);
 }
 

src/tallies/filter_meshsurface.cpp

@@ -22,6 +22,7 @@ MeshSurfaceFilter::text_label(int bin) const
   int n_dim = mesh.n_dimension_;
 
   // Get flattend mesh index and surface index.
+  //TODO: off-by-one
   int i_mesh = (bin - 1) / (4 * n_dim) + 1;
   int i_surf = ((bin - 1) % (4 * n_dim)) + 1;
 

src/tallies/filter_mu.cpp

@@ -57,6 +57,7 @@ std::string
 MuFilter::text_label(int bin) const
 {
   std::stringstream out;
+  //TODO: off-by-one
   out << "Change-in-Angle [" << bins_[bin-1] << ", " << bins_[bin] << ")";
   return out.str();
 }

src/tallies/filter_particle.cpp

@@ -0,0 +1,48 @@
+#include "openmc/tallies/filter_particle.h"
+
+#include "openmc/xml_interface.h"
+
+namespace openmc {
+
+void
+ParticleFilter::from_xml(pugi::xml_node node)
+{
+  particles_ = get_node_array<int>(node, "bins");
+  n_bins_ = particles_.size();
+}
+
+void
+ParticleFilter::get_all_bins(const Particle* p, int estimator,
+                             FilterMatch& match) const
+{
+  for (auto i = 0; i < particles_.size(); i++) {
+    if (particles_[i] == p->type) {
+      //TODO: off-by-one
+      match.bins_.push_back(i + 1);
+      match.weights_.push_back(1.0);
+    }
+  }
+}
+
+void
+ParticleFilter::to_statepoint(hid_t filter_group) const
+{
+  Filter::to_statepoint(filter_group);
+  write_dataset(filter_group, "bins", particles_);
+}
+
+std::string
+ParticleFilter::text_label(int bin) const
+{
+  //TODO: off-by-one
+  return "Particle " + std::to_string(particles_[bin-1]);
+}
+
+//==============================================================================
+// Fortran interoperability
+//==============================================================================
+
+extern "C" int particle_filter_particles(ParticleFilter* filt, int i)
+{return filt->particles_[i-1];}
+
+} // namespace openmc

src/tallies/filter_polar.cpp

@@ -65,6 +65,7 @@ std::string
 PolarFilter::text_label(int bin) const
 {
   std::stringstream out;
+  //TODO: off-by-one
   out << "Polar Angle [" << bins_[bin-1] << ", " << bins_[bin] << ")";
   return out.str();
 }

src/tallies/filter_sph_harm.cpp

@@ -78,6 +78,7 @@ SphericalHarmonicsFilter::text_label(int bin) const
   std::stringstream out;
   for (int n = 0; n < order_ + 1; n++) {
     if (bin <= (n + 1) * (n + 1)) {
+      //TODO: off-by-one
       int m = (bin - n*n - 1) - n;
       out << "Spherical harmonic expansion, Y" << n << "," << m;
       return out.str();

src/tallies/filter_sptl_legendre.cpp

@@ -88,6 +88,7 @@ SpatialLegendreFilter::text_label(int bin) const
   } else {
     out << "z";
   }
+  //TODO: off-by-one
   out << " axis, P" << std::to_string(bin - 1);
   return out.str();
 }

src/tallies/filter_surface.cpp

@@ -65,6 +65,7 @@ SurfaceFilter::to_statepoint(hid_t filter_group) const
 std::string
 SurfaceFilter::text_label(int bin) const
 {
+  //TODO: off-by-one
   return "Surface " + std::to_string(surfaces[surfaces_[bin-1]]->id_);
 }
 

src/tallies/filter_universe.cpp

@@ -63,6 +63,7 @@ UniverseFilter::to_statepoint(hid_t filter_group) const
 std::string
 UniverseFilter::text_label(int bin) const
 {
+  //TODO: off-by-one
   return "Universe " + std::to_string(universes[universes_[bin-1]]->id_);
 }
 

src/tallies/filter_zernike.cpp

@@ -62,6 +62,7 @@ ZernikeFilter::text_label(int bin) const
   std::stringstream out;
   for (int n = 0; n < order_+1; n++) {
     int last = (n + 1) * (n + 2) / 2;
+    //TODO: off-by-one
     if (bin <= last) {
       int first = last - n;
       int m = -n + (bin - first) * 2;
@@ -106,6 +107,7 @@ ZernikeRadialFilter::get_all_bins(const Particle* p, int estimator,
 std::string
 ZernikeRadialFilter::text_label(int bin) const
 {
+  //TODO: off-by-one
   return "Zernike expansion, Z" + std::to_string(2*(bin-1)) + ",0";
 }
 

src/tallies/tally.F90

@@ -2344,11 +2344,11 @@ contains
     integer :: d             ! delayed group
     integer :: g             ! another delayed group
     integer :: d_bin         ! delayed group bin index
-    integer :: n             ! number of energies on filter
     integer :: k             ! loop index for bank sites
     integer :: l             ! loop index for tally filters
     integer :: f             ! index in filters array
     integer :: b             ! index of filter bin
+    integer :: i_match       ! matching bin index on energyout filter
     integer :: i_bin         ! index of matching filter bin
     integer :: bin_energyout ! original outgoing energy bin
     integer :: i_filter      ! index for matching filter bin combination
@@ -2366,9 +2366,6 @@ contains
     select type(eo_filt => filters(i) % obj)
     type is (EnergyoutFilter)
 
-      ! Get number of energies on filter
-      n = size(eo_filt % bins)
-
       ! Since the creation of fission sites is weighted such that it is
       ! expected to create n_particles sites, we need to multiply the
       ! score by keff to get the true nu-fission rate. Otherwise, the sum
@@ -2397,7 +2394,7 @@ contains
 
           ! modify the value so that g_out = 1 corresponds to the highest
           ! energy bin
-          g_out = size(eo_filt % bins) - g_out
+          g_out = eo_filt % n_bins - g_out + 1
 
           ! change outgoing energy bin
           call filter_matches(i) % bins_set_data(i_bin, g_out)
@@ -2412,12 +2409,11 @@ contains
                  % E))
           end if
 
-          ! check if outgoing energy is within specified range on filter
-          if (E_out < eo_filt % bins(1) .or. E_out > eo_filt % bins(n)) cycle
-
-          ! change outgoing energy bin
-          call filter_matches(i) % bins_set_data(i_bin, &
-               binary_search(eo_filt % bins, n, E_out))
+          ! If this outgoing energy falls within the energyout filter's range,
+          ! set the appropriate filter_matches bin.
+          i_match = eo_filt % search(E_out)
+          if (i_match == -1) cycle
+          call filter_matches(i) % bins_set_data(i_bin, i_match)
 
         end if
 

src/tallies/tally_filter.F90

@@ -7,7 +7,6 @@ module tally_filter
   use tally_filter_header
 
   ! Inherit other filters
-  use tally_filter_cpp
   use tally_filter_delayedgroup
   use tally_filter_distribcell
   use tally_filter_energy
@@ -170,14 +169,11 @@ contains
           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
+        filters(index) % obj % ptr = allocate_filter(type)
+        if (.not. c_associated(filters(index) % obj % ptr)) then
+          err = E_UNASSIGNED
+          call set_errmsg("Could not allocate C++ tally filter")
+        end if
 
       end if
     else

src/tallies/tally_filter_cpp.F90

@@ -1,235 +0,0 @@
-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

@@ -2,13 +2,7 @@ module tally_filter_delayedgroup
 
   use, intrinsic :: ISO_C_BINDING
 
-  use constants,          only: ONE, MAX_LINE_LEN, MAX_DELAYED_GROUPS
-  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
@@ -20,67 +14,25 @@ module tally_filter_delayedgroup
 !===============================================================================
 
   type, public, extends(TallyFilter) :: DelayedGroupFilter
-    integer, allocatable :: groups(:)
   contains
-    procedure :: from_xml
-    procedure :: get_all_bins => get_all_bins_dg
-    procedure :: to_statepoint => to_statepoint_dg
-    procedure :: text_label => text_label_dg
+    procedure :: groups
   end type DelayedGroupFilter
 
 contains
 
-  subroutine from_xml(this, node)
-    class(DelayedGroupFilter), intent(inout) :: this
-    type(XMLNode), intent(in) :: node
-
-    integer :: i
-    integer :: n
-
-    n = node_word_count(node, "bins")
-
-    ! Allocate and store bins
-    this % n_bins = n
-    allocate(this % groups(n))
-    call get_node_array(node, "bins", this % groups)
-
-    ! Check that bins are all are between 1 and MAX_DELAYED_GROUPS
-    do i = 1, n
-      if (this % groups(i) < 1 .or. &
-           this % groups(i) > MAX_DELAYED_GROUPS) then
-        call fatal_error("Encountered delayedgroup bin with index " &
-             // trim(to_str(this % groups(i))) // " that is outside &
-             &the range of 1 to MAX_DELAYED_GROUPS ( " &
-             // trim(to_str(MAX_DELAYED_GROUPS)) // ")")
-      end if
-    end do
-  end subroutine from_xml
-
-  subroutine get_all_bins_dg(this, p, estimator, match)
-    class(DelayedGroupFilter), intent(in)  :: this
-    type(Particle),            intent(in)  :: p
-    integer,                   intent(in)  :: estimator
-    type(TallyFilterMatch),         intent(inout) :: match
-
-    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)
+  function groups(this, i) result(group)
     class(DelayedGroupFilter), intent(in) :: this
-    integer(HID_T),            intent(in) :: filter_group
-
-    call write_dataset(filter_group, "type", "delayedgroup")
-    call write_dataset(filter_group, "n_bins", this % n_bins)
-    call write_dataset(filter_group, "bins", this % groups)
-  end subroutine to_statepoint_dg
-
-  function text_label_dg(this, bin) result(label)
-    class(DelayedGroupFilter), intent(in) :: this
-    integer,                   intent(in) :: bin
-    character(MAX_LINE_LEN)               :: label
-
-    label = "Delayed Group " // to_str(this % groups(bin))
-  end function text_label_dg
+    integer,                   intent(in) :: i
+    integer                               :: group
+    interface
+      function delayedgroup_filter_groups(filt, i) result(group) bind(C)
+        import C_PTR, C_INT
+        type(C_PTR),    value :: filt
+        integer(C_INT), value :: i
+        integer(C_INT)        :: group
+      end function
+    end interface
+    group = delayedgroup_filter_groups(this % ptr, i)
+  end function groups
 
 end module tally_filter_delayedgroup

src/tallies/tally_filter_distribcell.F90

@@ -1,11 +1,11 @@
 module tally_filter_distribcell
 
-  use tally_filter_cpp
+  use tally_filter_header
 
   implicit none
   private
 
-  type, public, extends(CppTallyFilter) :: DistribcellFilter
+  type, public, extends(TallyFilter) :: DistribcellFilter
   contains
     procedure :: initialize => initialize_distribcell
   end type DistribcellFilter
@@ -14,7 +14,7 @@ contains
 
   subroutine initialize_distribcell(this)
     class(DistribcellFilter), intent(inout) :: this
-    call this % initialize_cpp_inner()
+    call this % initialize_cpp()
     this % n_bins = this % n_bins_cpp()
   end subroutine initialize_distribcell
 

src/tallies/tally_filter_energy.F90

@@ -2,14 +2,6 @@ module tally_filter_energy
 
   use, intrinsic :: ISO_C_BINDING
 
-  use algorithm,           only: binary_search
-  use constants
-  use error
-  use hdf5_interface
-  use mgxs_interface,      only: num_energy_groups, rev_energy_bins
-  use particle_header,     only: Particle
-  use settings,            only: run_CE
-  use string,              only: to_str
   use tally_filter_header
   use xml_interface
 
@@ -18,20 +10,33 @@ module tally_filter_energy
   public :: openmc_energy_filter_get_bins
   public :: openmc_energy_filter_set_bins
 
+  interface
+    function openmc_energy_filter_get_bins(index, energies, n) result(err) bind(C)
+      import C_INT32_T, C_PTR, C_INT
+      integer(C_INT32_T), value :: index
+      type(C_PTR), intent(out) :: energies
+      integer(C_INT32_T), intent(out) :: n
+      integer(C_INT) :: err
+    end function
+    function openmc_energy_filter_set_bins(index, n, energies) result(err) bind(C)
+      import C_INT32_T, C_DOUBLE, C_INT
+      integer(C_INT32_T), value, intent(in) :: index
+      integer(C_INT32_T), value, intent(in) :: n
+      real(C_DOUBLE), intent(in) :: energies(n)
+      integer(C_INT) :: err
+    end function
+  end interface
+
 !===============================================================================
 ! ENERGYFILTER bins the incident neutron energy.
 !===============================================================================
 
   type, public, extends(TallyFilter) :: EnergyFilter
-    real(8), allocatable :: bins(:)
-
     ! True if transport group number can be used directly to get bin number
     logical :: matches_transport_groups = .false.
   contains
     procedure :: from_xml => from_xml_energy
-    procedure :: get_all_bins => get_all_bins_energy
-    procedure :: to_statepoint => to_statepoint_energy
-    procedure :: text_label => text_label_energy
+    procedure :: search
   end type EnergyFilter
 
 !===============================================================================
@@ -41,11 +46,6 @@ module tally_filter_energy
 !===============================================================================
 
   type, public, extends(EnergyFilter) :: EnergyoutFilter
-  contains
-    ! Inherit from_xml from EnergyFilter
-    procedure :: get_all_bins => get_all_bins_energyout
-    procedure :: to_statepoint => to_statepoint_energyout
-    procedure :: text_label => text_label_energyout
   end type EnergyoutFilter
 
 contains
@@ -58,194 +58,34 @@ contains
     class(EnergyFilter), intent(inout) :: this
     type(XMLNode), intent(in) :: node
 
-    integer :: n
+    interface
+      function energy_filter_matches_transport_groups(filt) result(matches) &
+           bind(C)
+        import C_PTR, C_BOOL
+        type(C_PTR), value :: filt
+        logical(C_BOOL)    :: matches
+      end function
+    end interface
 
-    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)
-
-    ! We can save tallying time if we know that the tally bins match
-    ! the energy group structure.  In that case, the matching bin
-    ! index is simply the group (after flipping for the different
-    ! ordering of the library and tallying systems).
-    if (.not. run_CE) then
-      if (n == num_energy_groups + 1) then
-        if (all(this % bins == rev_energy_bins)) &
-             then
-          this % matches_transport_groups = .true.
-        end if
-      end if
-    end if
+    call this % from_xml_cpp(node)
+    this % n_bins = this % n_bins_cpp()
+    this % matches_transport_groups = &
+         energy_filter_matches_transport_groups(this % ptr)
   end subroutine from_xml_energy
 
-  subroutine get_all_bins_energy(this, p, estimator, match)
-    class(EnergyFilter), intent(in)  :: this
-    type(Particle),      intent(in)  :: p
-    integer,             intent(in)  :: estimator
-    type(TallyFilterMatch),   intent(inout) :: match
-
-    integer :: n
-    integer :: bin
-    real(8) :: E
-
-    n = this % n_bins
-
-    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)
-      else
-        call match % bins_push_back(num_energy_groups - p % last_g + 1)
-        call match % weights_push_back(ONE)
-      end if
-
-    else
-      ! Pre-collision energy of particle
-      E = p % last_E
-
-      ! 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)
-      end if
-    end if
-  end subroutine get_all_bins_energy
-
-  subroutine to_statepoint_energy(this, filter_group)
+  function search(this, val) result(bin)
     class(EnergyFilter), intent(in) :: this
-    integer(HID_T),      intent(in) :: filter_group
-
-    call write_dataset(filter_group, "type", "energy")
-    call write_dataset(filter_group, "n_bins", this % n_bins)
-    call write_dataset(filter_group, "bins", this % bins)
-  end subroutine to_statepoint_energy
-
-  function text_label_energy(this, bin) result(label)
-    class(EnergyFilter), 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 = "Incoming Energy [" // trim(to_str(E0)) // ", " &
-         // trim(to_str(E1)) // ")"
-  end function text_label_energy
-
-!===============================================================================
-! EnergyoutFilter methods
-!===============================================================================
-
-  subroutine get_all_bins_energyout(this, p, estimator, match)
-    class(EnergyoutFilter), 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
-
-    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)
-
-    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)
-      end if
-    end if
-  end subroutine get_all_bins_energyout
-
-  subroutine to_statepoint_energyout(this, filter_group)
-    class(EnergyoutFilter), intent(in) :: this
-    integer(HID_T),      intent(in) :: filter_group
-
-    call write_dataset(filter_group, "type", "energyout")
-    call write_dataset(filter_group, "n_bins", this % n_bins)
-    call write_dataset(filter_group, "bins", this % bins)
-  end subroutine to_statepoint_energyout
-
-  function text_label_energyout(this, bin) result(label)
-    class(EnergyoutFilter), 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 = "Outgoing Energy [" // trim(to_str(E0)) // ", " &
-         // trim(to_str(E1)) // ")"
-  end function text_label_energyout
-
-!===============================================================================
-!                               C API FUNCTIONS
-!===============================================================================
-
-  function openmc_energy_filter_get_bins(index, energies, n) result(err) bind(C)
-    ! Return the bounding energies for an energy filter
-    integer(C_INT32_T), value :: index
-    type(C_PTR), intent(out) :: energies
-    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 (EnergyFilter)
-        energies = C_LOC(f % bins)
-        n = size(f % bins)
-        err = 0
-      type is (EnergyoutFilter)
-        energies = C_LOC(f % bins)
-        n = size(f % bins)
-        err = 0
-        class default
-        err = E_INVALID_TYPE
-        call set_errmsg("Tried to get energy bins on a non-energy filter.")
-      end select
-    end if
-  end function openmc_energy_filter_get_bins
-
-
-  function openmc_energy_filter_set_bins(index, n, energies) result(err) bind(C)
-    ! Set the bounding energies for an energy filter
-    integer(C_INT32_T), value, intent(in) :: index
-    integer(C_INT32_T), value, intent(in) :: n
-    real(C_DOUBLE), intent(in) :: energies(n)
-    integer(C_INT) :: err
-
-    err = verify_filter(index)
-    if (err == 0) then
-      select type (f => filters(index) % obj)
-      type is (EnergyFilter)
-        f % n_bins = n - 1
-        if (allocated(f % bins)) deallocate(f % bins)
-        allocate(f % bins(n))
-        f % bins(:) = energies
-      type is (EnergyoutFilter)
-        f % n_bins = n - 1
-        if (allocated(f % bins)) deallocate(f % bins)
-        allocate(f % bins(n))
-        f % bins(:) = energies
-        class default
-        err = E_INVALID_TYPE
-        call set_errmsg("Tried to get energy bins on a non-energy filter.")
-      end select
-    end if
-  end function openmc_energy_filter_set_bins
+    real(8),             intent(in) :: val
+    integer                         :: bin
+    interface
+      function energy_filter_search(filt, val) result(bin) bind(C)
+        import C_PTR, C_DOUBLE, C_INT
+        type(C_PTR), value    :: filt
+        real(C_DOUBLE), value :: val
+        integer(C_INT)        :: bin
+      end function
+    end interface
+    bin = energy_filter_search(this % ptr, val)
+  end function search
 
 end module tally_filter_energy

src/tallies/tally_filter_header.F90

@@ -62,66 +62,62 @@ module tally_filter_header
   type, public, abstract :: TallyFilter
     integer :: id
     integer :: n_bins = 0
+    type(C_PTR) :: ptr
   contains
-    procedure(from_xml_),      deferred :: from_xml
-    procedure(get_all_bins_),  deferred :: get_all_bins
-    procedure(to_statepoint_), deferred :: to_statepoint
-    procedure(text_label_),    deferred :: text_label
-    procedure                           :: initialize => filter_initialize
+    procedure :: from_xml
+    procedure :: get_all_bins
+    procedure :: to_statepoint
+    procedure :: text_label
+    procedure :: initialize
+    procedure :: n_bins_cpp
+    procedure :: from_xml_cpp
+    procedure :: initialize_cpp
   end type TallyFilter
 
-  abstract interface
-
-    subroutine from_xml_(this, node)
-      import TallyFilter, XMLNode
-      class(TallyFilter), intent(inout) :: this
-      type(XMLNode), intent(in) :: node
-    end subroutine from_xml_
-
 !===============================================================================
-! GET_NEXT_BIN gives the index for the next valid filter bin and a weight that
-! will be applied to the flux.
-!
-! In principle, a filter can have multiple valid bins.  If current_bin =
-! NO_BIN_FOUND, then this method should give the first valid bin.  Providing the
-! first valid bin should then give the second valid bin, and so on.  When there
-! are no valid bins left, the next_bin should be NO_VALID_BIN.
-
-    subroutine get_all_bins_(this, p, estimator, match)
-      import TallyFilter
-      import Particle
-      import TallyFilterMatch
-      class(TallyFilter), intent(in)  :: this
-      type(Particle),     intent(in)  :: p
-      integer,            intent(in)  :: estimator
-      type(TallyFilterMatch), intent(inout) :: match
-    end subroutine get_all_bins_
-
+! Pure C++ filters
 !===============================================================================
-! TO_STATEPOINT writes all the information needed to reconstruct the filter to
-! the given filter_group.
 
-    subroutine to_statepoint_(this, filter_group)
-      import TallyFilter
-      import HID_T
-      class(TallyFilter), intent(in) :: this
-      integer(HID_T),     intent(in) :: filter_group
-    end subroutine to_statepoint_
+  type, public, extends(TallyFilter) :: AzimuthalFilter
+  end type
 
-!===============================================================================
-! TEXT_LABEL returns a string describing the given filter bin.  For example, an
-! energy filter might return the string "Incoming Energy [0.625E-6, 20.0)".
-! This is used to write the tallies.out file.
+  type, public, extends(TallyFilter) :: CellFilter
+  end type
 
-    function text_label_(this, bin) result(label)
-      import TallyFilter
-      import MAX_LINE_LEN
-      class(TallyFilter), intent(in) :: this
-      integer,            intent(in) :: bin
-      character(MAX_LINE_LEN)        :: label
-    end function text_label_
+  type, public, extends(TallyFilter) :: CellbornFilter
+  end type
 
-  end interface
+  type, public, extends(CellFilter) :: CellFromFilter
+  end type
+
+  type, public, extends(TallyFilter) :: EnergyFunctionFilter
+  end type
+
+  type, public, extends(TallyFilter) :: MaterialFilter
+  end type
+
+  type, public, extends(TallyFilter) :: MuFilter
+  end type
+
+  type, public, extends(TallyFilter) :: PolarFilter
+  end type
+
+  type, public, extends(TallyFilter) :: SpatialLegendreFilter
+  end type
+
+  type, public, extends(TallyFilter) :: SurfaceFilter
+    ! True if this filter is used for surface currents
+    logical :: current = .false.
+  end type
+
+  type, public, extends(TallyFilter) :: UniverseFilter
+  end type
+
+  type, public, extends(TallyFilter) :: ZernikeFilter
+  end type
+
+  type, public, extends(ZernikeFilter) :: ZernikeRadialFilter
+  end type
 
 !===============================================================================
 ! TALLYFILTERCONTAINER contains an allocatable TallyFilter object for arrays of
@@ -258,12 +254,128 @@ contains
   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.
+! TallyFilter implementation
+!===============================================================================
 
-  subroutine filter_initialize(this)
+  subroutine from_xml(this, node)
     class(TallyFilter), intent(inout) :: this
-  end subroutine filter_initialize
+    type(XMLNode),      intent(in)    :: node
+    call this % from_xml_cpp(node)
+    this % n_bins = this % n_bins_cpp()
+  end subroutine from_xml
+
+  subroutine get_all_bins(this, p, estimator, match)
+    class(TallyFilter),     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
+
+  subroutine to_statepoint(this, filter_group)
+    class(TallyFilter), 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
+
+  function text_label(this, bin) result(label)
+    class(TallyFilter), 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
+
+  subroutine initialize(this)
+    class(TallyFilter), intent(inout) :: this
+    call this % initialize_cpp()
+  end subroutine initialize
+
+  function n_bins_cpp(this) result(n_bins)
+    class(TallyFilter), intent(in) :: 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(this, node)
+    class(TallyFilter), 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
+
+  subroutine initialize_cpp(this)
+    class(TallyFilter), 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
+
+!===============================================================================
+! 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 = filters(index) % obj % ptr
+  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
+    filters(index) % obj % n_bins = filters(index) % obj % n_bins_cpp()
+  end subroutine
 
 !===============================================================================
 ! FREE_MEMORY_TALLY_FILTER deallocates global arrays defined in this module

src/tallies/tally_filter_legendre.F90

@@ -2,7 +2,7 @@ module tally_filter_legendre
 
   use, intrinsic :: ISO_C_BINDING
 
-  use tally_filter_cpp
+  use tally_filter_header
 
   implicit none
 
@@ -15,7 +15,7 @@ module tally_filter_legendre
     end function
   end interface
 
-  type, extends(CppTallyFilter) :: LegendreFilter
+  type, extends(TallyFilter) :: LegendreFilter
   end type LegendreFilter
 
 end module tally_filter_legendre

src/tallies/tally_filter_mesh.F90

@@ -2,7 +2,7 @@ module tally_filter_mesh
 
   use, intrinsic :: ISO_C_BINDING
 
-  use tally_filter_cpp
+  use tally_filter_header
 
   implicit none
 
@@ -15,7 +15,7 @@ module tally_filter_mesh
     end function
   end interface
 
-  type, public, extends(CppTallyFilter) :: MeshFilter
+  type, public, extends(TallyFilter) :: MeshFilter
   contains
     procedure :: mesh => get_mesh
   end type MeshFilter

src/tallies/tally_filter_meshsurface.F90

@@ -2,7 +2,7 @@ module tally_filter_meshsurface
 
   use, intrinsic :: ISO_C_BINDING
 
-  use tally_filter_cpp
+  use tally_filter_header
 
   implicit none
 
@@ -16,7 +16,7 @@ module tally_filter_meshsurface
     end function
   end interface
 
-  type, extends(CppTallyFilter) :: MeshSurfaceFilter
+  type, extends(TallyFilter) :: MeshSurfaceFilter
   end type MeshSurfaceFilter
 
 end module tally_filter_meshsurface

src/tallies/tally_filter_particle.F90

@@ -2,12 +2,7 @@ module tally_filter_particle
 
   use, intrinsic :: ISO_C_BINDING
 
-  use constants,          only: ONE, MAX_LINE_LEN
-  use hdf5_interface
-  use particle_header,    only: Particle
-  use string,             only: to_str
   use tally_filter_header
-  use xml_interface
 
   implicit none
   private
@@ -17,62 +12,25 @@ module tally_filter_particle
 !===============================================================================
 
   type, public, extends(TallyFilter) :: ParticleFilter
-    integer, allocatable :: particles(:)
   contains
-    procedure :: from_xml
-    procedure :: get_all_bins
-    procedure :: to_statepoint
-    procedure :: text_label
+    procedure :: particles
   end type ParticleFilter
 
 contains
 
-  subroutine from_xml(this, node)
-    class(ParticleFilter), 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 % particles(n))
-    call get_node_array(node, "bins", this % particles)
-  end subroutine from_xml
-
-  subroutine get_all_bins(this, p, estimator, match)
-    class(ParticleFilter), 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 (this % particles(i) == p % type) then
-        call match % bins_push_back(i)
-        call match % weights_push_back(ONE)
-      end if
-    end do
-  end subroutine get_all_bins
-
-  subroutine to_statepoint(this, filter_group)
+  function particles(this, i) result(ptype)
     class(ParticleFilter), intent(in) :: this
-    integer(HID_T),    intent(in) :: filter_group
-
-    call write_dataset(filter_group, "type", "particle")
-    call write_dataset(filter_group, "n_bins", this % n_bins)
-    call write_dataset(filter_group, "bins", this % particles)
-  end subroutine to_statepoint
-
-  function text_label(this, bin) result(label)
-    class(ParticleFilter), intent(in) :: this
-    integer,           intent(in) :: bin
-    character(MAX_LINE_LEN)       :: label
-
-    label = "Particle " // to_str(this % particles(bin))
-  end function text_label
+    integer,               intent(in) :: i
+    integer                           :: ptype
+    interface
+      function particle_filter_particles(filt, i) result(ptype) bind(C)
+        import C_PTR, C_INT
+        type(C_PTR),    value :: filt
+        integer(C_INT), value :: i
+        integer(C_INT)        :: ptype
+      end function
+    end interface
+    ptype = particle_filter_particles(this % ptr, i)
+  end function particles
 
 end module tally_filter_particle

src/tallies/tally_filter_sph_harm.F90

@@ -2,7 +2,7 @@ module tally_filter_sph_harm
 
   use, intrinsic :: ISO_C_BINDING
 
-  use tally_filter_cpp
+  use tally_filter_header
 
   implicit none
   private
@@ -10,7 +10,7 @@ module tally_filter_sph_harm
   integer, public, parameter :: COSINE_SCATTER = 1
   integer, public, parameter :: COSINE_PARTICLE = 2
 
-  type, public, extends(CppTallyFilter) :: SphericalHarmonicsFilter
+  type, public, extends(TallyFilter) :: SphericalHarmonicsFilter
   contains
     procedure :: cosine
   end type SphericalHarmonicsFilter