diff --git a/.gitignore b/.gitignore
index 3af67b1..4348ced 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,32 @@
-*.jl.*.cov
+
+
+# Files generated by invoking Julia with --code-coverage
 *.jl.cov
+*.jl.*.cov
+
+# Files generated by invoking Julia with --track-allocation
 *.jl.mem
-/docs/build/
+*.jl.*.mem
+
+# System-specific files and directories generated by the BinaryProvider and BinDeps packages
+# They contain absolute paths specific to the host computer, and so should not be committed
+deps/deps.jl
+deps/build.log
+deps/downloads/
+deps/usr/
+deps/src/
+
+# Build artifacts for creating documentation generated by the Documenter package
+docs/build/
+docs/site/
+/run/
+
+# File generated by Pkg, the package manager, based on a corresponding Project.toml
+# It records a fixed state of all packages used by the project. As such, it should not be
+# committed for packages, but should be committed for applications that require a static
+# environment.
+Manifest.toml
+
+*.DS_Store
+
+.ipynb_checkpoints
\ No newline at end of file
diff --git a/Manifest.toml b/Manifest.toml
deleted file mode 100644
index f45eecf..0000000
--- a/Manifest.toml
+++ /dev/null
@@ -1,2 +0,0 @@
-# This file is machine-generated - editing it directly is not advised
-
diff --git a/Project.toml b/Project.toml
index 00c8ad0..d217f50 100644
--- a/Project.toml
+++ b/Project.toml
@@ -3,6 +3,11 @@ uuid = "3d4565c3-192b-4227-89f7-5476f913ab49"
 authors = ["Tao Wang", "Xiansheng Cai"]
 version = "0.1.0"
 
+[deps]
+AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
 [compat]
 julia = "1.6"
 
