Move MaterialFilter to C++

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]
 
   //! \brief Default temperature for cells containing this material.

include/openmc/tallies/tally_filter_material.h

@@ -0,0 +1,81 @@
+#ifndef OPENMC_TALLY_FILTER_MATERIAL_H
+#define OPENMC_TALLY_FILTER_MATERIAL_H
+
+#include <sstream>
+#include <unordered_map>
+#include <vector>
+
+#include "openmc/error.h"
+#include "openmc/material.h"
+#include "openmc/tallies/tally_filter.h"
+
+
+namespace openmc {
+
+class MaterialFilter : public TallyFilter
+{
+public:
+  virtual std::string type() const override {return "material";}
+
+  virtual ~MaterialFilter() override = default;
+
+  virtual void
+  from_xml(pugi::xml_node node) override
+  {
+    materials_ = get_node_array<int32_t>(node, "bins");
+    n_bins_ = materials_.size();
+  }
+
+  virtual void
+  initialize() override
+  {
+    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);
+      }
+    }
+
+    for (int i = 0; i < materials_.size(); i++) {
+      map_[materials_[i]] = i;
+    }
+  }
+
+  virtual void
+  get_all_bins(Particle* p, int estimator, TallyFilterMatch& match)
+  const override
+  {
+    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);
+    }
+  }
+
+  virtual void
+  to_statepoint(hid_t filter_group) const override
+  {
+    TallyFilter::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);
+  }
+
+  virtual std::string
+  text_label(int bin) const override
+  {
+    return "Material " + std::to_string(materials[materials_[bin-1]]->id_);
+  }
+
+  std::vector<int32_t> materials_;
+  std::unordered_map<int32_t, int> map_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLY_FILTER_MATERIAL_H

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);
       }

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/tallies/tally_filter.cpp

@@ -8,6 +8,7 @@
 #include "openmc/tallies/tally_filter_cellborn.h"
 #include "openmc/tallies/tally_filter_cellfrom.h"
 #include "openmc/tallies/tally_filter_distribcell.h"
+#include "openmc/tallies/tally_filter_material.h"
 #include "openmc/tallies/tally_filter_mesh.h"
 #include "openmc/tallies/tally_filter_meshsurface.h"
 #include "openmc/tallies/tally_filter_mu.h"
@@ -93,6 +94,8 @@ extern "C" {
       tally_filters.push_back(new CellFromFilter());
     } else if (type_ == "distribcell") {
       tally_filters.push_back(new DistribcellFilter());
+    } else if (type_ == "material") {
+      tally_filters.push_back(new MaterialFilter());
     } else if (type_ == "mesh") {
       tally_filters.push_back(new MeshFilter());
     } else if (type_ == "meshsurface") {
@@ -140,6 +143,24 @@ extern "C" {
     *n = filt->cells_.size();
   }
 
+  void
+  material_filter_get_bins(MaterialFilter* filt, int32_t** bins, int32_t* n)
+  {
+    *bins = filt->materials_.data();
+    *n = filt->materials_.size();
+  }
+
+  void
+  material_filter_set_bins(MaterialFilter* filt, int32_t n, int32_t* bins)
+  {
+    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_[bins[i]] = i;
+  }
+
   int mesh_filter_get_mesh(MeshFilter* filt) {return filt->mesh_;}
 
   void

src/tallies/tally_filter_material.F90

@@ -2,117 +2,20 @@ 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
+  type, extends(CppTallyFilter) :: MaterialFilter
   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
 !===============================================================================
@@ -124,12 +27,20 @@ contains
     integer(C_INT32_T), intent(out) :: n
     integer(C_INT) :: err
 
+    interface
+      subroutine material_filter_get_bins(filt, bins, n) bind(C)
+        import C_PTR, C_INT32_T
+        type(C_PTR), value :: filt
+        type(C_PTR)        :: bins
+        integer(C_INT32_T) :: n
+      end subroutine material_filter_get_bins
+    end interface
+
     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)
+        call material_filter_get_bins(f % ptr, bins, n)
         err = 0
         class default
         err = E_INVALID_TYPE
@@ -147,22 +58,21 @@ contains
     integer(C_INT32_T), intent(in) :: bins(n)
     integer(C_INT) :: err
 
-    integer :: i
+    interface
+      subroutine material_filter_set_bins(filt, n, bins) bind(C)
+        import C_PTR, C_INT32_T
+        type(C_PTR),        value :: filt
+        integer(C_INT32_T), value :: n
+        integer(C_INT32_T)        :: bins(n)
+      end subroutine material_filter_set_bins
+    end interface
 
     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
+        call material_filter_set_bins(f % ptr, n, bins)
+        f % n_bins = f % n_bins_cpp()
 
         class default
         err = E_INVALID_TYPE