Type stable

Project.toml

@@ -1,30 +1,33 @@
 name = "BrillouinZoneMeshes"
 uuid = "9f696214-9961-49f3-89a2-cb1ea204eb6e"
 authors = ["Xiansheng Cai", "Tao Wang", "Kun Chen"]
-version = "0.1.0"
+version = "0.2.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
 CompositeGrids = "b5136c89-beeb-4521-9139-60d2cac8be56"
-FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
-QHull = "a8468747-bd6f-53ef-9e5c-744dbc5c59e7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
-Spglib = "f761d5c5-86db-4880-b97f-9680a7cccfb5"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 spglib_jll = "ac4a9f1e-bdb2-5204-990c-47c8b2f70d4e"
+Spglib = "f761d5c5-86db-4880-b97f-9680a7cccfb5"
 
 [compat]
 AbstractTrees = "0.3, 0.4"
 CompositeGrids = "0.0.6, 0.1"
 StaticArrays = "1"
+Roots = "2.0"
+Spglib = "0.6.0"
 julia = "1.6"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+QHull = "a8468747-bd6f-53ef-9e5c-744dbc5c59e7"
+
 
 [targets]
 test = ["Test"]

README.md

@@ -33,7 +33,7 @@ Setup with:
 
 and here are some examples of usage:
 
