Initial and reference state

examples/cantilever/cantilever.jl

@@ -70,9 +70,9 @@ function test(sol::AbstractSolution)
     @testset "Displacements at the right-most node" begin
         Izz = b * h^3 / 12
         A = b * h
-        @test displacements(sol, right_most_node, 2)[1]≈-Py * L^3 / (3 * E * Izz) rtol=RTOL
-        @test displacements(sol, right_most_node, 6)[1]≈-Py * L^2 / (2 * E * Izz) rtol=RTOL
-        @test displacements(sol, right_most_node, 1)[1]≈Px * L / (E * A) rtol=RTOL
+        @test displacements(sol, right_most_node, 2)[end]≈-Py * L^3 / (3 * E * Izz) rtol=RTOL
+        @test displacements(sol, right_most_node, 6)[end]≈-Py * L^2 / (2 * E * Izz) rtol=RTOL
+        @test displacements(sol, right_most_node, 1)[end]≈Px * L / (E * A) rtol=RTOL
     end
 end
 

examples/linear_extension/linear_extension.jl

@@ -139,7 +139,6 @@ end;
 function analytic_solution(sol::AbstractSolution, p1::Point{dim}, p2::Point{dim},
         e::AbstractElement{dim}) where {dim}
     (; E, ν, λ, G, tension_load) = parameters()
