diff --git a/nice_trees.png b/nice_trees.png
deleted file mode 100644
index de1d4b6..0000000
Binary files a/nice_trees.png and /dev/null differ
diff --git a/test/forest.jl b/test/forest.jl
index f6aaf67..05c10a1 100644
--- a/test/forest.jl
+++ b/test/forest.jl
@@ -207,7 +207,7 @@ And we can render the forest with the function `render` as in the binary tree
 example but passing the whole forest at once
 
 =#
-render(Scene(newforest))
+render(Mesh(newforest))
 #=
 
 If we iterate 4 more iterations we will start seeing the different individuals
@@ -215,7 +215,7 @@ diverging in size due to the differences in growth rates
 
 =#
 newforest = [simulate(tree, getInternode, 15) for tree in newforest];
-render(Scene(newforest))
+render(Mesh(newforest))
 #=
 
 ## Multithreaded simulation
@@ -233,7 +233,7 @@ newforest = deepcopy(forest)
 @threads for i in eachindex(forest)
     newforest[i] = simulate(forest[i], getInternode, 6)
 end
-render(Scene(newforest, parallel = true))
+render(Mesh(newforest, parallel = true))
 #=
 
 An alternative way to perform the simulation is to have an outer loop for each timestep and an internal loop over the different trees. Although this approach is not required for this simple model, most FSP models will probably need such a scheme as growth of each individual plant will depend on competition for resources with neighbouring plants. In this case, this approach would look as follows:
@@ -245,18 +245,18 @@ for step in 1:15
         newforest[i] = simulate(newforest[i], getInternode, 1)
     end
 end
-render(Scene(newforest, parallel = true))
+render(Mesh(newforest, parallel = true))
 #=
 
-# Customizing the scene
+# Customizing the mesh
 
-Here we are going to customize the scene of our simulation by adding a horizontal tile represting soil and
+Here we are going to customize the mesh of our simulation by adding a horizontal tile represting soil and
 tweaking the 3D representation. When we want to combine plants generated from graphs with any other
-geometric element it is best to combine all these geometries in a `GLScene` object. We can start the scene
+geometric element it is best to combine all these geometries in a `GLScene` object. We can start the mesh
 with the `newforest` generated in the above:
 
 =#
-scene = Scene(newforest);
+mesh = Mesh(newforest);
 #=
 
 We can create the soil tile directly, without having to create a graph. The simplest approach is two use
@@ -273,29 +273,29 @@ rotatey!(soil, pi/2)
 VirtualPlantLab.translate!(soil, Vec(0.0, 10.5, 0.0))
 #=
 
-We can now add the `soil` to the `scene` object with the `add!` function.
+We can now add the `soil` to the `mesh` object with the `add!` function.
 
 =#
-VirtualPlantLab.add!(scene, mesh = soil, colors = RGB(1,1,0))
+VirtualPlantLab.add!(mesh, soil, colors = RGB(1,1,0))
 #=
 
-We can now render the scene that combines the random forest of binary trees and a yellow soil. Notice that
+We can now render the mesh that combines the random forest of binary trees and a yellow soil. Notice that
 in all previous figures, a coordinate system with grids was being depicted. This is helpful for debugging
-your code but also to help setup the scene (e.g. if you are not sure how big the soil tile should be).
+your code but also to help setup the mesh (e.g. if you are not sure how big the soil tile should be).
 Howver, it may be distracting for the visualization. It turns out that we can turn that off with
 `show_axes = false`:
 
 =#
-render(scene, axes = false)
+render(mesh, axes = false)
 #=
 
-We may also want to save a screenshot of the scene. For this, we need to store the output of the `render` function.
-We can then resize the window rendering the scene, move around, zoom, etc. When we have a perspective that we like,
+We may also want to save a screenshot of the mesh. For this, we need to store the output of the `render` function.
+We can then resize the window rendering the mesh, move around, zoom, etc. When we have a perspective that we like,
 we can run the `save_scene` function on the object returned from `render`. The argument `resolution` can be adjusted in both
 `render` to increase the number of pixels in the final image. A helper function `calculate_resolution` is provided to
 compute the resolution from a physical width and height in cm and a dpi (e.g., useful for publications and posters):
 
 =#
 res = calculate_resolution(width = 16.0, height = 16.0, dpi = 1_000)
-output = render(scene, axes = false, size = res)
+output = render(mesh, axes = false, size = res)
 export_scene(scene = output, filename = "nice_trees.png")
