diff --git a/Project.toml b/Project.toml
index 1ace827..4e78701 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,7 +4,6 @@ authors = ["Alberto F. Martín <alberto.f.martin@anu.edu.au>", "Eric Neiva <enei
 version = "0.1.3"
 
 [deps]
-CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
@@ -13,13 +12,12 @@ Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 
 [compat]
-CairoMakie = "0.5, 0.6, 0.8"
 FileIO = "1"
-FillArrays = "0.11, 0.12, 0.13, 1"
-GLMakie = "0.3, 0.4, 0.6"
+FillArrays = "1"
+GLMakie = "0.10"
 GeometryBasics = "0.3, 0.4"
-Gridap = "0.16, 0.17, 0.18"
-Makie = "0.13, 0.15, 0.17"
+Gridap = "0.16, 0.17,0.18"
+Makie = "0.21"
 julia = "1.6"
 
 [extras]
diff --git a/README.md b/README.md
index 1821e7f..ec37416 100644
--- a/README.md
+++ b/README.md
@@ -135,7 +135,7 @@ Therefore, we may as well visualize such mesh
 model = DiscreteModelFromFile("models/model.json")
 Ω = Triangulation(model)
 ∂Ω = BoundaryTriangulation(model)
-fig = plot(Ω, shading=true)
+fig = plot(Ω)
 wireframe!(∂Ω, color=:black)
 save("images/3d_Fig1.png", fig)
 ````
@@ -146,7 +146,7 @@ save("images/3d_Fig1.png", fig)
 
 ````julia
 v(x) = sin(π*(x[1]+x[2]+x[3]))
-fig, ax, plt = plot(Ω, v, shading=true)
+fig, ax, plt = plot(Ω, v)
 Colorbar(fig[1,2], plt)
 save("images/3d_Fig2.png", fig)
 ````
@@ -159,7 +159,7 @@ we can even plot functions in certain subdomains, e.g.
 
 ````julia
 Γ = BoundaryTriangulation(model, tags=["square", "triangle", "circle"])
-fig = plot(Γ, v, colormap=:rainbow, shading=true)
+fig = plot(Γ, v, colormap=:rainbow)
 wireframe!(∂Ω, linewidth=0.5, color=:gray)
 save("images/3d_Fig3.png", fig)
 ````
@@ -178,7 +178,7 @@ t = Observable(0.0)
 u = lift(t) do t
     x->sin(π*(x[1]+x[2]+x[3]))*cos(π*t)
 end
-fig = plot(Ω, u, colormap=:rainbow, shading=true, colorrange=(-1,1))
+fig = plot(Ω, u, colormap=:rainbow, colorrange=(-1,1))
 wireframe!(∂Ω, color=:black, linewidth=0.5)
 framerate = 30
 timestamps = range(0, 2, step=1/framerate)
@@ -194,4 +194,3 @@ end
 ---
 
 *This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
-
diff --git a/src/recipes.jl b/src/recipes.jl
index 330c71c..8edf803 100644
--- a/src/recipes.jl
+++ b/src/recipes.jl
@@ -31,22 +31,24 @@ end
 mesh_theme = Makie.Theme(
         color      = :pink,
         colormap   = :bluesreds,
-        shading    = false,
+        shading    = Makie.NoShading,
         cycle      = nothing
 )
 
 # By merging the default theme for the Mesh type and some attributes whose values we want to impose from mesh_theme
 # (in this case: color, colormap, etc.), we may override the current default theme and use our settings.
 # Note: Makie.Theme() returns an Attributes type.
-function Makie.default_theme(scene,::Type{<:Makie.Mesh{<:Tuple{PlotGrid}}})
+@Makie.recipe(PlotGridMesh) do scene
     merge!(
         mesh_theme,
         Makie.default_theme(scene, Makie.Mesh)
     )
 end
 
+
 # The lift function is necessary when dealing with reactive attributes or Observables.
-function Makie.plot!(plot::Makie.Mesh{<:Tuple{PlotGrid}})
+#function Makie.plot!(plot::Makie.Mesh{<:Tuple{PlotGrid}})
+function Makie.plot!(plot::PlotGridMesh{<:Tuple{PlotGrid}})
     grid = Makie.lift(get_grid, plot[1])
     D = num_cell_dims(grid[])
     if D in (0,1,2)
@@ -66,20 +68,17 @@ function Makie.plot!(plot::Makie.Mesh{<:Tuple{PlotGrid}})
     # plot.attributes.attributes returns a Dict{Symbol, Observable} to be called by any function.
     # plot.attributes returns an Attributes type.
     if D in (2,3)
-        Makie.mesh!(plot, mesh;
-            plot.attributes.attributes...,
-            color = color,
-        )
+        valid_attributes = Makie.shared_attributes(plot,Makie.Mesh)
+        valid_attributes[:color] = color
+        Makie.mesh!(plot, valid_attributes, mesh)
     elseif D == 1
-        Makie.linesegments!(plot, mesh;
-            plot.attributes.attributes...,
-            color = color,
-        )
+        valid_attributes = Makie.shared_attributes(plot,Makie.LineSegments)
+        valid_attributes[:color] = color
+        Makie.linesegments!(plot, valid_attributes, mesh )
     elseif D == 0
-        Makie.scatter!(plot, mesh;
-            plot.attributes.attributes...,
-            color = color,
-        )
+        valid_attributes = Makie.shared_attributes(plot,Makie.Scatter)
+        valid_attributes[:color] = color
+        Makie.scatter!(plot, valid_attributes, mesh )
     else
         @unreachable
     end
@@ -99,7 +98,7 @@ function Makie.convert_arguments(::Type{<:Makie.Scatter}, pg::PlotGrid)
     (x, )
 end
 
-function Makie.convert_arguments(::Type{<:Makie.Mesh}, trian::Triangulation)
+function Makie.convert_arguments(::Type{PlotGridMesh}, trian::Triangulation)
     grid = to_grid(trian)
     (PlotGrid(grid), )
 end
@@ -110,7 +109,10 @@ function Makie.convert_arguments(t::Type{<:Union{Makie.Wireframe, Makie.Scatter}
 end
 
 # Set default plottype as mesh if argument is type Triangulation, i.e., mesh(Ω) == plot(Ω).
-Makie.plottype(::Triangulation) = Makie.Mesh
+Makie.plottype(::Triangulation) = PlotGridMesh
+Makie.args_preferred_axis(t::Triangulation)= num_point_dims(t)<=2 ? Makie.Axis : Makie.LScene
+Makie.plottype(::PlotGrid) = PlotGridMesh
+Makie.args_preferred_axis(pg::PlotGrid)= num_point_dims(pg.Grid)<=2 ? Makie.Axis : Makie.LScene
 
 @Makie.recipe(MeshField) do scene
     merge!(
@@ -131,7 +133,7 @@ function Makie.plot!(p::MeshField{<:Tuple{Triangulation, Any}})
     if p[:colorrange][] === Makie.automatic
         p[:colorrange] = Makie.lift(extrema, p[:color])
     end
-    Makie.mesh!(p, pg;
+    Makie.plot!(p, pg;
         p.attributes.attributes...
     )
 end
@@ -147,13 +149,13 @@ function Makie.plot!(p::MeshField{<:Tuple{CellField}})
     if p[:colorrange][] === Makie.automatic
         p[:colorrange] = Makie.lift(extrema, p[:color])
     end
-    Makie.mesh!(p, pg;
+    Makie.plot!(p, pg;
         p.attributes.attributes...
     )
 end
 
 Makie.plottype(::CellField) = MeshField
-
+Makie.args_preferred_axis(c::CellField)= num_point_dims(get_triangulation(c))<=2 ? Makie.Axis : Makie.LScene
 
 function Makie.point_iterator(pg::PlotGrid)
     UnstructuredGrid(pg.grid) |> to_dg_points
diff --git a/test/runtests.jl b/test/runtests.jl
index 194388e..a7d99eb 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -109,7 +109,7 @@ end
         fig
     end
     @test savefig("3d_Fig13") do
-      fig, _ , plt = plot(Ω, color=3*celldata, colormap=:heat)
+        fig, _ , plt = plot(Ω, color=3*celldata, colormap=:heat)
         Colorbar(fig[1,2], plt)
         fig
     end
@@ -119,7 +119,7 @@ end
         fig
     end
     @test savefig("3d_fig15") do 
-      fig, _ , plt = plot(uh, colormap=:Spectral, colorrange=(0,1))
+        fig, _ , plt = plot(uh, colormap=:Spectral, colorrange=(0,1))
         Colorbar(fig[1,2], plt)
         fig   
     end