diff --git a/apps/impact_io.cpp b/apps/impact_io.cpp index 347c16d..01d098e 100644 --- a/apps/impact_io.cpp +++ b/apps/impact_io.cpp @@ -34,11 +34,11 @@ int main(int argc, char *argv[]) { TCLAP::ValueArg arg_nb_bands("b", "nbands", "Number of bands to consider", false, 12, "int"); TCLAP::ValueArg arg_nb_threads("j", "nthreads", "number of threads to use.", false, 1, "int"); TCLAP::ValueArg arg_radius_BZ("R", - "radiusBZ", - "Max norm of G vector for which the BZ center in G is taken into account", - false, - 0, - "double"); + "radiusBZ", + "Max norm of G vector for which the BZ center in G is taken into account", + false, + 0, + "double"); TCLAP::SwitchArg plot_with_python("P", "plot", "Call a python script after the computation to plot the band structure.", false); cmd.add(plot_with_python); cmd.add(arg_mesh_file); @@ -71,14 +71,13 @@ int main(int argc, char *argv[]) { EmpiricalPseudopotential::Material current_material = materials.materials.at(arg_material.getValue()); - bz_mesh::DielectricMesh my_dielectric_mesh(current_material); - my_dielectric_mesh.read_dielectric_file(arg_dielectric_file.getValue()); + bz_mesh::ImpactIonization my_impact_ionization(current_material, arg_mesh_file.getValue()); + my_impact_ionization.read_dielectric_file(arg_dielectric_file.getValue()); - bz_mesh::ImpactIonization my_impact_ionization(current_material, arg_mesh_file.getValue()); my_impact_ionization.set_max_radius_G0_BZ(arg_radius_BZ.getValue()); my_impact_ionization.compute_eigenstates(); - + auto end = std::chrono::high_resolution_clock::now(); std::chrono::duration elapsed_seconds = end - start; diff --git a/src/BZ_MESH/CMakeLists.txt b/src/BZ_MESH/CMakeLists.txt index d056c0c..91a4e81 100644 --- a/src/BZ_MESH/CMakeLists.txt +++ b/src/BZ_MESH/CMakeLists.txt @@ -8,6 +8,8 @@ set(HEADER_FILES_LIBBZMESH electron_phonon.hpp impact_ionization.hpp dielectric_mesh.hpp + bbox_mesh.hpp + octree_bz.hpp ) @@ -18,6 +20,7 @@ set(SOURCE_FILES_LIBBZMESH electron_phonon.cpp impact_ionization.cpp dielectric_mesh.cpp + octree_bz.cpp ) add_library(lib_bzmesh STATIC ${SOURCE_FILES_LIBBZMESH} ${HEADER_FILES_LIBBZMESH}) diff --git a/src/BZ_MESH/bbox_mesh.hpp b/src/BZ_MESH/bbox_mesh.hpp new file mode 100644 index 0000000..958e9bf --- /dev/null +++ b/src/BZ_MESH/bbox_mesh.hpp @@ -0,0 +1,150 @@ +/** + * @file box_bz.hpp + * @author remzerrr (remi.helleboid@gmail.com) + * @brief + * @version 0.1 + * @date 2022-08-25 + * + * @copyright Copyright (c) 2022 + * + */ + +#pragma once + +#include +#include +#include + +#include "vector.hpp" + +namespace bz_mesh { + +/** + * @brief Class representing a box in 3D space. + * The box is of the form AABB (Axis-Aligned Bounding Box) + * + * + */ +class bbox_mesh { + private: + double m_x_min{0.0}; + double m_x_max{0.0}; + double m_y_min{0.0}; + double m_y_max{0.0}; + double m_z_min{0.0}; + double m_z_max{0.0}; + + public: + bbox_mesh(){}; + bbox_mesh(double x_min, double x_max, double y_min, double y_max, double z_min, double z_max) + : m_x_min(x_min), + m_x_max(x_max), + m_y_min(y_min), + m_y_max(y_max), + m_z_min(z_min), + m_z_max(z_max) { + if (!check_order()) { + throw std::invalid_argument("Order of bbox_mesh corner is bad."); + } + } + + bbox_mesh(const vector3 &bottom_left_front_corner, const vector3 &up_right_back_corner) + : m_x_min(bottom_left_front_corner.x()), + m_x_max(up_right_back_corner.x()), + m_y_min(bottom_left_front_corner.y()), + m_y_max(up_right_back_corner.y()), + m_z_min(bottom_left_front_corner.z()), + m_z_max(up_right_back_corner.z()) { + if (!check_order()) { + throw std::invalid_argument("Order of bbox_mesh corner is bad."); + } + } + + bool check_order() const { return (m_x_min <= m_x_max) && (m_y_min <= m_y_max) && (m_z_min <= m_z_max); } + + double get_x_min() const { return m_x_min; } + double get_x_max() const { return m_x_max; } + double get_y_min() const { return m_y_min; } + double get_y_max() const { return m_y_max; } + double get_z_min() const { return m_z_min; } + double get_z_max() const { return m_z_max; } + + double get_x_size() const { return m_x_max - m_x_min; } + double get_z_size() const { return m_y_max - m_y_min; } + double get_y_size() const { return m_z_max - m_z_min; } + + double get_diagonal_size() const { + return sqrt(get_x_size() * get_x_size() + get_y_size() * get_y_size() + get_z_size() * get_z_size()); + } + + double get_volume() const { + const double surface = fabs(m_x_max - m_x_min) * fabs(m_y_max - m_y_min) * fabs(m_z_max - m_z_min); + return surface; + } + + vector3 get_center() const { + constexpr double one_half = 1.0 / 2.0; + return {one_half * (m_x_min + m_x_max), one_half * (m_y_min + m_y_max), one_half * (m_z_min + m_z_max)}; + } + + bool is_inside(const vector3 &location) const { + return (location.x() > m_x_min) && (location.x() < m_x_max) && (location.y() > m_y_min) && (location.y() < m_y_max) && + (location.z() > m_z_min) && (location.z() < m_z_max); + } + + /** + * @brief Split the box into 8 equal sub-boxes. + * + * @return std::array + */ + std::array split_3d_box_in_octants() const { + std::array octants; + const vector3 box_center = get_center(); + + octants[0] = bbox_mesh(m_x_min, box_center.x(), m_y_min, box_center.y(), m_z_min, box_center.z()); + octants[1] = bbox_mesh(box_center.x(), m_x_max, m_y_min, box_center.y(), m_z_min, box_center.z()); + octants[2] = bbox_mesh(box_center.x(), m_x_max, box_center.y(), m_y_max, m_z_min, box_center.z()); + octants[3] = bbox_mesh(m_x_min, box_center.x(), box_center.y(), m_y_max, m_z_min, box_center.z()); + octants[4] = bbox_mesh(m_x_min, box_center.x(), m_y_min, box_center.y(), box_center.z(), m_z_max); + octants[5] = bbox_mesh(box_center.x(), m_x_max, m_y_min, box_center.y(), box_center.z(), m_z_max); + octants[6] = bbox_mesh(box_center.x(), m_x_max, box_center.y(), m_y_max, box_center.z(), m_z_max); + octants[7] = bbox_mesh(m_x_min, box_center.x(), box_center.y(), m_y_max, box_center.z(), m_z_max); + return octants; + } + + void dilate(double factor) { + m_x_min *= factor; + m_x_max *= factor; + m_y_min *= factor; + m_y_max *= factor; + m_z_min *= factor; + m_z_max *= factor; + } + + void translate(const vector3 &translation) { + m_x_min += translation.x(); + m_x_max += translation.x(); + m_y_min += translation.y(); + m_y_max += translation.y(); + m_z_min += translation.z(); + m_z_max += translation.z(); + } + + bool is_overlapping(const bbox_mesh &second_box) const { + const bool noOverlap = this->m_x_min > second_box.m_x_max || second_box.m_x_min > this->m_x_max || + this->m_y_min > second_box.m_y_max || second_box.m_y_min > this->m_y_max || + this->m_z_min > second_box.m_z_max || second_box.m_z_min > this->m_z_max; + return !noOverlap; + } + + vector3 draw_uniform_random_point_inside_box() const; + vector3 draw_uniform_random_point_inside_box(std::mt19937 &random_generator) const; + + friend std::ostream &operator<<(std::ostream &os, const bbox_mesh &my_box) { + os << "bbox_mesh: x_min = " << my_box.m_x_min << " x_max = " << my_box.m_x_max << " y_min = " << my_box.m_y_min + << " y_max = " << my_box.m_y_max << " z_min = " << my_box.m_z_min << " z_max = " << my_box.m_z_max; + return os; + } +}; + +} // namespace bz_mesh \ No newline at end of file diff --git a/src/BZ_MESH/bz_mesh.cpp b/src/BZ_MESH/bz_mesh.cpp index 038f333..f81e82a 100644 --- a/src/BZ_MESH/bz_mesh.cpp +++ b/src/BZ_MESH/bz_mesh.cpp @@ -24,7 +24,6 @@ namespace bz_mesh { - void MeshBZ::shift_bz_center(const vector3& center) { m_center = center; for (auto&& vtx : m_list_vertices) { @@ -40,7 +39,7 @@ void MeshBZ::shift_bz_center(const vector3& center) { * @param filename * @param lattice_constant */ -void MeshBZ::read_mesh_geometry_from_msh_file(const std::string& filename) { +void MeshBZ::read_mesh_geometry_from_msh_file(const std::string& filename, bool normalize_by_fourier_factor) { std::cout << "Opening file " << filename << std::endl; gmsh::initialize(); gmsh::option::setNumber("General.Verbosity", 1); @@ -62,9 +61,8 @@ void MeshBZ::read_mesh_geometry_from_msh_file(const std::string& filename) { double lattice_constant = m_material.get_lattice_constant_meter(); std::cout << "Lattice const: " << lattice_constant << std::endl; std::cout << "V: " << std::pow(2.0 * M_PI, 3) / std::pow(lattice_constant, 3.0) << std::endl; - - double normalization_factor = 2.0 * M_PI / lattice_constant; - // double normalization_factor = 1.0; + const double fourier_factor = 2.0 * M_PI / lattice_constant; + double normalization_factor = normalize_by_fourier_factor ? fourier_factor : 1.0; for (std::size_t index_vertex = 0; index_vertex < size_nodes_tags; ++index_vertex) { m_list_vertices.push_back(Vertex(index_vertex, normalization_factor * nodeCoords[3 * index_vertex], @@ -103,6 +101,41 @@ void MeshBZ::read_mesh_geometry_from_msh_file(const std::string& filename) { m_total_volume = compute_mesh_volume(); } +bbox_mesh MeshBZ::compute_bounding_box() const { + double x_min = std::numeric_limits::max(); + double y_min = std::numeric_limits::max(); + double z_min = std::numeric_limits::max(); + double x_max = std::numeric_limits::min(); + double y_max = std::numeric_limits::min(); + double z_max = std::numeric_limits::min(); + for (auto&& vtx : m_list_vertices) { + const vector3& position = vtx.get_position(); + x_min = std::min(x_min, position.x()); + y_min = std::min(y_min, position.y()); + z_min = std::min(z_min, position.z()); + x_max = std::max(x_max, position.x()); + y_max = std::max(y_max, position.y()); + z_max = std::max(z_max, position.z()); + } + vector3 min_corner(x_min, y_min, z_min); + vector3 max_corner(x_max, y_max, z_max); + return bbox_mesh(min_corner, max_corner); +} + +void MeshBZ::build_search_tree() { + bbox_mesh mesh_bbox = compute_bounding_box(); + std::cout << "Mesh bounding box: " << mesh_bbox << std::endl; + mesh_bbox.dilate(1.05); + // bbox_mesh.translate({0.0, 0.0, 0.0}); + auto start = std::chrono::high_resolution_clock::now(); + m_search_tree = std::make_unique(get_list_p_tetra(), mesh_bbox); + auto end = std::chrono::high_resolution_clock::now(); + auto total = std::chrono::duration_cast(end - start).count(); + std::cout << "Octree built in " << total / 1000.0 << "s" << std::endl; +} + +Tetra* MeshBZ::find_tetra_at_location(const vector3& location) const { return m_search_tree->find_tetra_at_location(location); } + /** * @brief Get the nearest k index object. * Brute force search of the nearest k-point index. :( @@ -337,7 +370,6 @@ vector3 MeshBZ::retrieve_k_inside_mesh_geometry(const vector3& k) const { } std::cout << "No k-point inside the mesh geometry found." << std::endl; throw std::runtime_error("No k-point inside the mesh geometry found."); - } } // namespace bz_mesh \ No newline at end of file diff --git a/src/BZ_MESH/bz_mesh.hpp b/src/BZ_MESH/bz_mesh.hpp index 5470d50..2e9e68d 100644 --- a/src/BZ_MESH/bz_mesh.hpp +++ b/src/BZ_MESH/bz_mesh.hpp @@ -15,6 +15,8 @@ #include "iso_triangle.hpp" #include "mesh_tetra.hpp" #include "mesh_vertex.hpp" +#include "octree_bz.hpp" +#include "vector.hpp" namespace bz_mesh { @@ -46,6 +48,12 @@ class MeshBZ { */ std::vector m_list_tetrahedra; + /** + * @brief Octree used to search for the tetrahedra that are overlapping with a given point. + * + */ + std::unique_ptr m_search_tree; + /** * @brief The indexes of the valence bands within the Vertex list of energies. * For example, if indices_valence_bands = [0, 1, 2] it means that @@ -77,13 +85,25 @@ class MeshBZ { double m_total_volume = 0.0; public: + MeshBZ() = default; MeshBZ(const EmpiricalPseudopotential::Material& material) : m_material(material){}; MeshBZ(const MeshBZ&) = default; - vector3 get_center() const { return m_center; } - void shift_bz_center(const vector3& shift); + vector3 get_center() const { return m_center; } + void shift_bz_center(const vector3& shift); + + + bbox_mesh compute_bounding_box() const; + void build_search_tree(); + std::vector get_list_p_tetra() { + std::vector tetra_pointers; + tetra_pointers.reserve(m_list_tetrahedra.size()); + std::transform(m_list_tetrahedra.begin(), m_list_tetrahedra.end(), std::back_inserter(tetra_pointers), [](Tetra& tetra) { return &tetra; }); + return tetra_pointers; + } + Tetra* find_tetra_at_location(const vector3& location) const; - void read_mesh_geometry_from_msh_file(const std::string& filename); + void read_mesh_geometry_from_msh_file(const std::string& filename, bool normalize_by_fourier_factor=true); void read_mesh_bands_from_msh_file(const std::string& filename); void add_new_band_energies_to_vertices(const std::vector& energies_at_vertices); void compute_min_max_energies_at_tetras(); diff --git a/src/BZ_MESH/dielectric_mesh.cpp b/src/BZ_MESH/dielectric_mesh.cpp index 1f4d20b..8a7f04e 100644 --- a/src/BZ_MESH/dielectric_mesh.cpp +++ b/src/BZ_MESH/dielectric_mesh.cpp @@ -11,9 +11,9 @@ #include "dielectric_mesh.hpp" -#include "gmsh.h" #include +#include "gmsh.h" namespace bz_mesh { @@ -24,15 +24,15 @@ void DielectricMesh::read_dielectric_file(const std::string& filename) { gmsh::open(filename); std::vector viewTags; gmsh::view::getTags(viewTags); - int nb_energies = viewTags.size(); + int nb_energies = viewTags.size(); if (nb_energies % 2 != 0) { // The number of viewTags should be even because we have the real and imaginary part of the dielectric function for each energy. std::cerr << "The number of energies is not even. The file is corrupted." << std::endl; return; } nb_energies /= 2; - int count_energy = 0; - std::vector energies; + int count_energy = 0; + std::vector energies; std::vector> real_dielectric; std::vector> imag_dielectric; for (auto&& tag : viewTags) { @@ -53,11 +53,11 @@ void DielectricMesh::read_dielectric_file(const std::string& filename) { // If the name of the view starts with "eps_r" we store the data in the real part of the dielectric function. // If the name of the view starts with "eps_i" we store the data in the imaginary part of the dielectric function. - std::regex re_real("eps_r"); - std::regex re_imag("eps_i"); + std::regex re_real("eps_r"); + std::regex re_imag("eps_i"); - std::smatch match_real; - std::smatch match_imag; + std::smatch match_real; + std::smatch match_imag; std::string type; std::vector tags; double time; @@ -86,10 +86,50 @@ void DielectricMesh::read_dielectric_file(const std::string& filename) { for (std::size_t idx_node = 0; idx_node < real_dielectric.size(); ++idx_node) { m_dielectric_function[idx_node].resize(energies.size()); for (std::size_t idx_energy = 0; idx_energy < energies.size(); ++idx_energy) { - m_dielectric_function[idx_node][idx_energy] = std::complex(real_dielectric[idx_node][idx_energy], imag_dielectric[idx_node][idx_energy]); + m_dielectric_function[idx_node][idx_energy] = + std::complex(real_dielectric[idx_node][idx_energy], imag_dielectric[idx_node][idx_energy]); } } - +} + +std::pair DielectricMesh::find_closest_energy(double energy) const { + if (m_energies.empty()) { + return std::make_pair(0, 0.0); + } + if (energy < m_energies.front()) { + return std::make_pair(0, 0.0); + } + if (energy > m_energies.back()) { + return std::make_pair(m_energies.size() - 1, 0.0); + } + auto it = std::lower_bound(m_energies.begin(), m_energies.end(), energy); + if (it == m_energies.begin()) { + return std::make_pair(0, 0.0); + } + if (it == m_energies.end()) { + return std::make_pair(m_energies.size() - 1, 0.0); + } + std::size_t idx = std::distance(m_energies.begin(), it); + double t = (energy - m_energies[idx - 1]) / (m_energies[idx] - m_energies[idx - 1]); + return std::make_pair(idx - 1, t); +} + +complex_d DielectricMesh::interpolate_dielectric_function(const vector3& k, double energy) const { + auto k_positive{vector3{std::abs(k.x()), std::abs(k.y()), std::abs(k.z())}}; + Tetra* p_tetra = find_tetra_at_location(k_positive); + if (p_tetra == nullptr) { + return 0.0; + } + const std::array barycentric_coordinates = p_tetra->compute_barycentric_coordinates(k); + + const std::size_t idx_node = p_tetra->get_index(); + const std::pair closest_energy = find_closest_energy(energy); + const std::size_t idx_energy = closest_energy.first; + const double t = closest_energy.second; + const complex_d dielectric_function_0 = m_dielectric_function[idx_node][idx_energy]; + const complex_d dielectric_function_1 = m_dielectric_function[idx_node][idx_energy + 1]; + complex_d dielectric_function = dielectric_function_0 * (1.0 - t) + dielectric_function_1 * t; + return dielectric_function; } } // namespace bz_mesh \ No newline at end of file diff --git a/src/BZ_MESH/dielectric_mesh.hpp b/src/BZ_MESH/dielectric_mesh.hpp index 94c2560..11084bf 100644 --- a/src/BZ_MESH/dielectric_mesh.hpp +++ b/src/BZ_MESH/dielectric_mesh.hpp @@ -43,6 +43,7 @@ class DielectricMesh : public MeshBZ { std::vector> m_dielectric_function; public: + DielectricMesh() = default; DielectricMesh(const EmpiricalPseudopotential::Material& material) : MeshBZ(material) {} /** @@ -51,6 +52,27 @@ class DielectricMesh : public MeshBZ { * @param filename Path to the file containing the dielectric function. */ void read_dielectric_file(const std::string& filename); + + /** + * @brief Find the closest energy in the list of energies. + * Return the index (idx) of the stored energy directly below the given energy and the fraction (t) of the distance between the two + * closest energies. The dielectric function at the given energy can be interpolated as: m_dielectric_function[idx_node][idx] * (1 - t) + * + m_dielectric_function[idx_node][idx + 1] * t + * + * + * @param energy + * @return std::pair + */ + std::pair find_closest_energy(double energy) const; + + /** + * @brief Interpolate the dielectric function at a given k-point and energy. + * + * @param k + * @param energy + * @return complex_d + */ + complex_d interpolate_dielectric_function(const vector3& k, double energy) const; }; } // namespace bz_mesh \ No newline at end of file diff --git a/src/BZ_MESH/impact_ionization.cpp b/src/BZ_MESH/impact_ionization.cpp index 5c46447..95c6f12 100644 --- a/src/BZ_MESH/impact_ionization.cpp +++ b/src/BZ_MESH/impact_ionization.cpp @@ -41,6 +41,13 @@ ImpactIonization::ImpactIonization(const EmpiricalPseudopotential::Material& mat m_initial_mesh_path = initial_mesh_path; } +void ImpactIonization::read_dielectric_file(const std::string& filename) { + bool normalize_by_fourier_factor = false; + m_dielectric_mesh.read_mesh_geometry_from_msh_file(filename, normalize_by_fourier_factor); + m_dielectric_mesh.build_search_tree(); + m_dielectric_mesh.read_dielectric_file(filename); +} + void ImpactIonization::compute_eigenstates() { int nb_bands_to_use = 4; bz_mesh::BZ_States my_bz_mesh(m_material); @@ -75,7 +82,10 @@ void ImpactIonization::compute_eigenstates() { continue; } std::cout << "G_BZ: " << G_BZ << std::endl; - auto ptr_BZ_states = std::make_unique(my_bz_mesh); + auto ptr_BZ_states = std::make_unique(m_material); + ptr_BZ_states->set_nb_bands(nb_bands_to_use); + ptr_BZ_states->set_basis_vectors(basis); + ptr_BZ_states->read_mesh_geometry_from_msh_file(m_initial_mesh_path); ptr_BZ_states->shift_bz_center(G_BZ); ptr_BZ_states->compute_eigenstates(nb_threads); m_list_BZ_states.push_back(std::move(ptr_BZ_states)); diff --git a/src/BZ_MESH/impact_ionization.hpp b/src/BZ_MESH/impact_ionization.hpp index 94a806c..2cc9aa8 100644 --- a/src/BZ_MESH/impact_ionization.hpp +++ b/src/BZ_MESH/impact_ionization.hpp @@ -18,6 +18,7 @@ #include "Material.h" #include "bz_mesh.hpp" #include "bz_states.hpp" +#include "dielectric_mesh.hpp" namespace bz_mesh { @@ -51,6 +52,12 @@ class ImpactIonization { */ EmpiricalPseudopotential::Material m_material; + /** + * @brief Mesh of the dielectric function. + * + */ + DielectricMesh m_dielectric_mesh; + /** * @brief List of the results of the impact ionization rate. * @@ -59,6 +66,7 @@ class ImpactIonization { public: ImpactIonization(const EmpiricalPseudopotential::Material& material, const std::string& initial_mesh_path); + void read_dielectric_file(const std::string& filename); double get_max_radius_G0_BZ() const { return m_max_radius_G0_BZ; } void set_max_radius_G0_BZ(double max_radius_G0_BZ) { m_max_radius_G0_BZ = max_radius_G0_BZ; } diff --git a/src/BZ_MESH/mesh_tetra.cpp b/src/BZ_MESH/mesh_tetra.cpp index 2325869..22f2bf5 100644 --- a/src/BZ_MESH/mesh_tetra.cpp +++ b/src/BZ_MESH/mesh_tetra.cpp @@ -17,13 +17,30 @@ #include #include "Constants.hpp" +#include "bbox_mesh.hpp" #include "iso_triangle.hpp" - namespace bz_mesh { std::vector Tetra::ms_case_stats = {0, 0, 0, 0, 0}; +bbox_mesh Tetra::compute_bounding_box() const { + std::array coordinates_x; + std::array coordinates_y; + std::array coordinates_z; + for (int i = 0; i < 4; ++i) { + coordinates_x[i] = m_list_vertices[i]->get_position().x(); + coordinates_y[i] = m_list_vertices[i]->get_position().y(); + coordinates_z[i] = m_list_vertices[i]->get_position().z(); + } + auto min_max_x = std::minmax_element(coordinates_x.begin(), coordinates_x.end()); + auto min_max_y = std::minmax_element(coordinates_y.begin(), coordinates_y.end()); + auto min_max_z = std::minmax_element(coordinates_z.begin(), coordinates_z.end()); + return bbox_mesh(*min_max_x.first, *min_max_x.second, *min_max_y.first, *min_max_y.second, *min_max_z.first, *min_max_z.second); +} + +const bbox_mesh& Tetra::get_bounding_box() const { return m_bbox; } + /** * @brief Construct a new Tetra by passing directly the array of the four pointers to the vertices. * @@ -40,6 +57,7 @@ Tetra::Tetra(std::size_t index, const std::array& list_vertices) m_list_edges[4] = compute_edge(3, 1); m_list_edges[5] = compute_edge(3, 2); m_signed_volume = compute_signed_volume(); + m_bbox = compute_bounding_box(); } vector3 Tetra::compute_barycenter() const { diff --git a/src/BZ_MESH/mesh_tetra.hpp b/src/BZ_MESH/mesh_tetra.hpp index 647a8d9..928a5e1 100644 --- a/src/BZ_MESH/mesh_tetra.hpp +++ b/src/BZ_MESH/mesh_tetra.hpp @@ -17,6 +17,7 @@ #include #include +#include "bbox_mesh.hpp" #include "mesh_vertex.hpp" namespace bz_mesh { @@ -43,6 +44,12 @@ class Tetra { */ std::array m_list_edges{}; + /** + * @brief Bounding box of the tetrahedra. + * + */ + bbox_mesh m_bbox; + /** * @brief Signed volume of the tetrahedra. * The sign depends on the "orientation" of the tetrahedra. @@ -87,7 +94,10 @@ class Tetra { Tetra() = delete; Tetra(std::size_t index, const std::array& list_vertices); - void compute_min_max_energies_at_bands(); + + const bbox_mesh& get_bounding_box() const; + bbox_mesh compute_bounding_box() const; + void compute_min_max_energies_at_bands(); std::size_t get_index() const { return m_index; } const std::array& get_list_vertices() const { return m_list_vertices; } diff --git a/src/BZ_MESH/octree_bz.cpp b/src/BZ_MESH/octree_bz.cpp new file mode 100644 index 0000000..eb890d9 --- /dev/null +++ b/src/BZ_MESH/octree_bz.cpp @@ -0,0 +1,80 @@ +/** + * @file Octree_mesh.cpp + * @author remzerrr (remi.helleboid@gmail.com) + * @brief + * @version 0.1 + * @date 2022-08-25 + * + * @copyright Copyright (c) 2022 + * + */ + +#include "octree_bz.hpp" + + +#include +#include +#include +#include + +namespace bz_mesh { + +Octree_mesh::Octree_mesh(const std::vector &list_tetras, const bbox_mesh &bounding_box) { + if (list_tetras.size() <= max_number_of_elements_per_node || bounding_box.get_diagonal_size() < min_size_of_a_node) { + m_is_leaf = true; + m_list_tetras = list_tetras; + m_node_box = bounding_box; + return; + } + + m_is_leaf = false; + m_node_box = bounding_box; + std::array sub_boxes = m_node_box.split_3d_box_in_octants(); + + for (int i = 0; i < 8; i++) { + m_list_sub_nodes.push_back( + std::make_unique(Octree_mesh::find_overlapping_tetras(list_tetras, sub_boxes[i]), sub_boxes[i])); + } +} + +/** + * @brief Find the tetra from the list of tetras that are overlapping with the bounding box. + * + * @note The list of pointers to tetras is passed by copy because it is modified during the search. + * Perf should be compared with a list of pointers to tetras that is passed by reference and copied within the function. + * + * @param list_p_tetras + * @param bounding_box + * @return std::vector + */ +std::vector Octree_mesh::find_overlapping_tetras(const std::vector &list_p_tetras, const bbox_mesh &bounding_box) { + std::vector list_overlapping_tetras; + for (auto &p_tetra : list_p_tetras) { + if (bounding_box.is_overlapping(p_tetra->get_bounding_box())) { + list_overlapping_tetras.push_back(p_tetra); + } + } + return list_overlapping_tetras; +} + +/** + * @brief Find the tetra that contains the location. + * If none is found, return nullptr, else return a pointer to the tetra. + * + * @param location + * @return Tetra* + */ +Tetra *Octree_mesh::find_tetra_at_location(const vector3 &location) const { + if (!m_is_leaf) { + for (auto &&p_sub_node : m_list_sub_nodes) { + if (p_sub_node->is_inside(location)) { + return p_sub_node->find_tetra_at_location(location); + } + } + } + const auto it_tetra = + std::find_if(m_list_tetras.begin(), m_list_tetras.end(), [&](const Tetra *p_tetra) { return p_tetra->is_location_inside(location); }); + return it_tetra != m_list_tetras.end() ? *it_tetra : nullptr; +} + +} // namespace bz_mesh diff --git a/src/BZ_MESH/octree_bz.hpp b/src/BZ_MESH/octree_bz.hpp new file mode 100644 index 0000000..59ba6e5 --- /dev/null +++ b/src/BZ_MESH/octree_bz.hpp @@ -0,0 +1,70 @@ +/** + * @file Octree_mesh.hpp + * @author remzerrr (remi.helleboid@gmail.com) + * @brief + * @version 0.1 + * @date 2022-08-25 + * + * @copyright Copyright (c) 2022 + * + */ + +#pragma once + +#include +#include + +#include "mesh_tetra.hpp" +#include "vector.hpp" +#include "bbox_mesh.hpp" + +namespace bz_mesh { + + + +class Octree_mesh { + protected: + /** + * @brief Maximum number of elements in a leaf node. + * If the number of elements in a node is greater than this value, the node is split into 8 children. + */ + static constexpr int max_number_of_elements_per_node = 32; + + /** + * @brief Minimum size of a node. + * If the size of a node is smaller than this value, the node is not split anymore. + */ + static constexpr double min_size_of_a_node = 1e-2; + + /** + * @brief Bounding box of the node. + */ + bbox_mesh m_node_box; + + /** + * @brief A leaf node is a node that has no children. + * + */ + bool m_is_leaf = false; + + /** + * @brief List of pointers to elements overlapping in the node. + */ + std::vector m_list_tetras; + + /** + * @brief List of children of the node. + */ + std::vector> m_list_sub_nodes; + + std::vector find_overlapping_tetras(const std::vector &list_p_tetras, const bbox_mesh &bounding_box); + bool is_inside(const vector3 &location) const { return m_node_box.is_inside(location); } + + public: + Octree_mesh(){}; + Octree_mesh(const std::vector &list_tetras, const bbox_mesh &bounding_box); + + Tetra *find_tetra_at_location(const vector3 &location) const; +}; + +} // namespace bz_mesh