diff --git a/test/growthforest.jl b/test/growthforest.jl
index 3812fc4..3920e26 100644
--- a/test/growthforest.jl
+++ b/test/growthforest.jl
@@ -11,7 +11,7 @@ Centre for Crop Systems Analysis - Wageningen University
 > - Growth rules, based on information stored in organs (dimensions, carbon assimilation)
 > - Update dimensions in function of assimilation
 > - Compute sink strength
-> - Merge Scenes
+> - Merge meshes
 > - Generate forest on grid and retrieve canopy-level data (e.g., LAI)
 >
 
@@ -395,8 +395,8 @@ end
 
 As in the previous example, it makes sense to visualize the forest with a soil
 tile beneath it. Unlike in the previous example, we will construct the soil tile
-using a dedicated graph and generate a `Scene` object which can later be
-merged with the rest of scene generated in daily step:
+using a dedicated graph and generate a `Mesh` object which can later be
+merged with the rest of mesh generated in daily step:
 
 =#
 Base.@kwdef struct Soil <: VirtualPlantLab.Node
@@ -408,19 +408,19 @@ function VirtualPlantLab.feed!(turtle::Turtle, s::Soil, vars)
 end
 soil_graph = RA(-90.0) + T(Vec(0.0, 10.0, 0.0)) + ## Moves into position
              Soil(length = 20.0, width = 20.0) ## Draws the soil tile
-soil = Scene(Graph(axiom = soil_graph));
+soil = Mesh(Graph(axiom = soil_graph));
 render(soil, axes = false)
 #=
 
