diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7c3c8d8f..efc5d6fe 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,6 +1,9 @@
 # CHANGES
 
-## v1.6.0
+## v1.6.0 April 28, 2026
+  - unicode plotting extension marked as deprecated (will be removed in next major version), they can be realized now with GridVisualize 1.9 and backend Plotter = UnicdePlots
+  - added scalarplot, vectorplot and broken_scalarplot defined for FEVectorBlocks (moved from ExtendableFEM)
+  - updated plotting code in all examples
   - added Curl2D and Curl3D evaluation for Hdiv
   
 
diff --git a/Project.toml b/Project.toml
index aea5af29..a09366c5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -10,6 +10,7 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 ExtendableGrids = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
 ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+GridVisualize = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
@@ -25,12 +26,16 @@ UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
 ExtendableFEMBaseUnicodePlotsExt = ["UnicodePlots"]
 
 [compat]
+Aqua = "0.8"
 DiffResults = "1"
 DifferentiationInterface = "0.7.10"
 DocStringExtensions = "0.8,0.9"
+ExampleJuggler = "2.2.1"
+ExplicitImports = "1"
 ExtendableGrids = "1.13.0"
 ExtendableSparse = "1.5.1, 2"
 ForwardDiff = "0.10.35, 1"
+GridVisualize = "1.20"
 LinearAlgebra = "1.9"
 Polynomials = "2.0.21, 3, 4"
 Printf = "1.9"
@@ -38,5 +43,19 @@ SparseArrays = "1.9"
 SparseConnectivityTracer = "1.1.2"
 SparseMatrixColorings = "0.4.22"
 SpecialPolynomials = "0.4.9, 0.5"
+Term = "2.0.8"
+Test = "1"
 UnicodePlots = "3.6"
 julia = "1.9"