diff --git a/example/testabstracttree.jl b/example/testabstracttree.jl
new file mode 100644
index 0000000..4cd00d0
--- /dev/null
+++ b/example/testabstracttree.jl
@@ -0,0 +1,23 @@
+using AbstractTrees
+
+struct MyNode
+    val::Int
+    children::Vector{MyNode}
+end
+
+AbstractTrees.children(node::MyNode) = node.children
+
+function MyNode(val::Int)
+    val = val
+    children = Vector{MyNode}()
+    for i in 1:val
+        push!(children, MyNode(val-1))
+    end
+
+    return MyNode(val, children)
+end
+
+tree = MyNode(3)
+
+print_tree(tree)
+
diff --git a/example/tree.jl b/example/tree.jl
new file mode 100644
index 0000000..6fcb81d
--- /dev/null
+++ b/example/tree.jl
@@ -0,0 +1,153 @@
+
+#--- test tree divided BZ ---
+
+using Random
+using LinearAlgebra
+using Plots
+
+struct SquareNod
+    gridpoints::Matrix{Float64}
+    sons::Vector{SquareNod}
+end
+
+function minarea(sn::SquareNod)
+    return (isempty(sn.sons)) ? area(sn) : min([minarea(ssn) for ssn in sn.sons]...)
+end
+
+function SquareNod!(gridpoints, weightcap::Float64, gpcontainer; treelevelcap = 5)
+
+    K00, K01, K10, K11 = Matrix(gridpoints[1,:]'), Matrix(gridpoints[2,:]'), Matrix(gridpoints[3,:]'), Matrix(gridpoints[4,:]')
+    K0h = 0.5 .* (K00 .+ K01)
+    Kh0 = 0.5 .* (K00 .+ K10)
+    K1h = 0.5 .* (K10 .+ K11)
+    Kh1 = 0.5 .* (K01 .+ K11)
+    Khh = 0.5 .* (K00 .+ K11)
+
+    push!(gpcontainer, K0h)
+    push!(gpcontainer, Kh0)
+    # push!(gpcontainer, K1h)
+    # push!(gpcontainer, Kh1)
+    push!(gpcontainer, Khh)
+
+    sons = []
+    if treelevelcap > 1 && weight(gridpoints) > weightcap
+        sn00 = SquareNod!(Matrix([K00; K0h; Kh0; Khh]), weightcap, gpcontainer;treelevelcap= treelevelcap-1)
+        sn01 = SquareNod!(Matrix([K0h; K01; Khh; Kh1]), weightcap, gpcontainer;treelevelcap= treelevelcap-1)
+        sn10 = SquareNod!(Matrix([Kh0; Khh; K10; K1h]), weightcap, gpcontainer;treelevelcap= treelevelcap-1)
+        sn11 = SquareNod!(Matrix([Khh; Kh1; K1h; K11]), weightcap, gpcontainer;treelevelcap= treelevelcap-1)
+
+        push!(sons, sn00)
+        push!(sons, sn01)
+        push!(sons, sn10)
+        push!(sons, sn11)
+    end
+
+    return SquareNod(gridpoints, sons)
+end
+
+function SquareNod(gridpoints, weightcap::Float64; treelevelcap = 5)
+
+    K00, K01, K10, K11 = Matrix(gridpoints[1,:]'), Matrix(gridpoints[2,:]'), Matrix(gridpoints[3,:]'), Matrix(gridpoints[4,:]')
+    K0h = 0.5 .* (K00 .+ K01)
+    Kh0 = 0.5 .* (K00 .+ K10)
+    K1h = 0.5 .* (K10 .+ K11)
+    Kh1 = 0.5 .* (K01 .+ K11)
+    Khh = 0.5 .* (K00 .+ K11)
+
+    sons = []
+    if treelevelcap > 1 && weight(gridpoints) > weightcap
+        sn00 = SquareNod(Matrix([K00; K0h; Kh0; Khh]), weightcap;treelevelcap= treelevelcap-1)
+        sn01 = SquareNod(Matrix([K0h; K01; Khh; Kh1]), weightcap;treelevelcap= treelevelcap-1)
+        sn10 = SquareNod(Matrix([Kh0; Khh; K10; K1h]), weightcap;treelevelcap= treelevelcap-1)
+        sn11 = SquareNod(Matrix([Khh; Kh1; K1h; K11]), weightcap;treelevelcap= treelevelcap-1)
+
+        push!(sons, sn00)
+        push!(sons, sn01)
+        push!(sons, sn10)
+        push!(sons, sn11)
+    end
+
+    return SquareNod(gridpoints, sons)
+end
+
+@inline function density(K::AbstractVector)
+    me = 0.5
+    T = 0.01
+
+    # μ = 1.0
+    # k = norm(K)
+    # ϵ = k^2 / (2me)
+    # dϵdk = k / me
+
+    μ = 0.5
+    ϵ = -(cos(π*K[1]) + cos(π*K[2])) / (2me)
+    dϵdk = sqrt(sin(π*K[1])^2 + sin(π*K[2])^2) / (2me)
+
+    return (exp((ϵ - μ) / T) / T)/(exp((ϵ - μ) / T) + 1.0)^2 * dϵdk
+end
+
+function area(gridpoints)
+    a = norm(gridpoints[2,:]-gridpoints[1,:])
+    b = norm(gridpoints[3,:]-gridpoints[1,:])
+    return a*b
+end
+
+area(sn::SquareNod) = area(sn.gridpoints)
+
+function weight(gridpoints; N=400)
+    result = 0.0
+    a = area(gridpoints)
+
+    K1, K2 = gridpoints[2,:]-gridpoints[1,:], gridpoints[3,:]-gridpoints[1,:]
+    for i in 1:N
+        K = gridpoints[1, :] .+ rand() .* K1 .+ rand() .* K2
+        result += density(K)
+    end
+
+    return result * a
+end
+
+weight(sn::SquareNod; N=100) = weight(sn.gridpoints; N=N)
+
+gridpoints = Matrix([-1.0 -1.0 ; -1 1 ; 1 -1 ; 1 1])
+K00, K01, K10, K11 = Matrix(gridpoints[1,:]'), Matrix(gridpoints[2,:]'), Matrix(gridpoints[3,:]'), Matrix(gridpoints[4,:]')
+gpcontainer = [K00,]
+
+sn = SquareNod!(gridpoints, weight(gridpoints)/128^2, gpcontainer; treelevelcap = 12)
+
+println(area(sn))
+println(weight(sn))
+
+function Base.print(sn::SquareNod; treelevel = 1)
+    for i in 1:treelevel
+        print(" ")
+    end
+    print("area=$(area(sn)),")
+    print("weight=$(weight(sn)),")
+    print("gridpoints:",sn.gridpoints)
+    print("\n")
+
+    for ssn in sn.sons
+        print(ssn, ;treelevel = treelevel + 1)
+    end
+end
+
+# print(sn)
+
+println("minarea = $(minarea(sn))")
+# print(gpcontainer)
+println("length:",length(gpcontainer))
+
+println("compress:", minarea(sn)*length(gpcontainer)/4)
+
+X, Y = zeros(Float64, length(gpcontainer)), zeros(Float64, length(gpcontainer))
+for (Ki, K) in enumerate(gpcontainer)
+    X[Ki], Y[Ki] = K[1], K[2]
+end
+
+# p = plot(legend = false, size = (1024, 1024))
+# scatter!(p, X, Y, marker = :cross, markersize = 2)
+# savefig(p, "run/grid.pdf")
+# display(p)
+# readline()
+
diff --git a/example/treegrid.jl b/example/treegrid.jl
new file mode 100644
index 0000000..704203a
--- /dev/null
+++ b/example/treegrid.jl
@@ -0,0 +1,98 @@
+#--- test tree divided BZ ---
+
+using Random
+using LinearAlgebra
+using Plots
+
+struct MeshNod
+    depth::Int
+    pos::Vector{Int}
+    sons::Vector{MeshNod}
+end
+
+struct TreeMesh
+    dim::Int
+    root::MeshNod
+    latticevectors::Matrix{Float64}
+    gridpoints::Matrix{Float64}
+end
+
+function TreeMesh(latticevectors::Matrix, weightcap::Float64; depthcap = 5)
+    
+end
+
+@inline function density(K::AbstractVector)
+    me = 0.5
+    T = 0.01
+
+    μ = 1.0
+    k = norm(K)
+    ϵ = k^2 / (2me)
+    dϵdk = k / me
+
+    # μ = 0.5
+    # ϵ = -(cos(π*K[1]) + cos(π*K[2])) / (2me)
+    # dϵdk = sqrt(sin(π*K[1])^2 + sin(π*K[2])^2) / (2me)
+
+    return (exp((ϵ - μ) / T) / T)/(exp((ϵ - μ) / T) + 1.0)^2 * dϵdk
+end
+
+function area(gridpoints)
+    a = norm(gridpoints[2,:]-gridpoints[1,:])
+    b = norm(gridpoints[3,:]-gridpoints[1,:])
+    return a*b
+end
+
+area(sn::SquareNod) = area(sn.gridpoints)
+
+function weight(gridpoints; N=100)
+    result = 0.0
+    a = area(gridpoints)
+
+    K1, K2 = gridpoints[2,:]-gridpoints[1,:], gridpoints[3,:]-gridpoints[1,:]
+    for i in 1:N
+        K = gridpoints[1, :] .+ rand() .* K1 .+ rand() .* K2
+        result += density(K)
+    end
+
+    return result * a
+end
+
+weight(sn::SquareNod; N=100) = weight(sn.gridpoints; N=N)
+
+gridpoints = Matrix([-1.0 -1.0 ; -1 1 ; 1 -1 ; 1 1])
+K00, K01, K10, K11 = Matrix(gridpoints[1,:]'), Matrix(gridpoints[2,:]'), Matrix(gridpoints[3,:]'), Matrix(gridpoints[4,:]')
+gpcontainer = [K00, K01, K10, K11]
+
+sn = SquareNod!(gridpoints, weight(gridpoints)/500, gpcontainer; treelevelcap = 10)
+
+println(area(sn))
+println(weight(sn))
+
+function Base.print(sn::SquareNod; treelevel = 1)
+    for i in 1:treelevel
+        print(" ")
+    end
+    print("area=$(area(sn)),")
+    print("weight=$(weight(sn)),")
+    print("gridpoints:",sn.gridpoints)
+    print("\n")
+
+    for ssn in sn.sons
+        print(ssn, ;treelevel = treelevel + 1)
+    end
+end
+
+print(sn)
+# print(gpcontainer)
+println("length:",length(gpcontainer))
+
+X, Y = zeros(Float64, length(gpcontainer)), zeros(Float64, length(gpcontainer))
+for (Ki, K) in enumerate(gpcontainer)
+    X[Ki], Y[Ki] = K[1], K[2]
+end
+
+p = scatter(X, Y)
+display(p)
+readline()
+
diff --git a/example/treenode.jl b/example/treenode.jl
new file mode 100644
index 0000000..7a206d8
--- /dev/null
+++ b/example/treenode.jl
@@ -0,0 +1,66 @@
+using SpaceGrid
+using SpaceGrid.AbstractTrees, SpaceGrid.GridTree, SpaceGrid.BaseMesh
+using Plots
+using LinearAlgebra
+
+mesh = UniformMesh{2, 4}([0.0, 0.0], [0 1 ; 1 0])
+show(mesh)
+
+naiveisfine(depth, pos) = depth >= 2
+
+# tree = GridNode{2}(naiveisfine)
+# print_tree(tree)
+
+# for node in PostOrderDFS(tree)
+#     if isempty(node.children)
+#         println(node.pos)
+#     end
+# end
+
+function dispersion(k)
+    me = 0.5
+    μ = 1.0
+    return norm(k)^2/(2me)-μ
+
+    # t = 0.5
+    # μ = 0.0
+    # return sum([t * cos(π*p) for p in k]) + μ
+ end
+
+function density(K, latvec)
+    DIM = length(K)
+
+    T = 0.001
+
+    kpoints = [K, ]
+    for i in 1:DIM
+        push!(kpoints, K .+ latvec[i, :])
+        push!(kpoints, K .- latvec[i, :])
+    end
+
+    ϵ = minimum([dispersion(k) for k in kpoints])
+
+    return 1 / (exp((ϵ) / T) + 1.0)
+    # return 1 / ((π * T)^2 + ϵ^2)
+end
+
+latvec = [2 0; 1 sqrt(3)]
+# latvec = [2 0; 0 2]
+# tg = uniformtreegrid(naiveisfine, latvec; N = 3)
+# tg = treegridfromdensity(density, latvec; rtol = 1e-4, maxdepth = 5)
+tg = treegridfromdensity(k->density(k, latvec), latvec; rtol = 1e-1, maxdepth = 6)
+
+X, Y = zeros(Float64, size(tg)), zeros(Float64, size(tg))
+for (pi, p) in enumerate(tg)
+    X[pi], Y[pi] = p[1], p[2]
+    # println(p)
+end
+
+println("size:$(size(tg))")
+println("length:$(length(tg))")
+println("efficiency:$(efficiency(tg))")
+
+# p = plot(legend = false, size = (1024, 1024), xlim = (-1, 1), ylim = (-1, 1))
+p = plot(legend = false, size = (1024, 1024), xlim = (-1.5, 1.5), ylim = (-1.5, 1.5))
+scatter!(p, X, Y, marker = :cross, markersize = 2)
+savefig(p, "run/tg.pdf")
diff --git a/src/SpaceGrid.jl b/src/SpaceGrid.jl
index ef87a54..1a620ab 100644
--- a/src/SpaceGrid.jl
+++ b/src/SpaceGrid.jl
@@ -1,5 +1,15 @@
 module SpaceGrid
 