-And the following function renders the entire scene (notice that we need to
-use `display()` to force the rendering of the scene when called within a loop
+And the following function renders the entire mesh (notice that we need to
+use `display()` to force the rendering of the mesh when called within a loop
 or a function):
 
 =#
 function render_forest(forest, soil)
-    scene = Scene(vec(forest)) ## create scene from forest
-    scene = Scene([scene, soil]) ## merges the two scenes
-    render(scene)
+    mesh = Mesh(vec(forest)) ## create mesh from forest
+    mesh = Mesh([mesh, soil]) ## merges the two scenes
+    render(mesh)
 end
 #=
 
diff --git a/test/raytracedforest.jl b/test/raytracedforest.jl
index 6e8dae5..78a1c82 100644
--- a/test/raytracedforest.jl
+++ b/test/raytracedforest.jl
@@ -206,14 +206,14 @@ end
 
 As growth is now dependent on intercepted PAR via RUE, we now need to simulate
 light interception by the trees. We will use a ray-tracing approach to do so.
-The first step is to create a scene with the trees and the light sources. As for
-rendering, the scene can be created from the `forest` object by simply calling
-`Scene(forest)` that will generate the 3D meshes and connect them to their
+The first step is to create a mesh with the trees and the light sources. As for
+rendering, the mesh can be created from the `forest` object by simply calling
+`Mesh(forest)` that will generate the 3D meshes and connect them to their
 optical properties.
 
 However, we also want to add the soil surface as this will affect the light
-distribution within the scene due to reflection from the soil surface. This is
-similar to the customized scene that we created before for rendering, but now
+distribution within the mesh due to reflection from the soil surface. This is
+similar to the customized mesh that we created before for rendering, but now
 for the light simulation.
 
 =#
@@ -225,17 +225,17 @@ function create_soil()
 end
 function create_scene(forest)
     # These are the trees
-    scene = Scene(vec(forest))
+    mesh = Mesh(vec(forest))
     # Add a soil surface
     soil = create_soil()
     soil_material = Lambertian(τ = 0.0, ρ = 0.21)
-    add!(scene, mesh = soil, materials = soil_material)
-    # Return the scene
-    return scene
+    add!(mesh, soil, materials = soil_material)
+    # Return the mesh
+    return mesh
 end
 #=
 
-Given the scene, we can create the light sources that can approximate the solar
+Given the mesh, we can create the light sources that can approximate the solar
 irradiance on a given day, location and time of the day using the functions from
 the  package (see package documentation for details). Given the latitude,
 day of year and fraction of the day (`f = 0` being sunrise and `f = 1` being sunset),
@@ -243,11 +243,11 @@ the function `clear_sky()` computes the direct and diffuse solar radiation assum
 a clear sky. These values may be converted to different wavebands and units using
 `waveband_conversion()`. Finally, the collection of light sources approximating
 the solar irradiance distribution over the sky hemisphere is constructed with the
-function `sky()` (this last step requires the 3D scene as input in order to place
+function `sky()` (this last step requires the 3D mesh as input in order to place
 the light sources adequately).
 
 =#
-function create_sky(;scene, lat = 52.0*π/180.0, DOY = 182)
+function create_sky(;mesh, lat = 52.0*π/180.0, DOY = 182)
     # Fraction of the day and day length
     fs = collect(0.1:0.1:0.9)
     dec = declination(DOY)
@@ -264,7 +264,7 @@ function create_sky(;scene, lat = 52.0*π/180.0, DOY = 182)
     Idir_PAR = f_dir.*Idir
     Idif_PAR = f_dif.*Idif
     # Create the dome of diffuse light
-    dome = sky(scene,
+    dome = sky(mesh,
                   Idir = 0.0, ## No direct solar radiation
                   Idif = sum(Idir_PAR)/10*DL, ## Daily Diffuse solar radiation
                   nrays_dif = 1_000_000, ## Total number of rays for diffuse solar radiation
@@ -274,20 +274,20 @@ function create_sky(;scene, lat = 52.0*π/180.0, DOY = 182)
                   nphi = 12) ## Number of discretization steps in the azimuth angle
     # Add direct sources for different times of the day
     for I in Idir_PAR
-        push!(dome, sky(scene, Idir = I/10*DL, nrays_dir = 100_000, Idif = 0.0)[1])
+        push!(dome, sky(mesh, Idir = I/10*DL, nrays_dir = 100_000, Idif = 0.0)[1])
     end
     return dome
 end
 #=
 
-The 3D scene and the light sources are then combined into a `RayTracer` object,
+The 3D mesh and the light sources are then combined into a `RayTracer` object,
 together with general settings for the ray tracing simulation chosen via `RTSettings()`.
 The most important settings refer to the Russian roulette system and the grid
 cloner (see section on Ray Tracing for details). The settings for the Russian
 roulette system include the number of times a ray will be traced
 deterministically (`maxiter`) and the probability that a ray that exceeds `maxiter`
 is terminated (`pkill`). The grid cloner is used to approximate an infinite canopy
-by replicating the scene in the different directions (`nx` and `ny` being the
+by replicating the mesh in the different directions (`nx` and `ny` being the
 number of replicates in each direction along the x and y axes, respectively). It
 is also possible to turn on parallelization of the ray tracing simulation by
 setting `parallel = true` (currently this uses Julia's builtin multithreading
@@ -296,26 +296,26 @@ capabilities).
 In addition `RTSettings()`, an acceleration structure and a splitting rule can
 be defined when creating the `RayTracer` object (see ray tracing documentation
 for details). The acceleration structure allows speeding up the ray tracing
-by avoiding testing all rays against all objects in the scene.
+by avoiding testing all rays against all objects in the mesh.
 
 =#
-function create_raytracer(scene, sources)
+function create_raytracer(mesh, sources)
     settings = RTSettings(pkill = 0.9, maxiter = 4, nx = 5, ny = 5, parallel = true)
-    RayTracer(scene, sources, settings = settings, acceleration = BVH,
+    RayTracer(mesh, sources, settings = settings, acceleration = BVH,
                      rule = SAH{3}(5, 10));
 end
 #=
 
 The actual ray tracing simulation is performed by calling the `trace!()` method
 on the ray tracing object. This will trace all rays from all light sources and
-update the radiant power absorbed by the different surfaces in the scene inside
+update the radiant power absorbed by the different surfaces in the mesh inside
 the `Material` objects (see `feed!()` above):
 
 =#
 function run_raytracer!(forest; DOY = 182)
-    scene   = create_scene(forest)
-    sources = create_sky(scene = scene, DOY = DOY)
-    rtobj   = create_raytracer(scene, sources)
+    mesh   = create_scene(forest)
+    sources = create_sky(mesh = mesh, DOY = DOY)
+    rtobj   = create_raytracer(mesh, sources)
     trace!(rtobj)
     return nothing
 end
@@ -544,8 +544,8 @@ end
 
 As in the previous example, it makes sense to visualize the forest with a soil
 tile beneath it. Unlike in the previous example, we will construct the soil tile
-using a dedicated graph and generate a `Scene` object which can later be
-merged with the rest of scene generated in daily step:
+using a dedicated graph and generate a `Mesh` object which can later be
+merged with the rest of mesh generated in daily step:
 
 =#
 Base.@kwdef struct Soil <: VirtualPlantLab.Node
@@ -553,23 +553,24 @@ Base.@kwdef struct Soil <: VirtualPlantLab.Node
     width::Float64
 end
 function VirtualPlantLab.feed!(turtle::Turtle, s::Soil, data)
-    Rectangle!(turtle, length = s.length, width = s.width, colors = RGB(255/255, 236/255, 179/255))
+    Rectangle!(turtle, length = s.length, width = s.width, colors = RGB(255/255, 236/255, 179/255),
+               materials = Lambertian(τ = 0.0, ρ = 0.21))
 end
 soil_graph = RA(-90.0) + T(Vec(0.0, 10.0, 0.0)) + ## Moves into position
              Soil(length = 20.0, width = 20.0) ## Draws the soil tile
-soil = Scene(Graph(axiom = soil_graph));
+soil = Mesh(Graph(axiom = soil_graph));
 render(soil, axes = false)
 #=
 
-And the following function renders the entire scene (notice that we need to
-use `display()` to force the rendering of the scene when called within a loop
+And the following function renders the entire mesh (notice that we need to
+use `display()` to force the rendering of the mesh when called within a loop
 or a function):
 
 =#
 function render_forest(forest, soil)
-    scene = Scene(vec(forest)) ## create scene from forest
-    scene = Scene([scene, soil]) ## merges the two scenes
-    render(scene)
+    mesh = Mesh(vec(forest)) ## create mesh from forest
+    mesh = Mesh([mesh, soil]) ## merges the two scenes
+    render(mesh)
 end
 #=
 
diff --git a/test/snowflakes.jl b/test/snowflakes.jl
index abf6ab0..161c8c4 100644
--- a/test/snowflakes.jl
+++ b/test/snowflakes.jl
@@ -48,7 +48,7 @@ class to implement the edges of the snowflake. This can be achieved as follows:
 
 =#
 using VirtualPlantLab
-import GLMakie ## Import rather than "using" to avoid masking Scene
+import GLMakie ## Import rather than "using" to avoid masking Mesh
 using ColorTypes ## To define colors for the rendering
 module sn
     import VirtualPlantLab
@@ -128,7 +128,7 @@ After defining the method, we can now call the function render on the graph to
 generate a 3D interactive image of the Koch snowflake in the current state
 
 =#
-sc = Scene(Koch)
+sc = Mesh(Koch)
 render(sc, axes = false)
 #=
 
@@ -138,7 +138,7 @@ snowflake. Let's execute the rules once to verify that we get the 2nd iteration
 
 =#
 rewrite!(Koch)
-render(Scene(Koch), axes = false)
+render(Mesh(Koch), axes = false)
 #=
 
 And two more times
@@ -147,7 +147,7 @@ And two more times
 for i in 1:3
     rewrite!(Koch)
 end
-render(Scene(Koch), axes = false)
+render(Mesh(Koch), axes = false)
 #=
 
 # Other snowflake fractals
@@ -172,13 +172,13 @@ iterations and render the results
 =#
 ## First iteration
 rewrite!(Koch2)
-render(Scene(Koch2), axes = false)
+render(Mesh(Koch2), axes = false)
 ## Second iteration
 rewrite!(Koch2)
-render(Scene(Koch2), axes = false)
+render(Mesh(Koch2), axes = false)
 ## Third iteration
 rewrite!(Koch2)
-render(Scene(Koch2), axes = false)
+render(Mesh(Koch2), axes = false)
 #=
 
 This is know as [Koch
@@ -190,7 +190,7 @@ axiom:
 =#
 axiomCesaro = sn.E(L) + RU(90.0) + sn.E(L) + RU(90.0) + sn.E(L) + RU(90.0) + sn.E(L)
 Cesaro = Graph(axiom = axiomCesaro, rules = (rule2,))
-render(Scene(Cesaro), axes = false)
+render(Mesh(Cesaro), axes = false)
 #=
 
 And, as before, let's go through the first three iterations
@@ -198,10 +198,10 @@ And, as before, let's go through the first three iterations
 =#
 ## First iteration
 rewrite!(Cesaro)
-render(Scene(Cesaro), axes = false)
+render(Mesh(Cesaro), axes = false)
 ## Second iteration
 rewrite!(Cesaro)
-render(Scene(Cesaro), axes = false)
+render(Mesh(Cesaro), axes = false)
 ## Third iteration
 rewrite!(Cesaro)
-render(Scene(Cesaro), axes = false)
+render(Mesh(Cesaro), axes = false)
diff --git a/test/tree.jl b/test/tree.jl
index 3bec99b..6b627c7 100644
--- a/test/tree.jl
+++ b/test/tree.jl
@@ -186,11 +186,11 @@ newtree = simulate(tree, getInternode, 2)
 The binary tree after two iterations has two branches, as expected:
 
 =#
-render(Scene(newtree))
+render(Mesh(newtree))
 #=
 
 Notice how the lengths of the prisms representing internodes decreases as the branching order increases, as the internodes are younger (i.e. were generated fewer generations ago). Further steps will generate a structure that is more tree-like.
 
 =#
 newtree = simulate(newtree, getInternode, 15)
-render(Scene(newtree))
+render(Mesh(newtree))