diff --git a/examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl b/examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl
new file mode 100644
index 00000000000..67c664471d5
--- /dev/null
+++ b/examples/t8code_3d_fv/elixir_advection_basic_hybrid.jl
@@ -0,0 +1,43 @@
+using OrdinaryDiffEq
+using Trixi
+
+advection_velocity = (0.2, -0.7, 0.5)
+equations = LinearScalarAdvectionEquation3D(advection_velocity)
+
+initial_condition = initial_condition_convergence_test
+
+solver = FV(order = 2, extended_reconstruction_stencil = false,
+            surface_flux = flux_lax_friedrichs)
+
+# cmesh = Trixi.cmesh_new_periodic_hybrid()
+cmesh = Trixi.cmesh_quad_3d(periodicity = (true, true, true))
+# cmesh = Trixi.cmesh_new_periodic_tri()
+mesh = T8codeMesh(cmesh, solver,
+                  initial_refinement_level = 5)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+ode = semidiscretize(semi, (0.0, 10.0));
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 1,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 0.3)
+
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback,
+                        stepsize_callback, save_solution)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+summary_callback()
diff --git a/src/meshes/t8code_mesh.jl b/src/meshes/t8code_mesh.jl
index 8f6ec2a93da..0dafdf0af9e 100644
--- a/src/meshes/t8code_mesh.jl
+++ b/src/meshes/t8code_mesh.jl
@@ -1579,6 +1579,57 @@ function cmesh_quad(; comm = mpi_comm(), periodicity = (true, true))::t8_cmesh_t
     return cmesh
 end
 
+function cmesh_quad_3d(; comm = mpi_comm(), periodicity = (true, true, true))::t8_cmesh_t
+    n_dims = 3
+    vertices = [ # Just all vertices of all trees. partly duplicated
+        -1.0, -1.0, -1.0, # tree 0, quad
+        1.0, -1.0, -1.0,
+        -1.0, 1.0, -1.0,
+        1.0, 1.0, -1.0,
+        -1.0, -1.0, 1.0,
+        1.0, -1.0, 1.0,
+        -1.0, 1.0, 1.0,
+        1.0, 1.0, 1.0,
+    ]
+
+    # Generally, one can define other geometries. But besides linear the other
+    # geometries in t8code do not have C interface yet.
+    linear_geom = t8_geometry_linear_new(n_dims)
+
+    # This is how the cmesh looks like. The numbers are the tree numbers:
+    #
+    #   +---+
+    #   |   |
+    #   | 0 |
+    #   |   |
+    #   +---+
+    #
+
+    cmesh_ref = Ref(t8_cmesh_t())
+    t8_cmesh_init(cmesh_ref)
+    cmesh = cmesh_ref[]
+
+    # Use linear geometry
+    t8_cmesh_register_geometry(cmesh, linear_geom)
+    t8_cmesh_set_tree_class(cmesh, 0, T8_ECLASS_HEX)
+
+    t8_cmesh_set_tree_vertices(cmesh, 0, @views(vertices[(1 + 0):end]), 8)
+
+    if periodicity[1]
+        t8_cmesh_set_join(cmesh, 0, 0, 0, 1, 0)
+    end
+    if periodicity[2]
+        t8_cmesh_set_join(cmesh, 0, 0, 2, 3, 0)
+    end
+    if periodicity[3]
+        t8_cmesh_set_join(cmesh, 0, 0, 4, 5, 0)
+    end
+
+    t8_cmesh_commit(cmesh, comm)
+
+    return cmesh
+end
+
 function cmesh_new_periodic_tri(; comm = mpi_comm())::t8_cmesh_t
     n_dims = 2
     vertices = [ # Just all vertices of all trees. partly duplicated
diff --git a/src/solvers/fv_t8code/containers.jl b/src/solvers/fv_t8code/containers.jl
index 5694e6e5dd2..44dcd19a5e5 100644
--- a/src/solvers/fv_t8code/containers.jl
+++ b/src/solvers/fv_t8code/containers.jl
@@ -248,7 +248,7 @@ function init_elements!(elements, mesh::T8codeMesh, solver::FV)
     return nothing
 end
 
-@inline function init_extended_reconstruction_stencil!(corners, elements, solver)
+@inline function init_extended_reconstruction_stencil!(corners, elements, solver::FV)
     if solver.order != 2
         return nothing
     end
@@ -327,7 +327,7 @@ end
 
 function init_reconstruction_stencil!(elements, interfaces, boundaries,
                                       communication_data,
-                                      mesh, equations, solver)
+                                      mesh, equations, solver::FV)
     if solver.order != 2
         return nothing
     end
