From 5aed3962c349ac1d0bd03a8664a1be04a692abe3 Mon Sep 17 00:00:00 2001
From: Pablo San-Jose <lekand@gmail.com>
Date: Tue, 16 Jun 2020 14:02:18 +0200
Subject: [PATCH 1/6] draft of linecut API

fix tests

MeshSpec

fix anyoldmatrix

fix tests

degeneracy fix + closed kwarg

fix finite-difference Parametric + BZpoints

fix test

docstring

BZpoint typo
---
 src/Quantica.jl            |   3 +-
 src/bandstructure.jl       | 110 ++++++++++++-------
 src/diagonalizer.jl        |  61 ++++++-----
 src/mesh.jl                | 209 ++++++++++++++++++++++++++++++++++---
 src/parametric.jl          |   2 +-
 src/plot_vegalite.jl       |   2 +-
 src/tools.jl               |  82 ++-------------
 test/test_bandstructure.jl |  48 +++++----
 test/test_mesh.jl          |   4 +-
 9 files changed, 341 insertions(+), 180 deletions(-)

diff --git a/src/Quantica.jl b/src/Quantica.jl
index aba811a6..8b19b112 100644
--- a/src/Quantica.jl
+++ b/src/Quantica.jl
@@ -21,7 +21,8 @@ export sublat, bravais, lattice, dims, sites, supercell, unitcell,
        hopping, onsite, @onsite!, @hopping!, parameters,
        hamiltonian, parametric, bloch, bloch!, optimize!, similarmatrix,
        flatten, wrap, transform!, combine,
-       spectrum, bandstructure, marchingmesh, defaultmethod, bands, vertices,
+       spectrum, bandstructure, marchingmesh, linearmesh, buildmesh, buildlift,
+       defaultmethod, bands, vertices,
        energies, states,
        momentaKPM, dosKPM, averageKPM, densityKPM, bandrangeKPM
 
diff --git a/src/bandstructure.jl b/src/bandstructure.jl
index 63532f6a..2fa29fb4 100644
--- a/src/bandstructure.jl
+++ b/src/bandstructure.jl
@@ -139,48 +139,59 @@ end
 # bandstructure
 #######################################################################
 """
-bandstructure(h::Hamiltonian, mesh::Mesh; cut = missing, minprojection = 0.5, method = defaultmethod(h), transform = missing)
+    bandstructure(h::Hamiltonian; resolution = 13, kw...)
 
-Compute the bandstructure of Bloch Hamiltonian `bloch(h, ϕs)`, with `ϕs` evaluated on the
-vertices of `mesh`.
+Compute the bandstructure of `h` on a mesh over `h`'s full Brillouin zone, with `resolution`
+points along each axis, spanning the interval [-π,π].
 
-    bandstructure(h::Hamiltonian; resolution = 13, kw...)
+    bandstructure(h::Hamiltonian, spec::MeshSpec; lift = missing, kw...)
+
+Call `bandstructure(h, mesh; lift = lift, kw...)` with `mesh = buildmesh(spec, h)` and `lift
+= buildlift(spec, h)` if not provided. See `MeshSpec` for available mesh specs. If the `lift
+= missing` and the dimensions of the mesh do not match the Hamiltonian's, a `lift` function
+is used that lifts the mesh onto the dimensions `h` by appending vertex coordinates with
+zeros.
 
-Same as above with a uniform `mesh` of marching tetrahedra (generalized to the lattice
-dimensions of the Hamiltonian), with points `range(-π, π, length = resolution)` along each
-Bravais axis. Note that `resolution` denotes the number of points along each Bloch axis,
-including endpoints (can be a tuple for axis-dependent points).
+    bandstructure(h::Hamiltonian, mesh::Mesh; lift = missing, kw...)
 
-    bandstructure(ph::ParametricHamiltonian, mesh; kw...)
+Compute the bandstructure `bandstructure(h, mesh; kw...)` of Bloch Hamiltonian `bloch(h,
+ϕ)`, with `ϕ = v` taken on each vertex `v` of `mesh` (or `ϕ = lift(v...)` if a `lift`
+function is provided).
+
+    bandstructure(ph::ParametricHamiltonian, ...; kw...)
 
 Compute the bandstructure of a `ph` with `i` parameters (see `parameters(ph)`), where `mesh`
 is interpreted as a discretization of parameter space ⊗ Brillouin zone, so that each vertex
 reads `v = (p₁,..., pᵢ, ϕ₁,..., ϕⱼ)`, with `p` the values assigned to `parameters(ph)` and
-`ϕ` the Bloch phases.
+`ϕᵢ` the Bloch phases.
 
     bandstructure(matrixf::Function, mesh::Mesh; kw...)
 
-Compute the bandstructure of the Hamiltonian matrix `m = matrixf(ϕs)`, with `ϕs` evaluated
-on the vertices of `mesh`. No `cut` option is allowed in this case. If needed, include it in
-the definition of `matrixf`. Note that `ϕs` is either a `Tuple` or an `SVector`, but it is
+Compute the bandstructure of the Hamiltonian matrix `m = matrixf(ϕ)`, with `ϕ` evaluated on
+the vertices `v` of `mesh`. No `lift` option is allowed in this case. If needed, include it
+in the definition of `matrixf`. Note that `ϕ` is either a `Tuple` or an `SVector`, but it is
 not splatted into `matrixf`, i.e. `matrixf(x) = ...` or `matrixf(x, y) = ...` will not work,
 use `matrixf((x,)) = ...` or `matrixf((x, y)) = ...` instead.
 
 # Options
 
-The option `cut`, when not `missing`, is a function `cut = (mesh_coordinates...) -> φs`,
-where `φs` are Bloch phases if sampling a `h::Hamiltonian`, or `(params..., φs...)` if
-sampling a `ph::ParametricHamiltonian`, and `params` are values for `parameters(ph)`. This
+The default options are
+
+    (lift = missing, minprojection = 0.5, method = defaultmethod(h), transform = missing)
+
+`lift`: when not `missing`, `lift` is a function `lift = (vs...) -> ϕ`, where `vs` are the
+coordinates of a mesh vertex and `ϕ` are Bloch phases if sampling a `h::Hamiltonian`, or
+`(paramsⱼ..., ϕᵢ...)` if sampling a `ph::ParametricHamiltonian`, and `params` are values for
+`parameters(ph)`. It represents a mapping from a mesh and a Brillouin/parameter space. This
 allows to compute a bandstructure along a cut in the Brillouin zone/parameter space, see
 below for examples.
 
-The option `minprojection` determines the minimum projection between eigenstates to connect
+`minprojection`: determines the minimum projection between eigenstates to connect
 them into a common subband.
 
-The option `method` is chosen automatically if unspecified, and
-can be one of the following
+`method`: it is chosen automatically if unspecified, and can be one of the following
 
-    method                    diagonalization function
+    method                     diagonalization function
     --------------------------------------------------------------
     LinearAlgebraPackage()     LinearAlgebra.eigen!
     ArpackPackage()            Arpack.eigs (must be `using Arpack`)
@@ -189,46 +200,67 @@ Options passed to the `method` will be forwarded to the diagonalization function
 `method = ArpackPackage(nev = 8, sigma = 1im)` will use `Arpack.eigs(matrix; nev = 8,
 sigma = 1im)` to compute the bandstructure.
 
-The option `transform = ε -> f(ε)` allows to transform eigenvalues by `f` in the returned
+`transform`: the option `transform = ε -> f(ε)` allows to transform eigenvalues by `f` in the returned
 bandstructure (useful for performing shifts or other postprocessing).
 
 # Examples
 ```
-julia> h = LatticePresets.honeycomb() |> unitcell(3) |> hamiltonian(hopping(-1, range = 1/√3));
+julia> h = LatticePresets.honeycomb() |> hamiltonian(hopping(-1, range = 1/√3)) |> unitcell(3);
 
 julia> bandstructure(h; resolution = 25, method = LinearAlgebraPackage())
-Bandstructure: bands for a 2D hamiltonian
+Bandstructure{2}: collection of 2D bands
   Bands        : 8
   Element type : scalar (Complex{Float64})
   Mesh{2}: mesh of a 2-dimensional manifold
     Vertices   : 625
-    Edges      : 3552
+    Edges      : 1776
 
-julia> bandstructure(h, marchingmesh(range(0, 2π, length = 11)); cut = φ -> (φ, 0))
-       ## Computes a `Γ-X` cut of the bandstructure
+julia> bandstructure(h, linearmesh(:Γ, :X, :Y, :Γ))
+Bandstructure{1}: collection of 1D bands
+  Bands        : 17
+  Element type : scalar (Complex{Float64})
+  Mesh{1}: mesh of a 1-dimensional manifold
+    Vertices   : 37
+    Edges      : 36
+
+julia> bandstructure(h, marchingmesh((0, 2π); resolution = 25); lift = φ -> (φ, 0))
+       # Equivalent to bandstructure(h, linearmesh(:Γ, :X; resolution = 11))
 Bandstructure{1}: collection of 1D bands
   Bands        : 18
   Element type : scalar (Complex{Float64})
   Mesh{1}: mesh of a 1-dimensional manifold
-    Vertices   : 11
-    Edges      : 10
+    Vertices   : 25
+    Edges      : 24
 ```
 
 # See also
-    `marchingmesh`
+    `marchingmesh`, `linearmesh`
 """
