Fix BUG-004/005: mixed-materials load accumulation and plane-stress postprocessing

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,
@@ -82,12 +83,21 @@ buffer = setup_domainbuffers(domains);
 # [MaterialModelsBase#12](https://github.com/KnutAM/MaterialModelsBase.jl/issues/12).
 
 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, σ)
+    ϵ = symmetric(∇u) # Already the reduced (in-plane) strain
+    return calculate_stress(m.stress_state, m.material, ϵ, qp_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);
+function calculate_stress(stress_state, m::LinearElastic, ϵ, qp_state)
+    σ, _, _ = material_response(stress_state, m, ϵ, qp_state)
+    return σ
+end
+function calculate_stress(stress_state, m::Plastic, ϵ, qp_state)
+    ## `qp_state.ϵp` is the full 3d converged plastic strain, whose out-of-plane
+    ## component already accounts for the plane-stress constraint. Using it here,
+    ## rather than re-running `Plastic`'s own `material_response`, avoids advancing
+    ## the (already converged) state a second time during postprocessing.
+    ϵₑ = ϵ - MaterialModelsBase.reduce_tensordim(stress_state, qp_state.ϵp)
+    return calculate_stress(stress_state, m.elastic, ϵₑ, qp_state)
+end;
 
 # And then we create the QuadPointEvaluator including this function
 qe = QuadPointEvaluator{SymmetricTensor{2,2,Float64,3}}(buffer, calculate_stress);
@@ -111,7 +121,11 @@ function solve_nonlinear_timehistory(buffer, dh, ch, lh, l2_proj, qp_evaluator;
         ## Update and apply the Dirichlet boundary conditions
         update!(ch, t)
         apply!(a, ch)
+        fill!(fext, 0)
         apply!(fext, lh, t)
+        fext_check = zeros(length(fext)) #src
+        apply!(fext_check, lh, t)        #src
+        @test fext ≈ fext_check          #src
         for i in 1:maxiter
             ## Assemble the system
             assembler = start_assemble(K, r)
@@ -145,6 +159,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.