Collision-based is working. Initialization is a nightmare right now.

include/openmc/mesh.h

@@ -462,6 +462,10 @@ public:
 
   int get_bin(Position r) const;
 
+  int n_bins() const;
+
+  void get_indices(Position r, int* ijk, bool* in_mesh) const;
+
   std::pair<double, libMesh::Elem*>
   locate_boundary_element(const Position& r0,
                           const Position& r1) const;
@@ -485,8 +489,23 @@ public:
   double first(const libMesh::Node& a,
                const libMesh::Node& b) const;
 
+  int get_bin_from_indices(const int* ijk) const;
+
+  void get_indices_from_bin(int bin, int* ijk) const;
+
+
   bool intersects(Position& r0, Position r1, int* ijk) const;
 
+  int n_surface_bins() const;
+
+  void surface_bins_crossed(const Particle* p,
+                             std::vector<int>& bins) const;
+
+  std::pair<std::vector<double>, std::vector<double>> plot(Position plot_ll,
+                                                           Position plot_ur) const;
+
+  void to_hdf5(hid_t group) const;
+
 private:
   std::unique_ptr<libMesh::Mesh> m_;
   std::unique_ptr<libMesh::PointLocatorBase> point_locator_;

include/openmc/settings.h

@@ -6,6 +6,7 @@
 
 #include <array>
 #include <cstdint>
+#include <memory>
 #include <string>
 #include <unordered_set>
 #include <vector>
@@ -14,6 +15,10 @@
 
 #include "openmc/constants.h"
 
+#ifdef LIBMESH
+#include "libmesh/libmesh.h"
+#endif
+
 namespace openmc {
 
 //==============================================================================
@@ -22,6 +27,10 @@ namespace openmc {
 
 namespace settings {
 
+#ifdef LIBMESH
+extern std::unique_ptr<libMesh::LibMeshInit> LMI;
+#endif
+
 // Boolean flags
 extern bool assume_separate;          //!< assume tallies are spatially separate?
 extern bool check_overlaps;           //!< check overlaps in geometry?

src/initialize.cpp

@@ -3,6 +3,7 @@
 #include <cstddef>
 #include <cstdlib> // for getenv
 #include <cstring>
+#include <memory>
 #include <string>
 #include <vector>
 
@@ -32,6 +33,10 @@
 #include "openmc/thermal.h"
 #include "openmc/timer.h"
 
+#ifdef LIBMESH
+#include "libmesh/libmesh.h"
+#endif
+
 
 int openmc_init(int argc, char* argv[], const void* intracomm)
 {
@@ -50,6 +55,11 @@ int openmc_init(int argc, char* argv[], const void* intracomm)
   initialize_mpi(comm);
 #endif
 
+#ifdef LIBMESH
+  settings::LMI = std::unique_ptr<libMesh::LibMeshInit>(new libMesh::LibMeshInit(argc, argv));
+#endif
+
+
   // Parse command-line arguments
   int err = parse_command_line(argc, argv);
   if (err) return err;

src/mesh.cpp

@@ -2048,12 +2048,16 @@ UnstructuredMeshBase::UnstructuredMeshBase(pugi::xml_node node) : Mesh(node) {
 
 #ifdef LIBMESH
 LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMeshBase(node) {
-  libMesh::LibMeshInit init(0, NULL);
 
-  m_ = std::unique_ptr<libMesh::Mesh>(new libMesh::Mesh(init.comm(), 3));
+  // always 3 for unstructured meshes
+  n_dimension_ = 3;
+
+  m_ = std::unique_ptr<libMesh::Mesh>(new libMesh::Mesh(settings::LMI->comm(), 3));
   m_->read(filename_);
 
   point_locator_ = m_->sub_point_locator();
+  point_locator_->enable_out_of_mesh_mode();
+
   m_->find_neighbors();
 
   auto e = *m_->elements_begin();
@@ -2071,6 +2075,49 @@ LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMeshBase(node) {
 
 }
 
+void
+LibMesh::get_indices(Position r, int* ijk, bool* in_mesh) const {
+  int bin = get_bin(r);
+  *in_mesh = bin != -1;
+  ijk[0] = bin;
+  return;
+}
+
+int LibMesh::n_bins() const {
+  return m_->n_elem();
+}
+
+int LibMesh::n_surface_bins() const {
+  return 0;
+}
+
+void
+LibMesh::surface_bins_crossed(const Particle* p,
+                               std::vector<int>& bins) const
+{}
+
+std::pair<std::vector<double>, std::vector<double>>
+LibMesh::plot(Position plot_ll,
+              Position plot_ur) const { return {}; }
+
+
+int
+LibMesh::get_bin_from_indices(const int* ijk) const {
+  if (ijk[0] >= n_bins()) {
+    std::stringstream s;
+    s << "Invalid bin: " << ijk[0];
+    fatal_error(s);
+  }
+  int bin = first_element_->id() + ijk[0];
+  return bin;
+}
+
+void
+LibMesh::get_indices_from_bin(int bin, int* ijk) const {
+  ijk[0] = bin;
+}
+
+
 void
 LibMesh::bins_crossed(const Particle* p,
                       std::vector<int>& bins,
@@ -2305,6 +2352,22 @@ LibMesh::first(const libMesh::Node& a,
   }
 }
 
+void LibMesh::to_hdf5(hid_t group) const
+{
+  hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_));
+
+  write_dataset(mesh_group, "type", "unstructured");
+  write_dataset(mesh_group, "filename", filename_);
+
+  // write volume of each tet
+  std::vector<double> tet_vols;
+  for (int i = 0; i < m_->n_elem(); i++) {
+    tet_vols.emplace_back(m_->elem_ref(i).volume());
+  }
+  write_dataset(mesh_group, "volumes", tet_vols);
+
+  close_group(mesh_group);
+}
 
 double
 LibMesh::plucker_edge_test(const libMesh::Node& vertexa,
@@ -2401,17 +2464,31 @@ void read_meshes(pugi::xml_node root)
       mesh_type = "regular";
     }
 
+    std::string mesh_lib;
+    if (check_for_node(node, "library")) {
+      mesh_lib = get_node_value(node, "library", true, true);
+    } else {
+      mesh_lib = "moab";
+    }
+
     // Read mesh and add to vector
     if (mesh_type == "regular") {
       model::meshes.push_back(std::make_unique<RegularMesh>(node));
     } else if (mesh_type == "rectilinear") {
       model::meshes.push_back(std::make_unique<RectilinearMesh>(node));
 #ifdef DAGMC
-    } else if (mesh_type == "unstructured") {
+    }
+    else if (mesh_type == "unstructured" && mesh_lib == "moab") {
       model::meshes.push_back(std::make_unique<UnstructuredMesh>(node));
 #else
     } else if (mesh_type == "unstructured") {
       fatal_error("Unstructured mesh support is disabled.");
+#endif
+#ifdef LIBMESH
+    }
+    else if (mesh_type == "unstructured" && mesh_lib == "libmesh") {
+      std::cout << "Making libmesh mesh" << std::endl;
+      model::meshes.push_back(std::make_unique<LibMesh>(node));
 #endif
     } else {
       fatal_error("Invalid mesh type: " + mesh_type);

src/settings.cpp

@@ -37,6 +37,10 @@ namespace openmc {
 
 namespace settings {
 
+#ifdef LIBMESH
+std::unique_ptr<libMesh::LibMeshInit> LMI;
+#endif
+
 // Default values for boolean flags
 bool assume_separate         {false};
 bool check_overlaps          {false};