-    sa = analysis(sol)
     ## Displacements
     "Compute displacements νmeric solution ui, uj and uk for analytic validation."
     function u_ijk_analytic(λv::Vector{<:Real},
@@ -154,12 +153,12 @@ function analytic_solution(sol::AbstractSolution, p1::Point{dim}, p2::Point{dim}
         [[ui(t) for t in λv], [uj(t) for t in λv], [uk(t) for t in λv]]
     end
     # point 1
-    u_1 = u_ijk_analytic(load_factors(sa), p1[1], p1[2], p1[3])
+    u_1 = u_ijk_analytic(load_factors(analysis(sol)), p1[1], p1[2], p1[3])
     ui_1 = u_1[1]
     uj_1 = u_1[2]
     uk_1 = u_1[3]
     # point 2
-    u_2 = u_ijk_analytic(load_factors(sa), p2[1], p2[2], p2[3])
+    u_2 = u_ijk_analytic(load_factors(analysis(sol)), p2[1], p2[2], p2[3])
     ui_2 = u_2[1]
     uj_2 = u_2[2]
     uk_2 = u_2[3]
@@ -194,7 +193,7 @@ function analytic_solution(sol::AbstractSolution, p1::Point{dim}, p2::Point{dim}
     # point in the rand element selected
     p = rand(coordinates(e))
     # strain
-    λv = load_factors(sa)
+    λv = load_factors(analysis(sol))
     ϵ = ϵ_ijk_analytic(λv, p[1], p[2], p[3])
     ϵ_1 = ϵ[1]
     ϵ_2 = ϵ[2]

src/Interfaces/VTK.jl

@@ -199,7 +199,6 @@ function _extract_cell_data(sol::AbstractSolution, fields::Vector{Field},
     end
     cell_data
 end
-
 """
 Generate a VTK file from a solution structure, exporting simulation data such as displacements, strain, and stress.
 
@@ -246,6 +245,7 @@ end
 """
 Generate a time-series of VTK files for a solution struct and organize them into a ParaView collection.
 Each time step's VTK file is added to the collection for seamless visualization of the simulation over time.
+`states(sol)[1]` (initial state) is written at time 0.0; subsequent states at their respective times.
 Returns the path to the generated ParaView collection file.
 """
 function write_vtk(sol::AbstractSolution, base_filename::String;
@@ -254,16 +254,13 @@ function write_vtk(sol::AbstractSolution, base_filename::String;
     times_vector = times(analysis(sol))  # Extract the times from the solution analysis
     pvd = paraview_collection(base_filename)
 
-    # Loop over each time index and generate VTK files
-    for (time_index, time) in enumerate(times_vector)
-        vtk_filename = "$(base_filename)_timestep_$time_index.vtu"
-        vtk = write_vtk(sol, vtk_filename, time_index; fields, append)
+    for (step, time) in enumerate(times_vector)
+        vtk_filename = "$(base_filename)_timestep_$(step - 1).vtu"
+        vtk = write_vtk(sol, vtk_filename, step; fields, append)
         pvd[time] = vtk
     end
 
-    # Close the ParaView collection file
     close(pvd)
-
     @info "VTK file collection written to $base_filename.pvd"
     "$base_filename.pvd"
 end

src/StructuralAnalyses/LinearStaticAnalyses.jl

@@ -52,18 +52,19 @@ function LinearStaticAnalysis(s::S, λᵥ::LFV;
             LinearResidualsIterationStep),
         initial_step::Int = 1) where {S <: AbstractStructure,
         LFV <: Vector{<:Real}}
-    !(1 ≤ initial_step ≤ length(λᵥ)) &&
-        throw(ArgumentError("initial_step must be in [1, $(length(λᵥ))] but is: $initial_step."))
-    LinearStaticAnalysis(s, initial_state, λᵥ, initial_step)
+    λᵥ_full = iszero(first(λᵥ)) ? λᵥ : [zero(eltype(λᵥ)); λᵥ]
+    nsteps = length(λᵥ_full) - 1
+    !(1 ≤ initial_step ≤ nsteps) &&
+        throw(ArgumentError("initial_step must be in [1, $nsteps] but is: $initial_step."))
+    LinearStaticAnalysis(s, initial_state, λᵥ_full, initial_step)
 end
 
 "Constructor for linear analysis given a final time (or load factor) and the number of steps."
 function LinearStaticAnalysis(s::AbstractStructure, final_time::Real = 1.0; NSTEPS = 10,
         initial_state::FullStaticState = FullStaticState(s,
             LinearResidualsIterationStep),
         initial_step::Int = 1)
-    t₀ = final_time / NSTEPS
-    λᵥ = collect(LinRange(t₀, final_time, NSTEPS))
+    λᵥ = collect(LinRange(zero(final_time), final_time, NSTEPS + 1))
     LinearStaticAnalysis(s, λᵥ; initial_state, initial_step)
 end
 

src/StructuralAnalyses/NonLinearStaticAnalyses.jl

@@ -50,16 +50,17 @@ function NonLinearStaticAnalysis(s::S, λᵥ::LFV;
         initial_state::FullStaticState = FullStaticState(s),
         initial_step::Int = 1) where {S <: AbstractStructure,
         LFV <: AbstractVector{<:Real}}
-    !(1 ≤ initial_step ≤ length(λᵥ)) &&
-        throw(ArgumentError("initial_step must be in [1, $(length(λᵥ))] but is: $initial_step."))
-    NonLinearStaticAnalysis(s, initial_state, λᵥ, initial_step)
+    λᵥ_full = iszero(first(λᵥ)) ? λᵥ : [zero(eltype(λᵥ)); λᵥ]
+    nsteps = length(λᵥ_full) - 1
+    !(1 ≤ initial_step ≤ nsteps) &&
+        throw(ArgumentError("initial_step must be in [1, $nsteps] but is: $initial_step."))
+    NonLinearStaticAnalysis(s, initial_state, λᵥ_full, initial_step)
 end
 
 "Constructor for non linear static analysis given a final time (or load factor) and the number of steps."
 function NonLinearStaticAnalysis(
         s::AbstractStructure, t₁::Real = 1.0; NSTEPS = 10, initial_step::Int = 1)
-    t₀ = t₁ / NSTEPS
-    λᵥ = collect(LinRange(t₀, t₁, NSTEPS))
+    λᵥ = collect(LinRange(zero(t₁), t₁, NSTEPS + 1))
     NonLinearStaticAnalysis(s, λᵥ; initial_step)
 end
 

src/StructuralAnalyses/StaticAnalyses.jl

@@ -57,29 +57,29 @@ and extends the following methods:
 """
 abstract type AbstractStaticAnalysis <: AbstractStructuralAnalysis end
 
-"Return the initial load factor of an structural analysis."
-initial_time(sa::AbstractStaticAnalysis) = first(load_factors(sa))
+"Return the initial load factor of an structural analysis (always 0)."
+initial_time(sa::AbstractStaticAnalysis) = first(sa.λᵥ)
+
+"Return load factors for all N+1 states, starting at 0."
+load_factors(sa::AbstractStaticAnalysis) = sa.λᵥ
+times(sa::AbstractStaticAnalysis) = load_factors(sa)
 
 "Return the current load factor of an structural analysis."
-current_time(sa::AbstractStaticAnalysis) = load_factors(sa)[sa.current_step[]]
+current_time(sa::AbstractStaticAnalysis) = sa.λᵥ[sa.current_step[] + 1]
 
 "Return the final load factor of an structural analysis."
-final_time(sa::AbstractStaticAnalysis) = last(load_factors(sa))
+final_time(sa::AbstractStaticAnalysis) = last(sa.λᵥ)
 
 "Return true if the structural analysis is completed."
 function is_done(sa::AbstractStaticAnalysis)
-    is_done_bool = if sa.current_step[] > length(load_factors(sa))
+    if sa.current_step[] > length(sa.λᵥ) - 1
         sa.current_step -= 1
         true
     else
         false
     end
 end
 
-"Return the final load factor vector of an structural analysis."
-load_factors(sa::AbstractStaticAnalysis) = sa.λᵥ
-times(sa::AbstractStaticAnalysis) = load_factors(sa)
-
 "Return the current load factor of an structural analysis."
 current_load_factor(sa::AbstractStaticAnalysis) = current_time(sa)
 
@@ -155,7 +155,8 @@ function Base.push!(st_sol::Solution{<:FullStaticState}, c_state::FullStaticStat
 end
 
 function store!(sol::Solution{<:StaticState}, state::FullStaticState, step::Int)
-    solution_state = states(sol)[step]
+    # sol.states[1] is the initial state; load step k is stored at sol.states[k+1].
+    solution_state = states(sol)[step + 1]
     sol_Uᵏ = displacements(solution_state)
     sol_σᵏ = stress(solution_state)
     sol_ϵᵏ = strain(solution_state)

src/StructuralSolvers/Solutions.jl

@@ -27,7 +27,7 @@ using ..StructuralSolvers
                                       iterations
 @reexport import ..Entities: internal_forces, inertial_forces, strain, stress
 
-export AbstractSolution, Solution, analysis, solver, states,
+export AbstractSolution, Solution, analysis, solver, states, reference_state,
        displacements, external_forces, iteration_residuals, deformed_node_positions
 
 """
@@ -53,20 +53,34 @@ analysis(sol::AbstractSolution) = sol.analysis
 solver(sol::AbstractSolution) = sol.solver
 
 """
-Solution that stores all intermediate arrays during the analysis.
+Solution that stores all states of the analysis.
+
+`states[1]` is the initial (unloaded) state at time 0, capturing any applied boundary conditions
+before the first load step. `states[k+1]` is the converged state at load step `k`.
+This convention matches standard FEA solvers (FEBio, Abaqus, OpenSees): the reference state
+is always at t = 0, and each subsequent entry corresponds to a load/time increment.
 """
 struct Solution{ST <: AbstractStaticState,
     A <: AbstractStructuralAnalysis,
     SS <: Union{AbstractSolver, LinearSolver}} <: AbstractSolution
-    "Vector containing the converged structural states at each step."
+    """
+    Purely undeformed reference configuration (zero displacements, zero stress, zero strain).
+    Corresponds to the continuum mechanics reference configuration — independent of any applied BCs.
+    """
+    reference_state::ST
+    """
+    Vector of structural states. `states[1]` is the initial state at t = 0, capturing any applied
+    boundary conditions before the first load step. `states[k+1]` is the converged state at load
+    step `k`. Length is `N + 1` for `N` load steps.
+    """
     states::Vector{ST}
     "Analysis solved."
     analysis::A
     "Solver employed."
     solver::SS
 end
 
-"Constructor with empty `AbstractStructuralState`s `Vector` and type `S`."
+"Constructor pre-allocating states with an initial state at index 1."
 function Solution(analysis::A,
         solver::SS) where {A <: AbstractStructuralAnalysis,
         SS <: Union{AbstractSolver, LinearSolver}}
@@ -77,17 +91,33 @@ function Solution(analysis::A,
     ϵᵏ = strain(state)
     σᵏ = stress(state)
 
-    #TODO: Create a general way to pre-allocate with the number of analysis steps
-    states = Vector{StaticState}(undef, length(analysis.λᵥ))
-
-    for i in 1:length(analysis.λᵥ)
+    # λᵥ starts at 0; the N load steps are λᵥ[2:end].
+    nsteps = length(analysis.λᵥ) - 1
+    # states[1] = initial state; states[2..N+1] = one per load step.
+    states = Vector{StaticState}(undef, nsteps + 1)
+
+    # Reference configuration: purely undeformed (zero U, zero σ, zero ε).
+    ref_state = StaticState(
+        zero(Uᵏ),
+        dictionary([e => zero(ϵ) for (e, ϵ) in pairs(ϵᵏ)]),
+        dictionary([e => zero(σ) for (e, σ) in pairs(σᵏ)]))
+
+    # Capture the initial state from the analysis (reflects BCs and any initial conditions,
+    # e.g. prescribed displacements, pre-stress, thermal initial conditions — FEBio step 0).
+    states[1] = StaticState(
+        copy(Uᵏ),
+        dictionary([e => copy(ϵ) for (e, ϵ) in pairs(ϵᵏ)]),
+        dictionary([e => copy(σ) for (e, σ) in pairs(σᵏ)]))
+
+    # Pre-allocate the remaining states for each load step.
+    for i in 2:(nsteps + 1)
         sol_Uᵏ = similar(Uᵏ)
         sol_σᵏ = dictionary([e => similar(σ) for (e, σ) in pairs(σᵏ)])
         sol_ϵᵏ = dictionary([e => similar(ϵ) for (e, ϵ) in pairs(ϵᵏ)])
         states[i] = StaticState(sol_Uᵏ, sol_ϵᵏ, sol_σᵏ)
     end
 
-    Solution{StaticState, A, SS}(states, analysis, solver)
+    Solution{StaticState, A, SS}(ref_state, states, analysis, solver)
 end
 
 "Generic minimal show method for a generic `solution`"
@@ -107,14 +137,15 @@ function _print_table(solution::AbstractSolution)
     header = ["iter", "time", "||Uᵏ||", "||ΔUᵏ||/||Uᵏ||", "||ΔRᵏ||", "||ΔRᵏ||/||Fₑₓₜ||",
         "convergence criterion", "iterations"]
 
-    ΔU_rel = getindex.(displacement_tol(solution), 1)
-    ΔU = getindex.(displacement_tol(solution), 2)
-    ΔR_rel = getindex.(residual_forces_tol(solution), 1)
-    ΔR = getindex.(residual_forces_tol(solution), 2)
+    # states[1] is the initial state (not a load step); skip it so lengths match λᵥ.
+    ΔU_rel = getindex.(displacement_tol(solution), 1)[2:end]
+    ΔU = getindex.(displacement_tol(solution), 2)[2:end]
+    ΔR_rel = getindex.(residual_forces_tol(solution), 1)[2:end]
+    ΔR = getindex.(residual_forces_tol(solution), 2)[2:end]
     t = analysis(solution).λᵥ
     criterions = [string(criterion_step)[1:(end - 2)]
-                  for criterion_step in criterion(solution)]
-    iters = iterations(solution)
+                  for criterion_step in criterion(solution)][2:end]
+    iters = iterations(solution)[2:end]
     num_times = length(t)
     data = hcat(collect(1:num_times), t, ΔU, ΔU_rel, ΔR, ΔR_rel, criterions, iters)
 
@@ -133,7 +164,10 @@ function _print_table(solution::AbstractSolution)
         alignment = :c)
 end
 
-"Return the solved states."
+"Return the reference (undeformed) configuration: zero displacements, zero stress, zero strain."
+reference_state(sol::Solution) = sol.reference_state
+
+"Return all states: states[1] is the initial state, states[k+1] is load step k."
 states(sol::Solution) = sol.states
 
 for f in [:displacements, :internal_forces, :external_forces]

src/StructuralSolvers/SolversConfig.jl

@@ -44,9 +44,9 @@ end
 
 "Show convergence settings."
 function Base.show(io::IO, cs::ConvergenceSettings)
-    println("• ||ΔUᵏ||/||Uᵏ||| ≤ : $(residual_forces_tol(cs))")
-    println("• ||ΔRᵏ||/||Fₑₓₜ||| $(displacement_tol(cs))")
-    println("• iter k ≤: $(max_iter_tol(cs))")
+    println(io, "• ||ΔUᵏ||/||Uᵏ||  ≤ : $(displacement_tol(cs))")
+    println(io, "• ||ΔRᵏ||/||Fₑₓₜ|| ≤ : $(residual_forces_tol(cs))")
+    println(io, "• iter k            ≤ : $(max_iter_tol(cs))")
 end
 
 "Return residual forces tolerance set."
@@ -131,7 +131,7 @@ end
 "Sets the iteration step with nothing."
 function reset!(ri_step::ResidualsIterationStep{<:Nothing})
     ri_step.iter = 0
-    ri_step.criterion = ResidualForceCriterion()
+    ri_step.criterion = NotConvergedYet()
     ri_step
 end