diff --git a/Project.toml b/Project.toml
index a87bfbf9..0ea488a9 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "ExtendableGrids"
 uuid = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
 authors = ["Juergen Fuhrmann <juergen.fuhrmann@wias-berlin.de>, Christian Merdon  <christian.merdon@wias-berlin.de>"]
-version = "1.1.0"
+version = "1.2.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
diff --git a/docs/Project.toml b/docs/Project.toml
index 31a744c4..54e3d1f9 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,6 +1,7 @@
 [deps]
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+ExampleJuggler = "3bbe58f8-ed81-4c4e-a134-03e85fcf4a1a"
 ExtendableGrids = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
 Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb"
 GridVisualize = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8"
@@ -10,4 +11,4 @@ Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344"
 
 [compat]
 Documenter = "1"
-julia = "1.9"
\ No newline at end of file
+julia = "1.9"
diff --git a/docs/make.jl b/docs/make.jl
index c78897c8..fa22ad3f 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,29 +1,12 @@
-ENV["MPLBACKEND"] = "agg"
-using Documenter, ExtendableGrids, Literate, GridVisualize, SimplexGridFactory, Gmsh
-import CairoMakie, Triangulate
-
+using Documenter, ExtendableGrids, ExampleJuggler, Gmsh
+import CairoMakie
 CairoMakie.activate!(; type = "svg", visible = false)
-
-example_md_dir = joinpath(@__DIR__, "src", "examples")
-
-examples1d = joinpath(@__DIR__, "..", "examples", "examples1d.jl")
-include(examples1d)
-examples2d = joinpath(@__DIR__, "..", "examples", "examples2d.jl")
-include(examples2d)
-examples3d = joinpath(@__DIR__, "..", "examples", "examples3d.jl")
-include(examples3d)
-
-include("makeplots.jl")
+ExampleJuggler.verbose!(true)
 
 function mkdocs()