+
+[extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+ExampleJuggler = "3bbe58f8-ed81-4c4e-a134-03e85fcf4a1a"
+ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
+GridVisualize = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8"
+Term = "22787eb5-b846-44ae-b979-8e399b8463ab"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
+
+[targets]
+test = ["Aqua", "ExampleJuggler", "ExplicitImports", "GridVisualize", "Term", "Test", "UnicodePlots"]
diff --git a/docs/Project.toml b/docs/Project.toml
index 7a359d3f..6b6070d1 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -12,6 +12,7 @@ Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 PlutoStaticHTML = "359b1769-a58e-495b-9770-312e911026ad"
 PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
 UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
+Term = "22787eb5-b846-44ae-b979-8e399b8463ab"
 
 [compat]
 Literate = ">=0.2.7"
diff --git a/examples/Example200_LowLevelPoisson.jl b/examples/Example200_LowLevelPoisson.jl
index fb72582b..7d568d0e 100644
--- a/examples/Example200_LowLevelPoisson.jl
+++ b/examples/Example200_LowLevelPoisson.jl
@@ -27,13 +27,13 @@ using ExtendableFEMBase
 using ExtendableGrids
 using ExtendableSparse
 using GridVisualize
-using UnicodePlots
+using UnicodePlots, Term
 using Test #
 
 function main(;
         maxnref = 8,
         order = 2,
-        Plotter = nothing,
+        Plotter = UnicodePlots,
         mu = 1.0,
         rhs = x -> x[1] - x[2]
     )
@@ -47,7 +47,7 @@ function main(;
     solution, time_assembly, time_solve = solve_poisson_lowlevel(fe_space, mu, rhs)
 
     ## loop over uniform refinements + timings
-    plt = GridVisualizer(; Plotter = Plotter, layout = (1, 1), clear = true, resolution = (500, 500))
+    plt = GridVisualizer(; Plotter = Plotter, layout = (1, 3), clear = true, resolution = (1200, 400))
     loop_allocations = 0
     for level in 1:maxnref
         X = LinRange(0, 1, 2^level + 1)
@@ -55,22 +55,19 @@ function main(;
         time_facenodes = @elapsed xgrid[FaceNodes]
         fe_space = FESpace{FEType}(xgrid)
         println("\nLEVEL = $level, ndofs = $(fe_space.ndofs)\n")
-        if level < 4
-            println(stdout, unicode_gridplot(xgrid))
-        end
+
         time_dofmap = @elapsed fe_space[CellDofs]
         solution, time_assembly, time_solve = solve_poisson_lowlevel(fe_space, mu, rhs)
 
         ## plot statistics
         println(stdout, barplot(["Grid", "FaceNodes", "celldofs", "Assembly", "Solve"], [time_grid, time_facenodes, time_dofmap, time_assembly, time_solve], title = "Runtimes"))
 
-        ## plot
-        if Plotter !== nothing
-            scalarplot!(plt[1, 1], xgrid, view(solution.entries, 1:num_nodes(xgrid)), limits = (-0.0125, 0.0125))
-        else
-            solution_grad = continuify(solution[1], Gradient)
-            println(stdout, unicode_scalarplot(solution[1]))
-        end
+        # plot
+        scalarplot!(plt[1, 1], solution[1])
+        scalarplot!(plt[1, 2], solution[1], Gradient; abs = true, clear = true)
+        vectorplot!(plt[1, 2], solution[1], Gradient; clear = false)
+        gridplot!(plt[1, 3], xgrid; markersize = 0)
+        reveal(plt)
     end
 
     return solution, plt
diff --git a/examples/Example205_LowLevelSpaceTimePoisson.jl b/examples/Example205_LowLevelSpaceTimePoisson.jl
index 0f1cbc54..9d350f74 100644
--- a/examples/Example205_LowLevelSpaceTimePoisson.jl
+++ b/examples/Example205_LowLevelSpaceTimePoisson.jl
@@ -27,14 +27,14 @@ using ExtendableFEMBase
 using ExtendableGrids
 using ExtendableSparse
 using GridVisualize
-using UnicodePlots
+using UnicodePlots, Term
 using Test #
 
 ## data for Poisson problem
 const μ = (t) -> 1.0e-1 * t + 1 * max(0, (1 - 2 * t))
 const f = (x, t) -> sin(3 * pi * x[1]) * 4 * t - cos(3 * pi * x[2]) * 4 * (1 - t)
 
-function main(; dt = 0.01, Tfinal = 1, level = 5, order = 1, Plotter = nothing, produce_movie = false)
+function main(; dt = 0.01, Tfinal = 1, level = 5, order = 1, Plotter = UnicodePlots, produce_movie = false)
 
     ## Finite element type
     FEType_time = H1Pk{1, 1, order}
@@ -72,13 +72,14 @@ function main(; dt = 0.01, Tfinal = 1, level = 5, order = 1, Plotter = nothing,
     else
         @info "Plotting at five times..."
         plot_timesteps = [2, round(Int, length(T) / 4 + 0.25), round(Int, length(T) / 2 + 0.5), round(Int, length(T) - length(T) / 4), FES_time.ndofs]
-        plt = GridVisualizer(; Plotter = Plotter, layout = (1, length(plot_timesteps)), clear = true, resolution = (200 * length(plot_timesteps), 200))
+        plt = GridVisualizer(; Plotter = Plotter, layout = (1, length(plot_timesteps)), clear = true, resolution = (300 * length(plot_timesteps), 300))
         for tj in 1:length(plot_timesteps)
             t = plot_timesteps[tj]
             first = (t - 1) * FES_space.ndofs + 1
             last = t * FES_space.ndofs
             scalarplot!(plt[1, tj], grid_space, view(sol, first:last), title = "t = $(T[t])")
         end
+        reveal(plt)
         return sol, plt
     end
 
@@ -116,7 +117,7 @@ function solve_poisson_lowlevel(FES_time, FES_space, μ, f)
         ExtendableSparse.flush!(A)
     end
 
-    @info ".... spy plot of system matrix:\n$(UnicodePlots.spy(sparse(A.cscmatrix)))"
+    @info ".... spy plot of system matrix:\n$(UnicodePlots.spy(A.cscmatrix))"
 
     ## solve
     println("Solving linear system...")
diff --git a/examples/Example210_LowLevelNavierStokes.jl b/examples/Example210_LowLevelNavierStokes.jl
index 06b9c02c..160b778c 100644
--- a/examples/Example210_LowLevelNavierStokes.jl
+++ b/examples/Example210_LowLevelNavierStokes.jl
@@ -35,6 +35,7 @@ using ExtendableFEMBase
 using ExtendableGrids
 using ExtendableSparse
 using GridVisualize
+using UnicodePlots, Term
 using ForwardDiff
 using DiffResults
 
@@ -63,7 +64,7 @@ function p!(pval, qpinfo)
     return nothing
 end
 
-function main(; nref = 5, teval = 0, order = 2, Plotter = nothing)
+function main(; nref = 5, teval = 0, order = 2, Plotter = UnicodePlots)
 
     @assert order >= 2
 
@@ -98,9 +99,11 @@ function main(; nref = 5, teval = 0, order = 2, Plotter = nothing)
     println("l2 error pressure = $(error_p)")
 
     ## plot
-    plt = GridVisualizer(; Plotter = Plotter, layout = (1, 1), clear = true, resolution = (500, 500))
-    scalarplot!(plt[1, 1], xgrid, nodevalues(sol[1]; abs = true)[1, :]; title = "|u| + quiver", Plotter = Plotter)
-    vectorplot!(plt[1, 1], xgrid, eval_func_bary(PointEvaluator([(1, Identity)], sol)), clear = false)
+    plt = GridVisualizer(; Plotter = Plotter, layout = (1, 2), clear = true, resolution = (1200, 600))
+    scalarplot!(plt[1, 1], sol[1]; title = "|u| + quiver", abs = true)
+    vectorplot!(plt[1, 1], sol[1]; clear = false)
+    scalarplot!(plt[1, 2], sol[2]; title = "p")
+    reveal(plt)
 
     return sol, plt
 end
diff --git a/examples/Example280_BasisPlotter.jl b/examples/Example280_BasisPlotter.jl
index 73d0a9c7..22d00333 100644
--- a/examples/Example280_BasisPlotter.jl
+++ b/examples/Example280_BasisPlotter.jl
@@ -6,7 +6,7 @@
 This example plots all the basis functions of a H1 finite element on Edge1D or Triangle2D
 as unicode plots. This is the result with the default parameters (dim = 1, order = 3):
 
-![](https://github.com/chmerdon/ExtendableFEMBase.jl/blob/master/docs/src/assets/example280.png?raw=true")
+![](example280.png)
 
 =#
 
@@ -14,12 +14,15 @@ module Example280_BasisPlotter
 
 using ExtendableFEMBase
 using ExtendableGrids
+using GridVisualize
+using UnicodePlots, Term
 
 ## everything is wrapped in a main function
-function main(; dim = 1, order = 3)
+function main(; dim = 1, order = 3, Plotter = UnicodePlots)
 
     ## generate two grids
     @assert dim in [1, 2] "dim must be 1 or 2"
+    @assert order in 1:4 "order must be between 1 and 4"
     refgeom = dim == 1 ? Edge1D : Triangle2D
     xgrid = reference_domain(refgeom)
 
@@ -36,7 +39,41 @@ function main(; dim = 1, order = 3)
         FEFunc[1][j + coffsets[j]] = 1
     end
 
+    ## interpolate on finer grid
+    xgrid_plot = dim == 1 ? simplexgrid(0:0.01:1) : uniform_refine(xgrid, 4)
+    I = FEVector(FESpace{H1P1{ndofs}}(xgrid_plot))
+    lazy_interpolate!(I[1], FEFunc, [(1, Identity)])
+
     ## plot
-    return println(stdout, unicode_scalarplot(FEFunc[1]; title = "φ", ylim = (-0.5, 1), resolution = dim == 1 ? (40, 10) : (20, 15), nrows = order))
+    if dim == 1
+        layout = (1, 1) # everything is plotted into one plot
+        size = (600, 600)
+        colors = [:red, :green, :blue, :white, :yellow, :cyan, :magenta]
+    elseif dim == 2
+        l = round(Int, ceil(sqrt(ndofs)))
+        if l^2 - l >= ndofs
+            layout = (l - 1, l)
+        else
+            layout = (l, l)
+        end
+        size = (1000, 1000)
+        p = permutedims(reshape(1:prod(layout), layout))[:]
+    end
+    plt = GridVisualizer(; Plotter = Plotter, layout = layout, size = size)
+    for j in 1:ndofs
+        if dim == 1
+            ExtendableFEMBase.scalarplot!(plt[1], I[1], IdentityComponent{j}; Plotter = Plotter, clear = false, title = "dof $j", color = dim == 1 ? colors[j] : :white)
+        else
+            ExtendableFEMBase.scalarplot!(plt[p[j]], I[1], IdentityComponent{j}; Plotter = Plotter, clear = false, title = "dof $j", color = dim == 1 ? colors[j] : :white)
+        end
+    end
+    reveal(plt)
+    return FEFunc, plt
+end
+
+function generateplots(dir = pwd(); Plotter = nothing, kwargs...)
+    ~, plt = main(; Plotter = Plotter, kwargs...)
+    scene = GridVisualize.reveal(plt)
+    return GridVisualize.save(joinpath(dir, "example280.png"), scene; Plotter = Plotter)
 end
 end
diff --git a/examples/Example281_DiscontinuousPlot.jl b/examples/Example281_DiscontinuousPlot.jl
index 83bd0354..357b0154 100644
--- a/examples/Example281_DiscontinuousPlot.jl
+++ b/examples/Example281_DiscontinuousPlot.jl
@@ -16,6 +16,7 @@ module Example281_DiscontinuousPlot
 using ExtendableFEMBase
 using ExtendableGrids
 using GridVisualize
+using UnicodePlots, Term
 
 ## function to interpolate
 function u!(result, qpinfo)
@@ -30,7 +31,7 @@ function u!(result, qpinfo)
 end
 
 ## everything is wrapped in a main function
-function main(; broken = false, nrefs = 3, abs = false, Plotter = nothing)
+function main(; broken = false, nrefs = 3, abs = false, Plotter = UnicodePlots)
 
     ## generate two grids
     xgrid = grid_unitsquare(Triangle2D)
@@ -61,13 +62,14 @@ function main(; broken = false, nrefs = 3, abs = false, Plotter = nothing)
     nodevals4nodes2 = nodevalues(FEFunction[1], Identity; abs = abs, regions = [2], nodes = subnodes2)
 
     ## plot
-    if Plotter !== nothing
-        p = GridVisualizer(; Plotter = Plotter, layout = (2, 2), clear = true, resolution = (1000, 500))
-        gridplot!(p[1, 1], xgrid)
-        scalarplot!(p[1, 2], [subgrid1, subgrid2], xgrid, [view(nodevals4nodes1, :), view(nodevals4nodes2, :)], cellwise = false, levels = 11, title = "u")
-    end
-
-    return p
+    plt = GridVisualizer(; Plotter = Plotter, layout = (1, 5), clear = true, resolution = (1500, 300))
+    gridplot!(plt[1, 1], xgrid)
+    scalarplot!(plt[1, 2], subgrid1, view(nodevals4nodes1, :), cellwise = false, levels = 11, title = "u (region 1)")
+    scalarplot!(plt[1, 3], subgrid2, view(nodevals4nodes2, :), cellwise = false, levels = 11, title = "u (region 2)")
+    scalarplot!(plt[1, 4], [subgrid1, subgrid2], xgrid, [view(nodevals4nodes1, :), view(nodevals4nodes2, :)], cellwise = false, levels = 11, title = "u")
+    broken_scalarplot!(plt[1, 5], FEFunction[1])
+    reveal(plt)
+    return plt
 end
 
 function generateplots(dir = pwd(); Plotter = nothing, kwargs...)
diff --git a/examples/Example290_InterpolationBetweenMeshes.jl b/examples/Example290_InterpolationBetweenMeshes.jl
index bb3a3e48..a3a389b3 100644
--- a/examples/Example290_InterpolationBetweenMeshes.jl
+++ b/examples/Example290_InterpolationBetweenMeshes.jl
@@ -16,6 +16,7 @@ module Example290_InterpolationBetweenMeshes
 using ExtendableFEMBase
 using ExtendableGrids
 using GridVisualize
+using UnicodePlots, Term
 
 ## function to interpolate
 function u!(result, qpinfo)
@@ -25,7 +26,7 @@ function u!(result, qpinfo)
 end
 
 ## everything is wrapped in a main function
-function main(; ν = 1.0e-3, nrefs = 4, Plotter = nothing)
+function main(; ν = 1.0e-3, nrefs = 3, Plotter = UnicodePlots)
 
     ## generate two grids
     xgrid1 = uniform_refine(grid_unitsquare(Triangle2D), nrefs)
@@ -52,11 +53,13 @@ function main(; ν = 1.0e-3, nrefs = 4, Plotter = nothing)
     @time lazy_interpolate!(FEFunction2[1], FEFunction1; use_cellparents = true)
 
     ## plot
-    p = GridVisualizer(; Plotter = Plotter, layout = (1, 2), clear = true, resolution = (800, 400))
-    scalarplot!(p[1, 1], xgrid1, view(nodevalues(FEFunction1[1]), 1, :), levels = 11, title = "u_h ($FEType1, coarse grid)")
-    scalarplot!(p[1, 2], xgrid2, view(nodevalues(FEFunction2[1]), 1, :), levels = 11, title = "u_h ($FEType2, fine grid)")
-
-    return p
+    plt = GridVisualizer(; Plotter = Plotter, layout = (2, 2), clear = true, resolution = (1000, 1000))
+    scalarplot!(plt[1, 1], FEFunction1[1], levels = 11, title = "u_h ($FEType1, coarse grid)")
+    scalarplot!(plt[1, 2], FEFunction2[1], levels = 11, title = "u_h ($FEType2, fine grid)")
+    gridplot!(plt[2, 1], xgrid1, title = "coarse grid", markersize = 0)
+    gridplot!(plt[2, 2], xgrid2, title = "fine grid", markersize = 0)
+    reveal(plt)
+    return plt
 end
 
 function generateplots(dir = pwd(); Plotter = nothing, kwargs...)
diff --git a/ext/ExtendableFEMBaseUnicodePlotsExt.jl b/ext/ExtendableFEMBaseUnicodePlotsExt.jl
index a880f3cd..fb9e536c 100644
--- a/ext/ExtendableFEMBaseUnicodePlotsExt.jl
+++ b/ext/ExtendableFEMBaseUnicodePlotsExt.jl
@@ -18,6 +18,11 @@ function unicode_gridplot(
         plot_based = ON_CELLS,   # or ON_FACES/ON_EDGES
         kwargs...,
     )
+    Base.depwarn(
+        "The unicode extension of ExtendableFEMBase will be deprecated in a future release. A unicode gridplot can now be achieved
+         via GridVisualize and UnicodePlots as an backend, e.g., use gridplot(xgrid; Plotter = UnicodePlots) instead of unicode_gridplot(xgrid)",
+        :unicode_gridplot; force = true
+    )
     coords = xgrid[Coordinates]
     ex = extrema(view(coords, 1, :))
     ey = extrema(view(coords, 2, :))
@@ -85,6 +90,12 @@ function unicode_scalarplot(
         title = u.name,
         kwargs...
     )
+    Base.depwarn(
+        "The unicode extension of ExtendableFEMBase will be deprecated in a future release. A unicode scalarplot can now be achieved
+         via GridVisualize and UnicodePlots as an backend, e.g., use scalarplot(u; Plotter = UnicodePlots)
+         instead of unicode_scalarplot(u)",
+        :unicode_scalarplot; force = true
+    )
 
     xgrid = u.FES.xgrid
     coords = xgrid[Coordinates]
diff --git a/src/ExtendableFEMBase.jl b/src/ExtendableFEMBase.jl
index ed179357..6412c8f1 100644
--- a/src/ExtendableFEMBase.jl
+++ b/src/ExtendableFEMBase.jl
@@ -41,13 +41,15 @@ using ExtendableGrids: ExtendableGrids, AT_NODES, AbstractElementGeometry,
     local_celledgenodes, local_cellfacenodes, mapderiv!,
     max_num_targets_per_source, num_cells, num_edges,
     num_faces, num_nodes, num_sources, num_targets,
-    reference_domain, subgrid, unique,
+    reference_domain, subgrid, unique, simplexgrid,
     update_trafo!
 using ExtendableSparse: ExtendableSparse, ExtendableSparseMatrix, flush!,
     AbstractExtendableSparseMatrixCSC, ExtendableSparseMatrixCSC, MTExtendableSparseMatrixCSC,
     rawupdateindex!
 using DifferentiationInterface: AutoForwardDiff, AutoSparse, jacobian
 using ForwardDiff: ForwardDiff, DiffResults
+using GridVisualize: GridVisualize, GridVisualizer, reveal,
+    scalarplot, scalarplot!, vectorplot, vectorplot!, streamplot, streamplot!
 using LinearAlgebra: LinearAlgebra, convert, det, diagm, dot, eigen, ldiv!, lu,
     mul!, norm, transpose
 using Polynomials: Polynomials, Polynomial, coeffs
@@ -265,4 +267,11 @@ function unicode_scalarplot end
 export unicode_gridplot, unicode_scalarplot
 
 
+include("plots.jl")
+export scalarplot!, scalarplot
+export broken_scalarplot!, broken_scalarplot
+export vectorplot!, vectorplot
+export streamplot!, streamplot
+
+
 end # module ExtendableFEMBase.
diff --git a/src/plots.jl b/src/plots.jl
new file mode 100644
index 00000000..f9f5f413
--- /dev/null
+++ b/src/plots.jl
@@ -0,0 +1,93 @@
+"""
+    scalarplot!(vis, feVectorBlock::FEVectorBlock, operator = Identity; kwargs...)
+
+A standard scalarplot of (the operator evaluation of) a finite element vector into a existing visualizer.
+
+All kwargs of the calling method are transferred to the scalarplot.
+"""
+function GridVisualize.scalarplot!(vis::Union{Nothing, Dict{Symbol, Any}}, feVectorBlock::FEVectorBlock, operator = Identity; abs = false, component = 1, kwargs...)
+    return GridVisualize.scalarplot!(vis, feVectorBlock.FES.dofgrid, view(nodevalues(feVectorBlock, operator; abs = abs), component, :); kwargs...)
+end
+"""
+    scalarplot(feVectorBlock::FEVectorBlock, operator = Identity; kwargs...)
+
+A standard scalarplot of (the operator evaluation of) a finite element vector.
+
+All kwargs of the calling method are transferred to the scalarplot.
+"""
+function GridVisualize.scalarplot(feVectorBlock::FEVectorBlock, operator = Identity; abs = false, component = 1, kwargs...)
+    return GridVisualize.scalarplot(feVectorBlock.FES.dofgrid, view(nodevalues(feVectorBlock, operator; abs = abs), component, :); kwargs...)
+end
+
+
+"""
+    broken_scalarplot!(vis, feVectorBlock::FEVectorBlock, operator = Identity; kwargs...)
+
+A "broken" scalarplot (the operator evaluation of) a broken finite element vector into a existing visualizer.
+Instead of averaging the discontinuous values on the grid nodes, each grid cell is plotted
+independently. Thus, a discontinuous plot is generated.
+
+All kwargs of the calling method are transferred to the scalarplot.
+"""
+function broken_scalarplot!(vis, feVectorBlock::FEVectorBlock, operator = Identity; kwargs...)
+
+    dofgrid = feVectorBlock.FES.dofgrid
+    cell_nodes = dofgrid[CellNodes]
+    coords = dofgrid[Coordinates]
+
+    all_values = nodevalues(feVectorBlock, operator; cellwise = true) # cellwise evaluation of the FE
+    all_coords = @views coords[:, cell_nodes[:]]
+    all_cells = reshape(1:length(all_values), size(all_values))
+
+    return GridVisualize.scalarplot!(vis, simplexgrid(all_coords, all_cells, dofgrid[CellRegions]), view(all_values, :); kwargs...)
+end
+"""
+    broken_scalarplot(feVectorBlock::FEVectorBlock, operator = Identity; kwargs...)
+
+A "broken" scalarplot (the operator evaluation of) a broken finite element vector.
+Instead of averaging the discontinuous values on the grid nodes, each grid cell is plotted
+independently. Thus, a discontinuous plot is generated.
+
+All kwargs of the calling method are transferred to the scalarplot.
+"""
+function broken_scalarplot(feVectorBlock::FEVectorBlock, operator = Identity; kwargs...)
+    vis = GridVisualizer(; kwargs...)
+    broken_scalarplot!(vis, feVectorBlock, operator = operator; kwargs...)
+    return reveal(vis)
+end
+
+"""
+    vectorplot!(vis, feVectorBlock::FEVectorBlock, operator = Identity; kwargs...)
+
+A standard vectorplot of (the operator evaluation of) a finite element vector into a existing visualizer.
+
+All kwargs of the calling method are transferred to the vectorplot.
+"""
+function GridVisualize.vectorplot!(p, feVectorBlock::FEVectorBlock, operator = Identity; title = feVectorBlock.name, kwargs...)
+    return GridVisualize.vectorplot!(p, feVectorBlock.FES.dofgrid, eval_func_bary(PointEvaluator([(1, operator)], [feVectorBlock])); title = title, kwargs...)
+end
+"""
+    vectorplot(feVectorBlock::FEVectorBlock, operator = Identity; kwargs...)
+
+A standard vectorplot of (the operator evaluation of) a finite element vector.
+
+All kwargs of the calling method are transferred to the vectorplot.
+"""
+function GridVisualize.vectorplot(feVectorBlock::FEVectorBlock, operator = Identity; title = feVectorBlock.name, kwargs...)
+    return GridVisualize.vectorplot(feVectorBlock.FES.dofgrid, eval_func_bary(PointEvaluator([(1, operator)], [feVectorBlock])); title = title, kwargs...)
+end
+
+
+"""
+    vectorplot!(vis, feVectorBlock::FEVectorBlock, operator = Identity; kwargs...)
+
+A standard streamplot of (the operator evaluation of) a finite element vector.
+
+All kwargs of the calling method are transferred to the streamplot.
+"""
+function GridVisualize.streamplot!(p, feVectorBlock::FEVectorBlock, operator = Identity; title = feVectorBlock.name * " (streamlines)", kwargs...)
+    return streamplot!(p, feVectorBlock.FES.dofgrid, eval_func_bary(PointEvaluator([(1, operator)], [feVectorBlock])); title = title, kwargs...)
+end
+function GridVisualize.streamplot(feVectorBlock::FEVectorBlock, operator = Identity; title = feVectorBlock.name * " (streamlines)", kwargs...)
+    return streamplot(feVectorBlock.FES.dofgrid, eval_func_bary(PointEvaluator([(1, operator)], [feVectorBlock])); title = title, kwargs...)
+end
diff --git a/test/Project.toml b/test/Project.toml
deleted file mode 100644
index 7c4879be..00000000
--- a/test/Project.toml
+++ /dev/null
@@ -1,13 +0,0 @@
-[deps]
-Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
-ExampleJuggler = "3bbe58f8-ed81-4c4e-a134-03e85fcf4a1a"
-ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
-ExtendableGrids = "cfc395e8-590f-11e8-1f13-43a2532b2fa8"
-ExtendableSparse = "95c220a8-a1cf-11e9-0c77-dbfce5f500b3"
-GridVisualize = "5eed8a63-0fb0-45eb-886d-8d5a387d12b8"
-SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-UnicodePlots = "b8865327-cd53-5732-bb35-84acbb429228"
-
-[compat]
-Aqua = "0.8"