-```jldoctest
+```julia
 julia> using BrillouinZoneMeshes; DIM = 2; lattice = Matrix([1.0 0; 0 1]'); br = BZMeshes.Cell(lattice=lattice); umesh = BZMeshes.UniformBZMesh(cell=br, size=(4, 4)); mm = MeshMaps.MeshMap(umesh); rmesh = ReducedBZMesh(umesh, mm);
 
 julia> length(umesh), length(rmesh)
@@ -102,4 +102,4 @@ Monkhorst-Pack meshes, could be generated via `BZMeshes.UniformBZMesh` with
 different parameters. The default parameter of the constructor of 
 `BZMeshes.UniformBZMesh` generates Gamma-centered mesh, while two constructors
 for M-P mesh, `Monkhorst_Pack` and `DFTK_Monkhorst_Pack`, follow conventions
-from VASP and DFTK respectively.
+from VASP and DFTK respectively.

src/BaseMesh.jl

@@ -41,7 +41,8 @@ export fractional_coordinates, cartesian_coordinates
 
 function Base.getindex(mesh::AbstractUniformMesh{T,DIM}, inds...) where {T,DIM}
     n = SVector{DIM,Int}(inds)
-    return mesh.origin + lattice_vector(mesh) * ((n .- 1 .+ mesh.shift) ./ mesh.size)
+    mshift = SVector{DIM,T}(mesh.shift)
+    return mesh.origin + lattice_vector(mesh) * ((n .- 1 .+ mshift) ./ mesh.size)
 end
 
 function Base.getindex(mesh::AbstractUniformMesh, I::Int)
@@ -50,7 +51,8 @@ end
 
 function Base.getindex(mesh::AbstractUniformMesh{T,DIM}, ::Type{<:FracCoords}, I::Int) where {T,DIM}
     n = SVector{DIM,Int}(AbstractMeshes._ind2inds(mesh.size, I))
-    return inv_lattice_vector(mesh) * mesh.origin + (n .- 1 .+ mesh.shift) ./ mesh.size
+    mshift = SVector{DIM,T}(mesh.shift)
+    return inv_lattice_vector(mesh) * mesh.origin + (n .- 1 .+ mshift) ./ mesh.size
 end
 
 """

test/AbstractMeshes.jl

@@ -1,16 +1,6 @@
 @testset "AbstractMeshes" begin
     using BrillouinZoneMeshes.AbstractMeshes
 
-    function test_func_not_implemented(func, obj)
-        # if a func required is not implemented for obj
-        # an error occur
-        try
-            func(obj)
-        catch e
-            @test e isa ErrorException
-        end
-    end
-
     # create a random concrete mesh
     DIM = 2
     N1, N2 = 3, 5
@@ -38,24 +28,24 @@
     # test error thrown from funcs not implemented
     struct NotAMesh{T,DIM} <: AbstractMesh{T,DIM} end
     notamesh = NotAMesh{Float64,3}()
-    test_func_not_implemented(println, notamesh)
+    @test_throws ErrorException println(notamesh)
 
-    test_func_not_implemented(x -> getindex(x, 1), notamesh)
-    test_func_not_implemented(x -> getindex(x, 1, 2, 3), notamesh)
-    test_func_not_implemented(x -> getindex(x, FracCoords, 1), notamesh)
-    test_func_not_implemented(x -> getindex(x, FracCoords, 1, 2, 3), notamesh)
+    @test_throws ErrorException notamesh[1]
+    @test_throws ErrorException notamesh[1, 2, 3]
+    @test_throws ErrorException notamesh[FracCoords, 1]
+    @test_throws ErrorException notamesh[FracCoords, 1, 2, 3]
 
-    test_func_not_implemented(x -> locate(x, 1), notamesh)
-    test_func_not_implemented(x -> volume(x, 1), notamesh)
-    test_func_not_implemented(volume, notamesh)
+    @test_throws ErrorException locate(notamesh, 1)
+    @test_throws ErrorException volume(notamesh, 1)
+    @test_throws ErrorException volume(notamesh)
 
-    test_func_not_implemented(lattice_vector, notamesh)
-    test_func_not_implemented(inv_lattice_vector, notamesh)
-    test_func_not_implemented(cell_volume, notamesh)
+    @test_throws ErrorException lattice_vector(notamesh)
+    @test_throws ErrorException inv_lattice_vector(notamesh)
+    @test_throws ErrorException cell_volume(notamesh)
 
-    test_func_not_implemented(x -> integrate([1,], x), notamesh)
-    test_func_not_implemented(x -> interp([1,], x, 1), notamesh)
+    @test_throws ErrorException integrate([1,], notamesh)
+    @test_throws ErrorException interp([1,], notamesh, 1)
 
-    test_func_not_implemented(x -> interval(x, 1), notamesh)
+    @test_throws ErrorException interval(notamesh, 1)
 
 end

test/BZMeshes.jl

@@ -28,6 +28,10 @@
             println(bzmesh)
             display(bzmesh)
 
+            @inferred bzmesh[1]
+            @inferred bzmesh[1, 1]
+            @inferred bzmesh[AbstractMeshes.FracCoords, 1]
+
             for (pi, p) in enumerate(bzmesh)
                 @test bzmesh[pi] ≈ p # linear index
                 inds = AbstractMeshes._ind2inds(size(bzmesh), pi)

test/BaseMesh.jl

@@ -1,16 +1,6 @@
 @testset "Base Mesh" begin
     rng = MersenneTwister(1234)
 
-    function test_func_not_implemented(func, obj)
-        # if a func required is not implemented for obj
-        # an error occur
-        try
-            func(obj)
-        catch e
-            @test e isa ErrorException
-        end
-    end
-
     @testset "UMesh" begin
         DIM = 2
         N1, N2 = 3, 5
@@ -19,6 +9,10 @@
         cell = BZMeshes.Cell(lattice=lattice)
         mesh = BaseMesh.UMesh(br=cell, origin=ones(DIM) ./ 2, size=(N1, N2), shift=zeros(DIM))
 
+        @inferred mesh[1]
+        @inferred mesh[1, 1]
+        @inferred mesh[AbstractMeshes.FracCoords, 1]
+
         @test length(mesh) == N1 * N2
         @test size(mesh) == (N1, N2)
 
@@ -45,7 +39,7 @@
 
         # basics
         struct NotAPM{T,DIM} <: BaseMesh.AbstractProdMesh{T,DIM} end
-        test_func_not_implemented(x -> BaseMesh._getgrid(x, 1), NotAPM{Int,3}())
+        @test_throws ErrorException BaseMesh._getgrid(NotAPM{Int,3}(), 1)
 
         @testset "DirectProdMesh" begin
             N, M = 3, 2
@@ -62,6 +56,8 @@
             dpm = DirectProdMesh(r, theta, phi)
             println(size(dpm))
 
+            @inferred dpm[1]
+
             vol = 0.0
             for (pi, p) in enumerate(dpm)
                 i, j, k = AbstractMeshes._ind2inds(size(dpm), pi)
@@ -91,6 +87,9 @@
             cm = ProdMesh(grids, theta)
             println([cm.grids[i].panel[2] for i in 1:length(theta)])
             println(size(cm))
+
+            @inferred cm[1]
+
             for j in 1:length(cm.mesh)
                 for i in 1:length(cm.grids[j])
                     p = cm[i, j]
@@ -117,6 +116,8 @@
 
         cm2 = ChebMesh(origin, latvec, cm)
 
+        @inferred cm[1]
+
         vol = 0.0
         for (i, x) in enumerate(cm2)
             @test cm2[i] ≈ x

test/MeshMap.jl

@@ -1,22 +1,12 @@
 @testset "MeshMaps" begin
 
-    function test_func_not_implemented(func, obj)
-        # if a func required is not implemented for obj
-        # an error occur
-        try
-            func(obj)
-        catch e
-            @test e isa ErrorException
-        end
-    end
-
     locate, volume = MeshMaps.locate, MeshMaps.volume
 
     @testset "MeshMap" begin
 
         struct NotAMesh{T,DIM} <: AbstractMesh{T,DIM} end
         notamesh = NotAMesh{Float64,3}()
-        test_func_not_implemented(MeshMap, notamesh)
+        @test_throws ErrorException MeshMap(notamesh)
 
         # test MeshMap constructor
         map = [1, 2, 2, 1, 2, 6, 6, 2, 2, 6, 6, 2, 1, 2, 2, 1]

test/test_typestable.jl

@@ -0,0 +1,14 @@
+
+using BrillouinZoneMeshes
+using Test
+
+@testset "Type Stable?" begin
+    DIM = 2
+    N1, N2 = 3, 5
+    lattice = Matrix([1/N1/2 0; 0 1.0/N2/2]') .* 2π
+    # so that bzmesh[i,j] = (2i-1,2j-1)
+    cell = BZMeshes.Cell(lattice=lattice)
+    mesh = BaseMesh.UMesh(br=cell, origin=ones(DIM) ./ 2, size=(N1, N2), shift=zeros(DIM))
+
+    @inferred mesh[1]
+end