diff --git a/docs/Manifest.toml b/docs/Manifest.toml
index 592287a7..4d065df6 100644
--- a/docs/Manifest.toml
+++ b/docs/Manifest.toml
@@ -1198,7 +1198,7 @@ version = "0.1.3"
 
 [[deps.MaterialModelsBase]]
 deps = ["ForwardDiff", "StaticArrays", "Tensors"]
-git-tree-sha1 = "bc0d65963d330270c4d777d39cfba7f207cadb05"
+git-tree-sha1 = "cf91f9ea44df2be0f1c4d147038228421b4b660b"
 repo-rev = "main"
 repo-url = "https://github.com/KnutAM/MaterialModelsBase.jl.git"
 uuid = "af893363-701d-44dc-8b1e-d9a2c129bfc9"
@@ -1213,7 +1213,7 @@ version = "0.6.9"
 
 [[deps.MechanicalMaterialModels]]
 deps = ["ForwardDiff", "LinearAlgebra", "MaterialModelsBase", "Newton", "StaticArrays", "Tensors"]
-git-tree-sha1 = "ddfce3820a135ba597d356bf4cd877d4962670f2"
+git-tree-sha1 = "126b1b35473bfcd2f66cc7e42d6c90075f097c0d"
 repo-rev = "main"
 repo-url = "https://github.com/KnutAM/MechanicalMaterialModels.jl.git"
 uuid = "b3282f9b-607f-4337-ab95-e5488ab5652c"
diff --git a/docs/src/literate_tutorials/mixed_materials.jl b/docs/src/literate_tutorials/mixed_materials.jl
index 10c2e9d6..a75e1b5e 100644
--- a/docs/src/literate_tutorials/mixed_materials.jl
+++ b/docs/src/literate_tutorials/mixed_materials.jl
@@ -17,6 +17,7 @@
 using Ferrite, FerriteAssembly, FerriteMeshParser
 using MaterialModelsBase, MechanicalMaterialModels, WriteVTK
 using Downloads: download
+using Test #src
 
 # ## Setup Ferrite quantities
 # We start by the downloading and parsing the grid containing a central inclusion,
@@ -77,17 +78,15 @@ buffer = setup_domainbuffers(domains);
 # `Ferrite`'s `L2Projector`.
 # 
 # First, we define a function to calculate the stresses for each material.
-# Note that here we have to use some internals from `MechanicalMaterialModels.jl`,
-# but this should be solved with 
-# [MaterialModelsBase#12](https://github.com/KnutAM/MaterialModelsBase.jl/issues/12).
+# We use `MaterialModelsBase.stress_from_state`, which calculates the stress
+# conjugated to a given strain that is consistent with an already-converged
+# `state`, without invoking any local iteration that would advance history
+# variables. For a `ReducedStressState`, such as our plane-stress case, this
+# correctly accounts for the reduced dimensionality (e.g. plane stress).
 
 function calculate_stress(m::ReducedStressState, u, ∇u, qp_state)
-    ϵ = MaterialModelsBase.expand_tensordim(m.stress_state, symmetric(∇u))
-    σ = calculate_stress(m.material, ϵ, qp_state)
-    return MaterialModelsBase.reduce_tensordim(m.stress_state, σ)
-end
-calculate_stress(m::LinearElastic, ϵ, qp_state) = m.C ⊡ ϵ
-calculate_stress(m::Plastic, ϵ, qp_state) = calculate_stress(m.elastic, ϵ - qp_state.ϵp, qp_state);
+    return stress_from_state(m, symmetric(∇u), qp_state)
+end;
 
 # And then we create the QuadPointEvaluator including this function
 qe = QuadPointEvaluator{SymmetricTensor{2,2,Float64,3}}(buffer, calculate_stress);
@@ -105,13 +104,19 @@ function solve_nonlinear_timehistory(buffer, dh, ch, lh, l2_proj, qp_evaluator;
     r = zeros(ndofs(dh))
     fext = zeros(ndofs(dh))
     a = zeros(ndofs(dh))
+    fext_unit = zeros(ndofs(dh)) #src
+    apply!(fext_unit, lh, 1.0)   #src
     ## Prepare postprocessing
     pvd = paraview_collection("multiple_materials")
     for (n, t) in enumerate(time_history)
         ## Update and apply the Dirichlet boundary conditions
         update!(ch, t)
         apply!(a, ch)
+        fill!(fext, 0)
         apply!(fext, lh, t)
+        ## The applied traction is linear in `t`, so if `fext` accumulated loads #src
+        ## from previous steps instead of being reset, it would not match `t * fext_unit`. #src
+        @test fext ≈ t * fext_unit #src
         for i in 1:maxiter
             ## Assemble the system
             assembler = start_assemble(K, r)
@@ -145,6 +150,27 @@ function solve_nonlinear_timehistory(buffer, dh, ch, lh, l2_proj, qp_evaluator;
 end;
 solve_nonlinear_timehistory(buffer, dh, ch, lh, proj, qe; time_history=collect(range(0, 1, 20)));
 
+## Regression checks for `calculate_stress`'s plane-stress postprocessing (hidden from docs) #src
+## Analytical reference: for isotropic plane stress with E, ν and ϵ11=0.01, ϵ22=ϵ12=0,       #src
+## σ11 = E/(1-ν^2)*ϵ11 and σ22 = E/(1-ν^2)*ν*ϵ11.                                             #src
+let                                                                                           #src
+    E, ν = 210e3, 0.3                                                                         #src
+    ϵ11 = 0.01                                                                                #src
+    ∇u = Tensor{2,2}((ϵ11, 0.0, 0.0, 0.0))                                                    #src
+    qp_state = MaterialModelsBase.initial_material_state(elastic_material)                    #src
+    σ = calculate_stress(elastic_material, zero(Vec{2}), ∇u, qp_state)                        #src
+    σ11_ref = E / (1 - ν^2) * ϵ11                                                             #src
+    σ22_ref = E / (1 - ν^2) * ν * ϵ11                                                         #src
+    @test σ[1, 1] ≈ σ11_ref                                                                   #src
+    @test σ[2, 2] ≈ σ22_ref                                                                   #src
+    ## Check that the eliminated out-of-plane stress is indeed zero                           #src
+    _, _, _, ϵ_3d = material_response(elastic_material.stress_state, elastic_material.material, symmetric(∇u), qp_state) #src
+    σ_3d = elastic_material.material.C ⊡ ϵ_3d                                                 #src
+    @test σ_3d[3, 3] ≈ 0.0 atol = 1e-6 * abs(σ11_ref)                                         #src
+end                                                                                            #src
+## Pinned regression value for the full solve's postprocessed stresses #src
+@test norm(norm.(qe.data)) ≈ 62718.61437855114                         #src
+
 #md # ## [Plain program](@id mixed_materials_plain_program)
 #md #
 #md # Here follows a version of the program without any comments.