diff --git a/src/solvers/fv_t8code/fv.jl b/src/solvers/fv_t8code/fv.jl
index 08c82df056b..6445a2d7b23 100644
--- a/src/solvers/fv_t8code/fv.jl
+++ b/src/solvers/fv_t8code/fv.jl
@@ -216,7 +216,7 @@ function rhs!(du, u, t, mesh::T8codeMesh, equations,
     return nothing
 end
 
-function calc_gradient_reconstruction!(u, mesh, equations, solver, cache)
+function calc_gradient_reconstruction!(u, mesh::T8codeMesh{2}, equations, solver, cache)
     if solver.order == 1
         return nothing
     elseif solver.order > 2
@@ -237,7 +237,7 @@ function calc_gradient_reconstruction!(u, mesh, equations, solver, cache)
     # (A^T A)^-1 = determinant_factor * [a3 -a2; -a2, a1]
 
     # A^T b = [d1; d2]
-    # Note: A^T b depends on the variable v. Using structure [d/e, v]
+    # Note: A^T b depends on the variable v. Using structure [1/2, v]
     a = zeros(eltype(u), 3)             # [a1, a2, a3]
     d = zeros(eltype(u), size(u, 1), 2) # [d1(v), d2(v)]
 
@@ -322,7 +322,7 @@ end
 
 @inline function calc_gradient_simple!(stencil, n_stencil_neighbors, element, a, d,
                                        reconstruction_stencil, reconstruction_distance,
-                                       solution_data, equations)
+                                       solution_data, equations::AbstractEquations{2})
     for i in 1:n_stencil_neighbors
         # stencil=1 contains information from all neighbors
         # stencil=2,...,n_stencil_neighbors+1 is computed without (stencil+1)-th neighbor's information
@@ -397,6 +397,159 @@ end
     return nothing
 end
 
