Merge branch 'openmc-dev:develop' into infix-sense-evaluation

include/openmc/constants.h

@@ -24,7 +24,7 @@ using double_4dvec = vector<vector<vector<vector<double>>>>;
 constexpr int HDF5_VERSION[] {3, 0};
 
 // Version numbers for binary files
-constexpr array<int, 2> VERSION_STATEPOINT {17, 0};
+constexpr array<int, 2> VERSION_STATEPOINT {18, 0};
 constexpr array<int, 2> VERSION_PARTICLE_RESTART {2, 0};
 constexpr array<int, 2> VERSION_TRACK {3, 0};
 constexpr array<int, 2> VERSION_SUMMARY {6, 0};

include/openmc/endf.h

@@ -126,6 +126,26 @@ private:
     debye_waller_; //!< Debye-Waller integral divided by atomic mass in [eV^-1]
 };
 
+//==============================================================================
+//! Sum of multiple 1D functions
+//==============================================================================
+
+class Sum1D : public Function1D {
+public:
+  // Constructors
+  explicit Sum1D(hid_t group);
+
+  //! Evaluate each function and sum results
+  //! \param[in] x independent variable
+  //! \return Function evaluated at x
+  double operator()(double E) const override;
+
+  const unique_ptr<Function1D>& functions(int i) const { return functions_[i]; }
+
+private:
+  vector<unique_ptr<Function1D>> functions_; //!< individual functions
+};
+
 //! Read 1D function from HDF5 dataset
 //! \param[in] group HDF5 group containing dataset
 //! \param[in] name Name of dataset

include/openmc/hdf5_interface.h

@@ -61,6 +61,8 @@ void ensure_exists(hid_t obj_id, const char* name, bool attribute = false);
 vector<std::string> group_names(hid_t group_id);
 vector<hsize_t> object_shape(hid_t obj_id);
 std::string object_name(hid_t obj_id);
+hid_t open_object(hid_t group_id, const std::string& name);
+void close_object(hid_t obj_id);
 
 //==============================================================================
 // Fortran compatibility functions

include/openmc/mesh.h

@@ -7,13 +7,14 @@
 #include <unordered_map>
 
 #include "hdf5.h"
-#include "pugixml.hpp"
 #include "xtensor/xtensor.hpp"
+#include "pugixml.hpp"
 
 #include "openmc/memory.h" // for unique_ptr
 #include "openmc/particle.h"
 #include "openmc/position.h"
 #include "openmc/vector.h"
+#include "openmc/xml_interface.h"
 
 #ifdef DAGMC
 #include "moab/AdaptiveKDTree.hpp"