-function bandstructure(h::Hamiltonian{<:Any,L,M}; resolution = 13, kw...) where {L,M}
-    mesh = marchingmesh(filltuple(range(-π, π, length = resolution), Val(L))...)
-    return bandstructure(h, mesh; kw...)
+function bandstructure(h::Hamiltonian; resolution = 13, kw...)
+    meshspec = marchingmesh(filltuple((-π, π), Val(latdim(h)))...; resolution = resolution)
+    return bandstructure(h, meshspec; kw...)
+end
+
+function bandstructure(h::Union{Hamiltonian,ParametricHamiltonian}, spec::MeshSpec; lift = missing, kw...)
+    br = bravais_parameters(h)
+    mesh = buildmesh(spec, br)
+    lift´ = lift === missing ? buildlift(spec, br) : lift
+    return bandstructure(h, mesh; lift = lift´, kw...)
 end
 
+bravais_parameters(h::Hamiltonian) = bravais(h)
+bravais_parameters(ph::ParametricHamiltonian{P}) where {P} = _blockdiag(SMatrix{P,P}(I), bravais(ph))
+
+bandstructure(h::Function, spec::MeshSpec; kw...) =
+    bandstructure(h, buildmesh(spec); kw...)
+
 function bandstructure(h::Union{Hamiltonian,ParametricHamiltonian}, mesh::Mesh;
-                       method = defaultmethod(h), cut = missing, transform = missing, kw...)
+                       method = defaultmethod(h), lift = missing, transform = missing, kw...)
     # ishermitian(h) || throw(ArgumentError("Hamiltonian must be hermitian"))
     matrix = similarmatrix(h, method)
-    codiag = codiagonalizer(h, matrix, mesh, cut)
+    codiag = codiagonalizer(h, matrix, mesh, lift)
     d = DiagonalizeHelper(method, codiag; kw...)
-    matrixf(φs) = bloch!(matrix, h, applycut(cut, φs))
+    matrixf(φs) = bloch!(matrix, h, applylift(lift, φs))
     b = _bandstructure(matrixf, matrix, mesh, d)
     transform === missing || transform!(transform, b)
     return b
@@ -248,9 +280,9 @@ end
 
 _samplematrix(matrixf, mesh) = matrixf(Tuple(first(vertices(mesh))))
 
-@inline applycut(cut::Missing, ϕs) = toSVector(ϕs)
+@inline applylift(lift::Missing, ϕs) = toSVector(ϕs)
 
-@inline applycut(cut::Function, ϕs) = toSVector(cut(ϕs...))
+@inline applylift(lift::Function, ϕs) = toSVector(lift(ϕs...))
 
 function _bandstructure(matrixf::Function, matrix´::AbstractMatrix{M}, mesh::MD, d::DiagonalizeHelper) where {M,D,T,MD<:Mesh{D,T}}
     nϵ = 0                           # Temporary, to be reassigned