-    Literate.markdown(examples1d, example_md_dir; documenter = false, info = false)
-    Literate.markdown(examples2d, example_md_dir; documenter = false, info = false)
-    Literate.markdown(examples3d, example_md_dir; documenter = false, info = false)
-
-    generated_examples = joinpath.("examples", filter(x -> endswith(x, ".md"), readdir(example_md_dir)))
-
-    makeplots(example_md_dir; Plotter = CairoMakie)
-
+    cleanexamples()
+    generated_examples = @docscripts(joinpath(@__DIR__, "..", "examples"), ["examples1d.jl", "examples2d.jl", "examples3d.jl"],
+                                     Plotter=CairoMakie)
     makedocs(; sitename = "ExtendableGrids.jl",
              modules = [ExtendableGrids, Base.get_extension(ExtendableGrids, :ExtendableGridsGmshExt)],
              clean = false,
diff --git a/docs/makeplots.jl b/docs/makeplots.jl
deleted file mode 100644
index bd61573d..00000000
--- a/docs/makeplots.jl
+++ /dev/null
@@ -1,38 +0,0 @@
-
-function makeplot(func,picdir; Plotter=nothing)
-    f=getfield(Main,Symbol(func))
-    p=gridplot(f(); Plotter)
-    fname=joinpath(picdir,func*".svg")
-    save(fname,p;Plotter)
-    isfile(fname)
-end
-
-function makeplotx(func,picdir;Plotter=nothing)
-    f=getfield(Main,Symbol(func))
-    g,sg,sf=f()
-    p=GridVisualizer(;Plotter, layout=(1,3),resolution=(600,200))
-    gridplot!(p[1,1],g)
-    gridplot!(p[1,2],sg)
-    scalarplot!(p[1,3],sg,sf)
-    fname=joinpath(picdir,func*".svg")
-    save(fname,reveal(p);Plotter)
-    isfile(fname)
-end
-
-function makeplots(picdir;Plotter=nothing)
-
-    makeplot("interval_from_vector",picdir; Plotter)
-    makeplot("interval_localref",picdir; Plotter)
-    makeplot("interval_multiregion",picdir; Plotter)
-    makeplot("interval_subgrid",picdir; Plotter)
-    makeplot("rectangle",picdir; Plotter)
-    makeplot("rectangle_localref",picdir; Plotter)
-    makeplot("rectangle_multiregion",picdir; Plotter)
-    makeplot("rectangle_subgrid",picdir; Plotter)
-    makeplot("rect2d_bregion_function",picdir; Plotter)
-
-    makeplot("quadrilateral",picdir; Plotter)
-    makeplot("cross3d",picdir; Plotter)
-    makeplotx("sorted_subgrid",picdir; Plotter)        
-    
-end
diff --git a/examples/examples1d.jl b/examples/examples1d.jl
index a75376c0..e3cc8b03 100644
--- a/examples/examples1d.jl
+++ b/examples/examples1d.jl
@@ -5,33 +5,30 @@ using ExtendableGrids
 # ## Interval from vector
 #
 function interval_from_vector()
-    X=collect(0:0.05:1)
-    grid=simplexgrid(X)
+    X = collect(0:0.05:1)
+    grid = simplexgrid(X)
 end
 # ![](interval_from_vector.svg)
 
 # 
-# ##Interval with local refinement
+# ## Interval with local refinement
 # 
 function interval_localref()
-
-    XLeft=geomspace(0.0,0.5,0.1, 0.01)
-    XRight=geomspace(0.5,1.0,0.01,0.1)
-    X=glue(XLeft, XRight)
-    grid=simplexgrid(X)
+    XLeft = geomspace(0.0, 0.5, 0.1, 0.01)
+    XRight = geomspace(0.5, 1.0, 0.01, 0.1)
+    X = glue(XLeft, XRight)
+    grid = simplexgrid(X)
 end
 # ![](interval_localref.svg)
 
-
 #
 # ## Interval with  multiple regions
 #
 function interval_multiregion()
-
-    X=collect(0:0.05:1)
-    grid=simplexgrid(X)
-    cellmask!(grid,[0.0],[0.5],3)
-    bfacemask!(grid,[0.5], [0.5],4)
+    X = collect(0:0.05:1)
+    grid = simplexgrid(X)
+    cellmask!(grid, [0.0], [0.5], 3)
+    bfacemask!(grid, [0.5], [0.5], 4)
     grid
 end
 # ![](interval_multiregion.svg)
@@ -39,11 +36,34 @@ end
 # ## Multiple regions and subgrid
 #
 function interval_subgrid()
-    X=collect(0:0.01:1)
-    grid=simplexgrid(X)
-    bfacemask!(grid,[0.5],[0.5],3)
-    cellmask!(grid,[0.0],[0.25],2)
-    cellmask!(grid,[0.20],[0.5],3)
-    subgrid(grid,[2,3])
+    X = collect(0:0.01:1)
+    grid = simplexgrid(X)
+    bfacemask!(grid, [0.5], [0.5], 3)
+    cellmask!(grid, [0.0], [0.25], 2)
+    cellmask!(grid, [0.20], [0.5], 3)
+    subgrid(grid, [2, 3])
 end
 # ![](interval_subgrid.svg)
+# ## CI callbacks
+# Unit tests
+using Test
+
+function runtests()
+    @test numbers_match(interval_from_vector(), 21, 20, 2)
+    @test numbers_match(interval_localref(), 27, 26, 2)
+    @test numbers_match(interval_multiregion(), 21, 20, 3)
+    @test numbers_match(interval_subgrid(), 51, 50, 2)
+end
+
+# Plot generation
+using GridVisualize
+function generateplots(picdir; Plotter = nothing)
+    if isdefined(Plotter, :Makie)
+        size = (500, 200)
+        legend = :rt
+        Plotter.save(joinpath(picdir, "interval_from_vector.svg"), gridplot(interval_from_vector(); Plotter, size, legend))
+        Plotter.save(joinpath(picdir, "interval_localref.svg"), gridplot(interval_localref(); Plotter, size, legend))
+        Plotter.save(joinpath(picdir, "interval_multiregion.svg"), gridplot(interval_multiregion(); Plotter, size, legend))
+        Plotter.save(joinpath(picdir, "interval_subgrid.svg"), gridplot(interval_subgrid(); Plotter, size, legend))
+    end
+end
diff --git a/examples/examples2d.jl b/examples/examples2d.jl
index 0d6c5f13..74673f6e 100644
--- a/examples/examples2d.jl
+++ b/examples/examples2d.jl
@@ -1,104 +1,131 @@
 # 2D Grid examples
 # ===============
 # 
+using Triangulate, ExtendableGrids, SimplexGridFactory
 
 # ## Rectangle
 function rectangle()
-    X=collect(0:0.05:1)
-    Y=collect(0:0.05:1)
-    simplexgrid(X,X)
+    X = collect(0:0.05:1)
+    Y = collect(0:0.05:1)
+    simplexgrid(X, X)
 end
 # ![](rectangle.svg)
 # 
 # ## Rectangle with local refinement
 # 
 function rectangle_localref()
-    hmin=0.01
-    hmax=0.1
-    XLeft=geomspace(0.0,0.5,hmax,hmin)
-    XRight=geomspace(0.5,1.0,hmin,hmax)
-    X=glue(XLeft, XRight)
-    simplexgrid(X,X)
+    hmin = 0.01
+    hmax = 0.1
+    XLeft = geomspace(0.0, 0.5, hmax, hmin)
+    XRight = geomspace(0.5, 1.0, hmin, hmax)
+    X = glue(XLeft, XRight)
+    simplexgrid(X, X)
 end
 # ![](rectangle_localref.svg)
 
-
 # 
 # ## Rectangle with multiple regions
 # 
 function rectangle_multiregion()
-    X=collect(0:0.05:1)
-    Y=collect(0:0.05:1)
-    grid=simplexgrid(X,Y)
-    cellmask!(grid,[0.0,0.0],[1.0,0.5],3)
-    bfacemask!(grid,[0.0,0.0],[0.0,0.5],5)
-    bfacemask!(grid,[1.0,0.0],[1.0,0.5],6)
-    bfacemask!(grid,[0.0,0.5],[1.0,0.5],7)
+    X = collect(0:0.05:1)
+    Y = collect(0:0.05:1)
+    grid = simplexgrid(X, Y)
+    cellmask!(grid, [0.0, 0.0], [1.0, 0.5], 3)
+    bfacemask!(grid, [0.0, 0.0], [0.0, 0.5], 5)
+    bfacemask!(grid, [1.0, 0.0], [1.0, 0.5], 6)
+    bfacemask!(grid, [0.0, 0.5], [1.0, 0.5], 7)
 end
 # ![](rectangle_multiregion.svg)
 
-
-
 # 
 # ## Subgrid from rectangle
 # 
 function rectangle_subgrid()
-    X=collect(0:0.05:1)
-    Y=collect(0:0.05:1)
-    grid=simplexgrid(X,Y)
-    rect!(grid,[0.25,0.25],[0.75,0.75];region=2, bregion=5)
-    subgrid(grid,[1])
+    X = collect(0:0.05:1)
+    Y = collect(0:0.05:1)
+    grid = simplexgrid(X, Y)
+    rect!(grid, [0.25, 0.25], [0.75, 0.75]; region = 2, bregion = 5)
+    subgrid(grid, [1])
 end
 # ![](rectangle_subgrid.svg)
 
-
 # 
 # ## Rect2d with bregion function
 #
 # Here, we use function as bregion parameter - this allows to
 # have no bfaces at the interface between the two rects.
 function rect2d_bregion_function()
-    X=collect(0:0.5:10)
-    Y=collect(0:0.5:10)
-    grid=simplexgrid(X,Y)
-    rect!(grid,[5,4],[9,6];region=2, bregions=[5,5,5,5])
+    X = collect(0:0.5:10)
+    Y = collect(0:0.5:10)
+    grid = simplexgrid(X, Y)
+    rect!(grid, [5, 4], [9, 6]; region = 2, bregions = [5, 5, 5, 5])
+
+    rect!(grid, [4, 2], [5, 8]; region = 2, bregion = cur -> cur == 5 ? 0 : 8)
 
-    rect!(grid,[4,2],[5,8];region=2, bregion= cur-> cur == 5  ? 0 : 8   )
-    
-    subgrid(grid,[2])
-    
+    subgrid(grid, [2])
 end
 # ![](rect2d_bregion_function.svg)
 
-
-
-
-function sorted_subgrid(; maxvolume=0.01)
-    
-    builder=SimplexGridBuilder(Generator=Triangulate)
-    
-    p1=point!(builder,0,0)
-    p2=point!(builder,1,0)
-    p3=point!(builder,1,2)
-    p4=point!(builder,0,1)
-    p5=point!(builder,-1,2)
-
-    facetregion!(builder,1)
-    facet!(builder,p1,p2)
-    facetregion!(builder,2)
-    facet!(builder,p2,p3)
-    facetregion!(builder,3)
-    facet!(builder,p3,p4)
-    facetregion!(builder,4)
-    facet!(builder,p4,p5)
-    facetregion!(builder,5)
-    facet!(builder,p5,p1)
-
-    g=simplexgrid(builder;maxvolume)
-    sg=subgrid(g,[2],boundary=true,transform=(a,b)->a[1]=b[2])
-    f=map( (x,y)->sin(3x)*cos(3y),g)
-    sf=view(f,sg)
-    g,sg,sf
+function sorted_subgrid(; maxvolume = 0.01)
+    builder = SimplexGridBuilder(; Generator = Triangulate)
+
+    p1 = point!(builder, 0, 0)
+    p2 = point!(builder, 1, 0)
+    p3 = point!(builder, 1, 2)
+    p4 = point!(builder, 0, 1)
+    p5 = point!(builder, -1, 2)
+
+    facetregion!(builder, 1)
+    facet!(builder, p1, p2)
+    facetregion!(builder, 2)
+    facet!(builder, p2, p3)
+    facetregion!(builder, 3)
+    facet!(builder, p3, p4)
+    facetregion!(builder, 4)
+    facet!(builder, p4, p5)
+    facetregion!(builder, 5)
+    facet!(builder, p5, p1)
+
+    g = simplexgrid(builder; maxvolume)
+    sg = subgrid(g, [2]; boundary = true, transform = (a, b) -> a[1] = b[2])
+    f = map((x, y) -> sin(3x) * cos(3y), g)
+    sf = view(f, sg)
+    g, sg, sf
 end
 # ![](sorted_subgrid.svg)
+# ## CI callbacks
+# Unit tests
+using Test
+function runtests()
+    @test numbers_match(rectangle(), 441, 800, 80)
+    @test numbers_match(rectangle_localref(), 729, 1352, 104)
+    @test numbers_match(rectangle_multiregion(), 441, 800, 100)
+    @test numbers_match(rectangle_subgrid(), 360, 600, 120)
+    @test numbers_match(rect2d_bregion_function(), 79, 112, 44)
+
+    g, sg, sf = sorted_subgrid()
+    @test numbers_match(g, 187, 306, 66)
+    @test numbers_match(sg, 17, 16, 0)
+    @test issorted(view(sg[Coordinates], 1, :))
+end
 
+# Plot generation
+using GridVisualize
+function generateplots(picdir; Plotter = nothing)
+    if isdefined(Plotter, :Makie)
+        size = (300, 300)
+        Plotter.save(joinpath(picdir, "rectangle.svg"), gridplot(rectangle(); Plotter, size))
+        Plotter.save(joinpath(picdir, "rectangle_localref.svg"), gridplot(rectangle_localref(); Plotter, size))
+        Plotter.save(joinpath(picdir, "rectangle_multiregion.svg"), gridplot(rectangle_multiregion(); Plotter, size))
+        Plotter.save(joinpath(picdir, "rectangle_subgrid.svg"), gridplot(rectangle_subgrid(); Plotter, size))
+        Plotter.save(joinpath(picdir, "rect2d_bregion_function.svg"), gridplot(rect2d_bregion_function(); Plotter, size))
+
+        g, sg, sf = sorted_subgrid()
+        p = GridVisualizer(; Plotter, layout = (1, 3), size = (800, 300))
+        gridplot!(p[1, 1], g)
+        gridplot!(p[1, 2], sg)
+        scalarplot!(p[1, 3], sg, sf)
+        fname = joinpath(picdir, "sorted_subgrid.svg")
+        Plotter.save(fname, reveal(p))
+    end
+end
diff --git a/examples/examples3d.jl b/examples/examples3d.jl
index 08098dcf..929f0d0d 100644
--- a/examples/examples3d.jl
+++ b/examples/examples3d.jl
@@ -1,38 +1,36 @@
 # 3D Grid examples
 # ===============
 # 
+using ExtendableGrids
 
 # ## Quadrilateral
-function quadrilateral(;hx=0.25, hy=0.2, hz=0.1)
-    X=collect(0:hx:1)
-    Y=collect(0:hy:1)
-    Z=collect(0:hz:1)
-    simplexgrid(X,Y,Z)
+function quadrilateral(; hx = 0.25, hy = 0.2, hz = 0.1)
+    X = collect(0:hx:1)
+    Y = collect(0:hy:1)
+    Z = collect(0:hz:1)
+    simplexgrid(X, Y, Z)
 end
 # ![](quadrilateral.svg)
 
 # ## Cross3d
 function cross3d()
-    X=collect(0:1:10)
-    Y=collect(0:1:10)
-    Z=collect(0:1:10)
-    grid=simplexgrid(X,Y,Z)
+    X = collect(0:0.1:1)
+    Y = collect(0:0.1:1)
+    Z = collect(0:0.1:1)
+    grid = simplexgrid(X, Y, Z)
 
-    rect!(grid, (0,4,0), (10,6,2), region=2, bregions=[1,1,1,1,2,3])
+    rect!(grid, (0, 0.4, 0), (1, 0.6, 0.2); region = 2, bregions = [1, 1, 1, 1, 2, 3])
 
-    rect!(grid, (4,0,2), (6,10,4), region=2, bregions=[ 4,4,4,4, (cur)-> cur == 3 ? 0 : 5 , 6] )
-
-    subgrid(grid,[2])
+    rect!(grid, (0.4, 0, 0.2), (0.6, 1, 0.4); region = 2, bregions = [4, 4, 4, 4, (cur) -> cur == 3 ? 0 : 5, 6])
 
+    subgrid(grid, [2])
 end
 # ![](cross3d.svg)
-
-
-
-
+# ## CI callbacks
+# Unit tests
 
 function mask_bedges()
-    grid    = quadrilateral(hx=0.25, hy=0.25, hz=0.25)
+    grid = quadrilateral(; hx = 0.25, hy = 0.25, hz = 0.25)
 
     bedgemask!(grid, [0.0, 0.0, 0.0], [0.0, 0.0, 1.0], 1)
     bedgemask!(grid, [0.0, 0.0, 0.0], [0.0, 1.0, 0.0], 2)
@@ -43,4 +41,20 @@ function mask_bedges()
     true
 end
 
+using Test
+
+function runtests()
+    @test numbers_match(quadrilateral(), 330, 1200, 440)
+    @test mask_bedges()
+    @test numbers_match(cross3d(), 189, 480, 344)
+end
 
+# Plot generation
+using GridVisualize
+function generateplots(picdir; Plotter = nothing)
+    if isdefined(Plotter, :Makie)
+        size = (400, 400)
+        Plotter.save(joinpath(picdir, "quadrilateral.svg"), gridplot(quadrilateral(); Plotter, size))
+        Plotter.save(joinpath(picdir, "cross3d.svg"), gridplot(cross3d(); Plotter, size))
+    end
+end
diff --git a/src/ExtendableGrids.jl b/src/ExtendableGrids.jl
index 81075e01..02f46ce8 100644
--- a/src/ExtendableGrids.jl
+++ b/src/ExtendableGrids.jl
@@ -77,7 +77,7 @@ export dim_space, dim_grid
 export num_nodes, num_cells, num_bfaces, num_bedges
 export num_cellregions, num_bfaceregions, num_bedgeregions
 export gridcomponents
-export seemingly_equal
+export seemingly_equal, numbers_match
 
 include("subgrid.jl")
 export subgrid
diff --git a/src/extendablegrid.jl b/src/extendablegrid.jl
index 42373db7..55b9b783 100644
--- a/src/extendablegrid.jl
+++ b/src/extendablegrid.jl
@@ -562,6 +562,12 @@ seemingly_equal(x1::Type, x2::Type) = (x1 == x2)
 seemingly_equal(x1::Number, x2::Number) = (x1 ≈ x2)
 seemingly_equal(x1::Any, x2::Any) = (x1 == x2)
 
+function numbers_match(grid, nn, nc, nb)
+    num_nodes(grid) == nn &&
+        num_cells(grid) == nc &&
+        num_bfaces(grid) == nb
+end
+
 Base.extrema(grid::ExtendableGrid) = Base.extrema(grid[Coordinates]; dims = 2)
 
 function bbox(grid)
diff --git a/test/Project.toml b/test/Project.toml
index 4d53cf17..6dc69455 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,9 +1,11 @@
 [deps]
-CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+ExampleJuggler = "3bbe58f8-ed81-4c4e-a134-03e85fcf4a1a"
 Gmsh = "705231aa-382f-11e9-3f0c-b7cb4346fdeb"
 GridVisualize = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8"
 SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 SimplexGridFactory = "57bfcd06-606e-45d6-baf4-4ba06da0efd5"
-StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344"
+
+[compat]
+ExampleJuggler = "0.4"
diff --git a/test/runtests.jl b/test/runtests.jl
index 4ec308e7..86c49b29 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,15 +1,7 @@
-ENV["MPLBACKEND"] = "agg"
+using Test, ExampleJuggler
+using ExtendableGrids, SHA
 
-using Test, ExtendableGrids, GridVisualize, SHA, SimplexGridFactory, Triangulate
-import StatsBase
-using StatsBase: countmap
-
-import CairoMakie
-CairoMakie.activate!(; type = "svg", visible = false)
-
-function testgrid(grid, testdata)
-    (num_nodes(grid), num_cells(grid), num_bfaces(grid)) == testdata
-end
+ExampleJuggler.verbose!(true)
 
 @testset "Geomspace" begin
     function test_geomspace()
@@ -128,48 +120,32 @@ end
     @test testrw(simplexgrid(X, X, X), "sg")
 end
 
-examples1d = joinpath(@__DIR__, "..", "examples", "examples1d.jl")
-include(examples1d)
-examples2d = joinpath(@__DIR__, "..", "examples", "examples2d.jl")
-include(examples2d)
-examples3d = joinpath(@__DIR__, "..", "examples", "examples3d.jl")
-include(examples3d)
-
-@testset "1D" begin
+@testset "rectnd" begin
     function rect1d()
         X = collect(0:0.1:1)
         g = simplexgrid(X)
         rect!(g, [0.3], [0.6]; region = 2, bregions = [3, 4])
     end
 
-    @test testgrid(interval_from_vector(), (21, 20, 2))
-    @test testgrid(interval_localref(), (27, 26, 2))
-    @test testgrid(interval_multiregion(), (21, 20, 3))
-    @test testgrid(interval_subgrid(), (51, 50, 2))
-    @test testgrid(rect1d(), (11, 10, 4))
-end
-
-@testset "2D" begin
     function rect2d()
         X = collect(0:0.1:1)
         g = simplexgrid(X, X)
         rect!(g, [0.3, 0.3], [0.6, 0.6]; region = 2, bregions = [3, 4, 5, 6])
     end
 
-    @test testgrid(rectangle(), (441, 800, 80))
-    @test testgrid(rectangle_localref(), (729, 1352, 104))
-    @test testgrid(rectangle_multiregion(), (441, 800, 100))
-    @test testgrid(rectangle_subgrid(), (360, 600, 120))
-    @test testgrid(rect2d(), (121, 200, 52))
-    @test testgrid(rect2d_bregion_function(), (79, 112, 44))
-
-    g, sg, sf = sorted_subgrid()
-    @test testgrid(g, (187, 306, 66))
-    @test testgrid(sg, (17, 16, 0))
-    @test issorted(view(sg[Coordinates], 1, :))
+    function rect3d()
+        X = collect(0:0.1:1)
+        g = simplexgrid(X, X, X)
+        rect!(g, [0.3, 0.3, 0.4], [0.6, 0.6, 0.7]; region = 2, bregion = 8)
+        g
+    end
+
+    @test numbers_match(rect1d(), 11, 10, 4)
+    @test numbers_match(rect2d(), 121, 200, 52)
+    @test numbers_match(rect3d(), 1331, 6000, 1308)
 end
 
-@testset "3D" begin
+@testset "subgrid+extrusion" begin
     function subgen(; h = 0.2)
         X = collect(-1:h:2)
         Y = collect(0:h:1)
@@ -191,42 +167,14 @@ end
         # return subgrid of region 1
         subgrid(g, [1])
     end
-
-    function rect3d()
-        X = collect(0:0.1:1)
-        g = simplexgrid(X, X, X)
-        rect!(g, [0.3, 0.3, 0.4], [0.6, 0.6, 0.7]; region = 2, bregion = 8)
-        g
-    end
-
-    @test testgrid(quadrilateral(), (330, 1200, 440))
-    @test mask_bedges()
+    @test numbers_match(subgen(), 756, 3000, 950)
 
     X = collect(0:0.25:1)
     gxy = simplexgrid(X, X)
     gxyz = simplexgrid(gxy, X)
     g = simplexgrid(X, X, X)
-    @test testgrid(gxyz, (125, 384, 192))
+    @test numbers_match(gxyz, 125, 384, 192)
     @test g[Coordinates] ≈ gxyz[Coordinates]
-    @test testgrid(subgen(), (756, 3000, 950))
-    @test testgrid(rect3d(), (1331, 6000, 1308))
-    @test testgrid(cross3d(), (189, 480, 344))
-end
-
-@testset "plotting examples" begin
-    include("../docs/makeplots.jl")
-    picdir = mktempdir()
-
-    @test makeplot("interval_from_vector", picdir; Plotter = CairoMakie)
-    @test makeplot("interval_localref", picdir; Plotter = CairoMakie)
-    @test makeplot("interval_multiregion", picdir; Plotter = CairoMakie)
-    @test makeplot("interval_subgrid", picdir; Plotter = CairoMakie)
-    @test makeplot("rectangle", picdir; Plotter = CairoMakie)
-    @test makeplot("rectangle_localref", picdir; Plotter = CairoMakie)
-    @test makeplot("rectangle_multiregion", picdir; Plotter = CairoMakie)
-    @test makeplot("rectangle_subgrid", picdir; Plotter = CairoMakie)
-    @test makeplot("quadrilateral", picdir; Plotter = CairoMakie)
-    @test makeplotx("sorted_subgrid", picdir; Plotter = CairoMakie)
 end
 
 function tglue(; dim = 2, breg = 0)
@@ -257,12 +205,16 @@ function tglue(; dim = 2, breg = 0)
 end
 
 @testset "Glue" begin
-    @test testgrid(tglue(; dim = 1, breg = 0), (9, 8, 2))
-    @test testgrid(tglue(; dim = 1, breg = 1), (9, 8, 3))
-    @test testgrid(tglue(; dim = 2, breg = 0), (37, 48, 24))
-    @test testgrid(tglue(; dim = 2, breg = 1), (37, 48, 26))
-    @test testgrid(tglue(; dim = 3, breg = 0), (161, 480, 256))
-    @test testgrid(tglue(; dim = 3, breg = 1), (161, 480, 264))
+    @test numbers_match(tglue(; dim = 1, breg = 0), 9, 8, 2)
+    @test numbers_match(tglue(; dim = 1, breg = 1), 9, 8, 3)
+    @test numbers_match(tglue(; dim = 2, breg = 0), 37, 48, 24)
+    @test numbers_match(tglue(; dim = 2, breg = 1), 37, 48, 26)
+    @test numbers_match(tglue(; dim = 3, breg = 0), 161, 480, 256)
+    @test numbers_match(tglue(; dim = 3, breg = 1), 161, 480, 264)
+end
+
+@testset "Examples" begin
+    @testscripts(joinpath(@__DIR__, "..", "examples"), ["examples1d.jl", "examples2d.jl", "examples3d.jl"])
 end
 
 function voronoitest()