+using AbstractTrees
+using StaticArrays
+using Statistics
+
 # Write your package code here.
 
+include("basemesh.jl")
+using .BaseMesh
+
+include("tree.jl")
+using .GridTree
+
 end
diff --git a/src/basemesh.jl b/src/basemesh.jl
new file mode 100644
index 0000000..32fdb19
--- /dev/null
+++ b/src/basemesh.jl
@@ -0,0 +1,47 @@
+module BaseMesh
+
+using ..StaticArrays
+
+export UniformMesh
+
+abstract type AbstractMesh end
+
+struct UniformMesh{DIM, N} <: AbstractMesh
+    origin::SVector{DIM, Float64}
+    latvec::SMatrix{DIM, DIM, Float64}
+end
+
+Base.length(mesh::UniformMesh{DIM, N}) where {DIM, N} = N
+Base.size(mesh::UniformMesh{DIM, N}) where {DIM, N} = N^DIM
+function Base.show(io::IO, mesh::UniformMesh)
+    println("UniformMesh Grid:")
+    for (pi, p) in enumerate(mesh)
+        println(p)
+    end
+end
+
+# Base.view(mesh::UniformMesh, i::Int) = Base.view(mesh.mesh, :, i)
+
+# # set is not allowed for meshs
+function Base.getindex(mesh::UniformMesh{DIM, N}, inds...) where {DIM, N}
+    pos = Vector(mesh.origin)
+    for (ni, n) in enumerate(inds)
+        pos = pos .+ mesh.latvec[ni, :] .* (n - 1) ./ (length(mesh))
+    end
+    return pos
+end
+function Base.getindex(mesh::UniformMesh{DIM, N}, i::Int) where {DIM, N}
+    inds = digits(i-1, base = N)
+    pos = Vector(mesh.origin)
+    for (ni, n) in enumerate(inds)
+        pos = pos .+ mesh.latvec[ni, :] .* n ./ (N)
+    end
+    return pos
+end
+Base.firstindex(mesh::UniformMesh) = 1
+Base.lastindex(mesh::UniformMesh) = size(mesh)
+# # iterator
+Base.iterate(mesh::UniformMesh) = (mesh[1],1)
+Base.iterate(mesh::UniformMesh, state) = (state>=size(mesh)) ? nothing : (mesh[state+1],state+1)
+
+end
diff --git a/src/tree.jl b/src/tree.jl
new file mode 100644
index 0000000..f91a7aa
--- /dev/null
+++ b/src/tree.jl
@@ -0,0 +1,124 @@
+module GridTree
+
+using ..AbstractTrees
+using ..StaticArrays
+using ..BaseMesh
+using ..Statistics
+
+export GridNode, TreeGrid, uniformtreegrid, treegridfromdensity, efficiency
+
+struct GridNode{DIM}
+    depth::Int
+    pos::SVector{DIM, Int64}
+    children::Vector{GridNode{DIM}}
+end
+
+function GridNode{DIM}(isfine; depth = 0, pos = SVector{DIM, Int64}(zeros(Int64, DIM)), maxdepth = 10) where {DIM}
+    if isfine(depth, pos) || depth == maxdepth
+        return GridNode{DIM}(depth, pos, Vector{GridNode{DIM}}([]))
+    else
+        children = Vector{GridNode{DIM}}([])
+        for i in 0:2^DIM-1
+            bin = digits(i, base = 2, pad = DIM) |> reverse
+            childdepth = depth + 1
+            childpos = pos .* 2 .+ bin
+            push!(children, GridNode{DIM}(isfine; depth = childdepth, pos = childpos, maxdepth = 10))
+        end
+        return GridNode{DIM}(depth, pos, children)
+    end
+end
+
+AbstractTrees.children(node::GridNode) = node.children
+
+function efficiency(root::GridNode{DIM}) where {DIM}
+    np, depth = 0, 0
+    for node in PostOrderDFS(root)
+        if node.depth > depth
+            depth = node.depth
+        end
+        if isempty(node.children)
+            np = np + 1
+        end
+    end
+    return np / 2^(depth*DIM)
+end
+
+struct TreeGrid{DIM, SG}
+    root::GridNode{DIM}
+    latvec::SMatrix{DIM, DIM, Float64}
+    subgrids::Vector{SG}
+end
+
+efficiency(tg::TreeGrid) = efficiency(tg.root)
+
+function _calc_point(depth, pos, latvec)
+    DIM = length(pos)
+    ratio = pos ./ 2^depth .- 0.5
+    origin = zeros(size(ratio))
+
+    for i in 1:DIM
+        origin .+= ratio[i] .* latvec[i, :]
+    end
+
+    return origin
+end
+
+function _calc_origin(node::GridNode{DIM}, latvec) where {DIM}
+    return _calc_point(node.depth, node.pos, latvec)
+end
+
+function _calc_cornerpoints(depth, pos, latvec)
+    DIM = length(pos)
+    points = []
+    for i in 0:2^DIM-1
+        ii = digits(i, base = 2, pad = DIM)
+        push!(points, _calc_point(depth, pos .+ ii, latvec))
+    end
+    return points
+end
+
+function uniformtreegrid(isfine, latvec; maxdepth = 10, DIM = 2, N = 2) 
+    root = GridNode{DIM}(isfine; maxdepth = maxdepth)
+    subgrids = Vector{UniformMesh{DIM, N}}([])
+    for node in PostOrderDFS(root)
+        if isempty(node.children)
+            depth = node.depth
+            origin = _calc_origin(node, latvec )
+            mesh = UniformMesh{DIM, N}(origin, latvec ./ 2^depth)
+            push!(subgrids, mesh)
+        end
+    end
+    return TreeGrid{DIM, UniformMesh{DIM, N}}(root, latvec, subgrids)
+end
+
+Base.length(tg::TreeGrid{DIM, SG}) where {DIM, SG} = length(tg.subgrids)
+Base.size(tg::TreeGrid{DIM, SG}) where {DIM, SG} = length(tg) * size(tg.subgrids[1])
+# index and iterator
+function Base.getindex(tg::TreeGrid{DIM, SG}, i) where {DIM, SG}
+    sgsize = size(tg.subgrids[1])
+
+    tgi, sgi = (i - 1) ÷ sgsize + 1, (i - 1) % sgsize + 1
+
+    return getindex(tg.subgrids[tgi], sgi)
+end
+Base.firstindex(tg::TreeGrid) = 1
+Base.lastindex(tg::TreeGrid) = size(tg)
+
+Base.iterate(tg::TreeGrid) = (tg[1],1)
+Base.iterate(tg::TreeGrid, state) = (state>=size(tg)) ? nothing : (tg[state+1],state+1)
+
+function densityisfine(density, latvec, depth, pos, rtol)
+    cornerpoints = _calc_cornerpoints(depth, pos, latvec)
+    cpval = [density(p) for p in cornerpoints]
+    push!(cpval, density(_calc_point(depth, pos .+ 0.5, latvec)))
+    # return abs(sum(cpval) / length(cpval) - density(_calc_point(depth, pos .+ 0.5, latvec))) < rtol
+    return std(cpval) < rtol
+end
+
+function treegridfromdensity(density, latvec; rtol = 1e-4, maxdepth = 10, DIM = 2, N = 2)
+    isfine(depth, pos) = densityisfine(density, latvec, depth, pos, rtol)
+    return uniformtreegrid(isfine, latvec; maxdepth = maxdepth, DIM = DIM, N = N)
+end
+
+end
+