From 01c2a1590733c1fd8d182839dcb01b8166178fdd Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= <juergen-fuhrmann@web.de>
Date: Sun, 5 Feb 2023 17:16:30 +0100
Subject: [PATCH] Subgridplots (#16)

Handle plots of discontinuous functions in Makie,Pyplot, PlutoVista
---
 Project.toml           |   4 +-
 README.md              |   5 +-
 docs/makeplots.jl      |  49 ++-
 examples/plotting.jl   | 215 ++++++++--
 examples/plutovista.jl |  91 ++++-
 src/common.jl          |  19 +-
 src/dispatch.jl        |  70 +++-
 src/makie.jl           | 889 ++++++++++++++++++++++++++---------------
 src/meshcat.jl         | 156 ++++----
 src/plots.jl           | 264 +++++++-----
 src/plutovista.jl      | 406 ++++++++++++-------
 src/pyplot.jl          | 412 ++++++++++++-------
 src/vtkview.jl         |  26 +-
 13 files changed, 1792 insertions(+), 814 deletions(-)

diff --git a/Project.toml b/Project.toml
index ceb1a22..4c62af5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "GridVisualize"
 uuid = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8"
 authors = ["Juergen Fuhrmann <juergen.fuhrmann@wias-berlin.de>"]
-version = "0.6.4"
+version = "1.0.0"
 
 [deps]
 ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
@@ -26,7 +26,7 @@ DocStringExtensions = "0.8,0.9"
 ElasticArrays = "1"
 ExtendableGrids = "0.9"
 GeometryBasics = "0.4.1"
-GridVisualizeTools = "0.2.1"
+GridVisualizeTools = "0.3"
 HypertextLiteral = "0.9"
 Observables = "0.5"
 OrderedCollections = "1"
diff --git a/README.md b/README.md
index 2256646..d84732a 100644
--- a/README.md
+++ b/README.md
@@ -119,7 +119,7 @@ and all plotting functions will do nothing.
 - 'y': available
 - 'i': some level of interactive control
 - '(y)': availability only on rectangular grids
-- 'p':  planned
+- 'p':  planned (no schedule though)
 - 'n': probably not, also in the future
 
 |                | PyPlot | Makie   | PlutoVista | Plots | VTKView |
@@ -135,8 +135,7 @@ and all plotting functions will do nothing.
 | gridplot, 3D   | y      | y,i     | y,i        |  n    | y,i     |
 | vectorplot, 3D | p      | p       | p          |  n    |         |
 | streamplot, 3D |        | p       | p          |  n    |         |
-| movie          | p      | y       | n          |  y    |         |
-
+| movie          | n      | y       | n          |  y    |         |
 
 
 ## Sample output
diff --git a/docs/makeplots.jl b/docs/makeplots.jl
index 2cda5fe..fefd1cc 100644
--- a/docs/makeplots.jl
+++ b/docs/makeplots.jl
@@ -8,6 +8,19 @@ function makeplots(picdir; Plotter=GLMakie, extension="png")
     println("multiscene")
 
    
+
+    p=plotting_func1d(Plotter=Plotter)
+    fname=joinpath(picdir,"plotting_func1d."*extension)
+    save(fname,p,Plotter=Plotter)
+    @test isfile(fname)
+    println("func1d")
+
+    p=plotting_func2d(Plotter=Plotter)
+    fname=joinpath(picdir,"plotting_func2d."*extension)
+    save(fname,p,Plotter=Plotter)
+    @test isfile(fname)
+    println("func2d")
+
     p=plotting_func3d(Plotter=Plotter)
     fname=joinpath(picdir,"plotting_func3d."*extension)
     save(fname,p,Plotter=Plotter)
@@ -15,17 +28,23 @@ function makeplots(picdir; Plotter=GLMakie, extension="png")
     println("func3d")
 
 
-    p=plotting_func2d(Plotter=Plotter)
-    fname=joinpath(picdir,"plotting_func2d."*extension)
+    fname=joinpath(picdir,"plotting_jfunc1d."*"gif")
+    p=plotting_jfunc1d(Plotter=Plotter,filename=fname)
+    @test isfile(fname)
+    println("jfunc1d")
+
+    p=plotting_jfunc2d(Plotter=Plotter)
+    fname=joinpath(picdir,"plotting_jfunc2d."*extension)
     save(fname,p,Plotter=Plotter)
     @test isfile(fname)
-    println("func2d")
+    println("jfunc2d")
 
-    p=plotting_func1d(Plotter=Plotter)
-    fname=joinpath(picdir,"plotting_func1d."*extension)
+    p=plotting_jfunc3d(Plotter=Plotter)
+    fname=joinpath(picdir,"plotting_jfunc3d."*extension)
     save(fname,p,Plotter=Plotter)
     @test isfile(fname)
-    println("func1d")
+    println("jfunc3d")
+
 
     p=plotting_vec2d(Plotter=Plotter)
     fname=joinpath(picdir,"plotting_vec2d."*extension)
@@ -39,11 +58,11 @@ function makeplots(picdir; Plotter=GLMakie, extension="png")
     @test isfile(fname) || GridVisualize.plottername!="PyPlot"
     println("stream2d")
 
-    p=plotting_grid3d(Plotter=Plotter)
-    fname=joinpath(picdir,"plotting_grid3d."*extension)
+    p=plotting_grid1d(Plotter=Plotter)
+    fname=joinpath(picdir,"plotting_grid1d."*extension)
     save(fname,p,Plotter=Plotter)
     @test isfile(fname)
-    println("grid3d")
+    println("grid1d")
 
     p=plotting_grid2d(Plotter=Plotter)
     fname=joinpath(picdir,"plotting_grid2d."*extension)
@@ -51,17 +70,15 @@ function makeplots(picdir; Plotter=GLMakie, extension="png")
     @test isfile(fname)
     println("grid2d")
 
-
-    p=plotting_grid1d(Plotter=Plotter)
-    fname=joinpath(picdir,"plotting_grid1d."*extension)
+    p=plotting_grid3d(Plotter=Plotter)
+    fname=joinpath(picdir,"plotting_grid3d."*extension)
     save(fname,p,Plotter=Plotter)
     @test isfile(fname)
-    println("grid1d")
-
+    println("grid3d")
 
-    fname=joinpath(picdir,"plotting_video."*"gif")
+    fname=joinpath(picdir,"plotting_movie."*"gif")
     p=plotting_movie(;filename=fname,Plotter=Plotter)
     @test isfile(fname)
-    println("plotting_video")
+    println("plotting_movie")
     true
 end
diff --git a/examples/plotting.jl b/examples/plotting.jl
index e8ee6cf..34836c4 100644
--- a/examples/plotting.jl
+++ b/examples/plotting.jl
@@ -19,7 +19,7 @@ using GridVisualize
 
 # ### 1D grids
 function grid1d(; n = 50)
-    X = collect(0:(1 / n):1)
+    X = collect(0:(1/n):1)
     g = simplexgrid(X)
 end
 
@@ -31,7 +31,7 @@ end
 # ### 2D grids
 
 function grid2d(; n = 20)
-    X = collect(0:(1 / n):1)
+    X = collect(0:(1/n):1)
     g = simplexgrid(X, X)
 end
 
@@ -48,7 +48,7 @@ end
 # For Makie and VTKView, the cutplane values can be controlled interactively.
 
 function grid3d(; n = 15)
-    X = collect(0:(1 / n):1)
+    X = collect(0:(1/n):1)
     g = simplexgrid(X, X, X)
 end
 
@@ -98,10 +98,16 @@ function func3d(; n = 15)
     g, map((x, y, z) -> sinpi(2 * x) * sinpi(3.5 * y) * sinpi(1.5 * z), g)
 end
 
-function plotting_func3d(; Plotter = default_plotter(), zplane = 0.49, xplane = 0.49,
-                         levels = [0.25], kwargs...)
+function plotting_func3d(;
+    Plotter = default_plotter(),
+    zplane = 0.49,
+    xplane = 0.49,
+    levels = [0.25],
+    kwargs...,
+)
     g, f = func3d()
-    scalarplot(g, f; Plotter = Plotter, zplane, xplane, levels, kwargs...)
+    #    scalarplot(g, f; Plotter = Plotter, zplane, xplane, levels, kwargs...)
+    scalarplot(g, f; Plotter = Plotter)
 end
 # ![](plotting_func3d.svg)
 
@@ -110,8 +116,10 @@ end
 function vec2d(; n = 20)
     g = grid2d(; n = n)
     g,
-    vcat(map((x, y) -> sinpi(2 * x) * sinpi(3.5 * y), g)',
-         map((x, y) -> cospi(2 * x) * cospi(3.5 * y), g)')
+    vcat(
+        map((x, y) -> sinpi(2 * x) * sinpi(3.5 * y), g)',
+        map((x, y) -> cospi(2 * x) * cospi(3.5 * y), g)',
+    )
 end
 
 function plotting_vec2d(; Plotter = default_plotter(), kwargs...)
@@ -130,21 +138,39 @@ end
 
 # ### Movie
 # Movies can contain  any of the previous plots.
-function plotting_movie(;filename="plotting_video.gif", Plotter=default_plotter())
-    vis=GridVisualizer(Plotter=Plotter,size=(600,200),layout=(1,2))
-    X=0:0.2:10
-    grid=simplexgrid(X,X)
-    movie(vis,file=filename) do vis
-        for t=0:0.1:10
-	    f=map( (x,y)-> sin(x-t)*cos(y-t),grid)
-	    g=map( (x,y)-> sin(t)*sin(x)*cos(y),grid)
-	    scalarplot!(vis[1,1],grid,f,clear=true, title="t=$(t)",limits=(-1,1),levels=7,colormap=:hot)
-	    scalarplot!(vis[1,2],grid,g,clear=true, title="t=$(t)",limits=(-1,1),levels=7,colormap=:hot)
-	    reveal(vis)
+function plotting_movie(; filename = "plotting_video.gif", Plotter = default_plotter())
+    vis = GridVisualizer(Plotter = Plotter, size = (600, 200), layout = (1, 2))
+    X = 0:0.2:10
+    grid = simplexgrid(X, X)
+    movie(vis, file = filename) do vis
+        for t = 0:0.1:10
+            f = map((x, y) -> sin(x - t) * cos(y - t), grid)
+            g = map((x, y) -> sin(t) * sin(x) * cos(y), grid)
+            scalarplot!(
+                vis[1, 1],
+                grid,
+                f,
+                clear = true,
+                title = "t=$(t)",
+                limits = (-1, 1),
+                levels = 7,
+                colormap = :hot,
+            )
+            scalarplot!(
+                vis[1, 2],
+                grid,
+                g,
+                clear = true,
+                title = "t=$(t)",
+                limits = (-1, 1),
+                levels = 7,
+                colormap = :hot,
+            )
+            reveal(vis)
         end
     end
 end
-# ![](plotting_video.gif)
+# ![](plotting_movie.gif)
 
 
 # ## Multiscene plots
@@ -158,22 +184,155 @@ end
 # one subscene.
 function plotting_multiscene!(p)
     gridplot!(p[1, 1], grid1d(); title = "1D grid")
-    scalarplot!(p[2, 1], grid1d(), sin; title = "1D grid function", label = "sin",
-                markershape = :diamond, color = :red, legend = :rb)
-    scalarplot!(p[2, 1], grid1d(), cos; title = "1D grid function", label = "cos",
-                linestyle = :dash, markershape = :none, color = :green, clear = false)
+    scalarplot!(
+        p[2, 1],
+        grid1d(),
+        sin;
+        title = "1D grid function",
+        label = "sin",
+        markershape = :diamond,
+        color = :red,
+        legend = :rb,
+    )
+    scalarplot!(
+        p[2, 1],
+        grid1d(),
+        cos;
+        title = "1D grid function",
+        label = "cos",
+        linestyle = :dash,
+        markershape = :none,
+        color = :green,
+        clear = false,
+    )
     gridplot!(p[1, 2], grid2d(); title = "2D grid")
     scalarplot!(p[2, 2], func2d()...; colormap = :bamako, title = "2D grid function")
     gridplot!(p[1, 3], grid3d(); zplane = 0.49, title = "3D grid")
-    scalarplot!(p[2, 3], func3d()...; zplane = 0.49, flevel = 0.5, colormap = :bamako,
-                title = "3D grid function")
+    scalarplot!(
+        p[2, 3],
+        func3d()...;
+        zplane = 0.49,
+        flevel = 0.5,
+        colormap = :bamako,
+        title = "3D grid function",
+    )
     vectorplot!(p[1, 4], vec2d()...; title = "2D quiver")
     GridVisualize.streamplot!(p[2, 4], vec2d()...; title = "2D stream")
     reveal(p)
 end
 
 function plotting_multiscene(; Plotter = default_plotter(), resolution = (1000, 500))
-    plotting_multiscene!(GridVisualizer(; Plotter = Plotter, layout = (2, 4), clear = true,
-                                        resolution = resolution))
+    plotting_multiscene!(
+        GridVisualizer(;
+            Plotter = Plotter,
+            layout = (2, 4),
+            clear = true,
+            resolution = resolution,
+        ),
+    )
 end
 # ![](plotting_multiscene.svg)
+
+
+
+# ## Plots of functions on subgrids interfaces
+# We can jointly plot functions on different subgrids whic
+# e.g. model a particle density jumping at a heterointerface
+# Currently supported for PyPlot and Makie
+#
+# The general scheme is to pass a vector of subgrid, the parent grid and the corresponding
+# vector of functions on the respective subgrids
+#
+# ### 1D case
+function plotting_jfunc1d(; Plotter = default_plotter(), filename = "plotting_jfunc1d.gif")
+    X = 0:1:10
+    g = simplexgrid(X)
+    cellmask!(g, [0], [5], 2)
+    g1 = subgrid(g, [1])
+    g2 = subgrid(g, [2])
+
+
+    vis = GridVisualizer(; Plotter, color = :red)
+    movie(vis, file = filename) do vis
+        for t = 0:0.05:1
+            func1 = map((x) -> x - t, g1)
+            func2 = map((x) -> -x + t, g2)
+            func = map(x -> x^2 / 100 - t, g)
+            scalarplot!(
+                vis,
+                [g1, g2],
+                g,
+                [func1, func2];
+                Plotter,
+                elevation = 0.1,
+                clear = true,
+                color = :red,
+            )
+            scalarplot!(
+                vis,
+                g,
+                func;
+                Plotter,
+                elevation = 0.1,
+                clear = false,
+                color = :green,
+            )
+            reveal(vis)
+        end
+    end
+end
+
+# ![](plotting_jfunc1d.gif)
+
+
+# ### 2D case
+function plotting_jfunc2d(; Plotter = default_plotter(), kwargs...)
+    X = 0:0.7:10
+    g = simplexgrid(X, X)
+    cellmask!(g, [0, 0], [5, 5], 2)
+    g1 = subgrid(g, [1])
+    g2 = subgrid(g, [2])
+
+    func1 = map((x, y) -> x^2 + y, g1)
+    func2 = map((x, y) -> (x + y^2), g2)
+    scalarplot([g1, g2], g, [func1, func2]; Plotter, kwargs...)
+end
+
+# ![](plotting_jfunc2d.svg)
+
+# ### 3D case
+function plotting_jfunc3d(;
+    Plotter = default_plotter(),
+    levels = 0,
+    yplane = 0.25,
+    xplane = 0.25,
+    zplane = 0.25,
+    levelalpha = 1,
+    colormap = :hot,
+    kwargs...,
+)
+    X = 0:0.1:1
+    g = simplexgrid(X, X, X)
+    cellmask!(g, [0, 0, 0], [0.5, 0.5, 0.5], 2)
+    g1 = subgrid(g, [1])
+    g2 = subgrid(g, [2])
+    func1 = map((x, y, z) -> (x + y + z), g1)
+    func2 = map((x, y, z) -> (3 - x - y - z), g2)
+    scalarplot(
+        [g1, g2],
+        g,
+        [func1, func2];
+        Plotter,
+        levels,
+        xplane,
+        yplane,
+        zplane,
+        levelalpha,
+        colormap,
+        kwargs...,
+    )
+end
+
+# ![](plotting_jfunc3d.svg)
+
+
diff --git a/examples/plutovista.jl b/examples/plutovista.jl
index 25c5e5b..b9cd8eb 100644
--- a/examples/plutovista.jl
+++ b/examples/plutovista.jl
@@ -116,6 +116,44 @@ p = GridVisualizer(; resolution = (600, 200), dim = 1, legend = :lt);
 # ╔═╡ 661531f7-f740-4dd4-9a59-89ddff06ba5c
 scalarplot!(p, X2, f2(t4); show = true, clear = true, color = :red, label = "t=$(t4)")
 
+# ╔═╡ dda4599d-05a2-4131-899a-42a653a18b51
+md"""
+### Non-continuous functions
+"""
+
+# ╔═╡ 00b115d3-aa8e-43ef-be6b-3d9d7b42f8af
+let
+	X = 0:1:10
+    g = simplexgrid(X)
+    cellmask!(g, [0], [5], 2)
+    g1 = subgrid(g, [1])
+    g2 = subgrid(g, [2])
+
+
+    vis = GridVisualizer(color = :red,size=(600,200))
+    func1 = map((x) -> x, g1)
+    func2 = map((x) -> -x, g2)
+    func = map(x -> x^2 / 100, g)
+    scalarplot!(
+                vis,
+                [g1, g2],
+                g,
+                [func1, func2];
+                elevation = 0.1,
+                clear = true,
+                color = :red,
+            )
+    scalarplot!(
+                vis,
+                g,
+                func;
+                elevation = 0.1,
+                clear = false,
+                color = :green,
+            )
+            reveal(vis)
+end
+
 # ╔═╡ ed9b80e5-9678-4ba6-bb36-c2e0674ed9ba
 md"""
 ## 1D grid plot 
@@ -174,6 +212,24 @@ t= $(@bind t PlutoUI.Slider(-10:0.1:10, default=0, show_value=true))
 # ╔═╡ 412c905f-050c-4b78-a66f-0d03978e7edf
 scalarplot!(vis, grid, f(t); limits = (-π / 2, π / 2), show = true, levels = 5)
 
+# ╔═╡ e3c5a486-ec9b-4010-901b-07f7ac997355
+md"""
+### Non-continuous functions
+"""
+
+# ╔═╡ 9b41e77b-e6b9-43b2-89ce-14a9c0eb1242
+let
+	X = 0:0.7:10
+    g = simplexgrid(X, X)
+    cellmask!(g, [0, 0], [5, 5], 2)
+    g1 = subgrid(g, [1])
+    g2 = subgrid(g, [2])
+
+    func1 = map((x, y) -> x^2 + y, g1)
+    func2 = map((x, y) -> (x + y^2), g2)
+    scalarplot([g1, g2], g, [func1, func2])
+end
+
 # ╔═╡ e9bc2dae-c303-4063-9ea9-36f95f93371c
 md"""
 ## 2D Grid plot
@@ -237,6 +293,33 @@ end
 # ╔═╡ 57ed5eea-bc1c-45eb-b4d3-dc63088db21a
 scalarplot(g4, map((x, y, z) -> 0.01 * exp(-10 * (x^2 + y^2 + z^2)), g4); levels = 10)
 
+# ╔═╡ 943da8f0-d18f-40d5-8158-a3ab5793112f
+md"""
+### Non-continuous functions
+"""
+
+# ╔═╡ ef973737-5cc3-4a3c-8859-a86d9c12c976
+let
+	X = 0:0.1:1
+    g = simplexgrid(X, X, X)
+    cellmask!(g, [0, 0, 0], [0.5, 0.5, 0.5], 2)
+    g1 = subgrid(g, [1])
+    g2 = subgrid(g, [2])
+    func1 = map((x, y, z) -> (x + y + z), g1)
+    func2 = map((x, y, z) -> (3 - x - y - z), g2)
+    scalarplot(
+        [g1, g2],
+        g,
+        [func1, func2];
+        levels = 0,
+        yplane = 0.25,
+        xplane = 0.25,
+        zplane = 0.25,
+        levelalpha = 1,
+        colormap = :hot,
+    )
+end
+
 # ╔═╡ 4b9113d2-10bd-4f7a-a2b8-22092656c6b3
 md"""
 ## 3D grid plot
@@ -385,6 +468,8 @@ md"""
 # ╟─29fa4467-65ee-4dad-a660-5197864ddbdc
 # ╠═c1278fb2-3e75-445f-893a-b8b8a7e931d3
 # ╠═661531f7-f740-4dd4-9a59-89ddff06ba5c
+# ╟─dda4599d-05a2-4131-899a-42a653a18b51
+# ╠═00b115d3-aa8e-43ef-be6b-3d9d7b42f8af
 # ╟─ed9b80e5-9678-4ba6-bb36-c2e0674ed9ba
 # ╠═9ce4f63d-cd96-48d7-a637-07cb84fa88ab
 # ╠═d503ee1e-1e1f-4235-b286-dc3137a2c96a
@@ -396,6 +481,8 @@ md"""
 # ╠═faa59bbd-df1f-4c62-9a77-c4752c6a6df4
 # ╠═412c905f-050c-4b78-a66f-0d03978e7edf
 # ╟─6f1707ed-79ab-42dc-8ad8-d66a9e1a65b3
+# ╟─e3c5a486-ec9b-4010-901b-07f7ac997355
+# ╠═9b41e77b-e6b9-43b2-89ce-14a9c0eb1242
 # ╟─e9bc2dae-c303-4063-9ea9-36f95f93371c
 # ╠═2b3cb0f4-0656-4981-bec6-48785caf2994
 # ╠═1388c246-be49-4757-a2cc-a685642b6b37
@@ -408,6 +495,8 @@ md"""
 # ╠═ecd941a0-85b7-4bb7-a903-b19a500198e1
 # ╠═d924d90d-4102-4ae8-b8de-254a17a5d4df
 # ╠═57ed5eea-bc1c-45eb-b4d3-dc63088db21a
+# ╟─943da8f0-d18f-40d5-8158-a3ab5793112f
+# ╠═ef973737-5cc3-4a3c-8859-a86d9c12c976
 # ╟─4b9113d2-10bd-4f7a-a2b8-22092656c6b3
 # ╠═f78196ca-d972-4fa6-bdc2-e76eba7ca5a1
 # ╟─7dd92757-c100-4158-baa8-1d9218c39aa7
@@ -420,7 +509,7 @@ md"""
 # ╠═f6205299-d097-4e78-8488-b088475191f6
 # ╟─6915cc3e-ad9b-4721-9933-884cfc68a25a
 # ╟─49db8b25-50ce-4fb4-bea2-de8abfb53c56
-# ╠═15f4eeb3-c42e-449c-9161-f1df66de6cef
+# ╟─15f4eeb3-c42e-449c-9161-f1df66de6cef
 # ╟─75ffcd09-dfa8-42df-a3cd-a7e68786e73c
 # ╟─cf592b99-d596-4511-adbf-001145a59983
 # ╠═bc1c6d12-8d06-4f57-9044-8b5e86fd1c13
diff --git a/src/common.jl b/src/common.jl
index fc6e9a1..16ab8ea 100644
--- a/src/common.jl
+++ b/src/common.jl
@@ -56,6 +56,13 @@ function GridVisualizeTools.marching_tetrahedra(grid::ExtendableGrid, func, plan
     marching_tetrahedra(coord, cellnodes, func, planes, flevels; kwargs...)
 end
 
+function GridVisualizeTools.marching_tetrahedra(grids::Vector{ExtendableGrid{Tv,Ti}}, funcs, planes, flevels;
+                                                kwargs...) where {Tv,Ti}
+    coord = [grid[Coordinates] for grid in grids]
+    cellnodes = [grid[CellNodes] for grid in grids]
+    marching_tetrahedra(coord, cellnodes, funcs, planes, flevels; kwargs...)
+end
+
 ########################################################################################
 """
     $(SIGNATURES)
@@ -68,6 +75,14 @@ function GridVisualizeTools.marching_triangles(grid::ExtendableGrid, func, level
     marching_triangles(coord, cellnodes, func, levels)
 end
 
+
+function GridVisualizeTools.marching_triangles(grids::Vector{ExtendableGrid{Tv,Ti}}, funcs, levels) where {Tv,Ti}
+    coords = [grid[Coordinates] for grid in grids]
+    cellnodes = [ grid[CellNodes] for grid in grids]
+    marching_triangles(coords, cellnodes, funcs, levels)
+end
+
+
 ##############################################
 # Create meshes from grid data
 function regionmesh(grid, iregion)
@@ -280,8 +295,8 @@ function vectorsample(grid::ExtendableGrid{Tv, Ti}, v; offset = :default,
 end
 
 # Calculate isolevel values and function limits
-function isolevels(ctx, func)
-    makeisolevels(func,
+function isolevels(ctx, funcs)
+    makeisolevels(funcs,
                   ctx[:levels],
                   ctx[:limits] == :auto ? (1, -1) : ctx[:limits],
                   ctx[:colorbarticks] == :default ? nothing : ctx[:colorbarticks])
diff --git a/src/dispatch.jl b/src/dispatch.jl
index 0df4a42..8e54fc2 100644
--- a/src/dispatch.jl
+++ b/src/dispatch.jl
@@ -433,7 +433,38 @@ Keyword arguments: see [`available_kwargs`](@ref)
 function scalarplot!(ctx::SubVisualizer, grid::ExtendableGrid, func; kwargs...)
     _update_context!(ctx, Dict(:clear => true, :show => false, :reveal => false))
     _update_context!(ctx, kwargs)
-    scalarplot!(ctx, plottertype(ctx[:Plotter]), Val{dim_space(grid)}, grid, func)
+    scalarplot!(ctx, plottertype(ctx[:Plotter]), Val{dim_space(grid)}, [grid], grid, [func])
+end
+
+"""
+$(TYPEDSIGNATURES)
+
+Plot node vectors on subgrids of parent grid as P1 FEM function on the triangulation into subplot in the visualizer.
+If `[i,j]` is omitted, `[1,1]` is assumed.
+eyword arguments: see [`available_kwargs`](@ref)
+"""
+function scalarplot!(
+    ctx::SubVisualizer,
+    grids::Vector{ExtendableGrid{Tv,Ti}},
+    parentgrid::ExtendableGrid{Tv,Ti},
+    funcs::AbstractVector;
+    kwargs...,
+) where {Tv,Ti}
+    _update_context!(ctx, Dict(:clear => true, :show => false, :reveal => false))
+    _update_context!(ctx, kwargs)
+    if length(grids) != length(funcs)
+        error(
+            "number of subgrids: $(length(grids)) and number of functions: $(length(funcs)) not equal",
+        )
+    end
+    scalarplot!(
+        ctx,
+        plottertype(ctx[:Plotter]),
+        Val{dim_space(parentgrid)},
+        grids,
+        parentgrid,
+        funcs,
+    )
 end
 
 "$(TYPEDSIGNATURES)"
@@ -446,9 +477,23 @@ function scalarplot!(ctx::SubVisualizer, grid::ExtendableGrid, func::Function; k
     scalarplot!(ctx, grid, map(func, grid); kwargs...)
 end
 
+"$(TYPEDSIGNATURES)"
+function scalarplot!(
+    p::GridVisualizer,
+    grids::Vector{ExtendableGrid{Tv,Ti}},
+    parentgrid::ExtendableGrid{Tv,Ti},
+    funcs::AbstractVector;
+    kwargs...,
+) where {Tv,Ti}
+    scalarplot!(p[1, 1], grids, parentgrid, funcs; kwargs...)
+end
+
+
+
+
 "$(TYPEDSIGNATURES)"
 function scalarplot!(ctx::SubVisualizer, func::AbstractVector; kwargs...)
-    scalarplot!(ctx, simplexgrid(X), 1:length(func), func; kwargs...)
+    scalarplot!(ctx, simplexgrid(1:length(func)), func; kwargs...)
 end
 
 "$(TYPEDSIGNATURES)"
@@ -499,11 +544,30 @@ function scalarplot(grid::ExtendableGrid, func; Plotter = default_plotter(), kwa
     scalarplot!(GridVisualizer(; Plotter = Plotter, kwargs...), grid, func; show = true)
 end
 
+
+function scalarplot(
+    grids::Vector{ExtendableGrid{Tv,Ti}},
+    parentgrid::ExtendableGrid{Tv,Ti},
+    funcs::AbstractVector;
+    Plotter = default_plotter(),
+    kwargs...,
+) where {Tv,Ti}
+    scalarplot!(
+        GridVisualizer(; Plotter = Plotter, kwargs...),
+        grids,
+        parentgrid,
+        funcs;
+        show = true,
+    )
+end
+
+
 "$(TYPEDSIGNATURES)"
 scalarplot(func::AbstractVector; kwargs...) = scalarplot(1:length(func), func; kwargs...)
 
 "$(TYPEDSIGNATURES)"
-scalarplot(X::AbstractVector, func; kwargs...) = scalarplot(simplexgrid(X), func; kwargs...)
+scalarplot(X::AbstractVector{T}, func; kwargs...) where {T<:Number} =
+    scalarplot(simplexgrid(X), func; kwargs...)
 
 "$(TYPEDSIGNATURES)"
 function scalarplot(X::AbstractVector, Y::AbstractVector, func; kwargs...)
diff --git a/src/makie.jl b/src/makie.jl
index a00e221..b703ad4 100644
--- a/src/makie.jl
+++ b/src/makie.jl
@@ -12,11 +12,13 @@ function initialize!(p::GridVisualizer, ::Type{MakieType})
     version_installed = PkgVersion.Version(XMakie.Makie)
 
     if version_installed < version_min
-        error("Outdated version $(version_installed) of XMakie. Please upgrade to at least $(version_min)")
+        error(
+            "Outdated version $(version_installed) of Makie. Please upgrade to at least $(version_min)",
+        )
     end
 
     if version_installed > version_max
-        @warn("Possibly breaking version $(version_installed) of XMakie.")
+        @warn("Possibly breaking version $(version_installed) of Makie.")
     end
 
     # Prepare flippable layout
@@ -52,8 +54,8 @@ add_scene!(ctx, ax) = ctx[:flayout][ctx[:subplot]...] = ax
 
 # Revealing the  visualizer just returns the figure
 function reveal(p::GridVisualizer, ::Type{MakieType})
-    Plotter=p.context[:Plotter]
-    if haskey(p.context,:videostream)
+    Plotter = p.context[:Plotter]
+    if haskey(p.context, :videostream)
         Plotter.recordframe!(p.context[:videostream])
     else
         p.context[:figure]
@@ -76,20 +78,28 @@ end
 
 
 
-function movie(func, vis::GridVisualizer, ::Type{MakieType}; file=nothing, format="gif" ,kwargs...)
-    Plotter=vis.context[:Plotter]
+function movie(
+    func,
+    vis::GridVisualizer,
+    ::Type{MakieType};
+    file = nothing,
+    format = "gif",
+    kwargs...,
+)
+    Plotter = vis.context[:Plotter]
     if !isnothing(file)
         format = lstrip(splitext(file)[2], '.')
     end
-    
+
     if isdefined(Main, :PlutoRunner) || !isnothing(file)
-        vis.context[:videostream] = Plotter.VideoStream(vis.context[:figure]; format=format, kwargs...)
+        vis.context[:videostream] =
+            Plotter.VideoStream(vis.context[:figure]; format = format, kwargs...)
     end
-    
+
     func(vis)
-    
+
     if !isnothing(file)
-        Plotter.save(file,vis.context[:videostream])
+        Plotter.save(file, vis.context[:videostream])
     elseif isdefined(Main, :PlutoRunner)
         vis.context[:videostream]
     else
@@ -115,8 +125,12 @@ by triggering change values via up/down (± 1)  resp. page_up/page_down (±10) k
 
 The update_scene callback gets passed the change value and the symbol.
 """
-function scene_interaction(update_scene, scene, XMakie,
-                           switchkeys::Vector{Symbol} = [:nothing])
+function scene_interaction(
+    update_scene,
+    scene,
+    XMakie,
+    switchkeys::Vector{Symbol} = [:nothing],
+)
 
     # Check if pixel position pos sits within the scene
     function _inscene(scene, pos)
@@ -172,16 +186,17 @@ end
 
 # Standard kwargs for Makie scenes
 scenekwargs(ctx) = Dict(
-                        #:xticklabelsize => 0.5*ctx[:fontsize],
-                        #:yticklabelsize => 0.5*ctx[:fontsize],
-                        #:zticklabelsize => 0.5*ctx[:fontsize],
-                        #:xlabelsize => 0.5*ctx[:fontsize],
-                        #:ylabelsize => 0.5*ctx[:fontsize],
-                        #:zlabelsize => 0.5*ctx[:fontsize],
-                        #:xlabeloffset => 20,
-                        #:ylabeloffset => 20,
-                        #:zlabeloffset => 20,
-                        :titlesize => ctx[:fontsize])
+    #:xticklabelsize => 0.5*ctx[:fontsize],
+    #:yticklabelsize => 0.5*ctx[:fontsize],
+    #:zticklabelsize => 0.5*ctx[:fontsize],
+    #:xlabelsize => 0.5*ctx[:fontsize],
+    #:ylabelsize => 0.5*ctx[:fontsize],
+    #:zlabelsize => 0.5*ctx[:fontsize],
+    #:xlabeloffset => 20,
+    #:ylabeloffset => 20,
+    #:zlabeloffset => 20,
+    :titlesize => ctx[:fontsize],
+)
 
 #scenekwargs(ctx)=()
 
@@ -250,53 +265,65 @@ function gridplot!(ctx, TP::Type{MakieType}, ::Type{Val{1}}, grid)
     nbregions = num_bfaceregions(grid)
 
     if !haskey(ctx, :scene)
-        ctx[:scene] = XMakie.Axis(ctx[:figure];
-                                  yticklabelsvisible = false,
-                                  yticksvisible = false,
-                                  title = ctx[:title],
-                                  scenekwargs(ctx)...)
+        ctx[:scene] = XMakie.Axis(
+            ctx[:figure];
+            yticklabelsvisible = false,
+            yticksvisible = false,
+            title = ctx[:title],
+            scenekwargs(ctx)...,
+        )
 
         ctx[:grid] = Observable(grid)
         cmap = region_cmap(nregions)
         bcmap = bregion_cmap(nbregions)
 
         # Set scene size with invisible markers
-        XMakie.scatter!(ctx[:scene],
-                        map(g -> scenecorners1d(grid), ctx[:grid]);
-                        color = :white,
-                        markersize = 0.0,
-                        strokewidth = 0)
+        XMakie.scatter!(
+            ctx[:scene],
+            map(g -> scenecorners1d(grid), ctx[:grid]);
+            color = :white,
+            markersize = 0.0,
+            strokewidth = 0,
+        )
 
         # Draw node markers
-        XMakie.linesegments!(ctx[:scene],
-                             map(g -> basemesh1d(g), ctx[:grid]);
-                             color = :black)
+        XMakie.linesegments!(
+            ctx[:scene],
+            map(g -> basemesh1d(g), ctx[:grid]);
+            color = :black,
+        )
 
         # Colored cell regions
         for i = 1:nregions
-            XMakie.linesegments!(ctx[:scene],
-                                 map(g -> regionmesh1d(g, i), ctx[:grid]);
-                                 color = cmap[i],
-                                 linewidth = 4,
-                                 label = "c $(i)")
+            XMakie.linesegments!(
+                ctx[:scene],
+                map(g -> regionmesh1d(g, i), ctx[:grid]);
+                color = cmap[i],
+                linewidth = 4,
+                label = "c $(i)",
+            )
         end
 
         # Colored boundary grid
         for i = 1:nbregions
-            XMakie.linesegments!(ctx[:scene],
-                                 map(g -> bregionmesh1d(g, i), ctx[:grid]);
-                                 color = bcmap[i],
-                                 linewidth = 4,
-                                 label = "b$(i)")
+            XMakie.linesegments!(
+                ctx[:scene],
+                map(g -> bregionmesh1d(g, i), ctx[:grid]);
+                color = bcmap[i],
+                linewidth = 4,
+                label = "b$(i)",
+            )
         end
 
         # Legende
         if ctx[:legend] != :none
             pos = ctx[:legend] == :best ? :rt : ctx[:legend]
-            XMakie.axislegend(ctx[:scene];
-                              position = pos,
-                              labelsize = 0.5 * ctx[:fontsize],
-                              nbanks = 5)
+            XMakie.axislegend(
+                ctx[:scene];
+                position = pos,
+                labelsize = 0.5 * ctx[:fontsize],
+                nbanks = 5,
+            )
         end
         XMakie.reset_limits!(ctx[:scene])
         add_scene!(ctx, ctx[:scene])
@@ -309,9 +336,10 @@ end
 ########################################################################
 # 1D function
 
-function scalarplot!(ctx, TP::Type{MakieType}, ::Type{Val{1}}, grid, func)
+function scalarplot!(ctx, TP::Type{MakieType}, ::Type{Val{1}}, grids, parentgrid, funcs)
     XMakie = ctx[:Plotter]
 
+    nfuncs = length(funcs)
     if ctx[:title] == ""
         ctx[:title] = " "
     end
@@ -337,7 +365,7 @@ function scalarplot!(ctx, TP::Type{MakieType}, ::Type{Val{1}}, grid, func)
         points = [Point2f(coord[1, i], func[i]) for i = 1:length(func)]
     end
 
-    coord = grid[Coordinates]
+    coord = parentgrid[Coordinates]
     xlimits = ctx[:xlimits]
     ylimits = ctx[:limits]
     xmin = coord[1, 1]
@@ -347,50 +375,76 @@ function scalarplot!(ctx, TP::Type{MakieType}, ::Type{Val{1}}, grid, func)
         xmax = xlimits[2]
     end
 
-    (ymin, ymax) = extrema(func)
+
+
     if ylimits[1] < ylimits[2]
         ymin = ylimits[1]
         ymax = ylimits[2]
+    else
+        ext=extrema.(funcs)
+        (ymin, ymax)=(minimum(first.(ext)),maximum(last.(ext)))
     end
 
-    function update_lines(ctx)
+    function update_lines(ctx, newrange)
         if ctx[:markershape] == :none
             #line without marker
-            XMakie.lines!(ctx[:scene], map(a -> a, ctx[:lines][end]);
-                          linestyle = ctx[:linestyle],
-                          linewidth = ctx[:linewidth],
-                          color = RGB(ctx[:color]),
-                          label = ctx[:label])
+            for l in newrange
+                XMakie.lines!(
+                    ctx[:scene],
+                    map(a -> a, ctx[:lines][l]);
+                    linestyle = ctx[:linestyle],
+                    linewidth = ctx[:linewidth],
+                    color = RGB(ctx[:color]),
+                )
+            end
+            if ctx[:label] != ""
+                XMakie.scatterlines!(
+                    ctx[:scene],
+                    map(a -> a[1:1], ctx[:lines][newrange[begin]]);
+                    linestyle = ctx[:linestyle],
+                    linewidth = ctx[:linewidth],
+                    color = RGB(ctx[:color]),
+                    label = ctx[:label],
+                )
+            end
         else
             # line with markers separated by markevery
 
             # draw plain line without the label
-            XMakie.lines!(ctx[:scene], map(a -> a, ctx[:lines][end]);
-                          linestyle = ctx[:linestyle],
-                          color = RGB(ctx[:color]),
-                          linewidth = ctx[:linewidth])
-
-            # draw markers without label
-            XMakie.scatter!(ctx[:scene],
-                            map(a -> a[1:ctx[:markevery]:end], ctx[:lines][end]);
-                            color = RGB(ctx[:color]),
-                            marker = ctx[:markershape],
-                            markercolor = RGB(ctx[:color]),
-                            markersize = ctx[:markersize])
+            for l in newrange
+                XMakie.lines!(
+                    ctx[:scene],
+                    map(a -> a, ctx[:lines][l]);
+                    linestyle = ctx[:linestyle],
+                    color = RGB(ctx[:color]),
+                    linewidth = ctx[:linewidth],
+                )
+                # draw markers without label
+                XMakie.scatter!(
+                    ctx[:scene],
+                    map(a -> a[1:ctx[:markevery]:end], ctx[:lines][l]);
+                    color = RGB(ctx[:color]),
+                    marker = ctx[:markershape],
+                    markercolor = RGB(ctx[:color]),
+                    markersize = ctx[:markersize],
+                )
+            end
 
             # Draw  dummy line with marker on top ot the first
             # marker position already drawn in order to
             # get the proper legend entry
             if ctx[:label] != ""
-                XMakie.scatterlines!(ctx[:scene],
-                                     map(a -> a[1:1], ctx[:lines][end]);
-                                     linestyle = ctx[:linestyle],
-                                     linewidth = ctx[:linewidth],
-                                     marker = ctx[:markershape],
-                                     markersize = ctx[:markersize],
-                                     markercolor = RGB(ctx[:color]),
-                                     color = RGB(ctx[:color]),
-                                     label = ctx[:label])
+                XMakie.scatterlines!(
+                    ctx[:scene],
+                    map(a -> a[1:1], ctx[:lines][newrange[begin]]);
+                    linestyle = ctx[:linestyle],
+                    linewidth = ctx[:linewidth],
+                    marker = ctx[:markershape],
+                    markersize = ctx[:markersize],
+                    markercolor = RGB(ctx[:color]),
+                    color = RGB(ctx[:color]),
+                    label = ctx[:label],
+                )
             end
         end
 
@@ -405,47 +459,57 @@ function scalarplot!(ctx, TP::Type{MakieType}, ::Type{Val{1}}, grid, func)
         ctx[:xtitle] = Observable(ctx[:title])
 
         # Axis
-        ctx[:scene] = XMakie.Axis(ctx[:figure];
-                                  title = map(a -> a, ctx[:xtitle]),
-                                  xscale = ctx[:xscale] == :log ? log10 : identity,
-                                  yscale = ctx[:yscale] == :log ? log10 : identity,
-                                  scenekwargs(ctx)...)
+        ctx[:scene] = XMakie.Axis(
+            ctx[:figure];
+            title = map(a -> a, ctx[:xtitle]),
+            xscale = ctx[:xscale] == :log ? log10 : identity,
+            yscale = ctx[:yscale] == :log ? log10 : identity,
+            scenekwargs(ctx)...,
+        )
         # Plot size
-        XMakie.scatter!(ctx[:scene],
-                        [Point2f(xmin, ymin), Point2f(xmax, ymax)];
-                        color = :white,
-                        markersize = 0.0,
-                        strokewidth = 0)
+        XMakie.scatter!(
+            ctx[:scene],
+            [Point2f(xmin, ymin), Point2f(xmax, ymax)];
+            color = :white,
+            markersize = 0.0,
+            strokewidth = 0,
+        )
 
         # ctx[:lines]  is an array of lines to draw
         # Here, we start just with the first one.
-        ctx[:lines] = [Observable(polyline(grid, func))]
 
-        update_lines(ctx)
+        ctx[:lines] = [Observable(polyline(grids[i], funcs[i])) for i = 1:nfuncs]
+
+        update_lines(ctx, 1:nfuncs)
 
         XMakie.reset_limits!(ctx[:scene])
 
-        ctx[:nlines] = 1
+        ctx[:nlines] = nfuncs
 
         XMakie.reset_limits!(ctx[:scene])
         add_scene!(ctx, ctx[:scene])
 
     else
         if ctx[:clear]
-            ctx[:nlines] = 1
+            ctx[:nlines] = nfuncs
         else
-            ctx[:nlines] += 1
+            ctx[:nlines] += nfuncs
         end
 
-        p = polyline(grid, func)
         # Either update existing line, or
         # create new one. This works with repeating sequences of
         # updating lines.
         if ctx[:nlines] <= length(ctx[:lines])
-            ctx[:lines][ctx[:nlines]][] = p
+            for i = 1:nfuncs
+                ctx[:lines][ctx[:nlines]-nfuncs+i][] = polyline(grids[i], funcs[i])
+            end
         else
-            push!(ctx[:lines], Observable(p))
-            update_lines(ctx)
+            r0 = length(ctx[:lines])
+            for i = 1:nfuncs
+                push!(ctx[:lines], Observable(polyline(grids[i], funcs[i])))
+            end
+            r1 = length(ctx[:lines])
+            update_lines(ctx, r0+1:r1)
         end
 
         XMakie.reset_limits!(ctx[:scene])
@@ -466,25 +530,33 @@ function makescene_grid(ctx)
     nbcol = length(ctx[:cmap])
     # fontsize=0.5*ctx[:fontsize],ticklabelsize=0.5*ctx[:fontsize]
     if ctx[:colorbar] == :vertical
-        GL[1, 2] = XMakie.Colorbar(ctx[:figure];
-                                   colormap = XMakie.cgrad(ctx[:cmap]; categorical = true),
-                                   limits = (1, ncol),
-                                   width = 15)
-        GL[1, 3] = XMakie.Colorbar(ctx[:figure];
-                                   colormap = XMakie.cgrad(ctx[:bcmap]; categorical = true),
-                                   limits = (1, nbcol),
-                                   width = 15)
+        GL[1, 2] = XMakie.Colorbar(
+            ctx[:figure];
+            colormap = XMakie.cgrad(ctx[:cmap]; categorical = true),
+            limits = (1, ncol),
+            width = 15,
+        )
+        GL[1, 3] = XMakie.Colorbar(
+            ctx[:figure];
+            colormap = XMakie.cgrad(ctx[:bcmap]; categorical = true),
+            limits = (1, nbcol),
+            width = 15,
+        )
     elseif ctx[:colorbar] == :horizontal
-        GL[2, 1] = XMakie.Colorbar(ctx[:figure];
-                                   colormap = XMakie.cgrad(ctx[:cmap]; categorical = true),
-                                   limits = (1, ncol),
-                                   heigth = 15,
-                                   vertical = false)
-        GL[3, 1] = XMakie.Colorbar(ctx[:figure];
-                                   colormap = XMakie.cgrad(ctx[:bcmap]; categorical = true),
-                                   limits = (1, nbcol),
-                                   heigth = 15,
-                                   vertical = false)
+        GL[2, 1] = XMakie.Colorbar(
+            ctx[:figure];
+            colormap = XMakie.cgrad(ctx[:cmap]; categorical = true),
+            limits = (1, ncol),
+            heigth = 15,
+            vertical = false,
+        )
+        GL[3, 1] = XMakie.Colorbar(
+            ctx[:figure];
+            colormap = XMakie.cgrad(ctx[:bcmap]; categorical = true),
+            limits = (1, nbcol),
+            heigth = 15,
+            vertical = false,
+        )
     end
     GL
 end
@@ -509,31 +581,38 @@ function gridplot!(ctx, TP::Type{MakieType}, ::Type{Val{2}}, grid)
             ex = extrema(coord; dims = (2))
             data_aspect = (ex[2][2] - ex[2][1]) / (ex[1][2] - ex[1][1])
             makie_aspect = 1.0 / (data_aspect * ctx[:aspect])
-            ctx[:scene] = XMakie.Axis(ctx[:figure];
-                                      title = ctx[:title],
-                                      aspect = makie_aspect,
-                                      scenekwargs(ctx)...)
+            ctx[:scene] = XMakie.Axis(
+                ctx[:figure];
+                title = ctx[:title],
+                aspect = makie_aspect,
+                scenekwargs(ctx)...,
+            )
         end
 
         # Draw cells with region mark
         cmap = region_cmap(nregions)
         ctx[:cmap] = cmap
         for i = 1:nregions
-            XMakie.poly!(ctx[:scene], map(g -> regionmesh(g, i), ctx[:grid]);
-                         color = cmap[i],
-                         strokecolor = :black,
-                         strokewidth = ctx[:linewidth])
+            XMakie.poly!(
+                ctx[:scene],
+                map(g -> regionmesh(g, i), ctx[:grid]);
+                color = cmap[i],
+                strokecolor = :black,
+                strokewidth = ctx[:linewidth],
+            )
         end
 
         # Draw boundary lines
         bcmap = bregion_cmap(nbregions)
         ctx[:bcmap] = bcmap
         for i = 1:nbregions
-            lp = XMakie.linesegments!(ctx[:scene],
-                                      map(g -> bfacesegments(g, i), ctx[:grid]);
-                                      label = "$(i)",
-                                      color = bcmap[i],
-                                      linewidth = 4)
+            lp = XMakie.linesegments!(
+                ctx[:scene],
+                map(g -> bfacesegments(g, i), ctx[:grid]);
+                label = "$(i)",
+                color = bcmap[i],
+                linewidth = 4,
+            )
             XMakie.translate!(lp, 0, 0, 0.1)
         end
         XMakie.reset_limits!(ctx[:scene])
@@ -541,10 +620,12 @@ function gridplot!(ctx, TP::Type{MakieType}, ::Type{Val{2}}, grid)
         # Describe legend
         if ctx[:legend] != :none
             pos = ctx[:legend] == :best ? :rt : ctx[:legend]
-            XMakie.axislegend(ctx[:scene];
-                              position = pos,
-                              labelsize = 0.5 * ctx[:fontsize],
-                              backgroundcolor = :transparent)
+            XMakie.axislegend(
+                ctx[:scene];
+                position = pos,
+                labelsize = 0.5 * ctx[:fontsize],
+                backgroundcolor = :transparent,
+            )
         end
         add_scene!(ctx, makescene_grid(ctx))
     end
@@ -564,73 +645,138 @@ function makescene2d(ctx, key)
 
     # , fontsize=0.5*ctx[:fontsize],ticklabelsize=0.5*ctx[:fontsize]
     if ctx[:colorbar] == :vertical
-        GL[1, 2] = XMakie.Colorbar(ctx[:figure], ctx[key]; width = 10,
-                                   ticks = unique(ctx[:cbarticks]))
+        GL[1, 2] = XMakie.Colorbar(
+            ctx[:figure],
+            ctx[key];
+            width = 10,
+            ticks = unique(ctx[:cbarticks]),
+        )
     elseif ctx[:colorbar] == :horizontal
-        GL[2, 1] = XMakie.Colorbar(ctx[:figure], ctx[key]; height = 10,
-                                   ticks = ctx[:cbarticks], vertical = false)
+        GL[2, 1] = XMakie.Colorbar(
+            ctx[:figure],
+            ctx[key];
+            height = 10,
+            ticks = ctx[:cbarticks],
+            vertical = false,
+        )
     end
     GL
 end
 
 # 2D function
-function scalarplot!(ctx, TP::Type{MakieType}, ::Type{Val{2}}, grid, func)
+function scalarplot!(ctx, TP::Type{MakieType}, ::Type{Val{2}}, grids, parentgrid, funcs)
     XMakie = ctx[:Plotter]
 
     # Create GeometryBasics.mesh from grid data.
-    function make_mesh(grid::ExtendableGrid, func, elevation)
-        coord = grid[Coordinates]
-        npoints = num_nodes(grid)
-        cellnodes = grid[CellNodes]
+    function make_mesh(grids, funcs, elevation)
+        ngrids = length(grids)
+        coords = [grid[Coordinates] for grid in grids]
+        npoints = [num_nodes(grid) for grid in grids]
+        cellnodes = [grid[CellNodes] for grid in grids]
+        ncells = [num_cells(grid) for grid in grids]
+        offsets = zeros(Int, ngrids)
+
+        for i = 2:ngrids
+            offsets[i] = offsets[i-1] + npoints[i-1]
+        end
+
         if elevation ≈ 0.0
-            points = [Point2f(coord[1, i], coord[2, i]) for i = 1:npoints]
+            points = Vector{Point2f}(undef, sum(npoints))
+            k = 1
+            for j = 1:ngrids
+                for i = 1:npoints[j]
+                    points[k] = Point2f(coords[j][1, i], coords[j][2, i])
+                    k = k + 1
+                end
+            end
         else
-            points = [Point3f(coord[1, i], coord[2, i], func[i] * elevation)
-                      for i = 1:npoints]
+            points = Vector{Point3f}(undef, sum(npoints))
+            k = 1
+            for j = 1:ngrids
+                for i = 1:npoints[j]
+                    points[k] =
+                        Point3f(coords[j][1, i], coords[j][2, i], funcs[j][i] * elevation)
+                    k = k + 1
+                end
+            end
+        end
+        faces = Vector{TriangleFace{Int64}}(undef, sum(ncells))
+        k = 1
+        for j = 1:ngrids
+            for i = 1:ncells[j]
+                faces[k] = TriangleFace(
+                    cellnodes[j][1, i] + offsets[j],
+                    cellnodes[j][2, i] + offsets[j],
+                    cellnodes[j][3, i] + offsets[j],
+                )
+                k = k + 1
+            end
         end
-        faces = [TriangleFace(cellnodes[1, i], cellnodes[2, i], cellnodes[3, i])
-                 for i = 1:size(cellnodes, 2)]
         Mesh(points, faces)
     end
 
-    levels, crange, ctx[:cbarticks] = isolevels(ctx, func)
+    levels, crange, ctx[:cbarticks] = isolevels(ctx, funcs)
+
+
     eps = 1.0e-1
     if crange[1] == crange[2]
         crange = (crange[1] - eps, crange[1] + eps)
     end
 
-    set_plot_data!(ctx, :contourdata,
-                   (g = grid, f = func, e = ctx[:elevation], t = ctx[:title], l = levels,
-                    c = crange))
+    set_plot_data!(
+        ctx,
+        :contourdata,
+        (
+            g = grids,
+            f = funcs,
+            e = ctx[:elevation],
+            t = ctx[:title],
+            l = levels,
+            c = crange,
+        ),
+    )
 
     if !haskey(ctx, :contourplot)
         if !haskey(ctx, :scene)
-            coord = grid[Coordinates]
+            coord = parentgrid[Coordinates]
             ex = extrema(coord; dims = (2))
             data_aspect = (ex[2][2] - ex[2][1]) / (ex[1][2] - ex[1][1])
             makie_aspect = 1.0 / (data_aspect * ctx[:aspect])
-            # would need to switch to Axis3 for supporting elevation
-            ctx[:scene] = XMakie.Axis(ctx[:figure];
-                                      title = map(data -> data.t, ctx[:contourdata]),
-                                      aspect = makie_aspect,
-                                      scenekwargs(ctx)...)
+            if ctx[:elevation] ≈ 0
+                ctx[:scene] = XMakie.Axis(
+                    ctx[:figure];
+                    title = map(data -> data.t, ctx[:contourdata]),
+                    aspect = makie_aspect,
+                    scenekwargs(ctx)...,
+                )
+            else
+                ctx[:scene] = XMakie.Axis3(
+                    ctx[:figure];
+                    title = map(data -> data.t, ctx[:contourdata]),
+                    scenekwargs(ctx)...,
+                )
+            end
         end
 
         # Draw the mesh for the cells
-        ctx[:contourplot] = XMakie.poly!(ctx[:scene],
-                                         map(data -> make_mesh(data.g, data.f, data.e),
-                                             ctx[:contourdata]);
-                                         color = map(data -> data.f, ctx[:contourdata]),
-                                         colorrange = map(data -> data.c,
-                                                          ctx[:contourdata]),
-                                         colormap = ctx[:colormap])
+        ctx[:contourplot] = XMakie.poly!(
+            ctx[:scene],
+            map(data -> make_mesh(data.g, data.f, data.e), ctx[:contourdata]);
+            color = map(data -> vcat(data.f...), ctx[:contourdata]),
+            colorrange = map(data -> data.c, ctx[:contourdata]),
+            colormap = ctx[:colormap],
+        )
 
         # draw the isolines via marching triangles
-        XMakie.linesegments!(ctx[:scene],
-                             map(data -> marching_triangles(data.g, data.f, data.l),
-                                 ctx[:contourdata]);
-                             color = :black,
-                             linewidth = ctx[:linewidth])
+        if ctx[:elevation] ≈ 0
+            XMakie.linesegments!(
+                ctx[:scene],
+                map(data -> marching_triangles(data.g, data.f, data.l), ctx[:contourdata]);
+                color = :black,
+                linewidth = ctx[:linewidth],
+            )
+        end
+
 
         XMakie.reset_limits!(ctx[:scene])
         add_scene!(ctx, makescene2d(ctx, :contourplot))
@@ -649,20 +795,24 @@ function vectorplot!(ctx, TP::Type{MakieType}, ::Type{Val{2}}, grid, func)
 
     if !haskey(ctx, :arrowplot)
         if !haskey(ctx, :scene)
-            ctx[:scene] = XMakie.Axis(ctx[:figure];
-                                      title = ctx[:title],
-                                      aspect = XMakie.DataAspect(),
-                                      scenekwargs(ctx)...)
+            ctx[:scene] = XMakie.Axis(
+                ctx[:figure];
+                title = ctx[:title],
+                aspect = XMakie.DataAspect(),
+                scenekwargs(ctx)...,
+            )
             add_scene!(ctx, ctx[:scene])
         end
 
-        ctx[:arrowplot] = XMakie.arrows!(ctx[:scene],
-                                         map(data -> data.qc[1, :], ctx[:arrowdata]),
-                                         map(data -> data.qc[2, :], ctx[:arrowdata]),
-                                         map(data -> data.qv[1, :], ctx[:arrowdata]),
-                                         map(data -> data.qv[2, :], ctx[:arrowdata]);
-                                         color = :black,
-                                         linewidth = ctx[:linewidth])
+        ctx[:arrowplot] = XMakie.arrows!(
+            ctx[:scene],
+            map(data -> data.qc[1, :], ctx[:arrowdata]),
+            map(data -> data.qc[2, :], ctx[:arrowdata]),
+            map(data -> data.qv[1, :], ctx[:arrowdata]),
+            map(data -> data.qv[2, :], ctx[:arrowdata]);
+            color = :black,
+            linewidth = ctx[:linewidth],
+        )
         XMakie.reset_limits!(ctx[:scene])
     end
     reveal(ctx, TP)
@@ -700,14 +850,16 @@ function makeaxis3d(ctx)
         XMakie.LScene(ctx[:figure])
     else
         # "New" Axis3 with prospective new stuff by Julius.
-        XMakie.Axis3(ctx[:figure];
-                     aspect = :data,
-                     viewmode = :fit,
-                     elevation = ctx[:elev] * π / 180,
-                     azimuth = ctx[:azim] * π / 180,
-                     perspectiveness = ctx[:perspectiveness],
-                     title = map(data -> data.t, ctx[:data]),
-                     scenekwargs(ctx)...)
+        XMakie.Axis3(
+            ctx[:figure];
+            aspect = :data,
+            viewmode = :fit,
+            elevation = ctx[:elev] * π / 180,
+            azimuth = ctx[:azim] * π / 180,
+            perspectiveness = ctx[:perspectiveness],
+            title = map(data -> data.t, ctx[:data]),
+            scenekwargs(ctx)...,
+        )
     end
 end
 
@@ -722,39 +874,47 @@ function makescene3d(ctx)
     GL = XMakie.GridLayout(ctx[:figure]; default_rowgap = 0)
     if ctx[:scene3d] == "LScene"
         # LScene has no title, put the title into protrusion space on top  of the scene
-        GL[1, 1, XMakie.Top()] = XMakie.Label(ctx[:figure],
-                                              " $(map(data->data.t,ctx[:data])) ";
-                                              tellwidth = false,
-                                              height = 30,
-                                              fontsize = ctx[:fontsize])
+        GL[1, 1, XMakie.Top()] = XMakie.Label(
+            ctx[:figure],
+            " $(map(data->data.t,ctx[:data])) ";
+            tellwidth = false,
+            height = 30,
+            fontsize = ctx[:fontsize],
+        )
     end
     GL[1, 1] = ctx[:scene]
     # Horizontal or vertical colorbar
     if haskey(ctx, :crange)
         if ctx[:colorbar] == :vertical
-            GL[1, 2] = XMakie.Colorbar(ctx[:figure];
-                                       colormap = ctx[:colormap],
-                                       colorrange = ctx[:crange],
-                                       ticks = ctx[:levels],
-                                       width = 15,
-                                       ticklabelsize = 0.5 * ctx[:fontsize])
+            GL[1, 2] = XMakie.Colorbar(
+                ctx[:figure];
+                colormap = ctx[:colormap],
+                colorrange = ctx[:crange],
+                ticks = ctx[:levels],
+                width = 15,
+                ticklabelsize = 0.5 * ctx[:fontsize],
+            )
         elseif ctx[:colorbar] == :horizontal
-            GL[2, 1] = XMakie.Colorbar(ctx[:figure];
-                                       colormap = ctx[:colormap],
-                                       colorrange = ctx[:crange],
-                                       ticks = ctx[:levels],
-                                       height = 15,
-                                       ticklabelsize = 0.5 * ctx[:fontsize],
-                                       vertical = false)
+            GL[2, 1] = XMakie.Colorbar(
+                ctx[:figure];
+                colormap = ctx[:colormap],
+                colorrange = ctx[:crange],
+                ticks = ctx[:levels],
+                height = 15,
+                ticklabelsize = 0.5 * ctx[:fontsize],
+                vertical = false,
+            )
         end
     end
 
     # Put the status label into protrusion space on the bottom of the scene
-    GL[1, 1, XMakie.Bottom()] = XMakie.Label(ctx[:figure],
-                                             ctx[:status];
-                                             tellwidth = false,
-                                             height = 30,
-                                             fontsize = 0.5 * ctx[:fontsize])
+    GL[1, 1, XMakie.Bottom()] = XMakie.Label(
+        ctx[:figure],
+        ctx[:status];
+        tellwidth = false,
+        height = 30,
+        fontsize = 0.5 * ctx[:fontsize],
+    )
     GL
 end
 
@@ -788,8 +948,13 @@ function gridplot!(ctx, TP::Type{MakieType}, ::Type{Val{3}}, grid)
     adjust_planes()
 
     if !haskey(ctx, :scene)
-        ctx[:data] = Observable((g = grid, x = ctx[:xplanes][1], y = ctx[:yplanes][1],
-                                 z = ctx[:zplanes][1], t = ctx[:title]))
+        ctx[:data] = Observable((
+            g = grid,
+            x = ctx[:xplanes][1],
+            y = ctx[:yplanes][1],
+            z = ctx[:zplanes][1],
+            t = ctx[:title],
+        ))
 
         ctx[:scene] = makeaxis3d(ctx)
         cmap = region_cmap(nregions)
@@ -800,74 +965,100 @@ function gridplot!(ctx, TP::Type{MakieType}, ::Type{Val{3}}, grid)
         ############# Interior cuts
         # We draw a mesh for each color.
         if ctx[:interior]
-            ctx[:celldata] = map(d -> extract_visible_cells3D(d.g,
-                                                              (d.x, d.y, d.z);
-                                                              primepoints = hcat(xyzmin,
-                                                                                 xyzmax),
-                                                              Tp = Point3f,
-                                                              Tf = GLTriangleFace),
-                                 ctx[:data])
-
-            ctx[:cellmeshes] = map(d -> [make_mesh(d[1][i], d[2][i]) for i = 1:nregions],
-                                   ctx[:celldata])
+            ctx[:celldata] = map(
+                d -> extract_visible_cells3D(
+                    d.g,
+                    (d.x, d.y, d.z);
+                    primepoints = hcat(xyzmin, xyzmax),
+                    Tp = Point3f,
+                    Tf = GLTriangleFace,
+                ),
+                ctx[:data],
+            )
+
+            ctx[:cellmeshes] =
+                map(d -> [make_mesh(d[1][i], d[2][i]) for i = 1:nregions], ctx[:celldata])
 
             for i = 1:nregions
-                XMakie.mesh!(ctx[:scene], map(d -> d[i], ctx[:cellmeshes]);
-                             color = cmap[i],
-                             backlight = 1.0f0)
+                XMakie.mesh!(
+                    ctx[:scene],
+                    map(d -> d[i], ctx[:cellmeshes]);
+                    color = cmap[i],
+                    backlight = 1.0f0,
+                )
 
                 if ctx[:linewidth] > 0
-                    XMakie.wireframe!(ctx[:scene], map(d -> d[i], ctx[:cellmeshes]);
-                                      color = :black,
-                                      strokecolor = :black,
-                                      strokewidth = ctx[:linewidth],
-                                      linewidth = ctx[:linewidth])
+                    XMakie.wireframe!(
+                        ctx[:scene],
+                        map(d -> d[i], ctx[:cellmeshes]);
+                        color = :black,
+                        strokecolor = :black,
+                        strokewidth = ctx[:linewidth],
+                        linewidth = ctx[:linewidth],
+                    )
                 end
             end
         end
 
         ############# Visible boundary faces
-        ctx[:facedata] = map(d -> extract_visible_bfaces3D(d.g,
-                                                           (d.x, d.y, d.z);
-                                                           primepoints = hcat(xyzmin,
-                                                                              xyzmax),
-                                                           Tp = Point3f,
-                                                           Tf = GLTriangleFace),
-                             ctx[:data])
-
-        ctx[:facemeshes] = map(d -> [make_mesh(d[1][i], d[2][i]) for i = 1:nbregions],
-                               ctx[:facedata])
+        ctx[:facedata] = map(
+            d -> extract_visible_bfaces3D(
+                d.g,
+                (d.x, d.y, d.z);
+                primepoints = hcat(xyzmin, xyzmax),
+                Tp = Point3f,
+                Tf = GLTriangleFace,
+            ),
+            ctx[:data],
+        )
+
+        ctx[:facemeshes] =
+            map(d -> [make_mesh(d[1][i], d[2][i]) for i = 1:nbregions], ctx[:facedata])
 
         for i = 1:nbregions
-            XMakie.mesh!(ctx[:scene], map(d -> d[i], ctx[:facemeshes]);
-                         color = bcmap[i],
-                         backlight = 1.0f0)
+            XMakie.mesh!(
+                ctx[:scene],
+                map(d -> d[i], ctx[:facemeshes]);
+                color = bcmap[i],
+                backlight = 1.0f0,
+            )
             if ctx[:linewidth] > 0
-                XMakie.wireframe!(ctx[:scene], map(d -> d[i], ctx[:facemeshes]);
-                                  color = :black,
-                                  strokecolor = :black,
-                                  linewidth = ctx[:linewidth])
+                XMakie.wireframe!(
+                    ctx[:scene],
+                    map(d -> d[i], ctx[:facemeshes]);
+                    color = :black,
+                    strokecolor = :black,
+                    linewidth = ctx[:linewidth],
+                )
             end
         end
 
         ############# Transparent outline
 
         if ctx[:outlinealpha] > 0.0
-            ctx[:outlinedata] = map(d -> extract_visible_bfaces3D(d.g,
-                                                                  xyzmax;
-                                                                  primepoints = hcat(xyzmin,
-                                                                                     xyzmax),
-                                                                  Tp = Point3f,
-                                                                  Tf = GLTriangleFace),
-                                    ctx[:data])
-            ctx[:outlinemeshes] = map(d -> [make_mesh(d[1][i], d[2][i]) for i = 1:nbregions],
-                                      ctx[:outlinedata])
+            ctx[:outlinedata] = map(
+                d -> extract_visible_bfaces3D(
+                    d.g,
+                    xyzmax;
+                    primepoints = hcat(xyzmin, xyzmax),
+                    Tp = Point3f,
+                    Tf = GLTriangleFace,
+                ),
+                ctx[:data],
+            )
+            ctx[:outlinemeshes] = map(
+                d -> [make_mesh(d[1][i], d[2][i]) for i = 1:nbregions],
+                ctx[:outlinedata],
+            )
 
             for i = 1:nbregions
-                XMakie.mesh!(ctx[:scene], map(d -> d[i], ctx[:outlinemeshes]);
-                             color = (bcmap[i], ctx[:outlinealpha]),
-                             transparency = true,
-                             backlight = 1.0f0)
+                XMakie.mesh!(
+                    ctx[:scene],
+                    map(d -> d[i], ctx[:outlinemeshes]);
+                    color = (bcmap[i], ctx[:outlinealpha]),
+                    transparency = true,
+                    backlight = 1.0f0,
+                )
             end
         end
 
@@ -887,8 +1078,13 @@ function gridplot!(ctx, TP::Type{MakieType}, ::Type{Val{3}}, grid)
             end
             adjust_planes()
 
-            ctx[:data][] = (g = grid, x = ctx[:xplanes][1], y = ctx[:yplanes][1],
-                            z = ctx[:zplanes][1], t = ctx[:title])
+            ctx[:data][] = (
+                g = grid,
+                x = ctx[:xplanes][1],
+                y = ctx[:yplanes][1],
+                z = ctx[:zplanes][1],
+                t = ctx[:title],
+            )
         end
 
         ctx[:status] = Observable(" ")
@@ -896,16 +1092,21 @@ function gridplot!(ctx, TP::Type{MakieType}, ::Type{Val{3}}, grid)
         add_scene!(ctx, makescene_grid(ctx))
 
     else
-        ctx[:data][] = (g = grid, x = ctx[:xplanes][1], y = ctx[:yplanes][1],
-                        z = ctx[:zplanes][1], t = ctx[:title])
+        ctx[:data][] = (
+            g = grid,
+            x = ctx[:xplanes][1],
+            y = ctx[:yplanes][1],
+            z = ctx[:zplanes][1],
+            t = ctx[:title],
+        )
     end
 
     reveal(ctx, TP)
 end
 
 # 3d function
-function scalarplot!(ctx, TP::Type{MakieType}, ::Type{Val{3}}, grid, func)
-    levels, crange = isolevels(ctx, func)
+function scalarplot!(ctx, TP::Type{MakieType}, ::Type{Val{3}}, grids, parentgrid, funcs)
+    levels, crange = isolevels(ctx, funcs)
     ctx[:crange] = crange
 
     nan_replacement = 0.5 * (crange[1] + crange[2])
@@ -915,8 +1116,9 @@ function scalarplot!(ctx, TP::Type{MakieType}, ::Type{Val{3}}, grid, func)
         if length(fcs) > 0
             colors = XMakie.Makie.interpolated_getindex.((cmap,), vals, (crange,))
             if ctx[:levelalpha] > 0
-                colors = [RGBA(colors[i].r, colors[i].g, colors[i].b,
-                               Float32(ctx[:levelalpha])) for i = 1:length(colors)]
+                colors = [
+                RGBA(colors[i].r, colors[i].g, colors[i].b, Float32(ctx[:levelalpha])) for i = 1:length(colors)
+                     ]
             end
             GeometryBasics.Mesh(meta(pts; color = colors, normals = normals(pts, fcs)), fcs)
         else
@@ -924,12 +1126,12 @@ function scalarplot!(ctx, TP::Type{MakieType}, ::Type{Val{3}}, grid, func)
         end
     end
 
-    nregions = num_cellregions(grid)
-    nbregions = num_bfaceregions(grid)
+
+    nbregions = num_bfaceregions(parentgrid)
 
     XMakie = ctx[:Plotter]
     cmap = XMakie.to_colormap(ctx[:colormap])
-    xyzmin, xyzmax = xyzminmax(grid)
+    xyzmin, xyzmax = xyzminmax(parentgrid)
     xyzstep = (xyzmax - xyzmin) / 100
 
     fstep = (crange[2] - crange[1]) / 100
@@ -946,18 +1148,25 @@ function scalarplot!(ctx, TP::Type{MakieType}, ::Type{Val{3}}, grid, func)
 
     # adjust_planes()
 
+    ctx[:xplanes]=collect(ctx[:xplanes])
+    ctx[:yplanes]=collect(ctx[:yplanes])    
+    ctx[:zplanes]=collect(ctx[:zplanes])
+
     x = ctx[:xplanes]
     y = ctx[:yplanes]
     z = ctx[:zplanes]
 
     ε = 1.0e-5 * (xyzmax .- xyzmin)
 
-    ctx[:xplanes] = isa(x, Number) ?
-                    collect(range(xyzmin[1] + ε[1], xyzmax[1] - ε[1]; length = ceil(x))) : x
-    ctx[:yplanes] = isa(y, Number) ?
-                    collect(range(xyzmin[2] + ε[2], xyzmax[2] - ε[2]; length = ceil(y))) : y
-    ctx[:zplanes] = isa(z, Number) ?
-                    collect(range(xyzmin[3] + ε[3], xyzmax[3] - ε[3]; length = ceil(z))) : z
+    ctx[:xplanes] =
+        isa(x, Number) ?
+        collect(range(xyzmin[1] + ε[1], xyzmax[1] - ε[1]; length = ceil(x))) : x
+    ctx[:yplanes] =
+        isa(y, Number) ?
+        collect(range(xyzmin[2] + ε[2], xyzmax[2] - ε[2]; length = ceil(y))) : y
+    ctx[:zplanes] =
+        isa(z, Number) ?
+        collect(range(xyzmin[3] + ε[3], xyzmax[3] - ε[3]; length = ceil(z))) : z
 
     ctx[:xplanes][1] = min(xyzmax[1], ctx[:xplanes][1])
     ctx[:yplanes][1] = min(xyzmax[2], ctx[:yplanes][1])
@@ -966,78 +1175,114 @@ function scalarplot!(ctx, TP::Type{MakieType}, ::Type{Val{3}}, grid, func)
     ctx[:levels] = levels
 
     if !haskey(ctx, :scene)
-        ctx[:data] = Observable((g = grid, f = func, x = ctx[:xplanes], y = ctx[:yplanes],
-                                 z = ctx[:zplanes], l = ctx[:levels], t = ctx[:title]))
+        ctx[:data] = Observable((
+            g = grids,
+            p = parentgrid,
+            f = funcs,
+            x = ctx[:xplanes],
+            y = ctx[:yplanes],
+            z = ctx[:zplanes],
+            l = ctx[:levels],
+            t = ctx[:title],
+        ))
 
         ctx[:scene] = makeaxis3d(ctx)
 
-        #### Transparent outlne
+        #### Transparent outline
         if ctx[:outlinealpha] > 0.0
-            ctx[:outlinedata] = map(d -> extract_visible_bfaces3D(d.g,
-                                                                  xyzmax;
-                                                                  primepoints = hcat(xyzmin,
-                                                                                     xyzmax),
-                                                                  Tp = Point3f,
-                                                                  Tf = GLTriangleFace),
-                                    ctx[:data])
-            ctx[:facemeshes] = map(d -> [make_mesh(d[1][i], d[2][i]) for i = 1:nbregions],
-                                   ctx[:outlinedata])
+            ctx[:outlinedata] = map(
+                d -> extract_visible_bfaces3D(
+                    d.p,
+                    xyzmax;
+                    primepoints = hcat(xyzmin, xyzmax),
+                    Tp = Point3f,
+                    Tf = GLTriangleFace,
+                ),
+                ctx[:data],
+            )
+            ctx[:facemeshes] = map(
+                d -> [make_mesh(d[1][i], d[2][i]) for i = 1:nbregions],
+                ctx[:outlinedata],
+            )
             bcmap = bregion_cmap(nbregions)
             for i = 1:nbregions
-                XMakie.mesh!(ctx[:scene], map(d -> d[i], ctx[:facemeshes]);
-                             color = (bcmap[i], ctx[:outlinealpha]),
-                             transparency = true,
-                             backlight = 1.0f0)
+                XMakie.mesh!(
+                    ctx[:scene],
+                    map(d -> d[i], ctx[:facemeshes]);
+                    color = (bcmap[i], ctx[:outlinealpha]),
+                    transparency = true,
+                    backlight = 1.0f0,
+                )
             end
         end
 
-        f = d -> make_mesh(marching_tetrahedra(d.g,
-                                               d.f,
-                                               makeplanes(xyzmin, xyzmax, d.x, d.y, d.z),
-                                               d.l;
-                                               primepoints = hcat(xyzmin, xyzmax),
-                                               primevalues = crange,
-                                               tol = ctx[:tetxplane_tol],
-                                               Tp = Point3f,
-                                               Tf = GLTriangleFace,
-                                               Tv = Float32)...)
+        f =
+            d -> make_mesh(
+                marching_tetrahedra(
+                    d.g,
+                    d.f,
+                    makeplanes(xyzmin, xyzmax, d.x, d.y, d.z),
+                    d.l;
+                    primepoints = hcat(xyzmin, xyzmax),
+                    primevalues = crange,
+                    tol = ctx[:tetxplane_tol],
+                    Tp = Point3f,
+                    Tf = GLTriangleFace,
+                    Tv = Float32,
+                )...,
+            )
 
         #### Plane sections and isosurfaces
-        ctx[:mesh] = XMakie.mesh!(ctx[:scene], map(f, ctx[:data]); backlight = 1.0f0,
-                                  transparency = ctx[:levelalpha] < 1.0)
+        ctx[:mesh] = XMakie.mesh!(
+            ctx[:scene],
+            map(f, ctx[:data]);
+            backlight = 1.0f0,
+            transparency = ctx[:levelalpha] < 1.0,
+        )
         #### Interactions
         scene_interaction(ctx[:scene].scene, XMakie, [:z, :y, :x, :l, :q]) do delta, key
             if key == :x
                 ctx[:xplanes] .+= delta * xyzstep[1]
-                ctx[:status][] = "x=[" *
-                                 mapreduce(x -> @sprintf("%.3g,", x), *, ctx[:xplanes]) *
-                                 "]"
+                ctx[:status][] =
+                    "x=[" * mapreduce(x -> @sprintf("%.3g,", x), *, ctx[:xplanes]) * "]"
             elseif key == :y
                 ctx[:yplanes] .+= delta * xyzstep[2]
-                ctx[:status][] = "y=[" *
-                                 mapreduce(y -> @sprintf("%.3g,", y), *, ctx[:yplanes]) *
-                                 "]"
+                ctx[:status][] =
+                    "y=[" * mapreduce(y -> @sprintf("%.3g,", y), *, ctx[:yplanes]) * "]"
             elseif key == :z
                 ctx[:zplanes] .+= delta * xyzstep[3]
-                ctx[:status][] = "z=[" *
-                                 mapreduce(z -> @sprintf("%.3g,", z), *, ctx[:zplanes]) *
-                                 "]"
+                ctx[:status][] =
+                    "z=[" * mapreduce(z -> @sprintf("%.3g,", z), *, ctx[:zplanes]) * "]"
             elseif key == :l
                 ctx[:levels] .+= delta * fstep
-                ctx[:status][] = "l=[" *
-                                 mapreduce(l -> @sprintf("%.3g,", l), *, ctx[:levels]) * "]"
+                ctx[:status][] =
+                    "l=[" * mapreduce(l -> @sprintf("%.3g,", l), *, ctx[:levels]) * "]"
             elseif key == :q
                 ctx[:status][] = " "
             end
             #            adjust_planes()
-            ctx[:data][] = (g = grid, f = func, x = ctx[:xplanes], y = ctx[:yplanes],
-                            z = ctx[:zplanes], l = ctx[:levels], t = ctx[:title])
+            ctx[:data][] = (
+                g = grid,
+                f = func,
+                x = ctx[:xplanes],
+                y = ctx[:yplanes],
+                z = ctx[:zplanes],
+                l = ctx[:levels],
+                t = ctx[:title],
+            )
         end
         ctx[:status] = Observable(" ")
         add_scene!(ctx, makescene3d(ctx))
     else
-        ctx[:data][] = (g = grid, f = func, x = ctx[:xplanes], y = ctx[:yplanes],
-                        z = ctx[:zplanes], l = ctx[:levels], t = ctx[:title])
+        ctx[:data][] = (
+            g = grid,
+            f = func,
+            x = ctx[:xplanes],
+            y = ctx[:yplanes],
+            z = ctx[:zplanes],
+            l = ctx[:levels],
+            t = ctx[:title],
+        )
     end
     reveal(ctx, TP)
 end
diff --git a/src/meshcat.jl b/src/meshcat.jl
index 8d03e18..9ec8115 100644
--- a/src/meshcat.jl
+++ b/src/meshcat.jl
@@ -1,11 +1,13 @@
 function initialize!(p::GridVisualizer, ::Type{MeshCatType})
     MeshCat = p.context[:Plotter]
     layout = p.context[:layout]
-    @assert(layout==(1, 1))
+    @assert(layout == (1, 1))
     vis = MeshCat.Visualizer()
     MeshCat.send(vis.core, MeshCat.SetProperty(MeshCat.Path(["Grid"]), "visible", false))
-    MeshCat.send(vis.core,
-                 MeshCat.SetProperty(MeshCat.Path(["Background"]), "visible", false))
+    MeshCat.send(
+        vis.core,
+        MeshCat.SetProperty(MeshCat.Path(["Background"]), "visible", false),
+    )
     p.context[:scene] = vis
     for I in CartesianIndices(layout)
         ctx = p.subplots[I]
@@ -25,7 +27,7 @@ function reveal(ctx::SubVisualizer, TP::Type{MeshCatType})
 end
 
 gridplot!(ctx, TP::Type{MeshCatType}, ::Type{Val{1}}, grid) = nothing
-scalarplot!(ctx, TP::Type{MeshCatType}, ::Type{Val{1}}, grid, func) = nothing
+scalarplot!(ctx, TP::Type{MeshCatType}, ::Type{Val{1}}, grids, parentgrid, funcs) = nothing
 
 # 2D grid
 function gridplot!(ctx, TP::Type{MeshCatType}, ::Type{Val{2}}, grid)
@@ -39,15 +41,19 @@ function gridplot!(ctx, TP::Type{MeshCatType}, ::Type{Val{2}}, grid)
     bcmap = bregion_cmap(nbregions)
     for i = 1:nregions
         mesh = regionmesh(grid, i)
-        MeshCat.setobject!(vis["interior"]["r$(i)"],
-                           mesh,
-                           MeshCat.MeshLambertMaterial(;
-                                                       color = RGBA{Float32}(cmap[i], 1.0)))
-        MeshCat.setobject!(vis["interior"]["r$(i)_edges"], mesh,
-                           MeshCat.MeshPhongMaterial(;
-                                                     color = RGBA{Float32}(0.0, 0.0, 0.0,
-                                                                           1.0),
-                                                     wireframe = true))
+        MeshCat.setobject!(
+            vis["interior"]["r$(i)"],
+            mesh,
+            MeshCat.MeshLambertMaterial(; color = RGBA{Float32}(cmap[i], 1.0)),
+        )
+        MeshCat.setobject!(
+            vis["interior"]["r$(i)_edges"],
+            mesh,
+            MeshCat.MeshPhongMaterial(;
+                color = RGBA{Float32}(0.0, 0.0, 0.0, 1.0),
+                wireframe = true,
+            ),
+        )
     end
 
     for i = 1:nbregions
@@ -85,50 +91,63 @@ function gridplot!(ctx, TP::Type{MeshCatType}, ::Type{Val{3}}, grid)
     xyzcut = [ctx[:xplane], ctx[:yplane], ctx[:zplane]]
 
     if ctx[:interior]
-        pts, fcs = extract_visible_cells3D(grid,
-                                           xyzcut;
-                                           primepoints = hcat(xyzmin, xyzmax),
-                                           Tp = Point3f,
-                                           Tf = GLTriangleFace)
+        pts, fcs = extract_visible_cells3D(
+            grid,
+            xyzcut;
+            primepoints = hcat(xyzmin, xyzmax),
+            Tp = Point3f,
+            Tf = GLTriangleFace,
+        )
 
         for i = 1:nregions
             mesh = Mesh(pts[i], fcs[i])
-            MeshCat.setobject!(vis["r$(i)"],
-                               mesh,
-                               MeshCat.MeshLambertMaterial(;
-                                                           color = RGBA{Float32}(cmap[i],
-                                                                                 1.0)))
-            MeshCat.setobject!(vis["r$(i)_edges"], mesh,
-                               MeshCat.MeshPhongMaterial(;
-                                                         color = RGBA{Float32}(0.0, 0.0,
-                                                                               0.0, 1.0),
-                                                         wireframe = true))
+            MeshCat.setobject!(
+                vis["r$(i)"],
+                mesh,
+                MeshCat.MeshLambertMaterial(; color = RGBA{Float32}(cmap[i], 1.0)),
+            )
+            MeshCat.setobject!(
+                vis["r$(i)_edges"],
+                mesh,
+                MeshCat.MeshPhongMaterial(;
+                    color = RGBA{Float32}(0.0, 0.0, 0.0, 1.0),
+                    wireframe = true,
+                ),
+            )
         end
     end
 
-    pts, fcs = extract_visible_bfaces3D(grid,
-                                        xyzcut;
-                                        primepoints = hcat(xyzmin, xyzmax),
-                                        Tp = Point3f,
-                                        Tf = GLTriangleFace)
+    pts, fcs = extract_visible_bfaces3D(
+        grid,
+        xyzcut;
+        primepoints = hcat(xyzmin, xyzmax),
+        Tp = Point3f,
+        Tf = GLTriangleFace,
+    )
 
     for i = 1:nbregions
         mesh = Mesh(pts[i], fcs[i])
-        MeshCat.setobject!(vis["b$(i)"],
-                           mesh,
-                           MeshCat.MeshLambertMaterial(;
-                                                       color = RGBA{Float32}(bcmap[i], 1.0)))
-        MeshCat.setobject!(vis["b$(i)_edges"], mesh,
-                           MeshCat.MeshPhongMaterial(;
-                                                     color = RGBA{Float32}(0.0, 0.0, 0.0,
-                                                                           1.0),
-                                                     wireframe = true))
+        MeshCat.setobject!(
+            vis["b$(i)"],
+            mesh,
+            MeshCat.MeshLambertMaterial(; color = RGBA{Float32}(bcmap[i], 1.0)),
+        )
+        MeshCat.setobject!(
+            vis["b$(i)_edges"],
+            mesh,
+            MeshCat.MeshPhongMaterial(;
+                color = RGBA{Float32}(0.0, 0.0, 0.0, 1.0),
+                wireframe = true,
+            ),
+        )
     end
 
     reveal(ctx, TP)
 end
 
-function scalarplot!(ctx, TP::Type{MeshCatType}, ::Type{Val{3}}, grid, func)
+function scalarplot!(ctx, TP::Type{MeshCatType}, ::Type{Val{3}}, grids, parentgrid, funcs)
+    grid = parentgrid
+    func = funcs[1]
     MeshCat = ctx[:Plotter]
     vis = ctx[:figure]
 
@@ -150,20 +169,19 @@ function scalarplot!(ctx, TP::Type{MeshCatType}, ::Type{Val{3}}, grid, func)
     ctx[:zplane] = max(xyzmin[3], min(xyzmax[3], ctx[:zplane]))
     ctx[:flevel] = max(fminmax[1], min(fminmax[2], ctx[:flevel]))
 
-    makeplanes(x, y, z) = [[1, 0, 0, -x],
-        [0, 1, 0, -y],
-        [0, 0, 1, -z]]
-
-    ccoord, faces, values = marching_tetrahedra(grid,
-                                                func,
-                                                makeplanes(ctx[:xplane], ctx[:yplane],
-                                                           ctx[:zplane]),
-                                                [ctx[:flevel]];
-                                                primepoints = hcat(xyzmin, xyzmax),
-                                                primevalues = fminmax,
-                                                Tp = Point3f,
-                                                Tf = GLTriangleFace,
-                                                Tv = Float32)
+    makeplanes(x, y, z) = [[1, 0, 0, -x], [0, 1, 0, -y], [0, 0, 1, -z]]
+
+    ccoord, faces, values = marching_tetrahedra(
+        grid,
+        func,
+        makeplanes(ctx[:xplane], ctx[:yplane], ctx[:zplane]),
+        [ctx[:flevel]];
+        primepoints = hcat(xyzmin, xyzmax),
+        primevalues = fminmax,
+        Tp = Point3f,
+        Tf = GLTriangleFace,
+        Tv = Float32,
+    )
     mesh = Mesh(ccoord, faces)
 
     to01(v) = (v - fminmax[1]) / (fminmax[2] - fminmax[1])
@@ -176,18 +194,22 @@ function scalarplot!(ctx, TP::Type{MeshCatType}, ::Type{Val{3}}, grid, func)
     MeshCat.setobject!(vis[:marching_tets], mesh_meta, material)
 
     if ctx[:outline]
-        pts, fcs = extract_visible_bfaces3D(grid,
-                                            xyzmax;
-                                            primepoints = hcat(xyzmin, xyzmax),
-                                            Tp = Point3f,
-                                            Tf = GLTriangleFace)
+        pts, fcs = extract_visible_bfaces3D(
+            grid,
+            xyzmax;
+            primepoints = hcat(xyzmin, xyzmax),
+            Tp = Point3f,
+            Tf = GLTriangleFace,
+        )
         for i = 1:nbregions
             mesh = Mesh(pts[i], fcs[i])
-            MeshCat.setobject!(vis["b$(i)"],
-                               mesh,
-                               MeshCat.MeshLambertMaterial(;
-                                                           color = color = RGBA{Float32}(bcmap[i],
-                                                                                         0.35)))
+            MeshCat.setobject!(
+                vis["b$(i)"],
+                mesh,
+                MeshCat.MeshLambertMaterial(;
+                    color = color = RGBA{Float32}(bcmap[i], 0.35),
+                ),
+            )
         end
     end
     reveal(ctx, TP)
diff --git a/src/plots.jl b/src/plots.jl
index 340597a..c1b52ef 100644
--- a/src/plots.jl
+++ b/src/plots.jl
@@ -27,8 +27,8 @@ function reveal(p::GridVisualizer, ::Type{PlotsType})
         end
     end
     plt = Plots.plot(subplots...; layout = p.context[:layout], size = p.context[:size])
-    if haskey(p.context,:videostream)
-        Plots.frame(p.context[:videostream],plt)
+    if haskey(p.context, :videostream)
+        Plots.frame(p.context[:videostream], plt)
     else
         Plots.gui(plt)
         plt
@@ -47,8 +47,15 @@ function save(fname, scene, Plots, ::Type{PlotsType})
 end
 
 
-function movie(func, vis::GridVisualizer, ::Type{PlotsType}; file=nothing, format="gif" ,kwargs...)
-    Plotter=vis.context[:Plotter]
+function movie(
+    func,
+    vis::GridVisualizer,
+    ::Type{PlotsType};
+    file = nothing,
+    format = "gif",
+    kwargs...,
+)
+    Plotter = vis.context[:Plotter]
     if !isnothing(file)
         format = lstrip(splitext(file)[2], '.')
     end
@@ -57,15 +64,15 @@ function movie(func, vis::GridVisualizer, ::Type{PlotsType}; file=nothing, forma
     end
 
     func(vis)
-    
+
     if !isnothing(file)
-        if format=="mp4"
-            Plotter.mp4(vis.context[:videostream],file)
+        if format == "mp4"
+            Plotter.mp4(vis.context[:videostream], file)
         else
-            Plotter.gif(vis.context[:videostream],file)
+            Plotter.gif(vis.context[:videostream], file)
         end
     elseif isdefined(Main, :PlutoRunner)
-        if format=="mp4"
+        if format == "mp4"
             Plotter.mp4(vis.context[:videostream])
         else
             Plotter.gif(vis.context[:videostream])
@@ -102,16 +109,28 @@ function gridplot!(ctx, TP::Type{PlotsType}, ::Type{Val{1}}, grid)
         x2 = coord[1, cellnodes[2, icell]]
         Plots.plot!(p, [x1, x1], [-h, h]; linewidth = 0.5, color = :black, label = "")
         Plots.plot!(p, [x2, x2], [-h, h]; linewidth = 0.5, color = :black, label = "")
-        Plots.plot!(p, [x1, x2], [0, 0]; linewidth = 3.0, color = cmap[cellregions[icell]],
-                    label = "")
+        Plots.plot!(
+            p,
+            [x1, x2],
+            [0, 0];
+            linewidth = 3.0,
+            color = cmap[cellregions[icell]],
+            label = "",
+        )
     end
 
     cmap = bregion_cmap(nbfaceregions)
     for ibface = 1:num_bfaces(grid)
         if bfaceregions[ibface] > 0
             x1 = coord[1, bfacenodes[1, ibface]]
-            Plots.plot!(p, [x1, x1], [-2 * h, 2 * h]; linewidth = 3.0,
-                        color = cmap[bfaceregions[ibface]], label = "")
+            Plots.plot!(
+                p,
+                [x1, x1],
+                [-2 * h, 2 * h];
+                linewidth = 3.0,
+                color = cmap[bfaceregions[ibface]],
+                label = "",
+            )
         end
     end
     reveal(ctx, TP)
@@ -137,48 +156,73 @@ function gridplot!(ctx, TP::Type{PlotsType}, ::Type{Val{2}}, grid)
         inode2 = cellnodes[2, icell]
         inode3 = cellnodes[3, icell]
         # https://github.com/JuliaPlots/Plots.jl/issues/605
-        tri = Plots.Shape([coord[1, inode1], coord[1, inode2], coord[1, inode3]],
-                          [coord[2, inode1], coord[2, inode2], coord[2, inode3]])
+        tri = Plots.Shape(
+            [coord[1, inode1], coord[1, inode2], coord[1, inode3]],
+            [coord[2, inode1], coord[2, inode2], coord[2, inode3]],
+        )
         Plots.plot!(p, tri; color = cmap[cellregions[icell]], label = "")
     end
     for icell = 1:num_cells(grid)
         inode1 = cellnodes[1, icell]
         inode2 = cellnodes[2, icell]
         inode3 = cellnodes[3, icell]
-        Plots.plot!(p, [coord[1, inode1], coord[1, inode2]],
-                    [coord[2, inode1], coord[2, inode2]]; linewidth = 0.5, color = :black,
-                    label = "")
-        Plots.plot!(p, [coord[1, inode1], coord[1, inode3]],
-                    [coord[2, inode1], coord[2, inode3]]; linewidth = 0.5, color = :black,
-                    label = "")
-        Plots.plot!(p, [coord[1, inode2], coord[1, inode3]],
-                    [coord[2, inode2], coord[2, inode3]]; linewidth = 0.5, color = :black,
-                    label = "")
+        Plots.plot!(
+            p,
+            [coord[1, inode1], coord[1, inode2]],
+            [coord[2, inode1], coord[2, inode2]];
+            linewidth = 0.5,
+            color = :black,
+            label = "",
+        )
+        Plots.plot!(
+            p,
+            [coord[1, inode1], coord[1, inode3]],
+            [coord[2, inode1], coord[2, inode3]];
+            linewidth = 0.5,
+            color = :black,
+            label = "",
+        )
+        Plots.plot!(
+            p,
+            [coord[1, inode2], coord[1, inode3]],
+            [coord[2, inode2], coord[2, inode3]];
+            linewidth = 0.5,
+            color = :black,
+            label = "",
+        )
     end
 
     cmap = bregion_cmap(nbfaceregions)
     for ibface = 1:num_bfaces(grid)
         inode1 = bfacenodes[1, ibface]
         inode2 = bfacenodes[2, ibface]
-        Plots.plot!(p, [coord[1, inode1], coord[1, inode2]],
-                    [coord[2, inode1], coord[2, inode2]]; linewidth = 5,
-                    color = cmap[bfaceregions[ibface]],
-                    label = "")
+        Plots.plot!(
+            p,
+            [coord[1, inode1], coord[1, inode2]],
+            [coord[2, inode1], coord[2, inode2]];
+            linewidth = 5,
+            color = cmap[bfaceregions[ibface]],
+            label = "",
+        )
     end
     reveal(ctx, TP)
 end
 
-function scalarplot!(ctx, TP::Type{PlotsType}, ::Type{Val{1}}, grid, func)
-    legpos = Dict(:none => :none,
-                  :best => :best,
-                  :lt => :topleft,
-                  :ct => :topcenter,
-                  :rt => :topright,
-                  :lc => :centerleft,
-                  :rc => :centerright,
-                  :lb => :bottomleft,
-                  :cb => :bottomcenter,
-                  :rb => :bottomright)
+function scalarplot!(ctx, TP::Type{PlotsType}, ::Type{Val{1}}, grids, parentgrid, funcs)
+    grid = parentgrid
+    func = funcs[1]
+    legpos = Dict(
+        :none => :none,
+        :best => :best,
+        :lt => :topleft,
+        :ct => :topcenter,
+        :rt => :topright,
+        :lc => :centerleft,
+        :rc => :centerright,
+        :lb => :bottomleft,
+        :cb => :bottomcenter,
+        :rb => :bottomright,
+    )
 
     Plots = ctx[:Plotter]
     if !haskey(ctx, :ax)
@@ -205,11 +249,16 @@ function scalarplot!(ctx, TP::Type{PlotsType}, ::Type{Val{1}}, grid, func)
     ctx[:xscale] == :log ? ctx[:xscale] = :log10 : nothing
     ctx[:yscale] == :log ? ctx[:yscale] = :log10 : nothing
 
-    Plots.plot!(p, [xmin, xmax], [ymin, ymax]; seriestype = :scatter,
-                makersize = 0,
-                markercolor = :white,
-                markerstrokecolor = :white,
-                label = "")
+    Plots.plot!(
+        p,
+        [xmin, xmax],
+        [ymin, ymax];
+        seriestype = :scatter,
+        makersize = 0,
+        markercolor = :white,
+        markerstrokecolor = :white,
+        label = "",
+    )
 
     color = ctx[:color]
     if ctx[:cellwise] ## not checked
@@ -220,11 +269,21 @@ function scalarplot!(ctx, TP::Type{PlotsType}, ::Type{Val{1}}, grid, func)
             x1 = coord[1, i1]
             x2 = coord[1, i2]
             if icell == 1 && ctx[:label] != " "
-                Plots.plot!(p, [x1, x2], [func[i1], func[i2]];
-                            linecolor = Plots.RGB(color...), label = ctx[:label])
+                Plots.plot!(
+                    p,
+                    [x1, x2],
+                    [func[i1], func[i2]];
+                    linecolor = Plots.RGB(color...),
+                    label = ctx[:label],
+                )
             else
-                Plots.plot!(p, [x1, x2], [func[i1], func[i2]];
-                            linecolor = Plots.RGB(color...), label = "")
+                Plots.plot!(
+                    p,
+                    [x1, x2],
+                    [func[i1], func[i2]];
+                    linecolor = Plots.RGB(color...),
+                    label = "",
+                )
             end
         end
     else
@@ -232,43 +291,59 @@ function scalarplot!(ctx, TP::Type{PlotsType}, ::Type{Val{1}}, grid, func)
         markershape = ctx[:markershape]
         X = vec(grid[Coordinates])
         if markershape == :none
-            Plots.plot!(p, X, func;
-                        linecolor = Plots.RGB(color),
-                        linewidth = ctx[:linewidth],
-                        linestyle = ctx[:linestyle],
-                        legend = legpos[ctx[:legend]],
-                        xscale = ctx[:xscale],
-                        yscale = ctx[:yscale],
-                        label = ctx[:label])
+            Plots.plot!(
+                p,
+                X,
+                func;
+                linecolor = Plots.RGB(color),
+                linewidth = ctx[:linewidth],
+                linestyle = ctx[:linestyle],
+                legend = legpos[ctx[:legend]],
+                xscale = ctx[:xscale],
+                yscale = ctx[:yscale],
+                label = ctx[:label],
+            )
         else
             #Trick plots to use markers
-            Plots.plot!(p, X, func;
-                        linecolor = Plots.RGB(color),
-                        linewidth = ctx[:linewidth],
-                        linestyle = ctx[:linestyle],
-                        xscale = ctx[:xscale],
-                        yscale = ctx[:yscale],
-                        label = "")
-            Plots.plot!(p, [X[1]], [func[1]];
-                        markershape = markershape,
-                        label = ctx[:label],
-                        markersize = ctx[:markersize],
-                        linecolor = Plots.RGB(color),
-                        linewidth = ctx[:linewidth],
-                        xscale = ctx[:xscale],
-                        yscale = ctx[:yscale],
-                        legend = legpos[ctx[:legend]],
-                        linestyle = ctx[:linestyle],
-                        markercolor = Plots.RGB(color))
-            @views Plots.plot!(p, X[1:markevery:end], func[1:markevery:end],
-                               markercolor = Plots.RGB(color),
-                               label = "",
-                               linecolor = :white,
-                               xscale = ctx[:xscale],
-                               yscale = ctx[:yscale],
-                               markershape = markershape,
-                               markersize = ctx[:markersize],
-                               lines = false)
+            Plots.plot!(
+                p,
+                X,
+                func;
+                linecolor = Plots.RGB(color),
+                linewidth = ctx[:linewidth],
+                linestyle = ctx[:linestyle],
+                xscale = ctx[:xscale],
+                yscale = ctx[:yscale],
+                label = "",
+            )
+            Plots.plot!(
+                p,
+                [X[1]],
+                [func[1]];
+                markershape = markershape,
+                label = ctx[:label],
+                markersize = ctx[:markersize],
+                linecolor = Plots.RGB(color),
+                linewidth = ctx[:linewidth],
+                xscale = ctx[:xscale],
+                yscale = ctx[:yscale],
+                legend = legpos[ctx[:legend]],
+                linestyle = ctx[:linestyle],
+                markercolor = Plots.RGB(color),
+            )
+            @views Plots.plot!(
+                p,
+                X[1:markevery:end],
+                func[1:markevery:end],
+                markercolor = Plots.RGB(color),
+                label = "",
+                linecolor = :white,
+                xscale = ctx[:xscale],
+                yscale = ctx[:yscale],
+                markershape = markershape,
+                markersize = ctx[:markersize],
+                lines = false,
+            )
         end
     end
     reveal(ctx, TP)
@@ -290,7 +365,9 @@ function rectdata(grid, U)
     nothing
 end
 
-function scalarplot!(ctx, TP::Type{PlotsType}, ::Type{Val{2}}, grid, func)
+function scalarplot!(ctx, TP::Type{PlotsType}, ::Type{Val{2}}, grids, parentgrid, funcs)
+    grid = parentgrid
+    func = funcs[1]
     rdata = rectdata(grid, func)
     if rdata == nothing
         return nothing
@@ -302,10 +379,17 @@ function scalarplot!(ctx, TP::Type{PlotsType}, ::Type{Val{2}}, grid, func)
     p = ctx[:ax]
 
     levels, crange, colorbarticks = isolevels(ctx, func)
-    colorlevels = collect(crange[1]:((crange[2] - crange[1]) / (ctx[:colorlevels] - 1)):crange[2])
-
-    Plots.contourf!(p, rdata...; aspect_ratio = ctx[:aspect], fill = ctx[:colormap],
-                    linewidth = 0, levels = colorlevels, colorbar_ticks = colorbarticks)
+    colorlevels = collect(crange[1]:((crange[2]-crange[1])/(ctx[:colorlevels]-1)):crange[2])
+
+    Plots.contourf!(
+        p,
+        rdata...;
+        aspect_ratio = ctx[:aspect],
+        fill = ctx[:colormap],
+        linewidth = 0,
+        levels = colorlevels,
+        colorbar_ticks = colorbarticks,
+    )
     Plots.contour!(p, rdata...; aspect_ratio = ctx[:aspect], c = :black, levels = levels)
     reveal(ctx, TP)
 end
@@ -333,7 +417,9 @@ function gridplot!(ctx, TP::Type{PlotsType}, ::Type{Val{3}}, grid)
     reveal(ctx, TP)
 end
 
-function scalarplot!(ctx, TP::Type{PlotsType}, ::Type{Val{3}}, grid, func)
+function scalarplot!(ctx, TP::Type{PlotsType}, ::Type{Val{3}}, grids, parentgrid, funcs)
+    grid = parentgrid
+    func = funcs[1]
     Plots = ctx[:Plotter]
     if !haskey(ctx, :ax)
         ctx[:ax] = Plots.plot(; title = ctx[:title])
diff --git a/src/plutovista.jl b/src/plutovista.jl
index 569ec0a..46b7a6e 100644
--- a/src/plutovista.jl
+++ b/src/plutovista.jl
@@ -7,8 +7,11 @@ function initialize!(p::GridVisualizer, ::Type{PlutoVistaType})
     for I in CartesianIndices(layout)
         ctx = p.subplots[I]
         ctx[:figure] = PlutoVista.PlutoVistaPlot(; resolution = figres)
-        PlutoVista.backend!(ctx[:figure]; backend = p.context[:backend],
-                            datadim = p.context[:dim])
+        PlutoVista.backend!(
+            ctx[:figure];
+            backend = p.context[:backend],
+            datadim = p.context[:dim],
+        )
     end
 end
 
@@ -25,61 +28,75 @@ function reveal(p::GridVisualizer, ::Type{PlutoVistaType})
     if l == 1
         subplots[1][:figure]
     elseif l == 2
-        @htl("""<div><div style=" display: inline-block;">
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[1][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[2][:figure])</div>
-        </div>""")
+        @htl(
+            """<div><div style=" display: inline-block;">
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[1][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[2][:figure])</div>
+       </div>"""
+        )
     elseif l == 3
-        @htl("""<div><div style=" display: inline-block;">
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[1][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[2][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[3][:figure])</div>
-        </div>""")
+        @htl(
+            """<div><div style=" display: inline-block;">
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[1][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[2][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[3][:figure])</div>
+       </div>"""
+        )
     elseif l == 4
-        @htl("""<div><div style=" display: inline-block;">
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[1][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[2][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[3][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[4][:figure])</div>
-        </div>""")
+        @htl(
+            """<div><div style=" display: inline-block;">
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[1][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[2][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[3][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[4][:figure])</div>
+       </div>"""
+        )
     elseif l == 5
-        @htl("""<div><div style=" display: inline-block;">
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[1][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[2][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[3][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[4][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[5][:figure])</div>
-        </div>""")
+        @htl(
+            """<div><div style=" display: inline-block;">
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[1][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[2][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[3][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[4][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[5][:figure])</div>
+       </div>"""
+        )
     elseif l == 6
-        @htl("""<div><div style=" display: inline-block;">
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[1][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[2][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[3][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[4][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[5][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[6][:figure])</div>
-        </div>""")
+        @htl(
+            """<div><div style=" display: inline-block;">
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[1][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[2][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[3][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[4][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[5][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[6][:figure])</div>
+       </div>"""
+        )
     elseif l == 7
-        @htl("""<div><div style=" display: inline-block;">
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[1][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[2][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[3][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[4][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[5][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[6][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[7][:figure])</div>
-        </div>""")
+        @htl(
+            """<div><div style=" display: inline-block;">
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[1][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[2][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[3][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[4][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[5][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[6][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[7][:figure])</div>
+       </div>"""
+        )
     elseif l == 8
-        @htl("""<div><div style=" display: inline-block;">
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[1][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[2][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[3][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[4][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[5][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[6][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[7][:figure])</div>
-        <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[8][:figure])</div>
-        </div>""")
+        @htl(
+            """<div><div style=" display: inline-block;">
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[1][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[2][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[3][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[4][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[5][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[6][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[7][:figure])</div>
+       <div style= "width: $(figwidth)px;  display: inline-block;">$(subplots[8][:figure])</div>
+       </div>"""
+        )
     end
 end
 
@@ -124,12 +141,31 @@ function gridplot!(ctx, TP::Type{PlutoVistaType}, ::Type{Val{1}}, grid)
         x1 = coord[1, cellnodes[1, icell]]
         x2 = coord[1, cellnodes[2, icell]]
 
-        PlutoVista.plot!(ctx[:figure], [x1, x2], [0, 0]; clear = false, linewidth = 3.0,
-                         color = rgbtuple(cmap[cellregions[icell]]), label = label)
-        PlutoVista.plot!(ctx[:figure], [x1, x1], [-h, h]; clear = false,
-                         linewidth = ctx[:linewidth], color = :black)
-        PlutoVista.plot!(ctx[:figure], [x2, x2], [-h, h]; clear = false,
-                         linewidth = ctx[:linewidth], color = :black)
+        PlutoVista.plot!(
+            ctx[:figure],
+            [x1, x2],
+            [0, 0];
+            clear = false,
+            linewidth = 3.0,
+            color = rgbtuple(cmap[cellregions[icell]]),
+            label = label,
+        )
+        PlutoVista.plot!(
+            ctx[:figure],
+            [x1, x1],
+            [-h, h];
+            clear = false,
+            linewidth = ctx[:linewidth],
+            color = :black,
+        )
+        PlutoVista.plot!(
+            ctx[:figure],
+            [x2, x2],
+            [-h, h];
+            clear = false,
+            linewidth = ctx[:linewidth],
+            color = :black,
+        )
     end
 
     cmap = bregion_cmap(nbfaceregions)
@@ -139,37 +175,81 @@ function gridplot!(ctx, TP::Type{PlutoVistaType}, ::Type{Val{1}}, grid)
             label = brflag[ireg] ? "b$(ireg)" : ""
             brflag[ireg] = false
             x1 = coord[1, bfacenodes[1, ibface]]
-            PlutoVista.plot!(ctx[:figure], [x1, x1], [-2 * h, 2 * h]; clear = false,
-                             linewidth = 3.0, color = rgbtuple(cmap[ireg]), label = label,
-                             legend = leglocs[ctx[:legend]])
+            PlutoVista.plot!(
+                ctx[:figure],
+                [x1, x1],
+                [-2 * h, 2 * h];
+                clear = false,
+                linewidth = 3.0,
+                color = rgbtuple(cmap[ireg]),
+                label = label,
+                legend = leglocs[ctx[:legend]],
+            )
         end
     end
 
     reveal(ctx, TP)
 end
 
-function scalarplot!(ctx, TP::Type{PlutoVistaType}, ::Type{Val{1}}, grid, func)
+function scalarplot!(
+    ctx,
+    TP::Type{PlutoVistaType},
+    ::Type{Val{1}},
+    grids,
+    parentgrid,
+    funcs,
+)
+    nfuncs = length(funcs)
+
     PlutoVista = ctx[:Plotter]
-    coord = grid[Coordinates]
     PlutoVista.backend!(ctx[:figure]; backend = ctx[:backend], datadim = 1)
-    PlutoVista.plot!(ctx[:figure],
-                     coord[1, :],
-                     func;
-                     color = ctx[:color],
-                     markertype = ctx[:markershape],
-                     markercount = length(func) ÷ ctx[:markevery],
-                     linestyle = ctx[:linestyle],
-                     xlabel = ctx[:xlabel],
-                     ylabel = ctx[:ylabel],
-                     label = ctx[:label],
-                     linewidth = ctx[:linewidth],
-                     legend = ctx[:legend],
-                     xlimits = ctx[:xlimits],
-                     limits = ctx[:limits],
-                     clear = ctx[:clear],
-                     title = ctx[:title],
-                     xscale = ctx[:xscale],
-                     yscale = ctx[:yscale])
+
+    for ifunc = 1:nfuncs
+        func = funcs[ifunc]
+        grid = grids[ifunc]
+        coord = grid[Coordinates]
+
+        if ifunc == 1
+            PlutoVista.plot!(
+                ctx[:figure],
+                coord[1, :],
+                func;
+                color = ctx[:color],
+                markertype = ctx[:markershape],
+                markercount = length(func) ÷ ctx[:markevery],
+                linestyle = ctx[:linestyle],
+                xlabel = ctx[:xlabel],
+                ylabel = ctx[:ylabel],
+                label = ctx[:label],
+                linewidth = ctx[:linewidth],
+                legend = ctx[:legend],
+                xlimits = ctx[:xlimits],
+                limits = ctx[:limits],
+                clear = ctx[:clear],
+                title = ctx[:title],
+                xscale = ctx[:xscale],
+                yscale = ctx[:yscale],
+            )
+        else
+            PlutoVista.plot!(
+                ctx[:figure],
+                coord[1, :],
+                func;
+                color = ctx[:color],
+                markertype = ctx[:markershape],
+                markercount = length(func) ÷ ctx[:markevery],
+                linestyle = ctx[:linestyle],
+                xlabel = ctx[:xlabel],
+                ylabel = ctx[:ylabel],
+                linewidth = ctx[:linewidth],
+                xlimits = ctx[:xlimits],
+                limits = ctx[:limits],
+                clear = false,
+                xscale = ctx[:xscale],
+                yscale = ctx[:yscale],
+            )
+        end
+    end
     reveal(ctx, TP)
 end
 
@@ -187,35 +267,76 @@ function gridplot!(ctx, TP::Type{PlutoVistaType}, ::Type{Val{2}}, grid)
     edgemarkers = grid[BFaceRegions]
 
     PlutoVista.backend!(ctx[:figure]; backend = ctx[:backend], datadim = 2)
-    PlutoVista.trimesh!(ctx[:figure], pts, tris;
-                        zoom = ctx[:zoom],
-                        markers = markers, colormap = cmap, gridscale = ctx[:gridscale],
-                        edges = edges, edgemarkers = edgemarkers, edgecolormap = bcmap,
-                        aspect = ctx[:aspect],
-                        xlabel = ctx[:xlabel],
-                        ylabel = ctx[:ylabel],
-                        zlabel = ctx[:zlabel])
+    PlutoVista.trimesh!(
+        ctx[:figure],
+        pts,
+        tris;
+        zoom = ctx[:zoom],
+        markers = markers,
+        colormap = cmap,
+        gridscale = ctx[:gridscale],
+        edges = edges,
+        edgemarkers = edgemarkers,
+        edgecolormap = bcmap,
+        aspect = ctx[:aspect],
+        xlabel = ctx[:xlabel],
+        ylabel = ctx[:ylabel],
+        zlabel = ctx[:zlabel],
+    )
     reveal(ctx, TP)
 end
 
-function scalarplot!(ctx, TP::Type{PlutoVistaType}, ::Type{Val{2}}, grid, func)
+function scalarplot!(
+    ctx,
+    TP::Type{PlutoVistaType},
+    ::Type{Val{2}},
+    grids,
+    parentgrid,
+    funcs,
+)
+
+    ngrids = length(grids)
+    coords = [grid[Coordinates] for grid in grids]
+    npoints = [num_nodes(grid) for grid in grids]
+    cellnodes = [grid[CellNodes] for grid in grids]
+    ncells = [num_cells(grid) for grid in grids]
+    offsets = zeros(Int, ngrids)
+    for i = 2:ngrids
+        offsets[i] = offsets[i-1] + npoints[i-1]
+    end
+
+    allcoords = hcat(coords...)
+
+    allcellnodes = Matrix{Int}(undef, 3, sum(ncells))
+    k = 1
+    for j = 1:ngrids
+        for i = 1:ncells[j]
+            allcellnodes[1, k] = cellnodes[j][1, i] + offsets[j]
+            allcellnodes[2, k] = cellnodes[j][2, i] + offsets[j]
+            allcellnodes[3, k] = cellnodes[j][3, i] + offsets[j]
+            k = k + 1
+        end
+    end
+
     PlutoVista = ctx[:Plotter]
     PlutoVista.backend!(ctx[:figure]; backend = ctx[:backend], datadim = 2)
-    PlutoVista.tricontour!(ctx[:figure],
-                           grid[Coordinates],
-                           grid[CellNodes],
-                           func;
-                           colormap = ctx[:colormap],
-                           levels = ctx[:levels],
-                           colorbarticks = ctx[:colorbarticks],
-                           limits = ctx[:limits],
-                           backend = ctx[:backend],
-                           zoom = ctx[:zoom],
-                           gridscale = ctx[:gridscale],
-                           aspect = ctx[:aspect],
-                           xlabel = ctx[:xlabel],
-                           ylabel = ctx[:ylabel],
-                           zlabel = ctx[:zlabel])
+    PlutoVista.tricontour!(
+        ctx[:figure],
+        allcoords,
+        allcellnodes,
+        vcat(funcs...),
+        colormap = ctx[:colormap],
+        levels = ctx[:levels],
+        colorbarticks = ctx[:colorbarticks],
+        limits = ctx[:limits],
+        backend = ctx[:backend],
+        zoom = ctx[:zoom],
+        gridscale = ctx[:gridscale],
+        aspect = ctx[:aspect],
+        xlabel = ctx[:xlabel],
+        ylabel = ctx[:ylabel],
+        zlabel = ctx[:zlabel],
+    )
     reveal(ctx, TP)
 end
 
@@ -246,41 +367,58 @@ function gridplot!(ctx, TP::Type{PlutoVistaType}, ::Type{Val{3}}, grid)
     facemarkers = grid[BFaceRegions]
 
     PlutoVista.backend!(ctx[:figure]; backend = ctx[:backend], datadim = 3)
-    PlutoVista.tetmesh!(ctx[:figure], pts, tris;
-                        xplanes = ctx[:xplanes],
-                        yplanes = ctx[:yplanes],
-                        zplanes = ctx[:zplanes],
-                        zoom = ctx[:zoom],
-                        markers = markers, colormap = cmap,
-                        faces = faces, facemarkers = facemarkers, facecolormap = bcmap,
-                        outlinealpha = ctx[:outlinealpha])
+    PlutoVista.tetmesh!(
+        ctx[:figure],
+        pts,
+        tris;
+        xplanes = ctx[:xplanes],
+        yplanes = ctx[:yplanes],
+        zplanes = ctx[:zplanes],
+        zoom = ctx[:zoom],
+        markers = markers,
+        colormap = cmap,
+        faces = faces,
+        facemarkers = facemarkers,
+        facecolormap = bcmap,
+        outlinealpha = ctx[:outlinealpha],
+    )
     reveal(ctx, TP)
 end
 
-function scalarplot!(ctx, TP::Type{PlutoVistaType}, ::Type{Val{3}}, grid, func)
+function scalarplot!(
+    ctx,
+    TP::Type{PlutoVistaType},
+    ::Type{Val{3}},
+    grids,
+    parentgrid,
+    funcs,
+)
     PlutoVista = ctx[:Plotter]
-    nbregions = num_bfaceregions(grid)
+    nbregions = num_bfaceregions(parentgrid)
     bcmap = bregion_cmap(nbregions)
     PlutoVista.backend!(ctx[:figure]; backend = ctx[:backend], datadim = 3)
-    PlutoVista.tetcontour!(ctx[:figure],
-                           grid[Coordinates],
-                           grid[CellNodes],
-                           func;
-                           levels = ctx[:levels],
-                           colormap = ctx[:colormap],
-                           xplanes = ctx[:xplanes],
-                           yplanes = ctx[:yplanes],
-                           zplanes = ctx[:zplanes],
-                           limits = ctx[:limits],
-                           faces = grid[BFaceNodes],
-                           facemarkers = grid[BFaceRegions],
-                           facecolormap = bcmap,
-                           outlinealpha = ctx[:outlinealpha],
-                           levelalpha = ctx[:levelalpha],
-                           zoom = ctx[:zoom],
-                           aspect = ctx[:aspect],
-                           xlabel = ctx[:xlabel],
-                           ylabel = ctx[:ylabel],
-                           zlabel = ctx[:zlabel])
+    PlutoVista.tetcontour!(
+        ctx[:figure],
+        [grid[Coordinates] for grid in grids],
+        [grid[CellNodes] for grid in grids],
+        funcs;
+        parentpts = parentgrid[Coordinates],
+        levels = ctx[:levels],
+        colormap = ctx[:colormap],
+        xplanes = ctx[:xplanes],
+        yplanes = ctx[:yplanes],
+        zplanes = ctx[:zplanes],
+        limits = ctx[:limits],
+        faces = parentgrid[BFaceNodes],
+        facemarkers = parentgrid[BFaceRegions],
+        facecolormap = bcmap,
+        outlinealpha = ctx[:outlinealpha],
+        levelalpha = ctx[:levelalpha],
+        zoom = ctx[:zoom],
+        aspect = ctx[:aspect],
+        xlabel = ctx[:xlabel],
+        ylabel = ctx[:ylabel],
+        zlabel = ctx[:zlabel],
+    )
     reveal(ctx, TP)
 end
diff --git a/src/pyplot.jl b/src/pyplot.jl
index 6d182eb..384f6a9 100644
--- a/src/pyplot.jl
+++ b/src/pyplot.jl
@@ -42,49 +42,83 @@ function reveal(ctx::SubVisualizer, TP::Type{PyPlotType})
 end
 
 #translate Julia attribute symbols to pyplot-speak
-const mshapes = Dict(:dtriangle => "v",
-                     :utriangle => "^",
-                     :rtriangle => ">",
-                     :ltriangle => "^",
-                     :circle => "o",
-                     :square => "s",
-                     :cross => "+",
-                     :+ => "+",
-                     :xcross => "x",
-                     :x => "x",
-                     :diamond => "D",
-                     :star5 => "*",
-                     :pentagon => "p",
-                     :hexagon => "h")
-
-const lstyles = Dict(:solid => "-",
-                     :dot => "dotted",
-                     :dash => "--",
-                     :dashdot => "-.",
-                     :dashdotdot => (0, (3, 1, 1, 1)))
-
-const leglocs = Dict(:none => "",
-                     :best => "best",
-                     :lt => "upper left",
-                     :ct => "upper center",
-                     :rt => "upper right",
-                     :lc => "center left",
-                     :cc => "center center",
-                     :rc => "center right",
-                     :lb => "lower left",
-                     :cb => "lower center",
-                     :rb => "lower right")
+const mshapes = Dict(
+    :dtriangle => "v",
+    :utriangle => "^",
+    :rtriangle => ">",
+    :ltriangle => "^",
+    :circle => "o",
+    :square => "s",
+    :cross => "+",
+    :+ => "+",
+    :xcross => "x",
+    :x => "x",
+    :diamond => "D",
+    :star5 => "*",
+    :pentagon => "p",
+    :hexagon => "h",
+)
+
+const lstyles = Dict(
+    :solid => "-",
+    :dot => "dotted",
+    :dash => "--",
+    :dashdot => "-.",
+    :dashdotdot => (0, (3, 1, 1, 1)),
+)
+
+const leglocs = Dict(
+    :none => "",
+    :best => "best",
+    :lt => "upper left",
+    :ct => "upper center",
+    :rt => "upper right",
+    :lc => "center left",
+    :cc => "center center",
+    :rc => "center right",
+    :lb => "lower left",
+    :cb => "lower center",
+    :rb => "lower right",
+)
 
 """
 $(SIGNATURES)
 Return tridata to be splatted to PyPlot calls
 """
-function tridata(grid)
+function tridata(grid::ExtendableGrid)
     coord = grid[Coordinates]
     cellnodes = Matrix(grid[CellNodes])
     coord[1, :], coord[2, :], transpose(cellnodes .- 1)
 end
 
+function tridata(grids)
+    ngrids = length(grids)
+    coords = [grid[Coordinates] for grid in grids]
+    npoints = [num_nodes(grid) for grid in grids]
+    cellnodes = [grid[CellNodes] for grid in grids]
+    ncells = [num_cells(grid) for grid in grids]
+    offsets = zeros(Int, ngrids)
+    for i = 2:ngrids
+        offsets[i] = offsets[i-1] + npoints[i-1]
+    end
+
+    allcoords = hcat(coords...)
+
+    # transpose and subtract 1 !
+    allcellnodes = Matrix{Int}(undef, sum(ncells), 3)
+    k = 1
+    for j = 1:ngrids
+        for i = 1:ncells[j]
+            allcellnodes[k, 1] = cellnodes[j][1, i] + offsets[j] - 1
+            allcellnodes[k, 2] = cellnodes[j][2, i] + offsets[j] - 1
+            allcellnodes[k, 3] = cellnodes[j][3, i] + offsets[j] - 1
+            k = k + 1
+        end
+    end
+    allcoords[1, :], allcoords[2, :], allcellnodes
+end
+
+
 # Interfaces to Colors/Colorschemes
 plaincolormap(ctx) = colorschemes[ctx[:colormap]].colors
 
@@ -129,8 +163,13 @@ function gridplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{1}}, grid)
 
         x1 = coord[1, cellnodes[1, icell]]
         x2 = coord[1, cellnodes[2, icell]]
-        ax.plot([x1, x2], [0, 0]; linewidth = 3.0,
-                color = rgbtuple(cmap[cellregions[icell]]), label = label)
+        ax.plot(
+            [x1, x2],
+            [0, 0];
+            linewidth = 3.0,
+            color = rgbtuple(cmap[cellregions[icell]]),
+            label = label,
+        )
         ax.plot([x1, x1], [-h, h]; linewidth = ctx[:linewidth], color = "k", label = "")
         ax.plot([x2, x2], [-h, h]; linewidth = ctx[:linewidth], color = "k", label = "")
     end
@@ -142,8 +181,13 @@ function gridplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{1}}, grid)
             label = brflag[ireg] ? "b$(ireg)" : ""
             brflag[ireg] = false
             x1 = coord[1, bfacenodes[1, ibface]]
-            ax.plot([x1, x1], [-2 * h, 2 * h]; linewidth = 3.0,
-                    color = rgbtuple(cmap[ireg]), label = label)
+            ax.plot(
+                [x1, x1],
+                [-2 * h, 2 * h];
+                linewidth = 3.0,
+                color = rgbtuple(cmap[ireg]),
+                label = label,
+            )
         end
     end
     if ctx[:legend] != :none
@@ -189,19 +233,29 @@ function gridplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{2}}, grid)
     ax.set_aspect(ctx[:aspect])
     tridat = tridata(grid)
     cmap = region_cmap(ncellregions)
-    cdata = ax.tripcolor(tridat...;
-                         facecolors = grid[CellRegions],
-                         cmap = PyPlot.ColorMap(cmap, length(cmap)),
-                         vmin = 1.0,
-                         vmax = length(cmap))
+    cdata = ax.tripcolor(
+        tridat...;
+        facecolors = grid[CellRegions],
+        cmap = PyPlot.ColorMap(cmap, length(cmap)),
+        vmin = 1.0,
+        vmax = length(cmap),
+    )
     if ctx[:colorbar] == :horizontal
-        cbar = fig.colorbar(cdata; ax = ax, ticks = collect(1:length(cmap)),
-                            orientation = "horizontal")
+        cbar = fig.colorbar(
+            cdata;
+            ax = ax,
+            ticks = collect(1:length(cmap)),
+            orientation = "horizontal",
+        )
     end
 
     if ctx[:colorbar] == :vertical
-        cbar = fig.colorbar(cdata; ax = ax, ticks = collect(1:length(cmap)),
-                            orientation = "vertical")
+        cbar = fig.colorbar(
+            cdata;
+            ax = ax,
+            ticks = collect(1:length(cmap)),
+            orientation = "vertical",
+        )
     end
 
     ax.triplot(tridat...; color = "k", linewidth = ctx[:linewidth])
@@ -212,9 +266,13 @@ function gridplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{2}}, grid)
         c1 = [coord[:, bfacenodes[1, i]] for i = 1:num_sources(bfacenodes)]
         c2 = [coord[:, bfacenodes[2, i]] for i = 1:num_sources(bfacenodes)]
         rgb = [rgbtuple(cmap[bfaceregions[i]]) for i = 1:length(bfaceregions)]
-        ax.add_collection(PyPlot.matplotlib.collections.LineCollection(collect(zip(c1, c2));
-                                                                       colors = rgb,
-                                                                       linewidth = 3))
+        ax.add_collection(
+            PyPlot.matplotlib.collections.LineCollection(
+                collect(zip(c1, c2));
+                colors = rgb,
+                linewidth = 3,
+            ),
+        )
         for i = 1:nbfaceregions
             ax.plot(coord[1, 1:1], coord[2, 1:1]; label = "$(i)", color = rgbtuple(cmap[i]))
         end
@@ -264,37 +322,53 @@ function gridplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{3}}, grid)
     xyzcut = [ctx[:xplanes][1], ctx[:yplanes][1], ctx[:zplanes][1]]
 
     if ctx[:interior]
-        regpoints0, regfacets0 = extract_visible_cells3D(grid,
-                                                         xyzcut;
-                                                         primepoints = hcat(xyzmin, xyzmax))
-        regfacets = [reshape(reinterpret(Int32, regfacets0[i]), (3, length(regfacets0[i])))
-                     for i = 1:nregions]
-        regpoints = [reshape(reinterpret(Float32, regpoints0[i]),
-                             (3, length(regpoints0[i]))) for i = 1:nregions]
+        regpoints0, regfacets0 =
+            extract_visible_cells3D(grid, xyzcut; primepoints = hcat(xyzmin, xyzmax))
+        regfacets = [
+            reshape(reinterpret(Int32, regfacets0[i]), (3, length(regfacets0[i]))) for
+            i = 1:nregions
+        ]
+        regpoints = [
+            reshape(reinterpret(Float32, regpoints0[i]), (3, length(regpoints0[i]))) for
+            i = 1:nregions
+        ]
 
         for ireg = 1:nregions
             if size(regfacets[ireg], 2) > 0
-                ax.plot_trisurf(regpoints[ireg][1, :], regpoints[ireg][2, :],
-                                transpose(regfacets[ireg] .- 1), regpoints[ireg][3, :];
-                                color = rgbtuple(cmap[ireg]), edgecolors = :black,
-                                linewidth = 0.5)
+                ax.plot_trisurf(
+                    regpoints[ireg][1, :],
+                    regpoints[ireg][2, :],
+                    transpose(regfacets[ireg] .- 1),
+                    regpoints[ireg][3, :];
+                    color = rgbtuple(cmap[ireg]),
+                    edgecolors = :black,
+                    linewidth = 0.5,
+                )
             end
         end
     end
 
-    bregpoints0, bregfacets0 = extract_visible_bfaces3D(grid,
-                                                        xyzcut;
-                                                        primepoints = hcat(xyzmin, xyzmax))
-    bregfacets = [reshape(reinterpret(Int32, bregfacets0[i]), (3, length(bregfacets0[i])))
-                  for i = 1:nbregions]
-    bregpoints = [reshape(reinterpret(Float32, bregpoints0[i]),
-                          (3, length(bregpoints0[i]))) for i = 1:nbregions]
+    bregpoints0, bregfacets0 =
+        extract_visible_bfaces3D(grid, xyzcut; primepoints = hcat(xyzmin, xyzmax))
+    bregfacets = [
+        reshape(reinterpret(Int32, bregfacets0[i]), (3, length(bregfacets0[i]))) for
+        i = 1:nbregions
+    ]
+    bregpoints = [
+        reshape(reinterpret(Float32, bregpoints0[i]), (3, length(bregpoints0[i]))) for
+        i = 1:nbregions
+    ]
     for ireg = 1:nbregions
         if size(bregfacets[ireg], 2) > 0
-            ax.plot_trisurf(bregpoints[ireg][1, :], bregpoints[ireg][2, :],
-                            transpose(bregfacets[ireg] .- 1), bregpoints[ireg][3, :];
-                            color = rgbtuple(bcmap[ireg]), edgecolors = :black,
-                            linewidth = 0.5)
+            ax.plot_trisurf(
+                bregpoints[ireg][1, :],
+                bregpoints[ireg][2, :],
+                transpose(bregfacets[ireg] .- 1),
+                bregpoints[ireg][3, :];
+                color = rgbtuple(bcmap[ireg]),
+                edgecolors = :black,
+                linewidth = 0.5,
+            )
         end
     end
 
@@ -305,8 +379,9 @@ function gridplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{3}}, grid)
 end
 
 ### 1D Function
-function scalarplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{1}}, grid, func::AbstractVector)
+function scalarplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{1}}, grids, parentgrid, funcs)
     PyPlot = ctx[:Plotter]
+    nfuncs = length(funcs)
     if !haskey(ctx, :ax)
         ctx[:ax] = ctx[:figure].add_subplot(ctx[:layout]..., ctx[:iplot])
     end
@@ -328,9 +403,6 @@ function scalarplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{1}}, grid, func::Abst
     ax = ctx[:ax]
     fig = ctx[:figure]
 
-    cellnodes = grid[CellNodes]
-    coord = grid[Coordinates]
-
     pplot = ax.plot
     if ctx[:xscale] == :log
         if ctx[:yscale] == :log
@@ -347,51 +419,75 @@ function scalarplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{1}}, grid, func::Abst
         end
     end
 
-    if ctx[:cellwise] # not checked
+    if ctx[:cellwise] # not checked,  outdated
         for icell = 1:num_cells(grid)
             i1 = cellnodes[1, icell]
             i2 = cellnodes[2, icell]
             x1 = coord[1, i1]
             x2 = coord[1, i2]
             if icell == 1
-                ax.plot([x1, x2], [func[i1], func[i2]]; color = rgbtuple(ctx[:color]),
-                        label = ctx[:label])
+                ax.plot(
+                    [x1, x2],
+                    [func[i1], func[i2]];
+                    color = rgbtuple(ctx[:color]),
+                    label = ctx[:label],
+                )
             else
                 ax.plot([x1, x2], [func[i1], func[i2]]; color = rgbtuple(ctx[:color]))
             end
         end
     else
         if ctx[:markershape] == :none
-            if ctx[:label] !== ""
-                pplot(coord[1, :], func;
-                      linestyle = lstyles[ctx[:linestyle]],
-                      color = rgbtuple(ctx[:color]),
-                      linewidth = ctx[:linewidth],
-                      label = ctx[:label])
-            else
-                pplot(coord[1, :], func;
-                      linestyle = lstyles[ctx[:linestyle]],
-                      linewidth = ctx[:linewidth],
-                      color = rgbtuple(ctx[:color]))
+            for ifunc = 1:nfuncs
+                func = funcs[ifunc]
+                coord = grids[ifunc][Coordinates]
+                if ctx[:label] !== "" && ifunc == 1
+                    pplot(
+                        coord[1, :],
+                        func;
+                        linestyle = lstyles[ctx[:linestyle]],
+                        color = rgbtuple(ctx[:color]),
+                        linewidth = ctx[:linewidth],
+                        label = ctx[:label],
+                    )
+                else
+                    pplot(
+                        coord[1, :],
+                        func;
+                        linestyle = lstyles[ctx[:linestyle]],
+                        linewidth = ctx[:linewidth],
+                        color = rgbtuple(ctx[:color]),
+                    )
+                end
             end
         else
-            if ctx[:label] !== ""
-                pplot(coord[1, :], func;
-                      linestyle = lstyles[ctx[:linestyle]],
-                      color = rgbtuple(ctx[:color]),
-                      label = ctx[:label],
-                      marker = mshapes[ctx[:markershape]],
-                      markevery = ctx[:markevery],
-                      markersize = ctx[:markersize],
-                      linewidth = ctx[:linewidth])
-            else
-                pplot(coord[1, :], func;
-                      linestyle = lstyles[ctx[:linestyle]],
-                      color = rgbtuple(ctx[:color]),
-                      marker = mshapes[ctx[:markershape]],
-                      markevery = ctx[:markevery],
-                      markersize = ctx[:markersize],
-                      linewidth = ctx[:linewidth])
+            for ifunc = 1:nfuncs
+                func = funcs[ifunc]
+                coord = grids[ifunc][Coordinates]
+                if ctx[:label] !== "" && ifunc == 1
+                    pplot(
+                        coord[1, :],
+                        func;
+                        linestyle = lstyles[ctx[:linestyle]],
+                        color = rgbtuple(ctx[:color]),
+                        label = ctx[:label],
+                        marker = mshapes[ctx[:markershape]],
+                        markevery = ctx[:markevery],
+                        markersize = ctx[:markersize],
+                        linewidth = ctx[:linewidth],
+                    )
+                else
+                    pplot(
+                        coord[1, :],
+                        func;
+                        linestyle = lstyles[ctx[:linestyle]],
+                        color = rgbtuple(ctx[:color]),
+                        marker = mshapes[ctx[:markershape]],
+                        markevery = ctx[:markevery],
+                        markersize = ctx[:markersize],
+                        linewidth = ctx[:linewidth],
+                    )
+                end
             end
         end
         # points=[Point2f(coord[1,i],func[i]) for i=1:length(func)]
@@ -412,7 +508,8 @@ function scalarplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{1}}, grid, func::Abst
 end
 
 ### 2D Function
-function scalarplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{2}}, grid, func)
+function scalarplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{2}}, grids, parentgrid, funcs)
+
     PyPlot = ctx[:Plotter]
     if !haskey(ctx, :ax)
         ctx[:ax] = ctx[:figure].add_subplot(ctx[:layout]..., ctx[:iplot])
@@ -439,33 +536,53 @@ function scalarplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{2}}, grid, func)
     ax.set_aspect(ctx[:aspect])
     ax.set_title(ctx[:title])
 
-    levels, crange, colorbarticks = isolevels(ctx, func)
+    levels, crange, colorbarticks = isolevels(ctx, funcs)
     eps = 1.0e-5
     if crange[1] == crange[2]
         crange = (crange[1] - eps, crange[1] + eps)
-        colorlevels = collect(crange[1]:((crange[2] - crange[1]) / (1)):crange[2])
+        colorlevels = collect(crange[1]:((crange[2]-crange[1])/(1)):crange[2])
     else
-        colorlevels = collect(crange[1]:((crange[2] - crange[1]) / (ctx[:colorlevels] - 1)):crange[2])
+        colorlevels =
+            collect(crange[1]:((crange[2]-crange[1])/(ctx[:colorlevels]-1)):crange[2])
     end
 
-    if !haskey(ctx, :grid) || !seemingly_equal(ctx[:grid], grid)
-        ctx[:grid] = grid
-        ctx[:tridata] = tridata(grid)
-    end
-    cnt = ax.tricontourf(ctx[:tridata]..., func; levels = colorlevels,
-                         cmap = PyPlot.ColorMap(plaincolormap(ctx)))
+    #    if !haskey(ctx, :grid) || !seemingly_equal(ctx[:grid], grid)
+    #        ctx[:grid] = grids
+    #        ctx[:tridata] = tridata(grids)
+    #    end
+
+    tdat = tridata(grids)
+    func = vcat(funcs...)
+    cnt = ax.tricontourf(
+        tdat...,
+        func;
+        levels = colorlevels,
+        cmap = PyPlot.ColorMap(plaincolormap(ctx)),
+    )
+
     for c in cnt.collections
         c.set_edgecolor("face")
     end
-    ax.tricontour(ctx[:tridata]..., func; colors = "k", levels = levels)
+
+    ax.tricontour(tdat..., func; colors = "k", levels = levels)
 
     if ctx[:colorbar] == :horizontal
-        ctx[:cbar] = fig.colorbar(cnt; ax = ax, ticks = colorbarticks,
-                                  boundaries = colorlevels, orientation = "horizontal")
+        ctx[:cbar] = fig.colorbar(
+            cnt;
+            ax = ax,
+            ticks = colorbarticks,
+            boundaries = colorlevels,
+            orientation = "horizontal",
+        )
     end
     if ctx[:colorbar] == :vertical
-        ctx[:cbar] = fig.colorbar(cnt; ax = ax, ticks = colorbarticks,
-                                  boundaries = colorlevels, orientation = "vertical")
+        ctx[:cbar] = fig.colorbar(
+            cnt;
+            ax = ax,
+            ticks = colorbarticks,
+            boundaries = colorlevels,
+            orientation = "vertical",
+        )
     end
 
     ax.set_xlabel(ctx[:xlabel])
@@ -475,7 +592,7 @@ function scalarplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{2}}, grid, func)
     reveal(ctx, TP)
 end
 
-function scalarplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{3}}, grid, func)
+function scalarplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{3}}, grids, parentgrid, funcs)
     PyPlot = ctx[:Plotter]
     if !haskey(ctx, :ax)
         ctx[:ax] = ctx[:figure].add_subplot(ctx[:layout]..., ctx[:iplot]; projection = "3d")
@@ -483,29 +600,28 @@ function scalarplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{3}}, grid, func)
     ax = ctx[:ax]
     fig = ctx[:figure]
 
-    nregions = num_cellregions(grid)
-    nbregions = num_bfaceregions(grid)
-
     xyzmin = zeros(3)
     xyzmax = ones(3)
-    coord = grid[Coordinates]
+    coord = parentgrid[Coordinates]
     @views for idim = 1:3
         xyzmin[idim] = minimum(coord[idim, :])
         xyzmax[idim] = maximum(coord[idim, :])
     end
     xyzcut = [ctx[:xplanes], ctx[:yplanes], ctx[:zplanes]]
-    levels, crange, colorbarticks = isolevels(ctx, func)
+    levels, crange, colorbarticks = isolevels(ctx, funcs)
     eps = 1.0e-5
     if crange[1] == crange[2]
         crange = (crange[1] - eps, crange[1] + eps)
-        colorlevels = collect(crange[1]:((crange[2] - crange[1]) / (1)):crange[2])
+        colorlevels = collect(crange[1]:((crange[2]-crange[1])/(1)):crange[2])
     else
-        colorlevels = collect(crange[1]:((crange[2] - crange[1]) / (ctx[:colorlevels] - 1)):crange[2])
+        colorlevels =
+            collect(crange[1]:((crange[2]-crange[1])/(ctx[:colorlevels]-1)):crange[2])
     end
 
     planes = makeplanes(xyzmin, xyzmax, ctx[:xplanes], ctx[:yplanes], ctx[:zplanes])
 
-    ccoord0, faces0, values = marching_tetrahedra(grid, func, planes, levels, tol = ctx[:tetxplane_tol])
+    ccoord0, faces0, values =
+        marching_tetrahedra(grids, funcs, planes, levels, tol = ctx[:tetxplane_tol])
 
     faces = reshape(reinterpret(Int32, faces0), (3, length(faces0)))
     ccoord = reshape(reinterpret(Float32, ccoord0), (3, length(ccoord0)))
@@ -514,15 +630,19 @@ function scalarplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{3}}, grid, func)
     if nfaces > 0
         colors = zeros(nfaces)
         for i = 1:nfaces
-            colors[i] = (values[faces[1, i]] + values[faces[2, i]] + values[faces[3, i]]) /
-                        3
+            colors[i] =
+                (values[faces[1, i]] + values[faces[2, i]] + values[faces[3, i]]) / 3
         end
         # thx, https://stackoverflow.com/a/24229480/8922290 
-        collec = ctx[:ax].plot_trisurf(ccoord[1, :], ccoord[2, :], transpose(faces .- 1),
-                                       ccoord[3, :];
-                                       cmap = PyPlot.ColorMap(plaincolormap(ctx)),
-                                       vmin = crange[1],
-                                       vmax = crange[2])
+        collec = ctx[:ax].plot_trisurf(
+            ccoord[1, :],
+            ccoord[2, :],
+            transpose(faces .- 1),
+            ccoord[3, :];
+            cmap = PyPlot.ColorMap(plaincolormap(ctx)),
+            vmin = crange[1],
+            vmax = crange[2],
+        )
         collec.set_array(colors)
         collec.autoscale()
     end
@@ -533,12 +653,22 @@ function scalarplot!(ctx, TP::Type{PyPlotType}, ::Type{Val{3}}, grid, func)
     ax.view_init(ctx[:elev], ctx[:azim])
 
     if ctx[:colorbar] == :horizontal
-        ctx[:cbar] = fig.colorbar(collec; ax = ax, ticks = colorbarticks,
-                                  boundaries = colorlevels, orientation = "horizontal")
+        ctx[:cbar] = fig.colorbar(
+            collec;
+            ax = ax,
+            ticks = colorbarticks,
+            boundaries = colorlevels,
+            orientation = "horizontal",
+        )
     end
     if ctx[:colorbar] == :vertical
-        ctx[:cbar] = fig.colorbar(collec; ax = ax, ticks = colorbarticks,
-                                  boundaries = colorlevels, orientation = "vertical")
+        ctx[:cbar] = fig.colorbar(
+            collec;
+            ax = ax,
+            ticks = colorbarticks,
+            boundaries = colorlevels,
+            orientation = "vertical",
+        )
     end
     ax.set_title(ctx[:title])
 
diff --git a/src/vtkview.jl b/src/vtkview.jl
index ed743d2..5f63f30 100644
--- a/src/vtkview.jl
+++ b/src/vtkview.jl
@@ -68,7 +68,9 @@ end
 gridplot!(ctx, T::Type{VTKViewType}, ::Type{Val{2}}, grid) = gridplot!(ctx, T, grid)
 gridplot!(ctx, T::Type{VTKViewType}, ::Type{Val{3}}, grid) = gridplot!(ctx, T, grid)
 
-function scalarplot!(ctx, TP::Type{VTKViewType}, grid, func)
+function scalarplot!(ctx, TP::Type{VTKViewType}, grids, parentgrid, funcs)
+    grid = parentgrid
+    func = funcs[1]
     VTKView = ctx[:Plotter]
     frame = ctx[:frame]
     if !haskey(ctx, :dataset)
@@ -94,15 +96,17 @@ function scalarplot!(ctx, TP::Type{VTKViewType}, grid, func)
     reveal(ctx, TP)
 end
 
-function scalarplot!(ctx, T::Type{VTKViewType}, ::Type{Val{2}}, grid, func)
-    scalarplot!(ctx, T, grid, func)
+function scalarplot!(ctx, T::Type{VTKViewType}, ::Type{Val{2}}, grids, parentgrid, funcs)
+    scalarplot!(ctx, T, grids, parentgrid, funcs)
 end
-function scalarplot!(ctx, T::Type{VTKViewType}, ::Type{Val{3}}, grid, func)
-    scalarplot!(ctx, T, grid, func)
+function scalarplot!(ctx, T::Type{VTKViewType}, ::Type{Val{3}}, grids, parentgrid, funcs)
+    scalarplot!(ctx, T, grids, parentgrid, funcs)
 end
 gridplot!(ctx, T::Type{VTKViewType}, ::Type{Val{1}}, grid) = nothing
 
-function scalarplot!(ctx, TP::Type{VTKViewType}, ::Type{Val{1}}, grid, func)
+function scalarplot!(ctx, TP::Type{VTKViewType}, ::Type{Val{1}}, grids, parentgrid, funcs)
+    grid = parentgrid
+    func = funcs[1]
     VTKView = ctx[:Plotter]
     frame = ctx[:frame]
     if !haskey(ctx, :plot)
@@ -127,3 +131,13 @@ function scalarplot!(ctx, TP::Type{VTKViewType}, ::Type{Val{1}}, grid, func)
     VTKView.addplot!(plot, collect(grid[Coordinates][1, :]), collect(func))
     reveal(ctx, TP)
 end
+
+
+vectorplot!(ctx, TP::Type{VTKViewType}, ::Type{Val{1}}, grid, func) = nothing
+vectorplot!(ctx, TP::Type{VTKViewType}, ::Type{Val{2}}, grid, func) = nothing
+vectorplot!(ctx, TP::Type{VTKViewType}, ::Type{Val{3}}, grid, func) = nothing
+
+
+streamplot!(ctx, TP::Type{VTKViewType}, ::Type{Val{1}}, grid, func) = nothing
+streamplot!(ctx, TP::Type{VTKViewType}, ::Type{Val{2}}, grid, func) = nothing
+streamplot!(ctx, TP::Type{VTKViewType}, ::Type{Val{3}}, grid, func) = nothing