diff --git a/Project.toml b/Project.toml
index da323de..5dd99c8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,6 +5,7 @@ version = "0.1.1"
 
 [deps]
 BEAST = "bb4162c7-ba94-5a20-af32-d8ec4428bdd1"
+ClusterTrees = "5100927d-02aa-593a-b4f9-7235df19f0db"
 CompScienceMeshes = "3e66a162-7b8c-5da0-b8f8-124ecd2c3ae1"
 ExaFMMt = "add1291b-d076-47cc-8667-64576de04ee0"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/ext/NminClusterTrees/NminClusterTrees.jl b/ext/NminClusterTrees/NminClusterTrees.jl
new file mode 100644
index 0000000..8a4551b
--- /dev/null
+++ b/ext/NminClusterTrees/NminClusterTrees.jl
@@ -0,0 +1,9 @@
+module NminClusterTrees
+
+using StaticArrays
+using ClusterTrees
+import ClusterTrees.BlockTrees: BlockTree
+
+include("nmintree.jl")
+
+end
\ No newline at end of file
diff --git a/ext/NminClusterTrees/nmintree.jl b/ext/NminClusterTrees/nmintree.jl
new file mode 100644
index 0000000..42c1ebd
--- /dev/null
+++ b/ext/NminClusterTrees/nmintree.jl
@@ -0,0 +1,358 @@
+module NminTrees
+
+using StaticArrays
+using ParallelKMeans
+using LinearAlgebra
+using ClusterTrees
+
+struct Data{I, D, F, T}
+    sector::I
+    ct::SVector{D, F}
+    hs::F
+    values::Vector{T}
+end
+
+struct NminTree{D} <: ClusterTrees.PointerBasedTrees.APBTree
+    nodes::Vector{ClusterTrees.LevelledTrees.HNode{D}}
+    root::Int
+    num_elements::Int
+    levels::Vector{Int}
+end
+
+function NminTree(num_elements; center=SVector(0.0, 0.0, 0.0), radius=0.0)
+    root = ClusterTrees.LevelledTrees.HNode(
+        ClusterTrees.PointerBasedTrees.Node(Data(0, center, radius, Int[]), 0, 0, 0, 0), 0
+    )
+
+    return NminTree([root], 1, num_elements, Int[1])
+end
+
+ClusterTrees.root(tree::NminTree{D}) where {D} = tree.root
+ClusterTrees.data(tree::NminTree{D}, node) where {D} = tree.nodes[node].node.data
+ClusterTrees.parent(tree::NminTree{D}, node) where {D} = tree.nodes[node].node.parent
+ClusterTrees.PointerBasedTrees.nextsibling(tree::NminTree, node) =
+    tree.nodes[node].node.next_sibling
+ClusterTrees.PointerBasedTrees.firstchild(tree::NminTree, node) =
+    tree.nodes[node].node.first_child
+
+function value(tree, node::Int)
+    if !ClusterTrees.haschildren(tree, node)
+        return tree.nodes[node].node.data.values
+    else
+        values = Int[]
+        for leave in ClusterTrees.leaves(tree, node)
+            append!(values, tree.nodes[leave].node.data.values)
+        end
+        return values
+    end
+end
+
+# BoxTree
+struct BoxTreeOptions{I}
+    nchildren::I
+end
+
+function BoxTreeOptions(; dim=3)
+    return BoxTreeOptions(2^dim)
+end
+
+function boundingbox(points::Vector{SVector{D, F}}) where {D, F}
+    min_dim = Vector(points[1])
+    max_dim = Vector(points[1])
+
+    for i = 1:length(points)
+        for j = 1:D
+            min_dim[j] =  min_dim[j] < points[i][j] ? min_dim[j] : points[i][j]
+            max_dim[j] =  max_dim[j] > points[i][j] ? max_dim[j] : points[i][j]
+        end
+    end
+
+    center = @MVector zeros(D)
+
+    length_dim = zeros(F, D)
+    for j = 1:D
+        length_dim[j] = max_dim[j] - min_dim[j]
+        center[j] = (max_dim[j] + min_dim[j])/F(2.0)
+    end
+
+    halflength = maximum(length_dim)/F(2.0)
+    
+    return halflength, SVector(center)
+end
+
+function sector(pt::SVector{3, F}, ct::SVector{3, F}) where F
+    bl = pt .> ct
+    sc = sum(b ? 2^(i-1) : 0 for (i,b) in enumerate(bl))
+    return sc
+end
+
+function center(sibling::I, parent_ct::SVector{3, F}, hs::F) where {I, F}
+
+    ct = [
+        parent_ct + SVector(-hs, -hs, -hs),
+        parent_ct + SVector(hs, -hs, -hs),
+        parent_ct + SVector(-hs, hs, -hs),
+        parent_ct + SVector(hs, hs, -hs),
+        parent_ct + SVector(-hs, -hs, hs),
+        parent_ct + SVector(hs, -hs, hs),
+        parent_ct + SVector(-hs, hs, hs),
+        parent_ct + SVector(hs, hs, hs)
+    ]
+
+    return ct[sibling]
+end
+
+function clustering(
+    points::Vector{SVector{D, F}}, 
+    point_idcs::Vector{I},
+    parent_ct::SVector{D, F}, 
+    hs::F,
+    treeoptions::BoxTreeOptions{I}
+) where {I, D, F}
+    #computing true hs
+    hs = hs/sqrt(D)
+    sorted_point_idcs = zeros(Int, length(point_idcs) + 1, treeoptions.nchildren)
+    cluster = [sector(points[idc], parent_ct) for idc in point_idcs]
+    
+    for (index, value) in enumerate(cluster)
+        sorted_point_idcs[1, value+1] += 1
+        sorted_point_idcs[sorted_point_idcs[1, value+1] + 1, value+1] = point_idcs[index]
+    end
+
+    cts = [center(i, parent_ct, hs/2) for i = 1:treeoptions.nchildren]
+    #computing "radius" of cube
+    hss = sqrt(3)*[hs/2 for i = 1:treeoptions.nchildren]
+
+    return sorted_point_idcs, cts, hss
+
+end
+
+#KMeansTree
+struct KMeansSettings{I}
+    max_iters::I
+    n_threads::I
+end
+
+function KMeansSettings(; max_iters=100, n_threads=1)
+    return KMeansSettings(max_iters, n_threads)
+end
+
+struct KMeansTreeOptions{I}
+    nchildren::I
+    settings::KMeansSettings{I}
+end
+
+function KMeansTreeOptions(; nchildren=2, settings=KMeansSettings())
+    return KMeansTreeOptions(nchildren, settings)
+end
+
+function clustering(
+    points::Matrix{F}, 
+    point_idcs::Vector{I},
+    parent_ct::SVector{3, F}, 
+    hs::F,
+    treeoptions::KMeansTreeOptions{I}
+) where {I, F}
+
+    kmcluster = ParallelKMeans.kmeans(
+        points[:, point_idcs],
+        treeoptions.nchildren;
+        max_iters=treeoptions.settings.max_iters,
+        n_threads=treeoptions.settings.n_threads,
+    )
+
+    sorted_point_idcs = zeros(Int, length(point_idcs) + 1, treeoptions.nchildren)
+
+    for (index, value) in enumerate(kmcluster.assignments)
+        sorted_point_idcs[1, value] += 1
+        sorted_point_idcs[sorted_point_idcs[1, value] + 1, value] = point_idcs[index]
+    end
+
+    cts = [SVector{3,F}(kmcluster.centers[:, i]) for i = 1:treeoptions.nchildren]
+    hss = [maximum(
+        norm.(
+            eachcol(
+                points[:, sorted_point_idcs[2:(sorted_point_idcs[1, i] + 1), i]] .-
+                kmcluster.centers[:, i],
+            )
+        ),
+    ) for i = 1:treeoptions.nchildren]
+
+    return sorted_point_idcs, cts, hss
+
+end
+
+
+function child!(
+    tree::NminTree{D}, treeoptions, state, destination
+) where {D}
+    nmin, maxlevel = destination
+    parent_node_idx, ct, hs, level, sibling_idx, point_idcs, points = state
+    sorted_point_idcs, cts, hss = clustering(
+        points, 
+        point_idcs, 
+        ct,
+        hs,
+        treeoptions
+    )
+    ## check strong nmin
+    isnotnmin = (
+        sorted_point_idcs[1, 1] >= nmin || sorted_point_idcs[1, 1] == 0
+    ) && maximum(sorted_point_idcs[1, :]) != 0
+    for child in 2:treeoptions.nchildren
+        isnotnmin = isnotnmin && (
+            sorted_point_idcs[1, child] >= nmin || sorted_point_idcs[1, child] == 0
+        )
+    end
+    if !isnotnmin || level > maxlevel
+        append!(
+            tree.nodes[parent_node_idx].node.data.values,
+            point_idcs
+        )
+
+        return 0
+    end
+   
+    if sorted_point_idcs[1, 1] == 0
+        state_sibling = (
+            parent_node_idx,
+            cts,
+            hss,
+            level,
+            sibling_idx + 1,
+            sorted_point_idcs,
+            points
+        )
+        tree.nodes[parent_node_idx].node.first_child = sibling!(
+            tree, treeoptions, state_sibling, destination
+        )
+    else
+        push!(
+            tree.nodes,
+            ClusterTrees.LevelledTrees.HNode(
+                ClusterTrees.PointerBasedTrees.Node(
+                    Data(sibling_idx, cts[sibling_idx], hss[sibling_idx], Int[]), 
+                    treeoptions.nchildren, 
+                    0, 
+                    parent_node_idx, 
+                    0
+                ),
+                level,
+            ),
+        )
+
+        node_idx = length(tree.nodes)
+        level >= length(tree.levels) && resize!(tree.levels, level+1)
+        tree.levels[level+1] = node_idx
+        tree.nodes[parent_node_idx].node.first_child = node_idx
+        state_sibling = (
+            parent_node_idx,
+            cts,
+            hss,
+            level,
+            sibling_idx + 1,
+            sorted_point_idcs,
+            points
+        )
+        tree.nodes[node_idx].node.next_sibling = sibling!(tree, treeoptions, state_sibling, destination)
+
+        # Check if more than one node is left
+        if sorted_point_idcs[1, 1] == 1
+            append!(
+                tree.nodes[node_idx].node.data.values,
+                sorted_point_idcs[(2:(sorted_point_idcs[1, 1] + 1)), 1]
+            )
+
+            return node_idx
+        else
+            state_child = (
+                node_idx, 
+                cts[sibling_idx], 
+                hss[sibling_idx], 
+                level + 1, 
+                1, 
+                sorted_point_idcs[(2:(sorted_point_idcs[1, 1] + 1)), 1],
+                points 
+            )
+            tree.nodes[node_idx].node.first_child = child!(tree, treeoptions, state_child, destination)
+
+            return node_idx
+        end
+    end
+end
+
+function sibling!(tree::NminTree{D}, treeoptions, state, destination) where {D}
+    nmin, maxlevel = destination
+    parent_node_idx, cts, hss, level, sibling_idx, sorted_point_idcs, points = state
+    # Enough siblings?
+    sibling_idx > treeoptions.nchildren && return 0
+
+    if sorted_point_idcs[1, sibling_idx] == 0
+        
+        state_sibling = (
+            parent_node_idx,
+            cts, 
+            hss,
+            level,
+            sibling_idx + 1,
+            sorted_point_idcs,
+            points
+        )
+
+        return sibling!(tree, treeoptions, state_sibling, destination)
+    end
+
+    push!(
+        tree.nodes,
+        ClusterTrees.LevelledTrees.HNode(
+            ClusterTrees.PointerBasedTrees.Node(
+                Data(sibling_idx, cts[sibling_idx], hss[sibling_idx], Int[]),
+                treeoptions.nchildren,
+                0, 
+                parent_node_idx, 
+                0
+            ),
+            level,
+        ),
+    )
+
+    node_idx = length(tree.nodes)
+    level >= length(tree.levels) && resize!(tree.levels, level+1)
+    tree.levels[level+1] = node_idx
+
+    state_sibling = (
+        parent_node_idx,
+        cts, 
+        hss,
+        level,
+        sibling_idx + 1,
+        sorted_point_idcs,
+        points
+    )
+    tree.nodes[node_idx].node.next_sibling = sibling!(tree, treeoptions, state_sibling, destination)
+
+    # Check if more than one node is left
+    if sorted_point_idcs[1, sibling_idx] == 1
+        append!(
+            tree.nodes[node_idx].node.data.values,
+            sorted_point_idcs[(2:(sorted_point_idcs[1, sibling_idx] + 1)), sibling_idx]
+        )
+
+        return node_idx
+    else
+        state_child = (
+            node_idx, 
+            cts[sibling_idx], 
+            hss[sibling_idx], 
+            level + 1, 
+            1,
+            sorted_point_idcs[(2:(sorted_point_idcs[1, sibling_idx] + 1)), sibling_idx], 
+            points
+        )
+        tree.nodes[node_idx].node.first_child = child!(tree, treeoptions, state_child, destination)
+
+        return node_idx
+    end
+end
+
+end
\ No newline at end of file
diff --git a/src/FastBEAST.jl b/src/FastBEAST.jl
index 4c14f5f..440f42a 100644
--- a/src/FastBEAST.jl
+++ b/src/FastBEAST.jl
@@ -1,6 +1,11 @@
 module FastBEAST
 using LinearAlgebra
