diff --git a/example/compositemesh.jl b/example/compositemesh.jl index 54e44e5..6a316be 100644 --- a/example/compositemesh.jl +++ b/example/compositemesh.jl @@ -7,146 +7,57 @@ using BrillouinZoneMeshes.AbstractMeshes using BrillouinZoneMeshes.LinearAlgebra using BrillouinZoneMeshes.StaticArrays -using Test +lattice = Matrix([1 0; 0 1.0]') +atoms = [1] +positions = [zeros(2)] +br = BZMeshes.Cell(lattice=lattice, atoms=atoms, positions=positions) + +_dispersion(k) = -sum(cos.(k)) + 0.1 +# _dispersion(k) = dot(k, k) - π^2 / 4 +T = 0.001 + +function dispersion(k) + # constrain to 1st bz + _k = [(abs(p) < π) ? p : π for p in k] + return _dispersion(_k) +end -function AbstractMeshes.interval(mesh::ProdMesh, I::Int) - inds = AbstractMeshes._ind2inds(size(mesh), I) - J = Base._sub2ind(size(mesh)[2:end], inds[2:end]...) - x1, x2 = AbstractMeshes.interval(mesh.grids[J], inds[1]) - if length(inds) > 2 - return [(x1, x2), AbstractMeshes.interval(mesh.mesh, J)...] +function integrand(k) + #if dot(k, k) <= π^2 + if abs(k[1]) <= π && abs(k[2]) <= π + return 1 / (dispersion(k)^2 + T^2) else - return [(x1, x2), AbstractMeshes.interval(mesh.mesh, J)] + return 0.0 end end -struct CompositeMesh{T,PM,SM} <: AbstractMesh{T,1} - panelmesh::PM - submeshes::Vector{SM} -end - -function CompositeMesh(panelmesh::PM, N) where {PM} - T = eltype(panelmesh) - submeshes = [] - for (i, p) in enumerate(panelmesh) - intervals = AbstractMeshes.interval(panelmesh, i) - DIM = length(intervals) - origin = [intervals[j][1] for j in 1:DIM] - lattice = diagm(DIM, DIM, [intervals[j][2] - intervals[j][1] for j in 1:DIM]) - cm = ChebMesh(origin, lattice, DIM, N) - push!(submeshes, cm) +function int_mesh(mesh) + data = zeros(Float64, size(mesh)) + for i in 1:length(mesh) + data[i] = integrand(mesh[i]) end - - SM = typeof(submeshes[1]) - return CompositeMesh{T,PM,SM}(panelmesh, submeshes) + return integrate(data, mesh) end +N = 8 +bound = [-π, π] +theta = SimpleGrid.Uniform(bound, N; isperiodic=true) +# bzmesh = BZMeshes.PolarMesh(dispersion=dispersion, anglemesh=theta, cell=br, +# kmax=π, Nloggrid=5, Nbasegrid=4, minterval=0.001) +pm = BZMeshes.CompositePolarMesh(dispersion=dispersion, + anglemesh=theta, cell=br, basegridtype=:uniform, + kmax=π * sqrt(2.4), Nloggrid=8, Nbasegrid=3, minterval=0.1T, N=4) -@testset "CompositeMesh" begin - @testset "Constructor" begin - a, b = 0.8, 1.2 - - N, M = 3, 2 - # theta grid dense around 0 and π - theta = CompositeGrid.LogDensedGrid( - :cheb, - [-π, π], - [-π, 0, π], - N, - 0.1, - M - ) - println(theta) - grids = [CompositeGrid.LogDensedGrid(:cheb, [0.0, 2.0], [sqrt(a * cos(θ)^2 + b * sin(θ)^2),], N, 0.1, M) for θ in theta] - - pm = ProdMesh(theta, grids) - cm = CompositeMesh(pm, 3) - end - - # @testset "2D" begin - # # given dispersion function accept a k in cartesian - # # goal is to find k_F at direction specified by angle - # dispersion(k) = dot(k, k) - 1.0 - - # # 2d - # N = 10 - # bound = [-π, π] - # theta = SimpleGrid.Uniform(bound, N; isperiodic=true) - - # k_F_previous = 0.0 - # for θ in theta - # f(k) = dispersion(BZMeshes._polar2cart(Polar(k, θ))) - # k_F = find_zero(f, k_F_previous) - # @test k_F ≈ 1.0 - # @test find_kFermi(dispersion, θ; kinit=k_F_previous) ≈ 1.0 - # k_F_previous = k_F - # end - - # # grids = kF_densed_kgrids(dispersion=dispersion, anglemesh=theta, - # # bound=[0.0, 2.0]) - # # cm = CompositeMesh(theta, grids) - # DIM = 2 - # lattice = Matrix([1.0 0; 0 1]') - # br = BZMeshes.Cell(lattice=lattice) - - # pm = PolarMesh(dispersion=dispersion, anglemesh=theta, cell=br, kmax=2.0) - # @test AbstractMeshes.volume(pm) ≈ 4π - - # end - - # @testset "3D" begin - # dispersion(k) = dot(k, k) - 1.0 - - # N = 6 - # bound = [-π, π] - # phi = SimpleGrid.Uniform(bound, N; isperiodic=true) - - # N = 4 - # bound = [-π / 2, π / 2] - # theta = SimpleGrid.Uniform(bound, N; isperiodic=true) - - # am = CompositeMesh(phi, [theta for i in 1:length(phi)]) - - # k_F_previous = 0.0 - # for ap in am - # f(k) = dispersion(BZMeshes._spherical2cart(Spherical(k, ap...))) - # k_F = find_zero(f, k_F_previous) - # @test k_F ≈ 1.0 - # @test find_kFermi(dispersion, ap; kinit=k_F_previous) ≈ 1.0 - # k_F_previous = k_F - # end - - # DIM = 3 - # lattice = Matrix([1.0 1.0 0; 1 0 1; 0 1 1]') - # br = BZMeshes.Cell(lattice=lattice) - - # pm = PolarMesh(dispersion=dispersion, anglemesh=am, cell=br, kmax=2.0) - # @test AbstractMeshes.volume(pm) ≈ 32π / 3 - - # end - - # @testset "Radial rescale" begin - # @testset "RescaledGrid" begin - # rmax = 2.0 - # N = 10 - # rgrid = SimpleG.Uniform([0.0, rmax^2], N; gpbound=[rmax^2 / 2N, rmax^2 * (1 - 1 / 2N)]) - # rrg = radial_rescale(grid=rgrid, DIM=2) - # println(rrg.grid) - # bound = [-π, π] - # theta = SimpleGrid.Uniform(bound, N; gpbound=[π * (1 / N - 1), π * (1 - 1 / N)], isperiodic=true) +pm2 = BZMeshes.CompositePolarMesh(dispersion=dispersion, + anglemesh=theta, cell=br, basegridtype=:uniform, + kmax=π * sqrt(2.4), Nloggrid=10, Nbasegrid=6, minterval=0.01T, N=6) - # # for (i, p) in enumerate(theta) - # # println(volume(theta, i)) - # # end - # DIM = 2 - # lattice = Matrix([1.0 0; 0 1]') - # br = BZMeshes.Cell(lattice=lattice) - # bzmesh = PolarMesh(br, CompositeMesh(theta, [rrg for i in 1:length(theta)])) - # for (i, p) in enumerate(bzmesh) - # println(p, AbstractMeshes.volume(bzmesh, i)) - # end - # end +println("N=$(length(pm)), ", int_mesh(pm)) +println("N=$(length(pm2)), ", int_mesh(pm2)) - # end -end \ No newline at end of file +Nuni = 1000 +um = BZMeshes.UniformBZMesh(cell=br, size=(Nuni, Nuni)) +um2 = BZMeshes.UniformBZMesh(cell=br, size=(2Nuni, 2Nuni)) +println("N=$(length(um)), ", int_mesh(um)) +println("N=$(length(um2)), ", int_mesh(um2)) \ No newline at end of file diff --git a/src/BaseMesh.jl b/src/BaseMesh.jl index efebbe2..8c769a7 100644 --- a/src/BaseMesh.jl +++ b/src/BaseMesh.jl @@ -91,6 +91,16 @@ Here we assume periodic boundary condition so for all case it's the same. AbstractMeshes.volume(mesh::AbstractUniformMesh) = cell_volume(mesh) AbstractMeshes.volume(mesh::AbstractUniformMesh, i) = cell_volume(mesh) / length(mesh) +AbstractMeshes.interp(data, mesh::AbstractUniformMesh, x) = data[locate(mesh, x)] + +function AbstractMeshes.integrate(data, mesh::AbstractUniformMesh) + result = 0.0 + for i in 1:length(mesh) + result += data[i] * volume(mesh, i) + end + return result +end + struct HasLattice <: AbstractMeshes.LatticeStyle end # has mesh.lattice and mesh.inv_lattice AbstractMeshes.lattice_vector(::HasLattice, mesh) = mesh.lattice AbstractMeshes.inv_lattice_vector(::HasLattice, mesh) = mesh.inv_lattice diff --git a/src/PolarMeshes.jl b/src/PolarMeshes.jl index 586966d..67c8c40 100644 --- a/src/PolarMeshes.jl +++ b/src/PolarMeshes.jl @@ -162,7 +162,7 @@ function AbstractMeshes.integrate(data, mesh::PolarMesh{T,3,MT}) where {T,MT} weighteddata = similar(data) for pi in 1:length(mesh) r, θ, ϕ = _extract(mesh[AngularCoords, pi]) - weighteddata[pi] = r^2 * sin(θ) * data[pi] + weighteddata[pi] = r^2 * cos(θ) * data[pi] end return integrate(weighteddata, mesh.mesh) end @@ -270,7 +270,6 @@ function PolarMesh(; dispersion, anglemesh, cell, kmax, DIM = size(cell.lattice, 1) bound = [0.0, kmax] - println(typeof(dispersion)) grids = kF_densed_kgrids(; dispersion=dispersion, anglemesh=anglemesh, bound=bound, DIM=DIM, kwargs...) cm = ProdMesh(grids, anglemesh) pm = PolarMesh(cell, cm) @@ -282,7 +281,6 @@ function CompositePolarMesh(; dispersion, anglemesh, cell, kmax, N, DIM = size(cell.lattice, 1) bound = [0.0, kmax] - println(typeof(dispersion)) grids = kF_densed_kgrids(; dispersion=dispersion, anglemesh=anglemesh, bound=bound, DIM=DIM, kwargs...) prm = ProdMesh(grids, anglemesh) cm = CompositeMesh(prm, N) diff --git a/test/BZMeshes.jl b/test/BZMeshes.jl index 97015fc..36e5580 100644 --- a/test/BZMeshes.jl +++ b/test/BZMeshes.jl @@ -231,6 +231,44 @@ end + @testset "3D CompositePolarMesh" begin + dispersion(k) = dot(k, k) - 1.0 + + N = 6 + bound = [-π, π] + phi = SimpleGrid.Uniform(bound, N) + + N = 4 + bound = [-π / 2, π / 2] + theta = SimpleGrid.Uniform(bound, N) + + am = ProdMesh([theta for i in 1:length(phi)], phi) + + DIM = 3 + lattice = Matrix([1.0 1.0 0; 1 0 1; 0 1 1]') + br = BZMeshes.Cell(lattice=lattice) + + pm = CompositePolarMesh(dispersion=dispersion, anglemesh=am, cell=br, kmax=2.0, N=4) + @test AbstractMeshes.volume(pm) ≈ 32π / 3 + + data = zeros(size(pm)) + for (pi, p) in enumerate(pm) + data[pi] = dispersion(p) + end + + testN = 10 + for i in 1:testN + r, ϕ, θ = rand(rng) * 2.0, (rand(rng) * 2 - 1) * π, (rand(rng) * 2 - 1) * π / 2 + p = Spherical(r, θ, ϕ) + x = BZMeshes._cartesianize(p) + @test isapprox(AbstractMeshes.interp(data, pm, x), dispersion(x), rtol=1e-4) + @test isapprox(AbstractMeshes.interp(data, pm, p), dispersion(x), rtol=1e-4) + end + + @test isapprox(AbstractMeshes.integrate(data, pm), 4π * 56 / 15, rtol=1e-4) + + end + @testset "Radial rescale" begin @testset "RescaledGrid" begin rmax = 2.0