From 0111539c84d3398a10a5bc9f7031072ac22dafc2 Mon Sep 17 00:00:00 2001 From: ClaudeBot Date: Thu, 17 Sep 2026 08:57:52 -0400 Subject: [PATCH 1/4] Fix BUG-004/005: mixed-materials tutorial load accumulation and plane-stress postprocessing BUG-004: `fext` was allocated once outside the time loop and never cleared, so `apply!(fext, lh, t)` added each step's load onto the running total instead of replacing it, ending with 10x the intended final load after the 20-step history. Fixed by `fill!(fext, 0)` before each `apply!` call. BUG-005: `calculate_stress` computed plane-stress output by zero-padding the in-plane strain's out-of-plane components and applying the full 3d elastic stiffness tensor directly, which is only correct for plane strain, not plane stress (whose eliminated out-of-plane strain must satisfy sigma_33 = 0). Fixed by delegating to `MaterialModelsBase.material_response` on the actual stress state + material pair for the elastic case (which performs the same plane-stress-consistent solve used during assembly), and, for the plastic case, using the already-converged 3d plastic strain (reduced to in-plane components) to form the elastic strain directly rather than re-invoking `Plastic`'s own `material_response` a second time (which would otherwise double-advance the already-converged state during postprocessing). Added hidden (`#src`) regression checks: per-timestep verification that `fext` matches a freshly evaluated load vector, an elastic plane-stress analytical reference check (matches E/(1-nu^2) form) including an out-of-plane stress check, and a pinned end-to-end regression value for the full solve's postprocessed stresses. Verified numerically against MaterialModelsBase.material_response directly: for E=210e3, nu=0.3, eps11=0.01, the corrected formula gives (sigma11, sigma22) = (2307.69, 692.31), matching the documented correct reference (vs. the buggy (2826.92, 1211.54)). Test results: full Pkg.test() suite passes (3236 assertions across all testsets). All 5 unaffected tutorials, all 6 how-tos, and the fixed mixed_materials tutorial (with its new #src regression tests) run successfully. Full docs build succeeds with only the pre-existing benign viscoelasticity image-size-threshold warning. Independent Codex review (plan + final diff, same thread) raised one remaining medium finding both times: mixed_materials.png is a downloaded external asset (via FerriteAssembly.asset_url) rather than a file generated by the docs build, so this PR's fix does not regenerate the published tutorial image, which will keep showing results from the old, incorrect loading/postprocessing until that asset is separately replaced upstream. Co-Authored-By: Claude Sonnet 5 Claude-Session: https://claude.ai/code/session_01KySgpVgJ5fWR9AvPK1JQoU --- .../src/literate_tutorials/mixed_materials.jl | 45 ++++++++++++++++--- 1 file changed, 40 insertions(+), 5 deletions(-) diff --git a/docs/src/literate_tutorials/mixed_materials.jl b/docs/src/literate_tutorials/mixed_materials.jl index 10c2e9d6..d7d9e2b0 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, @@ -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. From 6f2893c97630b3b781244b80223bd8114a737795 Mon Sep 17 00:00:00 2001 From: ClaudeBot Date: Thu, 24 Sep 2026 07:45:28 -0400 Subject: [PATCH 2/4] Use MaterialModelsBase.stress_from_state for BUG-005 plane-stress postprocessing Replaces the hand-rolled calculate_stress dispatch (manual LinearElastic/Plastic branches calling material_response directly, with reduce_tensordim plumbing for the plastic elastic-strain computation) with MaterialModelsBase's new stress_from_state(m::ReducedStressState, strain, state) API (MaterialModelsBase.jl commit 03e818c, resolving the MaterialModelsBase#12 issue the tutorial's own comment previously referenced as a TODO). MechanicalMaterialModels.jl PR #13 (commit b90a3c9) adds the required stress_from_state implementations for LinearElastic and Plastic, so the tutorial's calculate_stress is now a single generic line relying on the public package API instead of tutorial-level internals. Updated docs/Manifest.toml to pick up the new MaterialModelsBase and MechanicalMaterialModels commits (both already tracked via [sources] at rev = "main"; no compat-bound changes needed). The root Manifest.toml is gitignored and was refreshed locally only. Verified the fix is behavior-preserving: the pinned end-to-end regression value (norm(norm.(qe.data)) from the previous BUG-005 fix) is bit-for-bit identical after this change, and the existing elastic analytical reference + out-of-plane checks still pass, confirming stress_from_state reproduces the same physics via the officially supported API instead of the manual workaround. Test results: full Pkg.test() suite passes (3236 assertions). All 5 unaffected tutorials, all 6 how-tos, and mixed_materials (with its #src regression tests) run successfully. Full docs build succeeds with only the pre-existing benign viscoelasticity image-size-threshold warning. Independent Codex review: plan review raised one finding (accepted) about not blindly updating the pinned regression value on any diff, since that could mask a real regression - addressed by only relying on the pinned value already unconditionally matching to the last digit, without touching it. Final diff review returned no findings. Co-Authored-By: Claude Sonnet 5 Claude-Session: https://claude.ai/code/session_01KySgpVgJ5fWR9AvPK1JQoU --- docs/Manifest.toml | 4 ++-- .../src/literate_tutorials/mixed_materials.jl | 23 +++++-------------- 2 files changed, 8 insertions(+), 19 deletions(-) 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 d7d9e2b0..163f282b 100644 --- a/docs/src/literate_tutorials/mixed_materials.jl +++ b/docs/src/literate_tutorials/mixed_materials.jl @@ -78,25 +78,14 @@ 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) - ϵ = symmetric(∇u) # Already the reduced (in-plane) strain - return calculate_stress(m.stress_state, m.material, ϵ, qp_state) -end -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) + return stress_from_state(m, symmetric(∇u), qp_state) end; # And then we create the QuadPointEvaluator including this function From 85e5d7b9a172ff351e4b09e8d92e336dcaf4c900 Mon Sep 17 00:00:00 2001 From: ClaudeBot Date: Thu, 24 Sep 2026 08:04:49 -0400 Subject: [PATCH 3/4] Address review: make the fext regression check independently meaningful MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit KnutAM noted the previous check ("This doesn't test anything?") compared fext against a duplicate of the same fill!+apply! computation performed immediately above it, which reads as a tautology rather than an independent verification. Replaced it with a check based on the load's known linearity in t: compute fext_unit = apply!(_, lh, 1.0) once outside the loop, then assert fext ≈ t * fext_unit at each step. This is an independent, closed-form reference rather than a re-run of the same code path, and was verified to actually catch the BUG-004 regression (manually re-removed the `fill!(fext, 0)` line in a scratch copy and confirmed the new test fails with a clear mismatch, then reverted). Re-ran full Pkg.test() (3236 assertions), all tutorials/how-tos, and the docs build - all pass; the pinned end-to-end regression value (norm(norm.(qe.data))) is unchanged. Co-Authored-By: Claude Sonnet 5 Claude-Session: https://claude.ai/code/session_01KySgpVgJ5fWR9AvPK1JQoU --- docs/src/literate_tutorials/mixed_materials.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/docs/src/literate_tutorials/mixed_materials.jl b/docs/src/literate_tutorials/mixed_materials.jl index 163f282b..725d5364 100644 --- a/docs/src/literate_tutorials/mixed_materials.jl +++ b/docs/src/literate_tutorials/mixed_materials.jl @@ -104,6 +104,8 @@ 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) @@ -112,9 +114,9 @@ function solve_nonlinear_timehistory(buffer, dh, ch, lh, l2_proj, qp_evaluator; 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 + ## The applied traction is linear in `t`, so if `fext` accumulated loads + ## 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) From 874e31e1517d1649e0d46af0110f0f05eb005b5a Mon Sep 17 00:00:00 2001 From: Knut Andreas Date: Thu, 24 Sep 2026 08:08:08 -0400 Subject: [PATCH 4/4] Hide test line --- docs/src/literate_tutorials/mixed_materials.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/literate_tutorials/mixed_materials.jl b/docs/src/literate_tutorials/mixed_materials.jl index 725d5364..a75e1b5e 100644 --- a/docs/src/literate_tutorials/mixed_materials.jl +++ b/docs/src/literate_tutorials/mixed_materials.jl @@ -114,7 +114,7 @@ function solve_nonlinear_timehistory(buffer, dh, ch, lh, l2_proj, qp_evaluator; apply!(a, ch) fill!(fext, 0) apply!(fext, lh, t) - ## The applied traction is linear in `t`, so if `fext` accumulated loads + ## 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