-include("tree/tree.jl")
+
+include("../ext/NminClusterTrees/NminClusterTrees.jl")
+
+include("clustertree/tree.jl")
+include("clustertree/kmeanstree.jl")
+include("clustertree/boxtree.jl")
 
 include("aca/aca_utils.jl")
 include("aca/pivoting.jl")
@@ -13,17 +18,11 @@ include("fmm.jl")
 include("beast.jl")
 include("fmm/operators/FMMoperator.jl")
 
-export BoundingBox
-export getboxframe
-export getchildbox
-export whichchildbox
-
-export BoxTreeNode
-export create_tree
-export BoxTreeOptions
-export ExaFMMOptions
 export KMeansTreeOptions
-export KMeansTreeNode
+export BoxTreeOptions
+export create_tree
+export computeinteractions
+export value
 
 export aca, allocate_aca_memory
 export LazyMatrix
diff --git a/src/beast.jl b/src/beast.jl
index a05c4b5..ad94bca 100644
--- a/src/beast.jl
+++ b/src/beast.jl
@@ -9,7 +9,8 @@ function hassemble(
     quadstratcbk=BEAST.defaultquadstrat(operator, test_functions, trial_functions),
     quadstratfbk=BEAST.defaultquadstrat(operator, test_functions, trial_functions),
     multithreading=false,
-    verbose=false
+    verbose=false,
+    η=1.5
 )
 
     @views blkasm = BEAST.blockassembler(
@@ -18,7 +19,6 @@ function hassemble(
         trial_functions,
         quadstrat=quadstratfbk
     )
-    
     @views function assembler(Z, tdata, sdata)
         @views store(v,m,n) = (Z[m,n] += v)
         blkasm(tdata,sdata,store)
@@ -30,26 +30,24 @@ function hassemble(
         trial_functions,
         quadstrat=quadstratcbk
     )
-    
     @views function farassembler(Z, tdata, sdata)
         @views store(v,m,n) = (Z[m,n] += v)
         farblkasm(tdata,sdata,store)
     end
 
-
     test_tree = create_tree(test_functions.pos, treeoptions)
     trial_tree = create_tree(trial_functions.pos, treeoptions)
 
     @time hmat = HMatrix(
         assembler,
-        test_tree,
-        trial_tree,
+        ClusterTrees.BlockTrees.BlockTree(test_tree, trial_tree),
         Int64,
         scalartype(operator),
         farmatrixassembler=farassembler,
         compressor=compressor,
         multithreading=multithreading,
-        verbose=verbose
+        verbose=verbose,
+        η=η
     )
     return hmat
 end
@@ -94,7 +92,7 @@ function fmmassemble(
         0,
         multithreading
     )
-    
+
     testpoints, testqp = meshtopoints(test_functions, quadstratcbk.outer_rule)
     trialpoints, trialqp = meshtopoints(trial_functions, quadstratcbk.outer_rule)
 
@@ -107,8 +105,8 @@ function fmmassemble(
 
     return FMMMatrix(
         operator,
-        test_functions, 
-        trial_functions, 
+        test_functions,
+        trial_functions,
         testqp,
         trialqp,
         fmm,
@@ -197,7 +195,7 @@ struct SafeDoubleQuadRule{P,Q}
     inner_quad_points::Q
 end
 
-function BEAST.quadrule(op, tref, bref, i ,τ, j, σ, qd, qs::SafeDoubleNumQStrat)
+function BEAST.quadrule(op, tref, bref, i, τ, j, σ, qd, qs::SafeDoubleNumQStrat)
     return SafeDoubleQuadRule(
         qd.tpoints[1,i],
         qd.bpoints[1,j])
@@ -212,10 +210,10 @@ function BEAST.quaddata(
     qs::SafeDoubleNumQStrat
 )
 
-    test_eval(x)  = test_refspace(x)
+    test_eval(x) = test_refspace(x)
     trial_eval(x) = trial_refspace(x)
 
-    tpoints = BEAST.quadpoints(test_eval,  test_elements,  (qs.outer_rule,))
+    tpoints = BEAST.quadpoints(test_eval, test_elements, (qs.outer_rule,))
     bpoints = BEAST.quadpoints(trial_eval, trial_elements, (qs.inner_rule,))
 
     return (;tpoints, bpoints)
@@ -230,27 +228,25 @@ function BEAST.quaddata(
     qs::SafeDoubleNumQStrat
 )
 
-    test_eval(x)  = test_refspace(x)
+    test_eval(x) = test_refspace(x)
     trial_eval(x) = trial_refspace(x)
 
-    tpoints = BEAST.quadpoints(test_eval,  test_elements,  (qs.outer_rule,))
+    tpoints = BEAST.quadpoints(test_eval, test_elements, (qs.outer_rule,))
     bpoints = BEAST.quadpoints(trial_eval, trial_elements, (qs.inner_rule,))
 
     return (;tpoints, bpoints)
 end
 
 function BEAST.momintegrals!(biop, tshs, bshs, tcell, bcell, z, strat::SafeDoubleQuadRule)
-    
+
     igd = BEAST.Integrand(biop, tshs, bshs, tcell, bcell)
     womps = strat.outer_quad_points
     wimps = strat.inner_quad_points
-    
     for womp in womps
         tgeo = womp.point
         tvals = womp.value
         M = length(tvals)
         jx = womp.weight
-        
         for wimp in wimps
             bgeo = wimp.point
             bvals = wimp.value
@@ -263,8 +259,10 @@ function BEAST.momintegrals!(biop, tshs, bshs, tcell, bcell, z, strat::SafeDoubl
                 z1 = j * igd(tgeo, bgeo, tvals, bvals)
                 for n in 1:N
                     for m in 1:M
-                        z[m,n] += z1[m,n]
-            end end end
+                        z[m, n] += z1[m, n]
+                    end
+                end
+            end
         end
     end
 
diff --git a/src/clustertree/boxtree.jl b/src/clustertree/boxtree.jl
new file mode 100644
index 0000000..c9af39f
--- /dev/null
+++ b/src/clustertree/boxtree.jl
@@ -0,0 +1,20 @@
+using ClusterTrees
+using LinearAlgebra
+
+function create_tree(
+    points::Vector{SVector{D, F}},
+    options::BoxTreeOptions{I}
+) where {D, I, F <: Real}
+
+    hs, ct = NminClusterTrees.NminTrees.boundingbox(points)
+    tree = NminClusterTrees.NminTrees.NminTree(length(points), center=ct, radius=F(sqrt(D))*hs)
+
+    D == 3 ? nchildren=8 : nchildren=4
+    treeoptions = NminClusterTrees.NminTrees.BoxTreeOptions(dim=D)
+    destination = (options.nmin, options.maxlevel)
+    state = (1, ct, F(sqrt(D))*hs, 1, 1, Vector(1:length(points)), points)
+
+    NminClusterTrees.NminTrees.child!(tree, treeoptions, state, destination)
+
+    return tree
+end
diff --git a/src/clustertree/kmeanstree.jl b/src/clustertree/kmeanstree.jl
new file mode 100644
index 0000000..7a74353
--- /dev/null
+++ b/src/clustertree/kmeanstree.jl
@@ -0,0 +1,23 @@
+using ClusterTrees
+using LinearAlgebra
+
+function create_tree(
+    points::Vector{SVector{D, F}},
+    options::KMeansTreeOptions{I}
+) where {I, D, F <: Real}
+    hs, ct = NminClusterTrees.NminTrees.boundingbox(points)
+    pointsM = reshape([point[i] for point in points for i in 1:3], (3, length(points)))
+
+    tree = NminClusterTrees.NminTrees.NminTree(length(points), center=ct, radius=F(sqrt(D))*hs)
+
+    treeoptions = NminClusterTrees.NminTrees.KMeansTreeOptions(
+        nchildren=options.nchildren,
+        settings=options.KMeansSettings
+    )
+    destination = (options.nmin, options.maxlevel)
+    state = (1, ct, F(sqrt(D))*hs, 1, 1, Vector(1:length(points)), pointsM)
+
+    NminClusterTrees.NminTrees.child!(tree, treeoptions, state, destination)
+
+    return tree
+end
diff --git a/src/clustertree/tree.jl b/src/clustertree/tree.jl
new file mode 100644
index 0000000..dc8cd9a
--- /dev/null
+++ b/src/clustertree/tree.jl
@@ -0,0 +1,90 @@
+using ClusterTrees
+using StaticArrays
+
+struct BoxTreeOptions{I}
+    nmin::I
+    maxlevel::I
+end
+
+function BoxTreeOptions(; nmin=1, maxlevel=10)
+
+    return BoxTreeOptions(nmin, maxlevel)
+end
+
+struct KMeansTreeOptions{I}
+    nmin::I
+    maxlevel::I
+    nchildren::I
+    KMeansSettings::NminClusterTrees.NminTrees.KMeansSettings{I}
+end
+
+function KMeansTreeOptions(; 
+    nmin=1, 
+    maxlevel=10,
+    nchildren=2,
+    KMeansSettings=NminClusterTrees.NminTrees.KMeansSettings()
+)
+
+    return KMeansTreeOptions(nmin, maxlevel, nchildren, KMeansSettings)
+end
+
+
+function isfar(
+    testnode::ClusterTrees.LevelledTrees.HNode{D}, 
+    trialnode::ClusterTrees.LevelledTrees.HNode{D};
+    η=2.0
+) where D
+    test_center = testnode.node.data.ct
+    test_radius =  testnode.node.data.hs
+    trial_center = trialnode.node.data.ct
+    trial_radius = trialnode.node.data.hs
+    center_dist = norm(test_center - trial_center)
+    dist = center_dist - (test_radius + trial_radius)
+    
+    if 2 * max(test_radius, trial_radius) < η * (dist) && dist != 0
+        return true
+    else
+        return false
+    end
+end
+
+
+function listnearfarinteractions(
+    block_tree::ClusterTrees.BlockTrees.BlockTree{T},
+    block,
+    state,
+    nears::Vector{Tuple{Int,Int}},
+    fars::Vector{Vector{Tuple{Int,Int}}},
+    level::Int;
+    η=1.5
+) where {T}
+    isfar(state[1], state[2], η=η) && (push!(fars[level], block); return nothing)
+    !ClusterTrees.haschildren(block_tree, block) && (push!(nears, block); return nothing)
+    for chd in ClusterTrees.children(block_tree, block)
+        children = ClusterTrees.data(block_tree, chd)
+        chd_state = (
+            block_tree.test_cluster.nodes[chd[1]],
+            block_tree.trial_cluster.nodes[chd[2]]
+        )
+        listnearfarinteractions(block_tree, chd, chd_state, nears, fars, level + 1, η=η)
+    end
+end
+
+
+function computeinteractions(tree::ClusterTrees.BlockTrees.BlockTree{T}; η=1.5) where {T}
+    nears = Tuple{Int,Int}[]
+    num_levels = length(tree.test_cluster.levels)
+    fars = [Tuple{Int,Int}[] for l in 1:num_levels]
+
+    root_state = (tree.test_cluster.nodes[1], tree.trial_cluster.nodes[1])
+    root_level = 1
+
+    listnearfarinteractions(
+        tree, ClusterTrees.root(tree), root_state, nears, fars, root_level, η=η
+    )
+    return nears, fars
+end
+
+function value(tree::NminClusterTrees.NminTrees.NminTree{D}, node::I) where {I, D}
+    return NminClusterTrees.NminTrees.value(tree, node)
+end
diff --git a/src/fmm.jl b/src/fmm.jl
index 5568674..4364752 100644
--- a/src/fmm.jl
+++ b/src/fmm.jl
@@ -126,18 +126,11 @@ function getfullrankblocks(
     correctionblocks = MBF[]
     correctionblocks_perthread = Vector{MBF}[]
 
-    test_tree = create_tree(test_functions.pos, treeoptions)
-    trial_tree = create_tree(trial_functions.pos, treeoptions)
-
-    fullinteractions = SVector{2}[]
-    compressableinteractions = SVector{2}[]
-
-    FastBEAST.computerinteractions!(
-        test_tree,
-        trial_tree,
-        fullinteractions,
-        compressableinteractions
+    tree = ClusterTrees.BlockTrees.BlockTree(
+        create_tree(test_functions.pos, treeoptions),
+        create_tree(trial_functions.pos, treeoptions)
     )
+    fullinteractions, compressableinteractions = computeinteractions(tree)
 
     if !multithreading
         for fullinteraction in fullinteractions
@@ -145,8 +138,8 @@ function getfullrankblocks(
                 fullrankblocks,
                 FastBEAST.getfullmatrixview(
                     assembler,
-                    fullinteraction[1],
-                    fullinteraction[2],
+                    value(tree.test_cluster, fullinteraction[1]),
+                    value(tree.trial_cluster, fullinteraction[2]),
                     Int,
                     scalartype(operator)
                 )
@@ -155,8 +148,8 @@ function getfullrankblocks(
                 correctionblocks,
                 FastBEAST.getfullmatrixview(
                     corassembler,
-                    fullinteraction[1],
-                    fullinteraction[2],
+                    value(tree.test_cluster, fullinteraction[1]),
+                    value(tree.trial_cluster, fullinteraction[2]),
                     Int,
                     scalartype(operator)
                 )
@@ -173,8 +166,8 @@ function getfullrankblocks(
                 fullrankblocks_perthread[Threads.threadid()],
                 FastBEAST.getfullmatrixview(
                     assembler,
-                    fullinteraction[1],
-                    fullinteraction[2],
+                    value(tree.test_cluster, fullinteraction[1]),
+                    value(tree.trial_cluster, fullinteraction[2]),
                     Int,
                     scalartype(operator)
                 )
@@ -183,8 +176,8 @@ function getfullrankblocks(
                 correctionblocks_perthread[Threads.threadid()],
                 FastBEAST.getfullmatrixview(
                     corassembler,
-                    fullinteraction[1],
-                    fullinteraction[2],
+                    value(tree.test_cluster, fullinteraction[1]),
+                    value(tree.trial_cluster, fullinteraction[2]),
                     Int,
                     scalartype(operator)
                 )
diff --git a/src/hmatrix.jl b/src/hmatrix.jl
index 3446e49..ebbc66b 100644
--- a/src/hmatrix.jl
+++ b/src/hmatrix.jl
@@ -202,33 +202,24 @@ end
 
 function HMatrix(
     matrixassembler::Function,
-    testtree::T,#::BoxTreeNode,
-    sourcetree::T,#::BoxTreeNode,
+    tree::ClusterTrees.BlockTrees.BlockTree{T},
     ::Type{I},
     ::Type{K};
     farmatrixassembler=matrixassembler,
     compressor=ACAOptions(),
     multithreading=false,
-    verbose=false
-) where {I, K, T <: AbstractNode}
+    verbose=false,
+    η=0.5
+) where {I, K, T}
     
-    fullinteractions = SVector{2}[]
-    compressableinteractions = SVector{2}[]
     fullcompressableinteractions = SVector{2}[]
-    
-    computerinteractions!(
-        testtree,
-        sourcetree,
-        fullinteractions,
-        compressableinteractions
-    )
-
+    fullinteractions, compressableinteractions = computeinteractions(tree, η=η)
     MBF = MatrixBlock{I, K, Matrix{K}}
     fullrankblocks_perthread = Vector{MBF}[]
     fullrankblocks = MBF[]
 
-    rowdim = numindices(testtree)
-    coldim = numindices(sourcetree)
+    rowdim = length(value(tree.test_cluster, 1))
+    coldim = length(value(tree.trial_cluster, 1))
 
     if !multithreading
         am = allocate_aca_memory(K, rowdim, coldim, maxrank=compressor.maxrank)
@@ -238,8 +229,7 @@ function HMatrix(
             push!(ams, allocate_aca_memory(K, rowdim, coldim, maxrank=compressor.maxrank))
         end
     end
-    nonzeros_perthread = I[]
-    nonzeros = 0
+
     verbose && println("Number of full interactions: ", length(fullinteractions))
     verbose && println(
         "Number of compressable interactions: ",
@@ -252,13 +242,12 @@ function HMatrix(
 
     if !multithreading
         for fullinteraction in fullinteractions
-            nonzeros += numindices(fullinteraction[1])*numindices(fullinteraction[2])
             push!(
                 fullrankblocks,
                 getfullmatrixview(
                     matrixassembler,
-                    fullinteraction[1],
-                    fullinteraction[2],
+                    value(tree.test_cluster, fullinteraction[1]),
+                    value(tree.trial_cluster, fullinteraction[2]),
                     I,
                     K
                 )
@@ -268,18 +257,15 @@ function HMatrix(
     elseif multithreading
         for i in 1:Threads.nthreads()
             push!(fullrankblocks_perthread, MBF[])
-            push!(nonzeros_perthread, 0)
         end
 
         Threads.@threads for fullinteraction in fullinteractions
-            nonzeros_perthread[Threads.threadid()] += 
-                numindices(fullinteraction[1])*numindices(fullinteraction[2])
             push!(
                 fullrankblocks_perthread[Threads.threadid()],
                 getfullmatrixview(
                     matrixassembler,
-                    fullinteraction[1],
-                    fullinteraction[2],
+                    value(tree.test_cluster, fullinteraction[1]),
+                    value(tree.trial_cluster, fullinteraction[2]),
                     I,
                     K
                 )
@@ -301,43 +287,46 @@ function HMatrix(
     end
 
     if !multithreading
-        for compressableinteraction in compressableinteractions
-            push!(
-                lowrankblocks, 
-                getcompressedmatrix(
-                    farmatrixassembler,
-                    compressableinteraction[1],
-                    compressableinteraction[2],
-                    I,
-                    K,
-                    am,
-                    compressor=compressor
+        for level in compressableinteractions
+            for compressableinteraction in level
+                push!(
+                    lowrankblocks, 
+                    getcompressedmatrix(
+                        farmatrixassembler,
+                        value(tree.test_cluster, compressableinteraction[1]),
+                        value(tree.trial_cluster, compressableinteraction[2]),
+                        I,
+                        K,
+                        am,
+                        compressor=compressor
+                    )
                 )
-            )
-            nonzeros += nnzs(lowrankblocks[end])
-            verbose && next!(p)
+                #nonzeros += nnzs(lowrankblocks[end])
+                verbose && next!(p)
+            end
         end
     elseif multithreading
         for i in 1:Threads.nthreads()
             push!(lowrankblocks_perthread, MBL[])
         end
-
-        Threads.@threads for compressableinteraction in compressableinteractions
-            push!(
-                lowrankblocks_perthread[Threads.threadid()],
-                getcompressedmatrix(
-                    farmatrixassembler,
-                    compressableinteraction[1],
-                    compressableinteraction[2],
-                    I,
-                    K,
-                    ams[Threads.threadid()],
-                    compressor=compressor
+        for level in compressableinteractions
+            Threads.@threads for compressableinteraction in level
+                push!(
+                    lowrankblocks_perthread[Threads.threadid()],
+                    getcompressedmatrix(
+                        farmatrixassembler,
+                        value(tree.test_cluster, compressableinteraction[1]),
+                        value(tree.trial_cluster, compressableinteraction[2]),
+                        I,
+                        K,
+                        ams[Threads.threadid()],
+                        compressor=compressor
+                    )
                 )
-            )
-            nonzeros_perthread[Threads.threadid()] += 
-                nnzs(lowrankblocks_perthread[Threads.threadid()][end])
-            verbose && next!(p)
+                #nonzeros_perthread[Threads.threadid()] += 
+                #    nnzs(lowrankblocks_perthread[Threads.threadid()][end])
+                verbose && next!(p)
+            end
         end
 
         for i in eachindex(lowrankblocks_perthread)
@@ -408,97 +397,27 @@ function HMatrix(
         lowrankblocks,
         rowdim,
         coldim,
-        nonzeros + sum(nonzeros_perthread),
+        0,
         compressor.maxrank,
         multithreading
     )
 end
 
-function computerinteractions!(
-    testnode::T,
-    sourcenode::T,
-    fullinteractions::Vector{SVector{2}},
-    compressableinteractions::Vector{SVector{2}}
-) where T <: AbstractNode
-    if iscompressable(sourcenode, testnode)
-        if level(sourcenode) == 0 && level(testnode) == 0
-            push!(compressableinteractions, SVector(testnode, sourcenode))
-            return
-        else
-            error("We do not expect this behavior")
-        end
-    end
-
-    if !haschildren(sourcenode) && !haschildren(testnode)
-        push!(fullinteractions, SVector(testnode, sourcenode))
-        return
-    else
-        if !haschildren(sourcenode)
-            schild = sourcenode
-            for tchild in children(testnode)
-                decide_compression(
-                    tchild, 
-                    schild, 
-                    fullinteractions, 
-                    compressableinteractions
-                )
-            end
-        elseif !haschildren(testnode)
-            tchild = testnode
-            for schild in children(sourcenode)
-                decide_compression(
-                    tchild, 
-                    schild, 
-                    fullinteractions, 
-                    compressableinteractions
-                )
-            end
-        else
-            for schild in children(sourcenode)
-                for tchild in children(testnode)
-                    decide_compression(
-                        tchild, 
-                        schild, 
-                        fullinteractions, 
-                        compressableinteractions
-                    )
-                end
-            end
-        end
-    end
-end
-
-function decide_compression(ttchild, sschild, fullinteractions, compressableinteractions)
-    if numindices(ttchild) > 0 && numindices(sschild) > 0
-        if iscompressable(ttchild, sschild)
-            push!(compressableinteractions, SVector(ttchild, sschild))
-
-            return
-        else
-            computerinteractions!(
-                ttchild,
-                sschild,
-                fullinteractions,
-                compressableinteractions
-            )
-        end
-    end
-end
 
 function getfullmatrixview(
     matrixassembler,
-    testnode,
-    sourcenode,
+    testidcs,
+    trialidcs,
     ::Type{I},
     ::Type{K};
 ) where {I, K}
-    matrix = zeros(K, numindices(testnode), numindices(sourcenode))
-    matrixassembler(matrix, indices(testnode), indices(sourcenode))
+    matrix = zeros(K, length(testidcs), length(trialidcs))
+    matrixassembler(matrix, testidcs, trialidcs)
 
     return MatrixBlock{I, K, Matrix{K}}(
         matrix,
-        indices(testnode),
-        indices(sourcenode)
+        testidcs,
+        trialidcs
     )
 end
 
@@ -521,15 +440,15 @@ end
 
 function getcompressedmatrix(
     matrixassembler::Function,
-    testnode,
-    sourcenode,
+    testidcs,
+    trialidcs,
     ::Type{I},
     ::Type{K},
     am;
     compressor=ACAOptions()
 ) where {I, K}
 
-        lm = LazyMatrix(matrixassembler, indices(testnode), indices(sourcenode), K)
+        lm = LazyMatrix(matrixassembler, testidcs, trialidcs, K)
 
         maxrank = compressor.maxrank
         maxrank == 0 && Int(round(length(lm.τ)*length(lm.σ)/(length(lm.τ)+length(lm.σ))))
@@ -547,13 +466,9 @@ function getcompressedmatrix(
 
         mbl = MatrixBlock{I, K, LowRankMatrix{K}}(
             LowRankMatrix(U, V),
-            indices(testnode),
-            indices(sourcenode)
+            testidcs,
+            trialidcs
         )
 
     return mbl
 end
-
-##
-
-real(ComplexF64)
\ No newline at end of file
diff --git a/src/tree/boundingbox.jl b/src/tree/boundingbox.jl
deleted file mode 100644
index 183943a..0000000
--- a/src/tree/boundingbox.jl
+++ /dev/null
@@ -1,104 +0,0 @@
-using StaticArrays
-
-struct BoundingBox{D, F} 
-    halflength::F
-    center::SVector{D, F}
-end
-
-function BoundingBox(points::Vector{SVector{D, F}}) where {D, F}
-    min_dim = Vector(points[1])
-    max_dim = Vector(points[1])
-
-    for i = 1:length(points)
-        for j = 1:D
-            min_dim[j] =  min_dim[j] < points[i][j] ? min_dim[j] : points[i][j]
-            max_dim[j] =  max_dim[j] > points[i][j] ? max_dim[j] : points[i][j]
-        end
-    end
-
-    center = @MVector zeros(D)
-
-    length_dim = zeros(F, D)
-    for j = 1:D
-        length_dim[j] = max_dim[j] - min_dim[j]
-        center[j] = (max_dim[j] + min_dim[j])/F(2.0)
-    end
-
-    halflength = maximum(length_dim)/F(2.0)
-    
-    return BoundingBox(halflength, SVector(center))
-end
-
-function getboxframe(bbox::BoundingBox{D, F}) where {D, F}
-    if D == 2
-        mat = zeros(5, 2)
-        sign_x = [1 -1 -1 1 1]
-        sign_y = [1 1 -1 -1 1]
-        for i = 1:5
-            mat[i,1] = bbox.center[1] + sign_x[i]*bbox.halflength 
-            mat[i,2] = bbox.center[2] + sign_y[i]*bbox.halflength 
-        end
-    elseif D == 3
-        sign_x = [1 -1 -1  1  1  1  1  1 NaN  1  1 -1 -1 NaN -1 -1 -1 NaN -1  1]
-        sign_y = [1  1 -1 -1  1  1 -1 -1 NaN -1 -1 -1 -1 NaN -1  1  1 NaN  1  1]
-        sign_z = [1  1  1  1  1 -1 -1  1 NaN -1 -1 -1  1 NaN -1 -1  1 NaN -1 -1]
-        mat = zeros(length(sign_x), 3)
-        for i = 1:length(sign_x)
-            mat[i,1] = bbox.center[1] + sign_x[i]*bbox.halflength 
-            mat[i,2] = bbox.center[2] + sign_y[i]*bbox.halflength 
-            mat[i,3] = bbox.center[3] + sign_z[i]*bbox.halflength 
-        end
-    else
-        error("Dimension must be 2 or 3")
-    end
-    return mat
-end
-
-"""
-    getchildbox(bbox, id)
-
-Returns a bounding box that this is the `id`th quadrants (2D) /
-octand (3D) of the bounding box `bbox`
-"""
-function getchildbox(bbox::BoundingBox{D, F}, id::Integer) where {D, F}
-
-    hl = bbox.halflength / F(2.0)
-
-    if D == 2
-        sign_x = [1 -1 -1  1]
-        sign_y = [1  1 -1 -1]
-
-        center = bbox.center + hl .* SVector(sign_x[id], sign_y[id])
-    elseif D == 3
-        sign_x = [1 -1 -1  1  1 -1 -1  1]
-        sign_y = [1  1 -1 -1  1  1 -1 -1]
-        sign_z = [1  1  1  1 -1 -1 -1 -1]
-        center = bbox.center + hl .* SVector(sign_x[id], sign_y[id],  sign_z[id])
-    else
-        error("We support only 2 or 3 dimensions")
-    end
-    return BoundingBox(hl, center)
-end
-
-function whichchildbox(bbox::BoundingBox{D, F}, point::SVector{D, F}) where {D, F}
-
-    if point[1] >= bbox.center[1]
-        if  point[2] >= bbox.center[2]
-            id = 1
-        else 
-            id = 4
-        end
-    else
-        if point[2] >= bbox.center[2]
-            id = 2
-        else 
-            id = 3
-        end
-    end
-
-    if D == 3 && point[3] < bbox.center[3]
-        id += 4
-    end
-
-    return id
-end
\ No newline at end of file
diff --git a/src/tree/boxtree.jl b/src/tree/boxtree.jl
deleted file mode 100644
index 5a23bf8..0000000
--- a/src/tree/boxtree.jl
+++ /dev/null
@@ -1,98 +0,0 @@
-struct BoxTreeOptions <: TreeOptions
-    nmin
-    maxlevel
-end
-
-function BoxTreeOptions(; nmin=1, maxlevel=64)
-    return BoxTreeOptions(nmin, maxlevel)
-end
-
-mutable struct BoxTreeNode{D, I, F, N <: NodeData{I, F}} <: AbstractNode{I, F, N}
-    parent::Union{BoxTreeNode{D, I, F, N}, Nothing}
-    children::Vector{BoxTreeNode{D, I, F, N}}
-    level::I
-    boundingbox::BoundingBox{D, F}
-    data::N
-end
-
-function BoxTreeNode(bbox::BoundingBox{D, F}, data::N) where {D, I, F, N <: NodeData{I, F}}
-    return BoxTreeNode(nothing, BoxTreeNode{D, I, F, N}[], I(0), bbox, data)
-end
-
-function BoxTreeNode(
-    parent::BoxTreeNode{D, I, F, N},
-    bbox::BoundingBox{D, F},
-    level::I,
-    data::N
-) where {D, I, F, N <: NodeData{I, F}}
-
-    return BoxTreeNode(parent, BoxTreeNode{D, I, F, N}[], level, bbox, data)
-end
-
-function add_child!(
-    parent::BoxTreeNode{D, I, F, N},
-    bbox::BoundingBox{D, F},
-    data::N
-) where {D, I, F, N <: NodeData{I, F}}
-
-    childnode = BoxTreeNode(parent, bbox, level(parent) + 1, data)
-    push!(parent.children, childnode)    
-end
-
-function create_tree(
-    points::Vector{SVector{D, F}},
-    treeoptions::BoxTreeOptions
-) where {D, F}
-
-    root = BoxTreeNode(
-        BoundingBox(points),
-        ACAData(Vector(1:length(points)), F)
-    )
-
-    fill_tree!(root, points, nmin=treeoptions.nmin, maxlevel=treeoptions.maxlevel)
-    
-    return root
-end
-
-function fill_tree!(
-    node::BoxTreeNode{D, I, F, N},
-    points::Vector{SVector{D, F}};
-    nmin=1,
-    maxlevel=-log2(eps(eltype(points)))
-) where {D, I, F, N  <: NodeData{I, F}}
-
-    if numindices(node) <= nmin || level(node) >= maxlevel - 1
-        return
-    end
-
-    nchildren = D == 2 ? 4 : 8 # Number of (possible) children
-    sorted_points = zeros(I, numindices(node) + 1, nchildren)
-
-    for pindex in indices(node)
-        id = whichchildbox(node.boundingbox, points[pindex])
-        sorted_points[1, id] +=  1
-        sorted_points[1 + sorted_points[1, id], id] = pindex
-    end
-
-    for i = 1:nchildren
-        if sorted_points[1, i] > 0
-            childpointindices = sorted_points[2:(sorted_points[1, i]+1), i]
-            add_child!(node, getchildbox(node.boundingbox, i), ACAData(childpointindices, F))
-            fill_tree!(lastchild(node), points, nmin=nmin, maxlevel=maxlevel)
-        end
-    end
-end
-
-function iscompressable(sourcenode::BoxTreeNode, testnode::BoxTreeNode)
-    mindistance = 
-        sqrt(length(sourcenode.boundingbox.center))*sourcenode.boundingbox.halflength 
-    mindistance += 
-        sqrt(length(testnode.boundingbox.center))*testnode.boundingbox.halflength 
-
-    distance_centers = norm(sourcenode.boundingbox.center - testnode.boundingbox.center)
-    if distance_centers > 1.1*mindistance
-        return true
-    else 
-        return false
-    end
-end
diff --git a/src/tree/kmeans.jl b/src/tree/kmeans.jl
deleted file mode 100644
index c0865cb..0000000
--- a/src/tree/kmeans.jl
+++ /dev/null
@@ -1,285 +0,0 @@
-using ParallelKMeans
-
-"""
-    KMeansTreeOptions <: TreeOptions
-
-Is the datatype that describes which tree the [`create_tree`](@ref) function creates.
-
-# Fields
-- `iterations`: number of iterations on each level, default is one iterations
-- `nchildren`: defines the number of children of each node, default is two
-- `nmin`: defines the minimum amount of datapoints which are needed in a 
-    cluster so that it gets split up in subclusters, default is 1
-- `maxlevel`: defines the maximum amount of levels, default is 100.
-- `algorithm`: defines which algorithm is used. The :naive approach is not recommended.
-    Default is the wrapped ParallelKMeans algorithm.
-"""
-struct KMeansTreeOptions <: TreeOptions
-    iterations
-    nchildren
-    nmin
-    maxlevel
-    algorithm
-end
-
-function KMeansTreeOptions(;
-    iterations=20,
-    nchildren=2,
-    nmin=1,
-    maxlevel=100,
-    algorithm=:ParallelKMeans
-)
-    return KMeansTreeOptions(iterations, nchildren, nmin, maxlevel, algorithm)
-end
-
-"""
-    KMeansTreeNode{T} <: AbstractNode
-
-Is the datatype of each node in the [K-Means Clustering Tree](@ref).
-
-# Fields
-- `parent::Union{KMeansTreeNode{T},Nothing}`: is the superordinate cluster of
-    of the represented cluster
-- `children::Union{Vector{KMeansTreeNode{T}}, Nothing}`: all directly 
-    subordinated clusters of the represented cluster
-- `level::Integer`: the level of the represented cluster
-- `center`: the center by which the cluster is defined
-- `radius`: the euclidian distance between the center and the farthest away
-    point
-- `data::T`: array containig the indices of the points in this cluster
-"""
-mutable struct KMeansTreeNode{D, I, F, N <: NodeData{I, F}} <: AbstractNode{I, F, N}
-    parent::Union{KMeansTreeNode{D, I, F, N}, Nothing}
-    children::Vector{KMeansTreeNode{D, I, F, N}}
-    level::I
-    center::SVector{D, F}
-    radius::F
-    data::N
-end
-
-function KMeansTreeNode(
-    center::SVector{D, F},
-    radius::F,
-    data::N
-) where {D, I, F, N <: NodeData{I, F}}
-
-    return KMeansTreeNode(nothing, KMeansTreeNode{D, I, F, N}[], I(0), center, radius, data)
-end
-
-function add_child!(
-    parent::KMeansTreeNode{D, I, F, N},
-    center::SVector{D, F},
-    radius::F,
-    data::N
-) where {D, I, F, N}
-    
-    childnode = KMeansTreeNode(
-        parent,
-        KMeansTreeNode{D, I, F, N}[],
-        level(parent) + 1,
-        center,
-        radius,
-        data
-    )
-
-    push!(parent.children, childnode)
-
-end
-
-"""
-    create_tree(points::Array{SVector{D, T}, 1}; treeoptions)
-
-Creates an algebraic tree for an givn set of datapoints. The returned 
-datastructure is the foundation for the algorithms in this package. 
-
-# Arguments
-- `points::Array{SVector{D, T}, 1}`: is an array of 
-    [SVector](https://juliaarrays.github.io/StaticArrays.jl/stable/pages/api/#SVector-1). 
-    Each 
-    [SVector](https://juliaarrays.github.io/StaticArrays.jl/stable/pages/api/#SVector-1)
-    contains in general two or three float values, which discribe the position 
-    in the space.
-
-# Keywords
-- `treeoptions::TreeOptions`: this keyword defines by which tree is build. 
-    `TreeOptions` is an abstract type which either can be `BoxTreeOptions` or
-    [`KMeansTreeOptions`](@ref). Default type is `BoxTreeOptions`.
-
-# Returns
-- `AbstractNode`: the root of the created tree. AbstractNode is an abstract type 
-    which either can be `BoxTreeNode` or [`KMeansTreeNode`](@ref), depending on the keyword.
-"""
-function create_tree(
-    points::Vector{SVector{D, F}},
-    treeoptions::KMeansTreeOptions
-) where {D, F}
-    root = KMeansTreeNode(zeros(SVector{D, F}), F(0.0), ACAData(Vector(1:length(points)), F))
-    if treeoptions.algorithm != :naive
-        pointsM = reshape(
-            [points[j][i] for j = 1:length(points) for i = 1:D], 
-            (D,length(points))
-        )
-        fill_tree!(
-            root,
-            pointsM,
-            nmin=treeoptions.nmin,
-            maxlevel=treeoptions.maxlevel,
-            iterations=treeoptions.iterations,
-            nchildren=treeoptions.nchildren
-        )
-    else
-        fill_tree!(
-            root,
-            points,
-            nmin=treeoptions.nmin,
-            maxlevel=treeoptions.maxlevel,
-            iterations=treeoptions.iterations,
-            nchildren=treeoptions.nchildren
-        )
-    end
-
-    return root
-end
-
-function whichcenter(center, point, nchildren)
-    centerindex = argmin(norm.([center[i]-point for i = 1:nchildren]))
-
-    return centerindex
-end
-
-function fill_tree!(
-    node::KMeansTreeNode{D, I, F, N},
-    points::Vector{SVector{D, F}};
-    nmin=1,
-    maxlevel=log2(eps(eltype(points))),
-    iterations, 
-    nchildren
-) where {D, I, F, N <: NodeData{I, F}}
-
-    if numindices(node) <= nmin|| level(node) >= maxlevel - 1 || numindices(node) <= nchildren
-        return
-    end
-
-    center = points[indices(node)[1:nchildren]]
-    sorted_points = zeros(I, numindices(node)+1, nchildren)
-
-    for _ = 1:iterations
-        sorted_points .= 0
-        for pindex in indices(node)
-            id = whichcenter(center, points[pindex], nchildren)
-            sorted_points[1, id] +=  1
-            sorted_points[1 + sorted_points[1, id], id] = pindex
-        end
-        for k = 1:nchildren
-            center[k] = sum(
-                points[sorted_points[2:(sorted_points[1, k]+1), k]]
-            ) ./ sorted_points[1, k]          
-        end
-    end
-    
-
-    for i = 1:nchildren
-        if sorted_points[1, i] > 0
-            childpointindices = sorted_points[2:(sorted_points[1, i]+1), i]
-            radius = maximum(norm.([points[childpointindices[j]].-center[i] for j = 1:length(childpointindices)]))
-            add_child!(node, center[i], radius, ACAData(childpointindices, F))
-            fill_tree!(
-                lastchild(node), 
-                points, 
-                nmin=nmin, 
-                maxlevel=maxlevel, 
-                iterations=iterations, 
-                nchildren=nchildren
-            )
-        end
-    end
-end
-
-function fill_tree!(
-    node::KMeansTreeNode{D, I, F, N},
-    pointsM::Matrix;
-    nmin=1,
-    maxlevel=log2(eps(eltype(points))),
-    iterations, 
-    nchildren
-) where {D, I, F, N <: NodeData{I, F}}
-
-    if numindices(node) <= nmin || level(node) >= maxlevel - 1 || numindices(node) <= nchildren
-        return
-    end
-    
-    sorted_points = zeros(I, length(indices(node))+1, nchildren)
-    
-    # A first (probably too naive?) heuristic
-    opt_nthreads = max(1, floor(I, log10(0.001*length(indices(node)))))
-    
-    kmcluster = ParallelKMeans.kmeans(
-        pointsM[:,indices(node)],
-        nchildren,
-        max_iters=iterations,
-        n_threads=opt_nthreads
-    )
-
-    for (index, value) in enumerate(kmcluster.assignments)
-        sorted_points[1, value] += 1
-        sorted_points[sorted_points[1,value]+1, value] = indices(node)[index] 
-    end
-
-    for i = 1:nchildren
-        if sorted_points[1, i] > 0
-            center = SVector{D, F}([kmcluster.centers[j, i] for j = 1:D])
-            radius = maximum(norm.(eachcol(
-                pointsM[:, sorted_points[2:(sorted_points[1, i]+1),i]].-kmcluster.centers[:,i]
-            )))
-            
-            add_child!(
-                node,
-                center,
-                radius,
-                ACAData(sorted_points[2:(sorted_points[1, i]+1), i], F)
-            )
-            fill_tree!(
-                lastchild(node), 
-                pointsM, 
-                nmin=nmin, 
-                maxlevel=maxlevel, 
-                iterations=iterations, 
-                nchildren=nchildren
-            )
-        end
-    end
-end
-
-"""
-    iscompressable(sourcenode::AbstractNode, testnode::AbstractNode)
-
-Determins whether two nodes of a tree are comressable. The criteria differs 
-between the [Box Tree](@ref) and the [K-Means Clustering Tree](@ref).
-For the [K-Means Clustering Tree](@ref) two nodes can be compressed, if the distance 
-between the centers of two clusters is greater than the sum of their radius 
-multiplied by a factor of 1.5.
-For the [Box Tree](@ref) two nodes can be compressed, if the distance between the centers 
-of the two boxes is greater than the sum of the distances of each box's center to 
-one of its corners multiplied by a factor of 1.1.
-
-# Arguments
-- `sourcenode::AbstractNode`: the node which is observed
-- `testnode::AbstractNode`: the node which is tested for compression
-
-# Returns
-- `true`: if the input nodes are compressable
-- `false`: if the input nodes are not compressable
-"""
-function iscompressable(sourcenode::KMeansTreeNode, testnode::KMeansTreeNode)
-    if level(sourcenode) > 0 && level(testnode) > 0
-        dist = norm(sourcenode.center - testnode.center)
-        factor = 1.5
-        if factor * (sourcenode.radius + testnode.radius) < dist
-            return true
-        else
-            return false
-        end
-    else
-        return false
-    end
-end
\ No newline at end of file
diff --git a/src/tree/tree.jl b/src/tree/tree.jl
deleted file mode 100644
index 45d1f71..0000000
--- a/src/tree/tree.jl
+++ /dev/null
@@ -1,35 +0,0 @@
-abstract type TreeOptions end
-abstract type AbstractNode{I, F, N} end
-
-children(node::AbstractNode{I, F, N}) where {I, F, N} = node.children
-numchildren(node::AbstractNode{I, F, N}) where {I, F, N} = length(node.children)
-haschildren(node::AbstractNode{I, F, N}) where {I, F, N} = numchildren(node) > 0 ? true : false
-lastchild(node::AbstractNode{I, F, N}) where {I, F, N} = last(node.children)
-level(node::AbstractNode{I, F, N}) where {I, F, N} = node.level
-
-abstract type NodeData{I, F} end
-
-indices(node::AbstractNode{I, F, N}) where {I, F, N <: NodeData{I, F}} = node.data.indices
-numindices(node::AbstractNode{I, F, N}) where {I, F, N <: NodeData{I, F}} = length(node.data.indices)
-
-struct ACAData{I, F} <: NodeData{I, F}
-    indices::Vector{I}
-end
-
-function ACAData(indices::Vector{I}, ::Type{F}) where {I, F}
-    return ACAData{I, F}(indices)
-end
-
-
-## Box tree
-# Implements quad- (2D) and oct-tree (3D) cluster strategies
-include("boundingbox.jl")
-include("boxtree.jl")
-
-## KMeans-Tree
-include("kmeans.jl")
-
-## Default functions
-function create_tree(points::Vector{SVector{D,T}}) where {D,T}
-    create_tree(points, BoxTreeOptions())
-end
\ No newline at end of file
diff --git a/test/geometries/ihplate.msh b/test/geometries/ihplate.msh
deleted file mode 100644
index 91011ac..0000000
--- a/test/geometries/ihplate.msh
+++ /dev/null
@@ -1,515 +0,0 @@
-$MeshFormat
-2.2 0 8
-$EndMeshFormat
-$Nodes
-168
-1 0 0 0
-2 1 0 0
-3 1 1 0
-4 0 1 0
-5 0.04956611949498604 0 0
-6 0.1065026421555667 0 0
-7 0.1719055292342407 0 0
-8 0.2470337154254791 0 0
-9 0.3333333432581525 0 0
-10 0.4324655797766708 0 0
-11 0.54633861455472 0 0
-12 0.677144398366437 0 0
-13 0.8274007767559597 0 0
-14 1 0.1999999999995579 0
-15 1 0.3999999999989749 0
-16 1 0.5999999999989468 0
-17 1 0.7999999999994734 0
-18 0.8274007767565476 1 0
-19 0.6771443983685665 1 0
-20 0.5463386145579807 1 0
-21 0.4324655797796574 1 0
-22 0.3333333432602792 1 0
-23 0.2470337154269195 1 0
-24 0.1719055292351374 1 0
-25 0.1065026421560601 1 0
-26 0.04956611949520884 1 0
-27 0 0.9499999999997918 0
-28 0 0.8999999999995836 0
-29 0 0.8499999999996529 0
-30 0 0.7999999999999998 0
-31 0 0.7500000000003466 0
-32 0 0.7000000000006934 0
-33 0 0.6500000000010401 0
-34 0 0.6000000000013869 0
-35 0 0.5500000000017335 0
-36 0 0.5000000000020587 0
-37 0 0.4500000000018723 0
-38 0 0.4000000000016644 0
-39 0 0.3500000000014564 0
-40 0 0.3000000000012483 0
-41 0 0.2500000000010403 0
-42 0 0.2000000000008322 0
-43 0 0.1500000000006241 0
-44 0 0.100000000000416 0
-45 0 0.05000000000020799 0
-46 0.04424531645915544 0.4686718627588514 0
-47 0.04330127018940205 0.3250000000013523 0
-48 0.0422897351513597 0.625583810083459 0
-49 0.4887358272837651 0.9051295415043563 0
-50 0.4784607674643949 0.08896347777231939 0
-51 0.04189793283862304 0.7759147895303854 0
-52 0.2912238380755942 0.9356178313362966 0
-53 0.04330127018940215 0.1750000000007281 0
-54 0.2930059336325294 0.07082355624494717 0
-55 0.8333783413509464 0.5118322150020913 0
-56 0.03996558777469131 0.875510486988799 0
-57 0.1360837279573675 0.05396368201333093 0
-58 0.1392040856955987 0.9470254816406143 0
-59 0.7557115201487716 0.1045938732441594 0
-60 0.7364118459299529 0.8790679684677253 0
-61 0.04330127018940205 0.07500000000031198 0
-62 0.8609210098255249 0.6930026921431607 0
-63 0.6994234490858564 0.615778155682118 0
-64 0.6818424299013865 0.4524595692773895 0
-65 0.5635812573927423 0.5234653702720717 0
-66 0.5517119899898626 0.397124843415637 0
-67 0.4515400986709663 0.4650497250818806 0
-68 0.4529480507983337 0.5778299369202551 0
-69 0.3597237330259118 0.5283011678329612 0
-70 0.3602113985519547 0.425753762343069 0
-71 0.3588678974009183 0.6269108412131941 0
-72 0.2780383790500585 0.5836448005311916 0
-73 0.2774537340910092 0.6698895990832181 0
-74 0.3534953313490144 0.7208639504566052 0
-75 0.205060713625232 0.6295790993276607 0
-76 0.2050157439595385 0.7077634061178897 0
-77 0.2021136896498163 0.5493706972086492 0
-78 0.1427193975906282 0.6724951497746396 0
-79 0.1340836028776653 0.7492336422904551 0
-80 0.2008062442648995 0.7824861419860103 0
-81 0.1402681356592492 0.8324309206173038 0
-82 0.2026605241480611 0.8446778655357992 0
-83 0.2624286115914775 0.8191118699488236 0
-84 0.7798137397627123 0.3548972052680547 0
-85 0.08860073456565437 0.7045362548748569 0
-86 0.1335663049530804 0.5957184005558058 0
-87 0.1350641092934091 0.5218078092501804 0
-88 0.1909291085987491 0.476928806139122 0
-89 0.1234195695080011 0.4522422104769451 0
-90 0.1818549073712395 0.4072908513813704 0
-91 0.1232996544844075 0.378132627750419 0
-92 0.1748740167785929 0.3347806195033934 0
-93 0.238343843278716 0.3656774363510955 0
-94 0.2289503470592066 0.2916245763452675 0
-95 0.1640543913459814 0.2624179216751565 0
-96 0.2185941549595059 0.2184423196861265 0
-97 0.2910973306994607 0.2464026215464323 0
-98 0.1439343197562555 0.1901441874353661 0
-99 0.208759900115411 0.1470225628847392 0
-100 0.0977302751706286 0.2414963134538471 0
-101 0.08576089561655775 0.553883759291011 0
-102 0.1060679400264971 0.305930902860195 0
-103 0.432891019070482 0.3386015821658118 0
-104 0.5373418972478688 0.2831898875646527 0
-105 0.4472214429184811 0.2448452370153111 0
-106 0.07861089378280442 0.9477420238235605 0
-107 0.113636469173417 0.8963213990840344 0
-108 0.5564226926753878 0.6501393057021035 0
-109 0.6534428550135551 0.7399692361786002 0
-110 0.5291045749904877 0.7766292690074454 0
-111 0.3762256035008104 0.07971911723821956 0
-112 0.3214185183386679 0.1561495855229389 0
-113 0.5993542324961327 0.09641741209352594 0
-114 0.6660150767042516 0.2227042126338564 0
-115 0.6051963526907179 0.8834660025263545 0
-116 0.4071313165495598 0.8076533628491167 0
-117 0.3828994615199683 0.9211341292579437 0
-118 0.2138754513088441 0.06873721729222695 0
-119 0.2097427093968656 0.9437840815618203 0
-120 0.1433707202859374 0.1226379480462509 0
-121 0.2749598545675043 0.4856695953782524 0
-122 0.07301523636956589 0.4171539205203737 0
-123 0.3452770353755762 0.3241498154393106 0
-124 0.4429949772939502 0.6933377776914532 0
-125 0.04533250365652295 0.1241068787287346 0
-126 0.09248691379084623 0.1501302828916216 0
-127 0.04466511650702323 0.222624620704423 0
-128 0.04263092300395721 0.825237546086409 0
-129 0.03820783060534923 0.3778589867546785 0
-130 0.04862743364892517 0.274175306170351 0
-131 0.04319995447402262 0.6756095315432964 0
-132 0.08976825278422541 0.09432973980370381 0
-133 0.3224431882742128 0.8607616259991144 0
-134 0.07748814530736782 0.04176641767706905 0
-135 0.2699714733081209 0.7473781599499345 0
-136 0.08477781873073356 0.8527689859446232 0
-137 0.09002209628229121 0.6500000000010404 0
-138 0.08660828638839649 0.7999999999999995 0
-139 0.04173946325444503 0.5249715179161629 0
-140 0.03800811574648807 0.9243911385568326 0
-141 0.2354979963755805 0.4270443052982039 0
-142 0.1722519356655662 0.8939955418949054 0
-143 0.04285544320060777 0.7269248855218664 0
-144 0.5208005219786414 0.1832737103706639 0
-145 0.3265494547400781 0.7892280007427425 0
-146 0.08512228430785181 0.5002998728138686 0
-147 0.6433450267212164 0.342075143631918 0
-148 0.855764977902764 0.2224206634809558 0
-149 0.8638883267373199 0.8567755009470851 0
-150 0.04047518416390562 0.5759454699752808 0
-151 0.07842284323343894 0.6002262879813194 0
-152 0.243691243788148 0.8852573043866973 0
-153 0.08204075114896542 0.3500000000014564 0
-154 0.9004733815328366 0.1012528677848339 0
-155 0.2582543456613133 0.1757032551994773 0
-156 0.07880919997418939 0.7513219144435127 0
-157 0.4072642788203681 0.1658680791419216 0
-158 0.08420480891422415 0.1954140322891254 0
-159 0.03323702580490027 0.964426632476037 0
-160 0.03591252214609406 0.03581053893296295 0
-161 0.03088972592416208 0.4250000000017683 0
-162 0.368212652747697 0.2507870194298796 0
-163 0.2590628298113531 0.1236872354288659 0
-164 0.08228605498010633 0.4595813028103135 0
-165 0.7596951329427534 0.765546205341723 0
-166 0.2790141648949306 0.325501857384657 0
-167 0.07095065558680878 0.9030629219344249 0
-168 0.2907467549366979 0.3977401727241164 0
-$EndNodes
-$Elements
-338
-1 15 2 0 1 1
-2 15 2 0 2 2
-3 15 2 0 3 3
-4 15 2 0 4 4
-5 1 2 0 1 1 5
-6 1 2 0 1 5 6
-7 1 2 0 1 6 7
-8 1 2 0 1 7 8
-9 1 2 0 1 8 9
-10 1 2 0 1 9 10
-11 1 2 0 1 10 11
-12 1 2 0 1 11 12
-13 1 2 0 1 12 13
-14 1 2 0 1 13 2
-15 1 2 0 2 2 14
-16 1 2 0 2 14 15
-17 1 2 0 2 15 16
-18 1 2 0 2 16 17
-19 1 2 0 2 17 3
-20 1 2 0 3 3 18
-21 1 2 0 3 18 19
-22 1 2 0 3 19 20
-23 1 2 0 3 20 21
-24 1 2 0 3 21 22
-25 1 2 0 3 22 23
-26 1 2 0 3 23 24
-27 1 2 0 3 24 25
-28 1 2 0 3 25 26
-29 1 2 0 3 26 4
-30 1 2 0 4 4 27
-31 1 2 0 4 27 28
-32 1 2 0 4 28 29
-33 1 2 0 4 29 30
-34 1 2 0 4 30 31
-35 1 2 0 4 31 32
-36 1 2 0 4 32 33
-37 1 2 0 4 33 34
-38 1 2 0 4 34 35
-39 1 2 0 4 35 36
-40 1 2 0 4 36 37
-41 1 2 0 4 37 38
-42 1 2 0 4 38 39
-43 1 2 0 4 39 40
-44 1 2 0 4 40 41
-45 1 2 0 4 41 42
-46 1 2 0 4 42 43
-47 1 2 0 4 43 44
-48 1 2 0 4 44 45
-49 1 2 0 4 45 1
-50 2 2 0 6 107 81 142
-51 2 2 0 6 81 107 136
-52 2 2 0 6 110 109 115
-53 2 2 0 6 57 118 120
-54 2 2 0 6 109 60 115
-55 2 2 0 6 116 74 124
-56 2 2 0 6 110 116 124
-57 2 2 0 6 49 110 115
-58 2 2 0 6 84 114 148
-59 2 2 0 6 110 49 116
-60 2 2 0 6 114 84 147
-61 2 2 0 6 114 59 148
-62 2 2 0 6 118 99 120
-63 2 2 0 6 15 84 148
-64 2 2 0 6 100 102 130
-65 2 2 0 6 113 114 144
-66 2 2 0 6 75 78 86
-67 2 2 0 6 102 47 130
-68 2 2 0 6 92 91 102
-69 2 2 0 6 54 9 111
-70 2 2 0 6 55 62 63
-71 2 2 0 6 86 78 137
-72 2 2 0 6 109 63 165
-73 2 2 0 6 116 49 117
-74 2 2 0 6 57 7 118
-75 2 2 0 6 64 63 65
-76 2 2 0 6 114 104 144
-77 2 2 0 6 108 110 124
-78 2 2 0 6 25 26 106
-79 2 2 0 6 67 70 103
-80 2 2 0 6 66 65 67
-81 2 2 0 6 67 68 69
-82 2 2 0 6 81 82 142
-83 2 2 0 6 9 10 111
-84 2 2 0 6 7 8 118
-85 2 2 0 6 69 71 72
-86 2 2 0 6 55 63 64
-87 2 2 0 6 60 109 165
-88 2 2 0 6 72 73 75
-89 2 2 0 6 102 91 153
-90 2 2 0 6 120 98 126
-91 2 2 0 6 64 65 66
-92 2 2 0 6 77 75 86
-93 2 2 0 6 75 76 78
-94 2 2 0 6 87 86 101
-95 2 2 0 6 78 79 85
-96 2 2 0 6 66 67 103
-97 2 2 0 6 79 138 156
-98 2 2 0 6 70 69 121
-99 2 2 0 6 54 111 112
-100 2 2 0 6 55 16 62
-101 2 2 0 6 79 81 138
-102 2 2 0 6 116 117 133
-103 2 2 0 6 77 88 121
-104 2 2 0 6 69 68 71
-105 2 2 0 6 103 70 123
-106 2 2 0 6 15 55 84
-107 2 2 0 6 126 98 158
-108 2 2 0 6 67 65 68
-109 2 2 0 6 98 95 100
-110 2 2 0 6 72 71 73
-111 2 2 0 6 67 69 70
-112 2 2 0 6 95 92 102
-113 2 2 0 6 75 73 76
-114 2 2 0 6 77 86 87
-115 2 2 0 6 78 76 79
-116 2 2 0 6 14 15 148
-117 2 2 0 6 72 75 77
-118 2 2 0 6 97 112 162
-119 2 2 0 6 74 116 145
-120 2 2 0 6 12 59 113
-121 2 2 0 6 72 77 121
-122 2 2 0 6 63 62 165
-123 2 2 0 6 108 109 110
-124 2 2 0 6 113 59 114
-125 2 2 0 6 15 16 55
-126 2 2 0 6 127 100 130
-127 2 2 0 6 58 25 106
-128 2 2 0 6 104 103 105
-129 2 2 0 6 62 149 165
-130 2 2 0 6 100 95 102
-131 2 2 0 6 122 129 153
-132 2 2 0 6 16 17 62
-133 2 2 0 6 65 63 108
-134 2 2 0 6 96 95 98
-135 2 2 0 6 70 121 168
-136 2 2 0 6 10 11 50
-137 2 2 0 6 112 157 162
-138 2 2 0 6 50 144 157
-139 2 2 0 6 8 9 54
-140 2 2 0 6 10 50 111
-141 2 2 0 6 99 98 120
-142 2 2 0 6 91 89 122
-143 2 2 0 6 105 103 162
-144 2 2 0 6 79 80 81
-145 2 2 0 6 20 21 49
-146 2 2 0 6 22 23 52
-147 2 2 0 6 49 21 117
-148 2 2 0 6 88 87 89
-149 2 2 0 6 90 89 91
-150 2 2 0 6 6 7 57
-151 2 2 0 6 8 54 118
-152 2 2 0 6 94 92 95
-153 2 2 0 6 52 23 119
-154 2 2 0 6 111 50 157
-155 2 2 0 6 24 25 58
-156 2 2 0 6 86 137 151
-157 2 2 0 6 22 52 117
-158 2 2 0 6 103 123 162
-159 2 2 0 6 106 26 159
-160 2 2 0 6 21 22 117
-161 2 2 0 6 62 17 149
-162 2 2 0 6 24 58 119
-163 2 2 0 6 23 24 119
-164 2 2 0 6 20 49 115
-165 2 2 0 6 69 72 121
-166 2 2 0 6 12 13 59
-167 2 2 0 6 18 19 60
-168 2 2 0 6 19 20 115
-169 2 2 0 6 57 132 134
-170 2 2 0 6 98 100 158
-171 2 2 0 6 66 103 104
-172 2 2 0 6 108 63 109
-173 2 2 0 6 60 19 115
-174 2 2 0 6 57 120 132
-175 2 2 0 6 68 108 124
-176 2 2 0 6 47 102 153
-177 2 2 0 6 91 122 153
-178 2 2 0 6 116 133 145
-179 2 2 0 6 55 64 84
-180 2 2 0 6 6 57 134
-181 2 2 0 6 58 106 107
-182 2 2 0 6 96 98 99
-183 2 2 0 6 83 133 152
-184 2 2 0 6 73 71 74
-185 2 2 0 6 125 61 132
-186 2 2 0 6 121 141 168
-187 2 2 0 6 94 95 96
-188 2 2 0 6 79 76 80
-189 2 2 0 6 68 65 108
-190 2 2 0 6 11 12 113
-191 2 2 0 6 89 87 146
-192 2 2 0 6 77 87 88
-193 2 2 0 6 84 64 147
-194 2 2 0 6 90 91 92
-195 2 2 0 6 81 80 82
-196 2 2 0 6 78 85 137
-197 2 2 0 6 88 89 90
-198 2 2 0 6 126 125 132
-199 2 2 0 6 54 112 163
-200 2 2 0 6 50 113 144
-201 2 2 0 6 73 74 135
-202 2 2 0 6 28 29 56
-203 2 2 0 6 104 114 147
-204 2 2 0 6 29 30 128
-205 2 2 0 6 44 45 61
-206 2 2 0 6 43 44 125
-207 2 2 0 6 42 43 53
-208 2 2 0 6 30 31 51
-209 2 2 0 6 2 14 154
-210 2 2 0 6 41 42 127
-211 2 2 0 6 40 41 130
-212 2 2 0 6 39 40 47
-213 2 2 0 6 33 34 48
-214 2 2 0 6 3 18 149
-215 2 2 0 6 122 46 161
-216 2 2 0 6 112 97 155
-217 2 2 0 6 32 33 131
-218 2 2 0 6 36 37 46
-219 2 2 0 6 38 39 129
-220 2 2 0 6 133 52 152
-221 2 2 0 6 5 6 134
-222 2 2 0 6 87 101 146
-223 2 2 0 6 93 92 94
-224 2 2 0 6 132 61 134
-225 2 2 0 6 71 68 124
-226 2 2 0 6 85 131 137
-227 2 2 0 6 85 79 156
-228 2 2 0 6 17 3 149
-229 2 2 0 6 80 76 135
-230 2 2 0 6 117 52 133
-231 2 2 0 6 119 142 152
-232 2 2 0 6 66 104 147
-233 2 2 0 6 56 29 128
-234 2 2 0 6 89 146 164
-235 2 2 0 6 50 11 113
-236 2 2 0 6 101 86 151
-237 2 2 0 6 76 73 135
-238 2 2 0 6 45 1 160
-239 2 2 0 6 142 82 152
-240 2 2 0 6 1 5 160
-241 2 2 0 6 44 61 125
-242 2 2 0 6 53 43 125
-243 2 2 0 6 94 96 97
-244 2 2 0 6 4 27 159
-245 2 2 0 6 59 13 154
-246 2 2 0 6 26 4 159
-247 2 2 0 6 30 51 128
-248 2 2 0 6 97 123 166
-249 2 2 0 6 82 80 83
-250 2 2 0 6 42 53 127
-251 2 2 0 6 144 105 157
-252 2 2 0 6 119 58 142
-253 2 2 0 6 90 92 93
-254 2 2 0 6 83 80 135
-255 2 2 0 6 74 71 124
-256 2 2 0 6 41 127 130
-257 2 2 0 6 47 40 130
-258 2 2 0 6 28 56 140
-259 2 2 0 6 35 36 139
-260 2 2 0 6 33 48 131
-261 2 2 0 6 121 88 141
-262 2 2 0 6 39 47 129
-263 2 2 0 6 120 126 132
-264 2 2 0 6 149 60 165
-265 2 2 0 6 53 125 126
-266 2 2 0 6 123 70 168
-267 2 2 0 6 46 122 164
-268 2 2 0 6 27 28 140
-269 2 2 0 6 36 46 139
-270 2 2 0 6 64 66 147
-271 2 2 0 6 51 31 143
-272 2 2 0 6 48 150 151
-273 2 2 0 6 13 2 154
-274 2 2 0 6 81 136 138
-275 2 2 0 6 32 131 143
-276 2 2 0 6 31 32 143
-277 2 2 0 6 131 48 137
-278 2 2 0 6 52 119 152
-279 2 2 0 6 104 105 144
-280 2 2 0 6 123 97 162
-281 2 2 0 6 99 118 163
-282 2 2 0 6 56 128 136
-283 2 2 0 6 48 34 150
-284 2 2 0 6 148 59 154
-285 2 2 0 6 166 123 168
-286 2 2 0 6 112 155 163
-287 2 2 0 6 46 37 161
-288 2 2 0 6 82 83 152
-289 2 2 0 6 38 129 161
-290 2 2 0 6 58 107 142
-291 2 2 0 6 107 106 167
-292 2 2 0 6 51 143 156
-293 2 2 0 6 135 74 145
-294 2 2 0 6 128 51 138
-295 2 2 0 6 97 96 155
-296 2 2 0 6 34 35 150
-297 2 2 0 6 35 139 150
-298 2 2 0 6 88 90 141
-299 2 2 0 6 56 136 167
-300 2 2 0 6 94 97 166
-301 2 2 0 6 90 93 141
-302 2 2 0 6 101 139 146
-303 2 2 0 6 18 60 149
-304 2 2 0 6 100 127 158
-305 2 2 0 6 131 85 143
-306 2 2 0 6 139 46 146
-307 2 2 0 6 136 128 138
-308 2 2 0 6 122 89 164
-309 2 2 0 6 112 111 157
-310 2 2 0 6 133 83 145
-311 2 2 0 6 157 105 162
-312 2 2 0 6 127 53 158
-313 2 2 0 6 150 101 151
-314 2 2 0 6 83 135 145
-315 2 2 0 6 129 122 161
-316 2 2 0 6 37 38 161
-317 2 2 0 6 96 99 155
-318 2 2 0 6 118 54 163
-319 2 2 0 6 136 107 167
-320 2 2 0 6 61 45 160
-321 2 2 0 6 143 85 156
-322 2 2 0 6 5 134 160
-323 2 2 0 6 139 101 150
-324 2 2 0 6 137 48 151
-325 2 2 0 6 53 126 158
-326 2 2 0 6 129 47 153
-327 2 2 0 6 140 106 159
-328 2 2 0 6 134 61 160
-329 2 2 0 6 93 94 166
-330 2 2 0 6 140 56 167
-331 2 2 0 6 27 140 159
-332 2 2 0 6 138 51 156
-333 2 2 0 6 93 166 168
-334 2 2 0 6 155 99 163
-335 2 2 0 6 106 140 167
-336 2 2 0 6 146 46 164
-337 2 2 0 6 141 93 168
-338 2 2 0 6 14 148 154
-$EndElements
diff --git a/test/test_beast.jl b/test/test_beast.jl
index ebf3b5c..f9c8c47 100644
--- a/test/test_beast.jl
+++ b/test/test_beast.jl
@@ -10,7 +10,6 @@ using Test
 k = 2 * π / λ
 
 Γs = meshrectangle(1.0,1.0, 0.3)
-Γt = translate(meshrectangle(1.0,1.0, 0.3), SVector(3.0, 0.0, 1.0))
 
 meshes = [
     (Γs, translate(Γs, SVector(3.0, 0.0, 1.0))),
@@ -27,7 +26,7 @@ for mesh in meshes
             SL, 
             Y,
             X;
-            treeoptions=BoxTreeOptions(nmin=100),
+            treeoptions=BoxTreeOptions(nmin=20),
             compressor=FastBEAST.ACAOptions(tol=1e-4),
             multithreading= multithreading
         )
diff --git a/test/test_fmm.jl b/test/test_fmm.jl
index 2e4aeb1..ee99728 100644
--- a/test/test_fmm.jl
+++ b/test/test_fmm.jl
@@ -62,13 +62,12 @@ for mesh in meshes
             ylF32 = matop(Ofl)*xF32
             ytF64 = matop(Oft)*xF64
             ylF64 = matop(Ofl)*xF64
-
             @test size(matop(Oft))[1] == size(matop(Oft), 1)
             @test size(matop(Oft))[2] == size(matop(Oft), 2)
             @test eltype(ytF32) == promote_type(eltype(xF32), eltype(Oft)) 
-            @test norm(ytF32 - ylF32)/norm(ylF32) ≈ 0 atol=1e-4
+            @test norm(ytF32 - ylF32)/norm(ylF32) ≈ 0 atol=2e-4
             @test eltype(ytF64) == promote_type(eltype(xF64), eltype(Oft)) 
-            @test norm(ytF64 - ylF64)/norm(ylF64) ≈ 0 atol=1e-4
+            @test norm(ytF64 - ylF64)/norm(ylF64) ≈ 0 atol=2e-4
         end
     end
 end
\ No newline at end of file
diff --git a/test/test_hmatrix.jl b/test/test_hmatrix.jl
index 7883ae6..a2ea7d8 100644
--- a/test/test_hmatrix.jl
+++ b/test/test_hmatrix.jl
@@ -1,3 +1,4 @@
+using ClusterTrees
 using FastBEAST
 using LinearAlgebra
 using Random
@@ -27,7 +28,7 @@ end
 
 
 function assembler(kernel, matrix, testpoints, sourcepoints)
-    for j in eachindex(sourcepoints)
+    for j in eachindex(sourcepoints) 
         for i in eachindex(testpoints)
             matrix[i,j] = kernel(testpoints[i], sourcepoints[j])
         end
@@ -37,26 +38,27 @@ end
 N = 2000
 
 spoints = [@SVector rand(3) for i = 1:N]
-tpoints = [@SVector rand(3) for i = 1:N]
+tpoints = [@SVector rand(3) for i = 1:N] 
 
 @views OneoverRkernelassembler(matrix, tdata, sdata) = assembler(
     OneoverRkernel, matrix, tpoints[tdata], spoints[sdata]
 )
 
-stree = create_tree(spoints, BoxTreeOptions(nmin=25))
-ttree = create_tree(tpoints, BoxTreeOptions(nmin=25))
+stree = create_tree(spoints, KMeansTreeOptions(nmin=20))
+ttree = create_tree(tpoints, KMeansTreeOptions(nmin=20))
+
 kmat = assembler(OneoverRkernel, tpoints, spoints)
 
 for multithreading in [true, false]
     hmat = HMatrix(
         OneoverRkernelassembler,
-        ttree,
-        stree,
+        ClusterTrees.BlockTrees.BlockTree(ttree,stree),
         Int64,
         Float64;
         compressor=FastBEAST.ACAOptions(tol=1e-4),
         verbose=true,
-        multithreading=multithreading
+        multithreading=multithreading, 
+        η=1.5
     )
 
     x = rand(N)
diff --git a/test/test_tree.jl b/test/test_tree.jl
index b5191cc..7b0824e 100644
--- a/test/test_tree.jl
+++ b/test/test_tree.jl
@@ -1,35 +1,28 @@
 using FastBEAST
-using LinearAlgebra
+using ClusterTrees
 using StaticArrays
 using Test
 
 #BoxTree
-points2D = [SVector(1.0, 1.0), #3
-            SVector(2.0, 2.0), #1
-            SVector(1.2, 1.7), #2
-            SVector(1.7, 1.2), #4
-            SVector(1.21, 1.7)] #2
+points = [SVector(1.0, 1.0, 1.0), SVector(0.5, 0.5, 0.5)]
 
-root = create_tree(points2D, BoxTreeOptions())
-
-@test length(root.children) == 4
-@test FastBEAST.indices(root.children[1])[1] == 2
-@test FastBEAST.indices(root.children[2])[1] == 3
-@test FastBEAST.indices(root.children[2])[2] == 5
-@test FastBEAST.indices(root.children[3])[1] == 1
-@test FastBEAST.indices(root.children[4])[1] == 4
+tree = create_tree(points, FastBEAST.KMeansTreeOptions())
+nears, fars = computeinteractions(ClusterTrees.BlockTrees.BlockTree(tree, tree))
+@test length(fars[2]) == 2
+##
 
 # KMeansTree
-points2D = [SVector(0.1, 0.1), #1
-            SVector(1.0, 1.0), #2
-            SVector(0.2, 0.2), #3
-            SVector(1.1, 1.1)] #4
-
-root = create_tree(points2D, KMeansTreeOptions(algorithm=:naive))
-
-@test length(root.children) == 2
-@test FastBEAST.indices(root.children[1])[1] == 1
-@test FastBEAST.indices(root.children[1])[2] == 3
-@test FastBEAST.indices(root.children[2])[1] == 2
-@test FastBEAST.indices(root.children[2])[2] == 4
-@test root.children[1].radius ≈ norm([0.05 0.05]) atol=1e-15
\ No newline at end of file
+points = [@SVector rand(3) for i = 1:1000]
+
+tree = create_tree(points, FastBEAST.KMeansTreeOptions(nmin=5))
+nears_strong, fars_strong = computeinteractions(ClusterTrees.BlockTrees.BlockTree(tree, tree), η=1.0)
+nears_weak, fars_weak = computeinteractions(ClusterTrees.BlockTrees.BlockTree(tree, tree), η=2.0)
+
+@test length(fars_weak[6]) > length(fars_strong[6])
+
+tree = create_tree(points, FastBEAST.BoxTreeOptions(nmin=5))
+nears_strong, fars_strong = computeinteractions(ClusterTrees.BlockTrees.BlockTree(tree, tree), η=1.0)
+nears_weak, fars_weak = computeinteractions(ClusterTrees.BlockTrees.BlockTree(tree, tree), η=2.0)
+@test length(fars_weak[3]) > length(fars_strong[3])
+
+