@@ -267,7 +299,7 @@ function _bandstructure(matrixf::Function, matrix´::AbstractMatrix{M}, mesh::MD
         matrix = matrixf(Tuple(ϕs))
         # (ϵk, ψk) = diagonalize(Hermitian(matrix), d)  ## This is faster (!)
         (ϵk, ψk) = diagonalize(matrix, d.method)
-        resolve_degeneracies!(ϵk, ψk, ϕs, matrix, d.codiag)
+        resolve_degeneracies!(ϵk, ψk, ϕs, d.codiag)
         if n == 1  # With first vertex can now know the number of eigenvalues... Reassign
             nϵ = size(ϵk, 1)
             ϵks = Matrix{T}(undef, nϵ, nk)
@@ -355,4 +387,4 @@ function findmostparallel(ψks::Array{M,3}, destk, srcb, srck) where {M}
         end
     end
     return maxproj, destb
-end
+end
\ No newline at end of file
diff --git a/src/diagonalizer.jl b/src/diagonalizer.jl
index 5474edfd..c9e5bc33 100644
--- a/src/diagonalizer.jl
+++ b/src/diagonalizer.jl
@@ -123,7 +123,7 @@ end
 # resolve_degeneracies
 #######################################################################
 # Tries to make states continuous at crossings. Here, ϵ needs to be sorted
-function resolve_degeneracies!(ϵ, ψ, ϕs, matrix, codiag)
+function resolve_degeneracies!(ϵ, ψ, ϕs, codiag)
     issorted(ϵ, by = real) || sorteigs!(codiag.perm, ϵ, ψ)
     hasapproxruns(ϵ, codiag.degtol) || return ϵ, ψ
     ranges, ranges´ = codiag.rangesA, codiag.rangesB
@@ -174,22 +174,22 @@ struct Codiagonalizer{T,F<:Function}
     perm::Vector{Int}               # Prealloc for sortperm!
 end
 
-function codiagonalizer(h::Union{Hamiltonian,ParametricHamiltonian}, matrix, mesh, cut; kw...)
+function codiagonalizer(h::Union{Hamiltonian,ParametricHamiltonian}, matrix, mesh, lift; kw...)
     veldirs = velocitydirections(parent(h); kw...)
-    meshdirs = meshdirections(mesh; kw...)
+    veldirs_with_params = padparams.(veldirs, Ref(h))
     nv = length(veldirs)
-    nm = length(meshdirs)
-    matrixindices = 1:(nv + nm + 1)
+    matrixindices = 1:(nv + nv + 1)
     degtol = sqrt(eps(real(blockeltype(h))))
-    delta = meshdelta(mesh)
+    delta = meshdelta(mesh, lift)
     delta = iszero(delta) ? degtol : delta
+    aom = anyoldmatrix(matrix)
     cmatrixf(meshϕs, n) =
         if n <= nv
-            bloch!(matrix, h, applycut(cut, meshϕs), dn -> im * veldirs[n]' * dn)
-        elseif n - nv <= nm # resort to finite differences
-            bloch!(matrix, h, applycut(cut, meshϕs + delta * meshdirs[n - nv]))
-        else # use a fixed random matrix
-            randomfill!(matrix)
+            bloch!(matrix, h, applylift(lift, meshϕs), dn -> im * veldirs[n]' * dn)
+        elseif n - nv <= nv # resort to finite differences
+            bloch!(matrix, h, applylift(lift, meshϕs) + delta * veldirs_with_params[n - nv])
+        else # use a fixed arbitrary matrix
+            aom
         end
     return Codiagonalizer(cmatrixf, matrixindices, degtol, UnitRange{Int}[], UnitRange{Int}[], Int[])
 end
@@ -201,18 +201,23 @@ function codiagonalizer(matrixf::Function, matrix::AbstractMatrix{T}, mesh; kw..
     degtol = sqrt(eps(real(eltype(T))))
     delta = meshdelta(mesh)
     delta = iszero(delta) ? degtol : delta
+    aom = anyoldmatrix(matrix)
     cmatrixf(meshϕs, n) =
         if n <= nm # finite differences
             matrixf(meshϕs + delta * meshdirs[n])
-        else # use a fixed random matrix
-            randomfill!(matrix)
+        else # use a fixed arbitrary matrix
+            aom
         end
     return Codiagonalizer(cmatrixf, matrixindices, degtol, UnitRange{Int}[], UnitRange{Int}[], Int[])
 end
 
 velocitydirections(::Hamiltonian{LA,L}; kw...) where {LA,L} = _directions(Val(L); kw...)
+
 meshdirections(::Mesh{L}; kw...) where {L} = _directions(Val(L); kw...)
 
+padparams(v, ::Hamiltonian) = v
+padparams(v::SVector{L,T}, ::ParametricHamiltonian{P}) where {L,T,P} = vcat(zero(SVector{P,T}), v)
+
 function _directions(::Val{L}; direlements = 0:1, onlypositive = true) where {L}
     directions = vec(SVector{L,Int}.(Iterators.product(ntuple(_ -> direlements, Val(L))...)))
     onlypositive && filter!(ispositive, directions)
@@ -221,20 +226,22 @@ function _directions(::Val{L}; direlements = 0:1, onlypositive = true) where {L}
     return directions
 end
 
-meshdelta(mesh::Mesh{<:Any,T}) where {T} = T(0.1) * first(minmax_edge_length(mesh))
+meshdelta(mesh::Mesh{<:Any,T}, lift = missing) where {T} =
+    T(0.1) * norm(applylift(lift, first(minmax_edge(mesh))))
 
-function randomfill!(matrix::SparseMatrixCSC, seed = 1234)
-    Random.seed!(seed)
-    rand!(nonzeros(matrix))  ## CAREFUL: this will be non-hermitian
-    return matrix
+# function anyoldmatrix(matrix::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti}
+#     n = size(matrix, 1)
+#     ri = one(Ti):Ti(n)
+#     rv = Tv.(im .* (1:n))
+#     s = sparse(ri, ri, rv, n, n)
+#     return s
+# end
+
+function anyoldmatrix(matrix::SparseMatrixCSC, rng = MersenneTwister(1))
+    s = copy(matrix)
+    rand!(rng, nonzeros(s))
+    return s
 end
 
-function randomfill!(matrix::AbstractArray{T}, seed = 1234) where {T}
-    Random.seed!(seed)
-    fill!(matrix, zero(T))
-    for i in 1:minimum(size(matrix))
-        r = rand(T)
-        @inbounds matrix[i, i] = r + r'
-    end
-    return matrix
-end
\ No newline at end of file
+# anyoldmatrix(m::M) where {T,M<:DenseArray{T}} = M(Diagonal(T.(im .* (1:size(m,1)))))
+anyoldmatrix(m::DenseArray, rng = MersenneTwister(1)) = rand!(rng, copy(m))
\ No newline at end of file
diff --git a/src/mesh.jl b/src/mesh.jl
index b04ef8ca..6491ec65 100644
--- a/src/mesh.jl
+++ b/src/mesh.jl
@@ -36,19 +36,21 @@ edgedest(m::Mesh, edge) = rowvals(m.adjmat)[edge]
 edgevertices(m::Mesh) =
     ((vsrc, m.vertices[edgedest(m, edge)]) for (i, vsrc) in enumerate(m.vertices) for edge in edges(m, i))
 
-function minmax_edge_length(m::Mesh{D,T}) where {D,T<:Real}
+function minmax_edge(m::Mesh{D,T}) where {D,T<:Real}
     minlen2 = typemax(T)
     maxlen2 = zero(T)
     verts = vertices(m)
+    minedge = zero(first(verts))
+    maxedge = zero(first(verts))
     for src in eachindex(verts), edge in edges(m, src)
         dest = edgedest(m, edge)
         dest > src || continue # Need only directed graph
         vec = verts[dest] - verts[src]
         norm2 = vec' * vec
-        norm2 < minlen2 && (minlen2 = norm2)
-        norm2 > maxlen2 && (maxlen2 = norm2)
+        norm2 < minlen2 && (minlen2 = norm2; minedge = vec)
+        norm2 > maxlen2 && (maxlen2 = norm2; maxedge = vec)
     end
-    return sqrt(minlen2), sqrt(maxlen2)
+    return minedge, maxedge
 end
 
 ######################################################################
@@ -110,24 +112,70 @@ elementvolume(mat::SMatrix{N,N}) where {N} = det(mat)
 switchlast(s::NTuple{N,T}) where {N,T} = ntuple(i -> i < N - 1 ? s[i] : s[2N - i - 1] , Val(N))
 
 ######################################################################
-# Special meshes
+# Mesh specifications
 ######################################################################
 """
-    marchingmesh(ranges::Vararg{AbstractRange,L}; axes = 1.0 * I)
+    MeshSpec
 
-Creates a L-dimensional marching-tetrahedra `Mesh` over a parallelepiped with axes given by
-the columns of `axes`. The points along axis `i` are given by `ranges[i]`.
+Parent type of mesh specifications, which are currently `MarchingMeshSpec` (constructed with
+`marchingmesh`) and `LinearMeshSpec` (constructed with `linearmesh`).
 
-# External links
+# See also
+    `marchingmesh`, `linearmesh`, `buildmesh`
+"""
+abstract type MeshSpec{L} end
+
+Base.show(io::IO, spec::MeshSpec{L}) where {L} =
+    print(io, "MeshSpec{$L} : specifications for building a $(L)D mesh.")
+
+#fallback
+buildlift(::MeshSpec, bravais) = missing
+
+
+#######################################################################
+# MarchingMeshSpec
+#######################################################################
+struct MarchingMeshSpec{L,R,T,M<:NTuple{L,Tuple{Number,Number}}} <: MeshSpec{L}
+    minmaxaxes::M
+    axes::SMatrix{L,L,T}
+    resolution::R
+end
+
+"""
+    marchingmesh(minmaxaxes...; axes = 1.0 * I, resolution = 13)
+
+Create a `MeshSpec` for a L-dimensional marching-tetrahedra `Mesh` over a parallelepiped
+with axes given by the columns of `axes`. The points along axis `i` are distributed between
+`first(minmaxaxes[i])` and `last(minmaxaxes[i])`. The number of points on each axis is given
+by `resolution`, or `resolution[i]` if several are given.
+
+# Examples
 
+```jldoctest
+julia> buildmesh(marchingmesh((-π, π), (0,2π); resolution = 25))
+Mesh{2}: mesh of a 2-dimensional manifold
+  Vertices   : 625
+  Edges      : 1776
+
+julia> buildmesh(marchingmesh((-π, π), (0,2π); resolution = (10,10)))
+Mesh{2}: mesh of a 2-dimensional manifold
+  Vertices   : 100
+  Edges      : 261
+```
+
+# See also
+    `linearmesh`, `buildmesh`
+
+# External links
 - Marching tetrahedra (https://en.wikipedia.org/wiki/Marching_tetrahedra) in Wikipedia
 """
-marchingmesh(ranges::Vararg{AbstractRange,L}; axes = 1.0 * I) where {L} =
-    _marchingmesh(ranges, SMatrix{L,L}(axes))
+marchingmesh(minmaxaxes::Vararg{Tuple,L}; axes = 1.0 * I, resolution = 13) where {L} =
+    MarchingMeshSpec(minmaxaxes, SMatrix{L,L}(axes), resolution)
 
-marchingmesh(; kw...) = throw(ArgumentError("Need a finite number of ranges to build a mesh"))
+marchingmesh(; kw...) = throw(ArgumentError("Need a finite number of ranges to define a marching mesh"))
 
-function _marchingmesh(ranges::NTuple{D,AbstractRange}, axes::SMatrix{D,D}) where {D}
+function buildmesh(m::MarchingMeshSpec{D}, bravais = missing) where {D}
+    ranges = ((b, r)->range(b...; length = r)).(m.minmaxaxes, m.resolution)
     npoints = length.(ranges)
     cs = CartesianIndices(ntuple(n -> 1:npoints[n], Val(D)))
     ls = LinearIndices(cs)
@@ -140,8 +188,8 @@ function _marchingmesh(ranges::NTuple{D,AbstractRange}, axes::SMatrix{D,D}) wher
             ntuple(i -> CartesianIndex(ntuple(j -> i == j ? 1 : 0, Val(D))), Val(D)))
     utets = [cumsum(pushfirst!(perm, zero(CartesianIndex{D}))) for perm in perms]
 
-    # We don't use generators because their non-inferreble eltype causes problems later
-    verts = [axes * SVector(getindex.(ranges, Tuple(c))) for c in cs]
+    # We don't use generators because their non-inferreble eltype causes problems elsewhere
+    verts = [m.axes * SVector(getindex.(ranges, Tuple(c))) for c in cs]
 
     s = SparseMatrixBuilder{Bool}(length(cs), length(cs))
     for c in cs
@@ -157,3 +205,134 @@ function _marchingmesh(ranges::NTuple{D,AbstractRange}, axes::SMatrix{D,D}) wher
 
     return Mesh(verts, adjmat)
 end
+
+buildlift(::MarchingMeshSpec{L}, bravais::SMatrix{E,L}) where {E,L} = missing
+buildlift(::MarchingMeshSpec{L´}, bravais::SMatrix{E,L}) where {E,L,L´} =
+    (pt...) -> padright(pt, Val(L))
+
+#######################################################################
+# LinearMeshSpec
+#######################################################################
+struct LinearMeshSpec{N,L,T,R} <: MeshSpec{1}
+    vertices::SVector{N,SVector{L,T}}
+    samelength::Bool
+    closed::Bool
+    resolution::R
+end
+
+"""
+    linearmesh(nodes...; resolution = 13, samelength = false, closed = false)
+
+Create a `MeshSpec` for a one-dimensional `Mesh` connecting the `nodes` with straight
+segments, each containing `resolution` points (endpoints included). If a different
+resolution for each of the `N` segments is required, use `resolution::NTuple{N,Int}`.
+If `samelength` each segment has equal length in mesh coordinates. If `closed` the last node
+is connected to the first node.
+
+# Examples
+
+```jldoctest
+julia> buildmesh(linearmesh(:Γ, :K, :M, :Γ; resolution = (101, 30, 30)))
+Mesh{1}: mesh of a 1-dimensional manifold
+  Vertices   : 159
+  Edges      : 158
+```
+
+# See also
+    `marchingmesh`, `buildmesh`
+"""
+linearmesh(nodes...; resolution = 13, samelength::Bool = false, closed::Bool = false) =
+    LinearMeshSpec(sanitize_BZpts(nodes), samelength, closed, resolution)
+
+function sanitize_BZpts(pts)
+    pts´ = parse_BZpoint.(pts)
+    dim = maximum(length.(pts´))
+    pts´´ = SVector(padright.(pts´, Val(dim)))
+    return pts´´
+end
+
+parse_BZpoint(p::SVector) = float.(p)
+parse_BZpoint(p::Tuple) = SVector(float.(p))
+function parse_BZpoint(p::Symbol)
+    pt = get(BZpoints, p, missing)
+    pt === missing && throw(ArgumentError("Unknown Brillouin zone point $p, use one of $(keys(BZpoints))"))
+    return SVector(float.(pt))
+end
+const BZpoints =
+    ( Γ  = (0,)
+    , X  = (pi,)
+    , Y  = (0, pi)
+    , Z  = (0, 0, pi)
+    , K  = (2pi/3, -2pi/3)
+    , Kp = (4pi/3, 2pi/3)
+    , M  = (pi, 0)
+    )
+
+linearmesh_nodes(l, br) = cumsum(SVector(0, segment_lengths(l, br)...))
+
+function segment_lengths(s::LinearMeshSpec{N,LS,TS}, br::SMatrix{E,LB,TB}) where {TS,TB,N,E,LS,LB}
+    T = promote_type(TS, TB)
+    verts = padright.(s.vertices, Val(LB))
+    dϕs = ntuple(i -> verts[i + 1] - verts[i], Val(N-1))
+    if s.samelength
+        ls = filltuple(T(1/(N-1)), Val(N-1))
+    else
+        ibr = pinverse(br)'
+        ls = (dϕ -> norm(ibr * dϕ)).(dϕs)
+        ls = ls ./ sum(ls)
+    end
+    return ls
+end
+
+function idx_to_node(s, br)
+    nodes = SVector.(linearmesh_nodes(s, br))
+    nmax = length(nodes)
+    nodefunc = nvec -> begin
+        n = only(nvec)
+        node = if n >= nmax
+            nodes[nmax]
+        else
+            nc = max(n, 1)
+            i = Int(floor(nc))
+            nodes[i] + rem(nc, 1) * (nodes[i+1] - nodes[i])
+        end
+        return node
+    end
+    return nodefunc
+end
+
+buildmesh(s::LinearMeshSpec{N,L,T}) where {N,L,T} = buildmesh(s, SMatrix{L,L,T}(I))
+
+function buildmesh(s::LinearMeshSpec{N}, br) where {N}
+    ranges = ((i, r) -> range(i, i+1, length = r)).(ntuple(identity, Val(N-1)), s.resolution)
+    verts = SVector.(first(ranges))
+    for r in Base.tail(ranges)
+        pop!(verts)
+        append!(verts, SVector.(r))
+    end
+    s.closed && pop!(verts)
+    nv = length(verts)
+    nodefunc = idx_to_node(s, br)
+    verts .= nodefunc.(verts)
+    adjmat = sparse(vcat(1:nv-1, 2:nv), vcat(2:nv, 1:nv-1), true, nv, nv)
+    s.closed && (adjmat[end, 1] = adjmat[1, end] = true)
+    return Mesh(verts, adjmat)
+end
+
+function buildlift(s::LinearMeshSpec, br::SMatrix{E,L}) where {E,L}
+    ls = segment_lengths(s, br)
+    nodes = linearmesh_nodes(s, br)
+    verts = padright.(s.vertices, Val(L))
+    l = sum(ls)
+    liftfunc = x -> begin
+        xc = clamp(only(x), 0, l)
+        for (i, node) in enumerate(nodes)
+            if node > xc
+                p = verts[i-1] + (xc - nodes[i-1])/ls[i-1] * (verts[i]-verts[i-1])
+                return p
+            end
+        end
+        return last(verts)
+    end
+    return liftfunc
+end
\ No newline at end of file
diff --git a/src/parametric.jl b/src/parametric.jl
index 0d3445ee..f15df362 100644
--- a/src/parametric.jl
+++ b/src/parametric.jl
@@ -192,7 +192,7 @@ matrixtype(ph::ParametricHamiltonian) = matrixtype(parent(ph))
 
 blockeltype(ph::ParametricHamiltonian) = blockeltype(parent(ph))
 
-bravais(ph::ParametricHamiltonian) = bravais(ph.hamiltonian.lattice)
+bravais(ph::ParametricHamiltonian) = bravais(ph.baseh.lattice)
 
 Base.parent(ph::ParametricHamiltonian) = ph.h
 
diff --git a/src/plot_vegalite.jl b/src/plot_vegalite.jl
index 282f130f..fcea79bb 100644
--- a/src/plot_vegalite.jl
+++ b/src/plot_vegalite.jl
@@ -16,7 +16,7 @@ Plots the the Hamiltonian lattice projected along `axes` using VegaLite.
     - `scaling`: `(scalex, scaley)` scalings for the x (Bloch phase) and y (energy) variables
     - `axes`: lattice axes to project onto the plot x-y plane
 """
-function VegaLite.vlplot(b::Bandstructure; labels = ("φ/2π", "ε"), scaling = (1/2π, 1), size = 640, points = false)
+function VegaLite.vlplot(b::Bandstructure; labels = ("φ", "ε"), scaling = (1, 1), size = 640, points = false)
     labelx, labely = labels
     table = bandtable(b, make_it_two(scaling))
     sizes = make_it_two(size)
diff --git a/src/tools.jl b/src/tools.jl
index e7f3475d..11057980 100644
--- a/src/tools.jl
+++ b/src/tools.jl
@@ -82,7 +82,7 @@ displayvectors(mat::SMatrix{E,L,<:AbstractFloat}; kw...) where {E,L} =
 displayvectors(mat::SMatrix{E,L,<:Integer}; kw...) where {E,L} =
     ntuple(l -> Tuple(mat[:,l]), Val(L))
 
-# pseudoinverse of s times an integer n, so that it is an integer matrix (for accuracy)
+# pseudoinverse of supercell s times an integer n, so that it is an integer matrix (for accuracy)
 pinvmultiple(s::SMatrix{L,0}) where {L} = (SMatrix{0,0,Int}(), 0)
 function pinvmultiple(s::SMatrix{L,L´}) where {L,L´}
     L < L´ && throw(DimensionMismatch("Supercell dimensions $(L´) cannot exceed lattice dimensions $L"))
@@ -94,6 +94,12 @@ function pinvmultiple(s::SMatrix{L,L´}) where {L,L´}
     return round.(Int, n * inv(qrfact.R) * qrfact.Q'), round(Int, n)
 end
 
+pinverse(m::SMatrix) = (f -> inv(f.R) * f.Q')(qr(m))
+
+_blockdiag(s1::SMatrix{M}, s2::SMatrix{N}) where {N,M} = hcat(
+    ntuple(j->vcat(s1[:,j], zero(s2[:,j])), Val(M))...,
+    ntuple(j->vcat(zero(s1[:,j]), s2[:,j]), Val(N))...)
+
 function isgrowing(vs::AbstractVector, i0 = 1)
     i0 > length(vs) && return true
     vprev = vs[i0]
@@ -299,80 +305,6 @@ end
 
 Base.factorial(n::Integer, k::Integer) = factorial(promote(n, k)...)
 
-# ######################################################################
-# # SparseMatrixIJV
-# ######################################################################
-
-# struct SparseMatrixIJV{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}
-#     I::Vector{Ti}
-#     J::Vector{Ti}
-#     V::Vector{Tv}
-#     m::Ti
-#     n::Ti
-#     klasttouch::Vector{Ti}
-#     csrrowptr::Vector{Ti}
-#     csrcolval::Vector{Ti}
-#     csrnzval::Vector{Tv}
-#     csccolptr::Vector{Ti}
-#     cscrowval::Vector{Ti}
-#     cscnzval::Vector{Tv}
-# end
-
-# SparseMatrixIJV{Tv}(m::Ti, n::Ti) where {Tv,Ti} = SparseMatrixIJV{Tv,Ti}(m,n)
-
-# function SparseMatrixIJV{Tv,Ti}(m::Integer, n::Integer; hintnnz = 0) where {Tv,Ti}
-#     I = Ti[]
-#     J = Ti[]
-#     V = Tv[]
-#     klasttouch = Vector{Ti}(undef, n)
-#     csrrowptr = Vector{Ti}(undef, m + 1)
-#     csrcolval = Vector{Ti}()
-#     csrnzval = Vector{Tv}()
-#     csccolptr = Vector{Ti}(undef, n + 1)
-#     cscrowval = Vector{Ti}()
-#     cscnzval = Vector{Tv}()
-
-#     if hintnnz > 0
-#         sizehint!(I, hintnnz)
-#         sizehint!(J, hintnnz)
-#         sizehint!(V, hintnnz)
-#         sizehint!(csrcolval, hintnnz)
-#         sizehint!(csrnzval, hintnnz)
-#         sizehint!(cscrowval, hintnnz)
-#         sizehint!(cscnzval, hintnnz)
-#     end
-
-#     return SparseMatrixIJV{Tv,Ti}(I, J, V, m, n, klasttouch, csrrowptr, csrcolval, csrnzval,
-#                                                              csccolptr, cscrowval, cscnzval)
-# end
-
-# Base.summary(::SparseMatrixIJV{Tv,Ti}) where {Tv,Ti} =
-#     "SparseMatrixIJV{$Tv,$Ti} : Sparse matrix builder using the IJV format"
-
-# function Base.show(io::IO, ::MIME"text/plain", s::SparseMatrixIJV)
-#     i = get(io, :indent, "")
-#     print(io, i, summary(s), "\n", "$i  Nonzero elements : $(length(s.I))")
-# end
-
-# function Base.push!(s::SparseMatrixIJV, (i, j, v))
-#     push!(s.I, i)
-#     push!(s.J, j)
-#     push!(s.V, v)
-#     return s
-# end
-
-# function SparseArrays.sparse(s::SparseMatrixIJV)
-#     numnz = length(s.I)
-#     resize!(s.csrcolval, numnz)
-#     resize!(s.csrnzval,  numnz)
-#     resize!(s.cscrowval, numnz)
-#     resize!(s.cscnzval,  numnz)
-#     return SparseArrays.sparse!(s.I, s.J, s.V, s.m, s.n, +, s.klasttouch,
-#         s.csrrowptr, s.csrcolval, s.csrnzval, s.csccolptr, s.cscrowval, s.cscnzval)
-# end
-
-# Base.size(s::SparseMatrixIJV) = (s.m, s.n)
-
 ############################################################################################
 ######## fast sparse copy #  Revise after #33589 is merged #################################
 ############################################################################################
diff --git a/test/test_bandstructure.jl b/test/test_bandstructure.jl
index 56a343d2..85a5eb63 100644
--- a/test/test_bandstructure.jl
+++ b/test/test_bandstructure.jl
@@ -6,19 +6,25 @@
     h = LatticePresets.honeycomb() |>
         hamiltonian(onsite(0.5, sublats = :A) + onsite(-0.5, sublats = :B) +
                     hopping(-1, range = 1/√3))
-    b = bandstructure(h, resolution = 13)
+    b = bandstructure(h, resolution = (13, 23))
     @test length(bands(b)) == 2
 
-    h = LatticePresets.square() |> hamiltonian(hopping(1)) |> unitcell(2)
-    b = bandstructure(h, resolution = 13)
-    @test length(bands(b)) == 4
+    h = LatticePresets.cubic() |> hamiltonian(hopping(1)) |> unitcell(2)
+    b = bandstructure(h, resolution = (5, 9, 5))
+    @test length(bands(b)) == 8
+
+    b = bandstructure(h, linearmesh(:Γ, :X, resolution = 4))
+    @test length(bands(b)) == 8
+
+    b = bandstructure(h, linearmesh(:Γ, :X, (0, π), :Z, :Γ, resolution = 4))
+    @test length(bands(b)) == 8
 end
 
 @testset "functional bandstructures" begin
     hc = LatticePresets.honeycomb() |> hamiltonian(hopping(-1, sublats = :A=>:B, plusadjoint = true))
     matrix = similarmatrix(hc, LinearAlgebraPackage())
     hf((x,)) = bloch!(matrix, hc, (x, -x))
-    mesh = marchingmesh(range(0, 1, length = 13))
+    mesh = marchingmesh((0, 1))
     b = bandstructure(hf, mesh)
     @test length(bands(b)) == 2
 
@@ -26,37 +32,41 @@ end
     hp2 = parametric(hc2, @hopping!((t; s) -> s*t))
     matrix2 = similarmatrix(hc2, LinearAlgebraPackage())
     hf2((s, x)) = bloch!(matrix2, hp2(s = s), (x, x))
-    mesh2 = marchingmesh(range(0, 1, length = 13), range(0, 1, length = 13))
+    mesh2 = marchingmesh((0, 1), (0, 1))
     b = bandstructure(hf2, mesh2)
     @test length(bands(b)) == 2
 end
 
-@testset "bandstructures cuts & transforms" begin
+@testset "bandstructures lifts & transforms" begin
     h = LatticePresets.honeycomb() |> hamiltonian(onsite(2) + hopping(-1, range = 1/√3))
-    mesh1D = marchingmesh(range(0, 2π, length = 13))
-    b = bandstructure(h, mesh1D, cut = φ -> (φ, -φ), transform = inv)
-    b´ = transform!(inv, bandstructure(h, mesh1D, cut = φ -> (φ, -φ)))
+    mesh1D = marchingmesh((0, 2π))
+    b = bandstructure(h, mesh1D, lift = φ -> (φ, -φ), transform = inv)
+    b´ = transform!(inv, bandstructure(h, mesh1D, lift = φ -> (φ, -φ)))
     @test length(bands(b)) == length(bands(b´)) == 2
     @test vertices(bands(b)[1]) == vertices(bands(b´)[1])
     h´ = unitcell(h)
     s1 = spectrum(h´, transform = inv)
     s2 = transform!(inv,spectrum(h´))
     @test energies(s1) == energies(s2)
+    # automatic lift from 2D to 3D
+    h = LatticePresets.cubic() |> hamiltonian(hopping(1)) |> unitcell(2)
+    b = bandstructure(h, marchingmesh((0, 2pi), (0, 2pi)))
+    @test length(bands(b)) == 8
 end
 
 @testset "parametric bandstructures" begin
     ph = LatticePresets.linear() |> hamiltonian(hopping(-I), orbitals = Val(2)) |> unitcell(2) |>
          parametric(@onsite!((o; k) -> o + k*I), @hopping!((t; k)-> t - k*I))
-    mesh2D = marchingmesh(range(0, 1, length = 13), range(0, 2π, length = 13))
-    mesh1D = marchingmesh(range(0, 2π, length = 13))
+    mesh2D = marchingmesh((0, 1), (0, 2π), resolution = 25)
+    mesh1D = marchingmesh((0, 2π), resolution = 25)
     b = bandstructure(ph, mesh2D)
+    @test length(bands(b)) == 2
+    b = bandstructure(ph, mesh2D, lift = (k, φ) -> (k + φ/2π, φ))
+    @test length(bands(b)) == 2
+    b = bandstructure(ph, mesh2D, lift = (k, φ) -> (k + φ/2π, φ))
+    @test length(bands(b)) == 2
+    b = bandstructure(ph, mesh1D, lift = k -> (k, 2π*k))
     @test length(bands(b)) == 4
-    b = bandstructure(ph, mesh2D, cut = (k, φ) -> (k + φ/2π, φ))
-    @test length(bands(b)) == 4
-    b = bandstructure(ph, mesh2D, cut = (k, φ) -> (k + φ/2π, φ))
-    @test length(bands(b)) == 4
-    b = bandstructure(ph, mesh1D, cut = k -> (k, 2π*k))
-    @test length(bands(b)) == 4
-    b = bandstructure(ph, mesh1D, cut = φ -> (φ, -φ))
+    b = bandstructure(ph, mesh1D, lift = φ -> (φ, -φ))
     @test length(bands(b)) == 4
 end
\ No newline at end of file
diff --git a/test/test_mesh.jl b/test/test_mesh.jl
index 48d95ef9..b43944df 100644
--- a/test/test_mesh.jl
+++ b/test/test_mesh.jl
@@ -1,6 +1,6 @@
 using Quantica: nvertices, nedges, nsimplices
 
 @test begin
-    mesh = marchingmesh(0:0.1:1, 1:0.1:2)
-    nvertices(mesh) == 121 && nedges(mesh) == 320 && nsimplices(mesh) == 200
+    mesh = buildmesh(marchingmesh((0,1), (0,2), resolution = (10, 20)))
+    nvertices(mesh) == 200 && nedges(mesh) == 541 && nsimplices(mesh) == 342
 end
\ No newline at end of file

From 332ed3298d8935dd515a2d68ee1dac7fb244704c Mon Sep 17 00:00:00 2001
From: Pablo San-Jose <lekand@gmail.com>
Date: Tue, 16 Jun 2020 18:21:52 +0200
Subject: [PATCH 2/6] resolution -> points

---
 src/bandstructure.jl       | 14 +++++++-------
 src/mesh.jl                | 32 ++++++++++++++++----------------
 test/test_bandstructure.jl | 14 +++++++-------
 test/test_mesh.jl          |  2 +-
 4 files changed, 31 insertions(+), 31 deletions(-)

diff --git a/src/bandstructure.jl b/src/bandstructure.jl
index 2fa29fb4..d8a7f3c9 100644
--- a/src/bandstructure.jl
+++ b/src/bandstructure.jl
@@ -139,9 +139,9 @@ end
 # bandstructure
 #######################################################################
 """
-    bandstructure(h::Hamiltonian; resolution = 13, kw...)
+    bandstructure(h::Hamiltonian; points = 13, kw...)
 
-Compute the bandstructure of `h` on a mesh over `h`'s full Brillouin zone, with `resolution`
+Compute the bandstructure of `h` on a mesh over `h`'s full Brillouin zone, with `points`
 points along each axis, spanning the interval [-π,π].
 
     bandstructure(h::Hamiltonian, spec::MeshSpec; lift = missing, kw...)
@@ -207,7 +207,7 @@ bandstructure (useful for performing shifts or other postprocessing).
 ```
 julia> h = LatticePresets.honeycomb() |> hamiltonian(hopping(-1, range = 1/√3)) |> unitcell(3);
 
-julia> bandstructure(h; resolution = 25, method = LinearAlgebraPackage())
+julia> bandstructure(h; points = 25, method = LinearAlgebraPackage())
 Bandstructure{2}: collection of 2D bands
   Bands        : 8
   Element type : scalar (Complex{Float64})
@@ -223,8 +223,8 @@ Bandstructure{1}: collection of 1D bands
     Vertices   : 37
     Edges      : 36
 
-julia> bandstructure(h, marchingmesh((0, 2π); resolution = 25); lift = φ -> (φ, 0))
-       # Equivalent to bandstructure(h, linearmesh(:Γ, :X; resolution = 11))
+julia> bandstructure(h, marchingmesh((0, 2π); points = 25); lift = φ -> (φ, 0))
+       # Equivalent to bandstructure(h, linearmesh(:Γ, :X; points = 11))
 Bandstructure{1}: collection of 1D bands
   Bands        : 18
   Element type : scalar (Complex{Float64})
@@ -236,8 +236,8 @@ Bandstructure{1}: collection of 1D bands
 # See also
     `marchingmesh`, `linearmesh`
 """
-function bandstructure(h::Hamiltonian; resolution = 13, kw...)
-    meshspec = marchingmesh(filltuple((-π, π), Val(latdim(h)))...; resolution = resolution)
+function bandstructure(h::Hamiltonian; points = 13, kw...)
+    meshspec = marchingmesh(filltuple((-π, π), Val(latdim(h)))...; points = points)
     return bandstructure(h, meshspec; kw...)
 end
 
diff --git a/src/mesh.jl b/src/mesh.jl
index 6491ec65..e6f0ceb3 100644
--- a/src/mesh.jl
+++ b/src/mesh.jl
@@ -138,26 +138,26 @@ buildlift(::MeshSpec, bravais) = missing
 struct MarchingMeshSpec{L,R,T,M<:NTuple{L,Tuple{Number,Number}}} <: MeshSpec{L}
     minmaxaxes::M
     axes::SMatrix{L,L,T}
-    resolution::R
+    points::R
 end
 
 """
-    marchingmesh(minmaxaxes...; axes = 1.0 * I, resolution = 13)
+    marchingmesh(minmaxaxes...; axes = 1.0 * I, points = 13)
 
 Create a `MeshSpec` for a L-dimensional marching-tetrahedra `Mesh` over a parallelepiped
 with axes given by the columns of `axes`. The points along axis `i` are distributed between
 `first(minmaxaxes[i])` and `last(minmaxaxes[i])`. The number of points on each axis is given
-by `resolution`, or `resolution[i]` if several are given.
+by `points`, or `points[i]` if several are given.
 
 # Examples
 
 ```jldoctest
-julia> buildmesh(marchingmesh((-π, π), (0,2π); resolution = 25))
+julia> buildmesh(marchingmesh((-π, π), (0,2π); points = 25))
 Mesh{2}: mesh of a 2-dimensional manifold
   Vertices   : 625
   Edges      : 1776
 
-julia> buildmesh(marchingmesh((-π, π), (0,2π); resolution = (10,10)))
+julia> buildmesh(marchingmesh((-π, π), (0,2π); points = (10,10)))
 Mesh{2}: mesh of a 2-dimensional manifold
   Vertices   : 100
   Edges      : 261
@@ -169,13 +169,13 @@ Mesh{2}: mesh of a 2-dimensional manifold
 # External links
 - Marching tetrahedra (https://en.wikipedia.org/wiki/Marching_tetrahedra) in Wikipedia
 """
-marchingmesh(minmaxaxes::Vararg{Tuple,L}; axes = 1.0 * I, resolution = 13) where {L} =
-    MarchingMeshSpec(minmaxaxes, SMatrix{L,L}(axes), resolution)
+marchingmesh(minmaxaxes::Vararg{Tuple,L}; axes = 1.0 * I, points = 13) where {L} =
+    MarchingMeshSpec(minmaxaxes, SMatrix{L,L}(axes), points)
 
 marchingmesh(; kw...) = throw(ArgumentError("Need a finite number of ranges to define a marching mesh"))
 
 function buildmesh(m::MarchingMeshSpec{D}, bravais = missing) where {D}
-    ranges = ((b, r)->range(b...; length = r)).(m.minmaxaxes, m.resolution)
+    ranges = ((b, r)->range(b...; length = r)).(m.minmaxaxes, m.points)
     npoints = length.(ranges)
     cs = CartesianIndices(ntuple(n -> 1:npoints[n], Val(D)))
     ls = LinearIndices(cs)
@@ -217,22 +217,22 @@ struct LinearMeshSpec{N,L,T,R} <: MeshSpec{1}
     vertices::SVector{N,SVector{L,T}}
     samelength::Bool
     closed::Bool
-    resolution::R
+    points::R
 end
 
 """
-    linearmesh(nodes...; resolution = 13, samelength = false, closed = false)
+    linearmesh(nodes...; points = 13, samelength = false, closed = false)
 
 Create a `MeshSpec` for a one-dimensional `Mesh` connecting the `nodes` with straight
-segments, each containing `resolution` points (endpoints included). If a different
-resolution for each of the `N` segments is required, use `resolution::NTuple{N,Int}`.
+segments, each containing a number `points` of points (endpoints included). If a different
+number of points for each of the `N` segments is required, use `points::NTuple{N,Int}`.
 If `samelength` each segment has equal length in mesh coordinates. If `closed` the last node
 is connected to the first node.
 
 # Examples
 
 ```jldoctest
-julia> buildmesh(linearmesh(:Γ, :K, :M, :Γ; resolution = (101, 30, 30)))
+julia> buildmesh(linearmesh(:Γ, :K, :M, :Γ; points = (101, 30, 30)))
 Mesh{1}: mesh of a 1-dimensional manifold
   Vertices   : 159
   Edges      : 158
@@ -241,8 +241,8 @@ Mesh{1}: mesh of a 1-dimensional manifold
 # See also
     `marchingmesh`, `buildmesh`
 """
-linearmesh(nodes...; resolution = 13, samelength::Bool = false, closed::Bool = false) =
-    LinearMeshSpec(sanitize_BZpts(nodes), samelength, closed, resolution)
+linearmesh(nodes...; points = 13, samelength::Bool = false, closed::Bool = false) =
+    LinearMeshSpec(sanitize_BZpts(nodes), samelength, closed, points)
 
 function sanitize_BZpts(pts)
     pts´ = parse_BZpoint.(pts)
@@ -304,7 +304,7 @@ end
 buildmesh(s::LinearMeshSpec{N,L,T}) where {N,L,T} = buildmesh(s, SMatrix{L,L,T}(I))
 
 function buildmesh(s::LinearMeshSpec{N}, br) where {N}
-    ranges = ((i, r) -> range(i, i+1, length = r)).(ntuple(identity, Val(N-1)), s.resolution)
+    ranges = ((i, r) -> range(i, i+1, length = r)).(ntuple(identity, Val(N-1)), s.points)
     verts = SVector.(first(ranges))
     for r in Base.tail(ranges)
         pop!(verts)
diff --git a/test/test_bandstructure.jl b/test/test_bandstructure.jl
index 85a5eb63..b0dfb371 100644
--- a/test/test_bandstructure.jl
+++ b/test/test_bandstructure.jl
@@ -1,22 +1,22 @@
 @testset "basic bandstructures" begin
     h = LatticePresets.honeycomb() |> hamiltonian(hopping(-1, range = 1/√3))
-    b = bandstructure(h, resolution = 13)
+    b = bandstructure(h, points = 13)
     @test length(bands(b)) == 1
 
     h = LatticePresets.honeycomb() |>
         hamiltonian(onsite(0.5, sublats = :A) + onsite(-0.5, sublats = :B) +
                     hopping(-1, range = 1/√3))
-    b = bandstructure(h, resolution = (13, 23))
+    b = bandstructure(h, points = (13, 23))
     @test length(bands(b)) == 2
 
     h = LatticePresets.cubic() |> hamiltonian(hopping(1)) |> unitcell(2)
-    b = bandstructure(h, resolution = (5, 9, 5))
+    b = bandstructure(h, points = (5, 9, 5))
     @test length(bands(b)) == 8
 
-    b = bandstructure(h, linearmesh(:Γ, :X, resolution = 4))
+    b = bandstructure(h, linearmesh(:Γ, :X, points = 4))
     @test length(bands(b)) == 8
 
-    b = bandstructure(h, linearmesh(:Γ, :X, (0, π), :Z, :Γ, resolution = 4))
+    b = bandstructure(h, linearmesh(:Γ, :X, (0, π), :Z, :Γ, points = 4))
     @test length(bands(b)) == 8
 end
 
@@ -57,8 +57,8 @@ end
 @testset "parametric bandstructures" begin
     ph = LatticePresets.linear() |> hamiltonian(hopping(-I), orbitals = Val(2)) |> unitcell(2) |>
          parametric(@onsite!((o; k) -> o + k*I), @hopping!((t; k)-> t - k*I))
-    mesh2D = marchingmesh((0, 1), (0, 2π), resolution = 25)
-    mesh1D = marchingmesh((0, 2π), resolution = 25)
+    mesh2D = marchingmesh((0, 1), (0, 2π), points = 25)
+    mesh1D = marchingmesh((0, 2π), points = 25)
     b = bandstructure(ph, mesh2D)
     @test length(bands(b)) == 2
     b = bandstructure(ph, mesh2D, lift = (k, φ) -> (k + φ/2π, φ))
diff --git a/test/test_mesh.jl b/test/test_mesh.jl
index b43944df..718878a0 100644
--- a/test/test_mesh.jl
+++ b/test/test_mesh.jl
@@ -1,6 +1,6 @@
 using Quantica: nvertices, nedges, nsimplices
 
 @test begin
-    mesh = buildmesh(marchingmesh((0,1), (0,2), resolution = (10, 20)))
+    mesh = buildmesh(marchingmesh((0,1), (0,2), points = (10, 20)))
     nvertices(mesh) == 200 && nedges(mesh) == 541 && nsimplices(mesh) == 342
 end
\ No newline at end of file

From 3cf86fe2e8a6ac86d699f5187f0381cf81ef1642 Mon Sep 17 00:00:00 2001
From: Pablo San-Jose <lekand@gmail.com>
Date: Tue, 16 Jun 2020 18:26:15 +0200
Subject: [PATCH 3/6] docstring

---
 src/bandstructure.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/bandstructure.jl b/src/bandstructure.jl
index d8a7f3c9..bdc76a58 100644
--- a/src/bandstructure.jl
+++ b/src/bandstructure.jl
@@ -142,7 +142,7 @@ end
     bandstructure(h::Hamiltonian; points = 13, kw...)
 
 Compute the bandstructure of `h` on a mesh over `h`'s full Brillouin zone, with `points`
-points along each axis, spanning the interval [-π,π].
+points along each axis, spanning the interval [-π,π] along each reciprocal axis.
 
     bandstructure(h::Hamiltonian, spec::MeshSpec; lift = missing, kw...)
 

From 4aa3647b66feb1ecb4602738ff42b8f408df1ec5 Mon Sep 17 00:00:00 2001
From: Pablo San-Jose <lekand@gmail.com>
Date: Thu, 18 Jun 2020 12:51:13 +0200
Subject: [PATCH 4/6] allow lifts in bandstructure(::Function, ::Mesh)

---
 src/bandstructure.jl | 25 ++++++++++++-------------
 1 file changed, 12 insertions(+), 13 deletions(-)

diff --git a/src/bandstructure.jl b/src/bandstructure.jl
index bdc76a58..391b1ecb 100644
--- a/src/bandstructure.jl
+++ b/src/bandstructure.jl
@@ -168,8 +168,7 @@ reads `v = (p₁,..., pᵢ, ϕ₁,..., ϕⱼ)`, with `p` the values assigned to
     bandstructure(matrixf::Function, mesh::Mesh; kw...)
 
 Compute the bandstructure of the Hamiltonian matrix `m = matrixf(ϕ)`, with `ϕ` evaluated on
-the vertices `v` of `mesh`. No `lift` option is allowed in this case. If needed, include it
-in the definition of `matrixf`. Note that `ϕ` is either a `Tuple` or an `SVector`, but it is
+the vertices `v` of the `mesh`. Note that `ϕ` is either a `Tuple` or an `SVector`, but it is
 not splatted into `matrixf`, i.e. `matrixf(x) = ...` or `matrixf(x, y) = ...` will not work,
 use `matrixf((x,)) = ...` or `matrixf((x, y)) = ...` instead.
 
@@ -251,35 +250,35 @@ end
 bravais_parameters(h::Hamiltonian) = bravais(h)
 bravais_parameters(ph::ParametricHamiltonian{P}) where {P} = _blockdiag(SMatrix{P,P}(I), bravais(ph))
 
-bandstructure(h::Function, spec::MeshSpec; kw...) =
-    bandstructure(h, buildmesh(spec); kw...)
-
 function bandstructure(h::Union{Hamiltonian,ParametricHamiltonian}, mesh::Mesh;
                        method = defaultmethod(h), lift = missing, transform = missing, kw...)
     # ishermitian(h) || throw(ArgumentError("Hamiltonian must be hermitian"))
     matrix = similarmatrix(h, method)
     codiag = codiagonalizer(h, matrix, mesh, lift)
     d = DiagonalizeHelper(method, codiag; kw...)
-    matrixf(φs) = bloch!(matrix, h, applylift(lift, φs))
+    matrixf(ϕs) = bloch!(matrix, h, applylift(lift, ϕs))
     b = _bandstructure(matrixf, matrix, mesh, d)
     transform === missing || transform!(transform, b)
     return b
 end
 
-# Should perhaps RFC to be merged with the above
 function bandstructure(matrixf::Function, mesh::Mesh;
-                       matrix = _samplematrix(matrixf, mesh),
-                       method = defaultmethod(matrix),
-                       minprojection = 0.5, transform = missing, kw...)
-    codiag = codiagonalizer(matrixf, matrix, mesh)
-    d = DiagonalizeHelper(method, codiag; kw...)
-    b = _bandstructure(matrixf, matrix, mesh, d)
+                       method = missing, lift = missing, minprojection = 0.5, transform = missing, kw...)
+    matrixf´ = _wraplift(matrixf, lift)
+    matrix = _samplematrix(matrixf´, mesh)
+    method´ = method === missing ? defaultmethod(matrix) : method
+    codiag = codiagonalizer(matrixf´, matrix, mesh)
+    d = DiagonalizeHelper(method´, codiag; kw...)
+    b = _bandstructure(matrixf´, matrix, mesh, d)
     transform === missing || transform!(transform, b)
     return b
 end
 
 _samplematrix(matrixf, mesh) = matrixf(Tuple(first(vertices(mesh))))
 
+_wraplift(matrixf, lift::Missing) = matrixf
+_wraplift(matrixf, lift) = ϕs -> matrixf(applylift(lift, ϕs))
+
 @inline applylift(lift::Missing, ϕs) = toSVector(ϕs)
 
 @inline applylift(lift::Function, ϕs) = toSVector(lift(ϕs...))

From 68566597531cd9bcb932e2013bf343b30b962e24 Mon Sep 17 00:00:00 2001
From: Pablo San-Jose <lekand@gmail.com>
Date: Thu, 18 Jun 2020 14:25:22 +0200
Subject: [PATCH 5/6] buildmesh(::MeshSpec, ::Hamiltonian) + fixes

---
 src/bandstructure.jl       | 10 +++-----
 src/mesh.jl                | 52 ++++++++++++++++++++++++++++++++------
 test/test_bandstructure.jl |  4 +--
 3 files changed, 49 insertions(+), 17 deletions(-)

diff --git a/src/bandstructure.jl b/src/bandstructure.jl
index 391b1ecb..7f6554b0 100644
--- a/src/bandstructure.jl
+++ b/src/bandstructure.jl
@@ -203,7 +203,7 @@ sigma = 1im)` to compute the bandstructure.
 bandstructure (useful for performing shifts or other postprocessing).
 
 # Examples
-```
+```jldoctest
 julia> h = LatticePresets.honeycomb() |> hamiltonian(hopping(-1, range = 1/√3)) |> unitcell(3);
 
 julia> bandstructure(h; points = 25, method = LinearAlgebraPackage())
@@ -241,15 +241,11 @@ function bandstructure(h::Hamiltonian; points = 13, kw...)
 end
 
 function bandstructure(h::Union{Hamiltonian,ParametricHamiltonian}, spec::MeshSpec; lift = missing, kw...)
-    br = bravais_parameters(h)
-    mesh = buildmesh(spec, br)
-    lift´ = lift === missing ? buildlift(spec, br) : lift
+    mesh = buildmesh(spec, h)
+    lift´ = lift === missing ? buildlift(spec, h) : lift
     return bandstructure(h, mesh; lift = lift´, kw...)
 end
 
-bravais_parameters(h::Hamiltonian) = bravais(h)
-bravais_parameters(ph::ParametricHamiltonian{P}) where {P} = _blockdiag(SMatrix{P,P}(I), bravais(ph))
-
 function bandstructure(h::Union{Hamiltonian,ParametricHamiltonian}, mesh::Mesh;
                        method = defaultmethod(h), lift = missing, transform = missing, kw...)
     # ishermitian(h) || throw(ArgumentError("Hamiltonian must be hermitian"))
diff --git a/src/mesh.jl b/src/mesh.jl
index e6f0ceb3..df619d77 100644
--- a/src/mesh.jl
+++ b/src/mesh.jl
@@ -128,9 +128,43 @@ abstract type MeshSpec{L} end
 Base.show(io::IO, spec::MeshSpec{L}) where {L} =
     print(io, "MeshSpec{$L} : specifications for building a $(L)D mesh.")
 
-#fallback
-buildlift(::MeshSpec, bravais) = missing
+"""
+    buildlift(s::MeshSpec{L}, h::Union{Hamiltonian,ParametricHamiltonian})
+
+Build a lift function that maps a `Mesh` built with `buildmesh(s, h)` to the
+Brillouin/parameter space of `h` (see `bandstructure` for details).
+
+# See also
+    `buildmesh`, `marchingmesh`, `linearmesh`
+"""
+buildlift(s::MeshSpec, h::Union{Hamiltonian,ParametricHamiltonian}) =
+    buildlift(s, bravais_parameters(h))
 
+"""
+    buildmesh(s::MeshSpec, h::Union{Hamiltonian,ParametricHamiltonian})
+
+Build a `Mesh` from the spec `s`, using properties of `h` as needed. The use of `h` depends
+on the spec. For a `LinearMeshSpec` with `samelength = false`, the Bravais matrix of `h` is
+needed to work out the length of each mesh segment in the Brillouin zone, while for other
+specs such as `MarchingMeshSpec`, `h` is not needed and may be omitted (see example).
+
+# Examples
+
+```jldoctest
+julia> buildmesh(marchingmesh((-π, π), (0, 2π), points = 10))
+Mesh{2}: mesh of a 2-dimensional manifold
+  Vertices   : 100
+  Edges      : 261
+```
+
+# See also
+    `buildlift`, `marchingmesh`, `linearmesh`
+"""
+buildmesh(s::MeshSpec, h::Union{Hamiltonian,ParametricHamiltonian}) =
+    buildmesh(s, bravais_parameters(h))
+
+bravais_parameters(h::Hamiltonian) = bravais(h)
+bravais_parameters(ph::ParametricHamiltonian{P}) where {P} = _blockdiag(SMatrix{P,P}(I), bravais(ph))
 
 #######################################################################
 # MarchingMeshSpec
@@ -174,7 +208,7 @@ marchingmesh(minmaxaxes::Vararg{Tuple,L}; axes = 1.0 * I, points = 13) where {L}
 
 marchingmesh(; kw...) = throw(ArgumentError("Need a finite number of ranges to define a marching mesh"))
 
-function buildmesh(m::MarchingMeshSpec{D}, bravais = missing) where {D}
+function buildmesh(m::MarchingMeshSpec{D}, bravais::SMatrix = SMatrix{D,D}(I)) where {D}
     ranges = ((b, r)->range(b...; length = r)).(m.minmaxaxes, m.points)
     npoints = length.(ranges)
     cs = CartesianIndices(ntuple(n -> 1:npoints[n], Val(D)))
@@ -242,10 +276,14 @@ Mesh{1}: mesh of a 1-dimensional manifold
     `marchingmesh`, `buildmesh`
 """
 linearmesh(nodes...; points = 13, samelength::Bool = false, closed::Bool = false) =
-    LinearMeshSpec(sanitize_BZpts(nodes), samelength, closed, points)
+    LinearMeshSpec(sanitize_BZpts(nodes, closed), samelength, closed, points)
 
-function sanitize_BZpts(pts)
+function sanitize_BZpts(pts, closed)
     pts´ = parse_BZpoint.(pts)
+    if closed
+        all(isapprox.(first(pts´), last(pts´))) ||
+            throw(ArgumentError("Closed linear meshes should have equal first and last nodes."))
+    end
     dim = maximum(length.(pts´))
     pts´´ = SVector(padright.(pts´, Val(dim)))
     return pts´´
@@ -301,9 +339,7 @@ function idx_to_node(s, br)
     return nodefunc
 end
 
-buildmesh(s::LinearMeshSpec{N,L,T}) where {N,L,T} = buildmesh(s, SMatrix{L,L,T}(I))
-
-function buildmesh(s::LinearMeshSpec{N}, br) where {N}
+function buildmesh(s::LinearMeshSpec{N}, br::SMatrix) where {N}
     ranges = ((i, r) -> range(i, i+1, length = r)).(ntuple(identity, Val(N-1)), s.points)
     verts = SVector.(first(ranges))
     for r in Base.tail(ranges)
diff --git a/test/test_bandstructure.jl b/test/test_bandstructure.jl
index b0dfb371..9b067e20 100644
--- a/test/test_bandstructure.jl
+++ b/test/test_bandstructure.jl
@@ -24,7 +24,7 @@ end
     hc = LatticePresets.honeycomb() |> hamiltonian(hopping(-1, sublats = :A=>:B, plusadjoint = true))
     matrix = similarmatrix(hc, LinearAlgebraPackage())
     hf((x,)) = bloch!(matrix, hc, (x, -x))
-    mesh = marchingmesh((0, 1))
+    mesh = buildmesh(marchingmesh((0, 1)))
     b = bandstructure(hf, mesh)
     @test length(bands(b)) == 2
 
@@ -32,7 +32,7 @@ end
     hp2 = parametric(hc2, @hopping!((t; s) -> s*t))
     matrix2 = similarmatrix(hc2, LinearAlgebraPackage())
     hf2((s, x)) = bloch!(matrix2, hp2(s = s), (x, x))
-    mesh2 = marchingmesh((0, 1), (0, 1))
+    mesh2 = buildmesh(marchingmesh((0, 1), (0, 1)))
     b = bandstructure(hf2, mesh2)
     @test length(bands(b)) == 2
 end

From cea2477826585af83e1b2ac7a80d5090fb19520e Mon Sep 17 00:00:00 2001
From: Pablo San-Jose <lekand@gmail.com>
Date: Thu, 18 Jun 2020 14:47:44 +0200
Subject: [PATCH 6/6] docfix

---
 src/bandstructure.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/bandstructure.jl b/src/bandstructure.jl
index 7f6554b0..106a1af4 100644
--- a/src/bandstructure.jl
+++ b/src/bandstructure.jl
@@ -168,9 +168,9 @@ reads `v = (p₁,..., pᵢ, ϕ₁,..., ϕⱼ)`, with `p` the values assigned to
     bandstructure(matrixf::Function, mesh::Mesh; kw...)
 
 Compute the bandstructure of the Hamiltonian matrix `m = matrixf(ϕ)`, with `ϕ` evaluated on
-the vertices `v` of the `mesh`. Note that `ϕ` is either a `Tuple` or an `SVector`, but it is
-not splatted into `matrixf`, i.e. `matrixf(x) = ...` or `matrixf(x, y) = ...` will not work,
-use `matrixf((x,)) = ...` or `matrixf((x, y)) = ...` instead.
+the vertices `v` of the `mesh`. Note that `ϕ` in `matrixf(ϕ)` is an unsplatted container.
+Hence, i.e. `matrixf(x) = ...` or `matrixf(x, y) = ...` will not work, use `matrixf((x,)) =
+...` or `matrixf((x, y)) = ...` instead.
 
 # Options