Start converting WMP class to C++

CMakeLists.txt

@@ -447,9 +447,9 @@ add_library(libopenmc SHARED
   src/tallies/tally.cpp
   src/timer.cpp
   src/thermal.cpp
+  src/wmp.cpp
   src/xml_interface.cpp
-  src/xsdata.cpp
-  src/tallies/tally.cpp)
+  src/xsdata.cpp)
 
 set_target_properties(libopenmc PROPERTIES
   OUTPUT_NAME openmc

include/openmc/hdf5_interface.h

@@ -253,15 +253,18 @@ void read_dataset(hid_t dset, xt::xarray<T>& arr, bool indep=false)
   std::size_t size = 1;
   for (const auto x : shape)
     size *= x;
-  std::vector<T> buffer(size);
+  T* buffer = new T[size];
 
   // Read data from attribute
-  read_dataset(dset, nullptr, H5TypeMap<T>::type_id, buffer.data(), indep);
+  read_dataset(dset, nullptr, H5TypeMap<T>::type_id, buffer, indep);
 
   // Adapt into xarray
-  arr = xt::adapt(buffer, shape);
+  arr = xt::adapt(buffer, size, xt::acquire_ownership(), shape);
 }
 
+template<>
+void read_dataset(hid_t dset, xt::xarray<std::complex<double>>& arr, bool indep);
+
 template <typename T>
 void read_dataset(hid_t obj_id, const char* name, xt::xarray<T>& arr, bool indep=false)
 {

include/openmc/wmp.h

@@ -0,0 +1,35 @@
+#ifndef OPENMC_WMP_H
+#define OPENMC_WMP_H
+
+#include "hdf5.h"
+#include "xtensor/xtensor.hpp"
+
+#include <string>
+
+namespace openmc {
+
+class WindowedMultipole {
+public:
+  // Types, aliases
+  using cdouble = std::complex<double>;
+
+  // Constructors, destructors
+  WindowedMultipole(hid_t group);
+
+  // Data members
+  std::string name_; //!< Name of nuclide
+  bool fissionable_; //!< Is the nuclide fissionable?
+  xt::xtensor<cdouble, 2> data_; //!< Poles and residues
+  double sqrt_awr_; //!< Square root of atomic weight ratio
+  double E_min_; //!< Minimum energy in [eV]
+  double E_max_; //!< Maximum energy in [eV]
+  double spacing_; //!< Spacing in sqrt(E) space
+  int fit_order_; //!< Order of the fit
+  xt::xtensor<int, 2> windows_; //!< Indices of pole at start/end of window
+  xt::xtensor<double, 2> curvefit_; //!< Fitting function (reaction, coeff index, window index)
+  xt::xtensor<bool, 1> broaden_poly_; //!< Whether to broaden curvefit
+};
+
+} // namespace openmc
+
+#endif // OPENMC_WMP_H

src/hdf5_interface.cpp

@@ -487,6 +487,25 @@ read_dataset(hid_t obj_id, const char* name, hid_t mem_type_id,
   if (name) H5Dclose(dset);
 }
 
+template<>
+void read_dataset(hid_t dset, xt::xarray<std::complex<double>>& arr, bool indep)
+{
+  // Get shape of dataset
+  std::vector<hsize_t> shape = object_shape(dset);
+
+  // Allocate new array to read data into
+  std::size_t size = 1;
+  for (const auto x : shape)
+    size *= x;
+  std::vector<std::complex<double>> buffer(size);
+
+  // Read data from attribute
+  read_complex(dset, nullptr, buffer.data(), indep);
+
+  // Adapt into xarray
+  arr = xt::adapt(buffer, shape);
+}
+
 
 void
 read_double(hid_t obj_id, const char* name, double* buffer, bool indep)
@@ -757,6 +776,8 @@ using_mpio_device(hid_t obj_id)
 
 // Specializations of the H5TypeMap template struct
 template<>
+const hid_t H5TypeMap<bool>::type_id = H5T_NATIVE_INT8;
+template<>
 const hid_t H5TypeMap<int>::type_id = H5T_NATIVE_INT;
 template<>
 const hid_t H5TypeMap<unsigned long>::type_id = H5T_NATIVE_ULONG;

src/wmp.cpp

@@ -0,0 +1,47 @@
+#include "openmc/wmp.h"
+
+#include "openmc/hdf5_interface.h"
+
+namespace openmc {
+
+WindowedMultipole::WindowedMultipole(hid_t group)
+{
+  // Get name of nuclide from group, removing leading '/'
+  name_ = object_name(group).substr(1);
+
+  // Read scalar values.
+  read_dataset(group, "spacing", spacing_);
+  read_dataset(group, "sqrtAWR", sqrt_awr_);
+  read_dataset(group, "E_min", E_min_);
+  read_dataset(group, "E_max", E_max_);
+
+  // Read the "data" array.  Use its shape to figure out the number of poles
+  // and residue types in this data.
+  read_dataset(group, "data", data_);
+  int n_residues = data_.shape()[1] - 1;
+
+  // Check to see if this data includes fission residues.
+  fissionable_ = (n_residues == 3);
+
+  // Read the "windows" array and use its shape to figure out the number of
+  // windows.
+  read_dataset(group, "windows", windows_);
+  int n_windows = windows_.shape()[0];
+
+  // Read the "broaden_poly" arrays.
+  read_dataset(group, "broaden_poly", broaden_poly_);
+  if (n_windows != broaden_poly_.shape()[0]) {
+    fatal_error("broaden_poly array shape is not consistent with the windows "
+      "array shape in WMP library for " + name_ + ".");
+  }
+
+  // Read the "curvefit" array.
+  read_dataset(group, "curvefit", curvefit_);
+  if (n_windows != broaden_poly_.shape()[0]) {
+    fatal_error("curvefit array shape is not consistent with the windows "
+      "array shape in WMP library for " + name_ + ".");
+  }
+  fit_order_ = curvefit_.shape()[1] - 1;
+}
+
+} // namespace openmc