diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 2591f7d2ce..20831fdc81 100644 --- a/include/openmc/mesh.h +++ b/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 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& bins) const; + + std::pair, std::vector> plot(Position plot_ll, + Position plot_ur) const; + + void to_hdf5(hid_t group) const; + private: std::unique_ptr m_; std::unique_ptr point_locator_; diff --git a/include/openmc/settings.h b/include/openmc/settings.h index 783c221266..ce2e7b10a3 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -6,6 +6,7 @@ #include #include +#include #include #include #include @@ -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 LMI; +#endif + // Boolean flags extern bool assume_separate; //!< assume tallies are spatially separate? extern bool check_overlaps; //!< check overlaps in geometry? diff --git a/src/initialize.cpp b/src/initialize.cpp index 2c645ace1b..edc1df5001 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -3,6 +3,7 @@ #include #include // for getenv #include +#include #include #include @@ -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(new libMesh::LibMeshInit(argc, argv)); +#endif + + // Parse command-line arguments int err = parse_command_line(argc, argv); if (err) return err; diff --git a/src/mesh.cpp b/src/mesh.cpp index 1c15b1e10c..3446308fe6 100644 --- a/src/mesh.cpp +++ b/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(new libMesh::Mesh(init.comm(), 3)); + // always 3 for unstructured meshes + n_dimension_ = 3; + + m_ = std::unique_ptr(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& bins) const +{} + +std::pair, std::vector> +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& 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 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(node)); } else if (mesh_type == "rectilinear") { model::meshes.push_back(std::make_unique(node)); #ifdef DAGMC - } else if (mesh_type == "unstructured") { + } + else if (mesh_type == "unstructured" && mesh_lib == "moab") { model::meshes.push_back(std::make_unique(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(node)); #endif } else { fatal_error("Invalid mesh type: " + mesh_type); diff --git a/src/settings.cpp b/src/settings.cpp index 41b18942c3..3326f0029a 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -37,6 +37,10 @@ namespace openmc { namespace settings { +#ifdef LIBMESH +std::unique_ptr LMI; +#endif + // Default values for boolean flags bool assume_separate {false}; bool check_overlaps {false};