+function calc_gradient_reconstruction!(u, mesh::T8codeMesh{3}, equations, solver, cache)
+    if solver.order == 1
+        return nothing
+    elseif solver.order > 2
+        error("Order $(solver.order) is not supported yet!")
+    end
+
+    (; elements, communication_data) = cache
+    (; reconstruction_stencil, reconstruction_distance, reconstruction_corner_elements, reconstruction_gradient, reconstruction_gradient_limited) = elements
+    (; solution_data) = communication_data
+
+    # A         N x 3 Matrix, where N is the number of stencil neighbors
+    # A^T A     3 x 3 Matrix
+    # b         N     Vector
+    # A^T b     3     Vector
+
+    # Matrix/vector notation
+    # A^T A = [a11 a12 a13; a12 a22 a23; a13 a23 a33]
+    # det(A^T A) = a11 * (a22*a33 - a23^2) - a12 * (a12*a33 - a13*a23) + a13 * (a12*a23 - a22*a13)
+    # (A^T A)^-1 = 1/det(A^T A) * [(a33*a22-a23^2) -(a33*a12-a23*a13) (a23*a12-a22*a13) ;
+    #                               ***             (a11*a33-a13^2)   -(a23*a11-a12*a13);
+    #                               ***            ***                (a11*a22-a12^2)   ]
+
+    # A^T b = [d1; d2; d3]
+    # Note: A^T b depends on the variable v. Using structure [1/2/3, v]
+    a = zeros(eltype(u), 6)             # [a11, a12, a13, a22, a23, a33]
+    d = zeros(eltype(u), size(u, 1), 3) # [d1(v), d2(v), d3(v)]
+
+    # Parameter for limiting weights
+    lambda = [0.0, 1.0]
+    # lambda = [1.0, 0.0] # No limiting
+    r = 4
+    epsilon = 1.0e-13
+
+    for element in eachelement(mesh, solver, cache)
+        if solver.extended_reconstruction_stencil
+            n_stencil_neighbors = elements.num_faces[element]
+            error("not yet supported!")
+        else
+            n_stencil_neighbors = length(reconstruction_stencil[element])
+        end
+
+        for stencil in 1:(n_stencil_neighbors + 1)
+            # Reset variables
+            a .= zero(eltype(u))
+            # A^T b = [d1(v), d2(v), d3(v)]
+            # Note: A^T b depends on the variable v. Using vectors with d[v, 1/2/3]
+            d .= zero(eltype(u))
+
+            if solver.extended_reconstruction_stencil
+                calc_gradient_extended!(stencil, n_stencil_neighbors, element, a, d,
+                                        reconstruction_stencil, reconstruction_distance,
+                                        reconstruction_corner_elements, solution_data,
+                                        equations, solver, cache)
+                error("Not supported yet")
+            else
+                calc_gradient_simple!(stencil, n_stencil_neighbors, element, a, d,
+                                      reconstruction_stencil, reconstruction_distance,
+                                      solution_data, equations)
+            end
+
+            # Divide by determinant
+            # Determinant    = a11  * (a22  * a33  -  a23^2 ) - a12  * (a12  * a33  - a13  * a23 ) + a13  * (a12  * a23  - a22  * a13 )
+            AT_A_determinant = a[1] * (a[4] * a[6] -  a[5]^2) - a[2] * (a[2] * a[6] - a[3] * a[5]) + a[3] * (a[2] * a[5] - a[4] * a[3])
+            if isapprox(AT_A_determinant, 0.0)
+                for v in eachvariable(equations)
+                    reconstruction_gradient[:, v, stencil, element] .= 0.0
+                end
+                continue
+            end
+
+            # det(A^T A) = a11 * (a22*a33 - a23^2) - a12 * (a12*a33 - a13*a23) + a13 * (a12*a23 - a22*a13)
+            # (A^T A)^-1 = 1/det(A^T A) * [(a33*a22-a23^2) -(a33*a12-a23*a13) (a23*a12-a22*a13) ;
+            #                               ***             (a11*a33-a13^2)   -(a23*a11-a12*a13);
+            #                               ***            ***                (a11*a22-a12^2)   ]
+            # = [m11 m12 m13; m12 m22 m23; m13 m23 m33]
+            m11 = a[6] * a[4] - a[5]^2
+            m12 = -(a[6] * a[2] - a[5] * a[3])
+            m13 = (a[5] * a[2] - a[4] * a[3])
+            m22 = (a[1] * a[6] - a[3]^2)
+            m23 = -(a[5] * a[1] - a[2] * a[3])
+            m33 = (a[1] * a[4] - a[2]^2)
+
+            # Solving least square problem
+            for v in eachvariable(equations)
+                reconstruction_gradient[1, v, stencil, element] = 1 / AT_A_determinant * (m11 * d[v, 1] + m12 * d[v, 2] + m13 * d[v, 3])
+                reconstruction_gradient[2, v, stencil, element] = 1 / AT_A_determinant * (m12 * d[v, 1] + m22 * d[v, 2] + m23 * d[v, 3])
+                reconstruction_gradient[3, v, stencil, element] = 1 / AT_A_determinant * (m13 * d[v, 1] + m23 * d[v, 2] + m33 * d[v, 3])
+            end
+        end
+
+        # Get limited reconstruction gradient by weighting the just computed
+        weight_sum = zero(eltype(reconstruction_gradient))
+        for v in eachvariable(equations)
+            reconstruction_gradient_limited[:, v, element] .= zero(eltype(reconstruction_gradient_limited))
+            weight_sum = zero(eltype(reconstruction_gradient))
+            for j in 1:(n_stencil_neighbors + 1)
+                gradient = get_node_coords(reconstruction_gradient, equations, solver,
+                                           v, j, element)
+                indicator = sum(gradient .^ 2)
+                lambda_j = (j == 1) ? lambda[1] : lambda[2]
+                weight = (lambda_j / (indicator + epsilon))^r
+                for x in axes(reconstruction_gradient_limited, 1)
+                    reconstruction_gradient_limited[x, v, element] += weight *
+                                                                      gradient[x]
+                end
+                weight_sum += weight
+            end
+            for x in axes(reconstruction_gradient_limited, 1)
+                reconstruction_gradient_limited[x, v, element] /= weight_sum
+            end
+        end
+    end
+
+    exchange_gradient_data!(reconstruction_gradient_limited, mesh, equations, solver,
+                            cache)
+
+    return nothing
+end
+
+@inline function calc_gradient_simple!(stencil, n_stencil_neighbors, element, a, d,
+                                       reconstruction_stencil, reconstruction_distance,
+                                       solution_data, equations::AbstractEquations{3})
+    for i in 1:n_stencil_neighbors
+        # stencil=1 contains information from all neighbors
+        # stencil=2,...,n_stencil_neighbors+1 is computed without (stencil+1)-th neighbor's information
+        if i + 1 != stencil
+            neighbor = reconstruction_stencil[element][i]
+            distance = reconstruction_distance[element][i]
+            a[1] += distance[1]^2 # = a11
+            a[2] += distance[1] * distance[2] # = a12
+            a[3] += distance[1] * distance[3]
+            a[4] += distance[2]^2
+            a[5] += distance[2] * distance[3]
+            a[6] += distance[3]^2
+
+            for v in eachvariable(equations)
+                d[v, 1] += distance[1] *
+                           (solution_data[neighbor].u[v] -
+                            solution_data[element].u[v])
+                d[v, 2] += distance[2] *
+                           (solution_data[neighbor].u[v] -
+                            solution_data[element].u[v])
+                d[v, 3] += distance[3] *
+                           (solution_data[neighbor].u[v] -
+                            solution_data[element].u[v])
+            end
+        end
+    end
+
+    return nothing
+end
+
 function prolong2interfaces!(cache, mesh::T8codeMesh, equations, solver::FV)
     (; interfaces, communication_data) = cache
     (; solution_data, domain_data, gradient_data) = communication_data