diff --git a/docs/src/literate/helmholtz.jl b/docs/src/literate/helmholtz.jl
index 1d72a60291..93c3b3fa87 100644
--- a/docs/src/literate/helmholtz.jl
+++ b/docs/src/literate/helmholtz.jl
@@ -93,14 +93,14 @@ function doassemble(cellvalues::CellScalarValues{dim}, facevalues::FaceScalarVal
     b = 1.0
     f = zeros(ndofs(dh))
     assembler = start_assemble(K, f)
-    
+
     n_basefuncs = getnbasefunctions(cellvalues)
     global_dofs = zeros(Int, ndofs_per_cell(dh))
 
     fe = zeros(n_basefuncs) # Local force vector
     Ke = zeros(n_basefuncs, n_basefuncs) # Local stiffness mastrix
 
-    @inbounds for (cellcount, cell) in enumerate(CellIterator(dh))
+    for (cellcount, cell) in enumerate(CellIterator(dh))
         fill!(Ke, 0)
         fill!(fe, 0)
         coords = getcoordinates(cell)
diff --git a/docs/src/literate/incompressible_elasticity.jl b/docs/src/literate/incompressible_elasticity.jl
index 9dd4f0c8bb..6a867c103c 100644
--- a/docs/src/literate/incompressible_elasticity.jl
+++ b/docs/src/literate/incompressible_elasticity.jl
@@ -127,7 +127,7 @@ function assemble_up!(Ke, fe, cell, cellvalues_u, cellvalues_p, facevalues_u, gr
     reinit!(cellvalues_p, cell)
 
     ## We only assemble lower half triangle of the stiffness matrix and then symmetrize it.
-    @inbounds for q_point in 1:getnquadpoints(cellvalues_u)
+    for q_point in 1:getnquadpoints(cellvalues_u)
         for i in 1:n_basefuncs_u
             ɛdev[i] = dev(symmetric(shape_gradient(cellvalues_u, q_point, i)))
         end
@@ -159,7 +159,7 @@ function assemble_up!(Ke, fe, cell, cellvalues_u, cellvalues_p, facevalues_u, gr
     ## We integrate the Neumann boundary using the facevalues.
     ## We loop over all the faces in the cell, then check if the face
     ## is in our `"traction"` faceset.
-    @inbounds for face in 1:nfaces(cell)
+    for face in 1:nfaces(cell)
         if onboundary(cell, face) && (cellid(cell), face) ∈ getfaceset(grid, "traction")
             reinit!(facevalues_u, cell, face)
             for q_point in 1:getnquadpoints(facevalues_u)
diff --git a/docs/src/literate/ns_vs_diffeq.jl b/docs/src/literate/ns_vs_diffeq.jl
index 9bbfa0ceeb..ef565eb247 100644
--- a/docs/src/literate/ns_vs_diffeq.jl
+++ b/docs/src/literate/ns_vs_diffeq.jl
@@ -222,7 +222,7 @@ function assemble_mass_matrix(cellvalues_v::CellVectorValues{dim}, cellvalues_p:
 
     ## It follows the assembly loop as explained in the basic tutorials.
     mass_assembler = start_assemble(M)
-    @inbounds for cell in CellIterator(dh)
+    for cell in CellIterator(dh)
         fill!(Mₑ, 0)
         Ferrite.reinit!(cellvalues_v, cell)
 
@@ -265,7 +265,7 @@ function assemble_stokes_matrix(cellvalues_v::CellVectorValues{dim}, cellvalues_
 
     ## Assembly loop
     stiffness_assembler = start_assemble(K)
-    @inbounds for cell in CellIterator(dh)
+    for cell in CellIterator(dh)
         ## Don't forget to initialize everything
         fill!(Kₑ, 0)
 
@@ -453,7 +453,7 @@ vtk_save(pvd);
 using Test                                                                  #hide
 function compute_divergence(dh, u, cellvalues_v)                            #hide
     divv = 0.0                                                              #hide
-    @inbounds for (i,cell) in enumerate(CellIterator(dh))                   #hide
+    for cell in CellIterator(dh)                                            #hide
         Ferrite.reinit!(cellvalues_v, cell)                                 #hide
         for q_point in 1:getnquadpoints(cellvalues_v)                       #hide
             dΩ = getdetJdV(cellvalues_v, q_point)                           #hide
diff --git a/docs/src/literate/quasi_incompressible_hyperelasticity.jl b/docs/src/literate/quasi_incompressible_hyperelasticity.jl
index 3f3cdde770..7a70fc7f7d 100644
--- a/docs/src/literate/quasi_incompressible_hyperelasticity.jl
+++ b/docs/src/literate/quasi_incompressible_hyperelasticity.jl
@@ -162,7 +162,7 @@ end;
 function calculate_element_volume(cell, cellvalues_u, ue)
     reinit!(cellvalues_u, cell)
     evol::Float64=0.0;
-    @inbounds for qp in 1:getnquadpoints(cellvalues_u)
+    for qp in 1:getnquadpoints(cellvalues_u)
         dΩ = getdetJdV(cellvalues_u, qp)
         ∇u = function_gradient(cellvalues_u, qp, ue)
         F = one(∇u) + ∇u
@@ -175,7 +175,7 @@ end;
 # and then assembled over all the cells (elements)
 function calculate_volume_deformed_mesh(w, dh::DofHandler, cellvalues_u)
     evol::Float64 = 0.0;
-    @inbounds for cell in CellIterator(dh)
+    for cell in CellIterator(dh)
         global_dofs = celldofs(cell)
         nu = getnbasefunctions(cellvalues_u)
         global_dofs_u = global_dofs[1:nu]
@@ -199,7 +199,7 @@ function assemble_element!(Ke, fe, cell, cellvalues_u, cellvalues_p, mp, ue, pe)
     n_basefuncs_u = getnbasefunctions(cellvalues_u)
     n_basefuncs_p = getnbasefunctions(cellvalues_p)
 
-    @inbounds for qp in 1:getnquadpoints(cellvalues_u)
+    for qp in 1:getnquadpoints(cellvalues_u)
         dΩ = getdetJdV(cellvalues_u, qp)
         ## Compute deformation gradient F
         ∇u = function_gradient(cellvalues_u, qp, ue)
diff --git a/docs/src/literate/transient_heat_equation.jl b/docs/src/literate/transient_heat_equation.jl
index 76e365ad76..c716b0af64 100644
--- a/docs/src/literate/transient_heat_equation.jl
+++ b/docs/src/literate/transient_heat_equation.jl
@@ -114,7 +114,7 @@ function doassemble_K!(K::SparseMatrixCSC, f::Vector, cellvalues::CellScalarValu
 
     assembler = start_assemble(K, f)
 
-    @inbounds for cell in CellIterator(dh)
+    for cell in CellIterator(dh)
 
         fill!(Ke, 0)
         fill!(fe, 0)
@@ -148,7 +148,7 @@ function doassemble_M!(M::SparseMatrixCSC, cellvalues::CellScalarValues{dim}, dh
 
     assembler = start_assemble(M)
 
-    @inbounds for cell in CellIterator(dh)
+    for cell in CellIterator(dh)
 
         fill!(Me, 0)