Continue implem of dielectricfunc for impact io

apps/impact_io.cpp

@@ -34,11 +34,11 @@ int main(int argc, char *argv[]) {
     TCLAP::ValueArg<int>         arg_nb_bands("b", "nbands", "Number of bands to consider", false, 12, "int");
     TCLAP::ValueArg<int>         arg_nb_threads("j", "nthreads", "number of threads to use.", false, 1, "int");
     TCLAP::ValueArg<double>      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<double> elapsed_seconds = end - start;

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})

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 <memory>
+#include <random>
+#include <vector>
+
+#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<bbox_mesh, 8>
+     */
+    std::array<bbox_mesh, 8> split_3d_box_in_octants() const {
+        std::array<bbox_mesh, 8> 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

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<double>::max();
+    double y_min = std::numeric_limits<double>::max();
+    double z_min = std::numeric_limits<double>::max();
+    double x_max = std::numeric_limits<double>::min();
+    double y_max = std::numeric_limits<double>::min();
+    double z_max = std::numeric_limits<double>::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<Octree_mesh>(get_list_p_tetra(), mesh_bbox);
+    auto end      = std::chrono::high_resolution_clock::now();
+    auto total    = std::chrono::duration_cast<std::chrono::milliseconds>(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

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<Tetra> m_list_tetrahedra;
 
+    /**
+     * @brief Octree used to search for the tetrahedra that are overlapping with a given point.
+     *
+     */
+    std::unique_ptr<Octree_mesh> 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<Tetra*> get_list_p_tetra()  {
+        std::vector<Tetra*> 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<double>& energies_at_vertices);
     void compute_min_max_energies_at_tetras();