@@ -36,6 +37,12 @@
 
 namespace openmc {
 
+//==============================================================================
+// Constants
+//==============================================================================
+
+enum class ElementType { UNSUPPORTED=-1, LINEAR_TET, LINEAR_HEX };
+
 //==============================================================================
 // Global variables
 //==============================================================================
@@ -53,7 +60,7 @@ extern vector<unique_ptr<Mesh>> meshes;
 
 #ifdef LIBMESH
 namespace settings {
-// used when creating new libMesh::Mesh instances
+// used when creating new libMesh::MeshBase instances
 extern unique_ptr<libMesh::LibMeshInit> libmesh_init;
 extern const libMesh::Parallel::Communicator* libmesh_comm;
 } // namespace settings
@@ -485,6 +492,23 @@ public:
   //! \return The centroid of the bin
   virtual Position centroid(int bin) const = 0;
 
+  //! Get the number of vertices in the mesh
+  //
+  //! \return Number of vertices
+  virtual int n_vertices() const = 0;
+
+  //! Retrieve a vertex of the mesh
+  //
+  //! \param[in] vertex ID
+  //! \return vertex coordinates
+  virtual Position vertex(int id) const = 0;
+
+  //! Retrieve connectivity of a mesh element
+  //
+  //! \param[in] element ID
+  //! \return element connectivity as IDs of the vertices
+  virtual std::vector<int> connectivity(int id) const = 0;
+
   //! Get the volume of a mesh bin
   //
   //! \param[in] bin Bin to return the volume for
@@ -557,6 +581,12 @@ public:
 
   Position centroid(int bin) const override;
 
+  int n_vertices() const override;
+
+  Position vertex(int id) const override;
+
+  std::vector<int> connectivity(int id) const override;
+
   double volume(int bin) const override;
 
 private:
@@ -621,6 +651,9 @@ private:
   //! \return MOAB EntityHandle of tet
   moab::EntityHandle get_ent_handle_from_bin(int bin) const;
 
+  //! Get a vertex index into the global range from a handle
+  int get_vert_idx_from_handle(moab::EntityHandle vert) const;
+
   //! Get the bin for a given mesh cell index
   //
   //! \param[in] idx Index of the mesh cell.
@@ -655,6 +688,7 @@ private:
 
   // Data members
   moab::Range ehs_; //!< Range of tetrahedra EntityHandle's in the mesh
+  moab::Range verts_; //!< Range of vertex EntityHandle's in the mesh
   moab::EntityHandle tetset_;      //!< EntitySet containing all tetrahedra
   moab::EntityHandle kdtree_root_; //!< Root of the MOAB KDTree
   std::shared_ptr<moab::Interface> mbi_;    //!< MOAB instance
@@ -671,7 +705,8 @@ class LibMesh : public UnstructuredMesh {
 public:
   // Constructors
   LibMesh(pugi::xml_node node);
-  LibMesh(const std::string& filename, double length_multiplier = 1.0);
+  LibMesh(const std::string & filename, double length_multiplier = 1.0);
+  LibMesh(libMesh::MeshBase & input_mesh, double length_multiplier = 1.0);
 
   static const std::string mesh_lib_type;
 
@@ -701,10 +736,19 @@ public:
 
   Position centroid(int bin) const override;
 
+  int n_vertices() const override;
+
+  Position vertex(int id) const override;
+
+  std::vector<int> connectivity(int id) const override;
+
   double volume(int bin) const override;
 
+  libMesh::MeshBase* mesh_ptr() const { return m_; };
+
 private:
   void initialize() override;
+  void set_mesh_pointer_from_filename(const std::string& filename);
 
   // Methods
 
@@ -715,7 +759,8 @@ private:
   int get_bin_from_element(const libMesh::Elem* elem) const;
 
   // Data members
-  unique_ptr<libMesh::Mesh> m_; //!< pointer to the libMesh mesh instance
+  unique_ptr<libMesh::MeshBase> unique_m_ = nullptr; //!< pointer to the libMesh MeshBase instance, only used if mesh is created inside OpenMC
+  libMesh::MeshBase* m_; //!< pointer to libMesh MeshBase instance, always set during intialization
   vector<unique_ptr<libMesh::PointLocatorBase>>
     pl_; //!< per-thread point locators
   unique_ptr<libMesh::EquationSystems>

include/openmc/secondary_thermal.h

@@ -150,6 +150,34 @@ private:
                                        //!< each incident energy
 };
 
+//==============================================================================
+//! Mixed coherent/incoherent elastic angle-energy distribution
+//==============================================================================
+
+class MixedElasticAE : public AngleEnergy {
+public:
+  //! Construct from HDF5 file
+  //
+  //! \param[in] group  HDF5 group
+  explicit MixedElasticAE(
+    hid_t group, const CoherentElasticXS& coh_xs, const Function1D& incoh_xs);
+
+  //! Sample distribution for an angle and energy
+  //! \param[in] E_in Incoming energy in [eV]
+  //! \param[out] E_out Outgoing energy in [eV]
+  //! \param[out] mu Outgoing cosine with respect to current direction
+  //! \param[inout] seed Pseudorandom number seed pointer
+  void sample(
+    double E_in, double& E_out, double& mu, uint64_t* seed) const override;
+
+private:
+  CoherentElasticAE coherent_dist_;         //!< Coherent distribution
+  unique_ptr<AngleEnergy> incoherent_dist_; //!< Incoherent distribution
+
+  const CoherentElasticXS& coherent_xs_; //!< Ref. to coherent XS
+  const Function1D& incoherent_xs_;      //!< Polymorphic ref. to incoherent XS
+};
+
 } // namespace openmc
 
 #endif // OPENMC_SECONDARY_THERMAL_H