Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix show method avoiding stats print #503

Merged
merged 1 commit into from
Jul 7, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
46 changes: 21 additions & 25 deletions examples/linear_extension/linear_extension.jl
Original file line number Diff line number Diff line change
Expand Up @@ -9,25 +9,25 @@ include("linear_extension_mesh.jl")

"Return problem parameters"
function parameters()
E = 2.0 # Young modulus in Pa
ν = 0.4 # Poisson's ratio
E = 2.0 # Young modulus in Pa
ν = 0.4 # Poisson's ratio
λ = E * ν / ((1 + ν) * (1 - 2 * ν)) # First Lamé parameter
G = E / (2 * (1 + ν)) # Second Lamé parameter
tension(t) = p * t # Tension load function
p = 3 # Tension load in Pa
Lx = 2.0 # Dimension in x of the box in m
Ly = 1.0 # Dimension in y of the box in m
Lz = 1.0 # Dimension in z of the box in m
RTOL = 1e-4 # Relative tolerance for tests
ATOL = 1e-6 # Absolute tolerance for tests
NSTEPS = 9 # νmber of steps for the test
ms = 0.5 # Refinement factor
(; Lx, Ly, Lz, E, ν, λ, G, tension, RTOL, ATOL, NSTEPS, ms)
G = E / (2 * (1 + ν)) # Second Lamé parameter
tension_load(t) = p * t # Tension load function
p = 3 # Tension load in Pa
Lx = 2.0 # Dimension in x of the box in m
Ly = 1.0 # Dimension in y of the box in m
Lz = 1.0 # Dimension in z of the box in m
RTOL = 1e-4 # Relative tolerance for tests
ATOL = 1e-6 # Absolute tolerance for tests
NSTEPS = 9 # number of steps for the test
ms = 0.5 # Refinement factor
(; Lx, Ly, Lz, E, ν, λ, G, tension_load, RTOL, ATOL, NSTEPS, ms)
end;

"Return the problem structural model"
function structure()
(; Lx, Ly, Lz, E, ν, tension, ms) = parameters()
(; Lx, Ly, Lz, E, ν, tension_load, ms) = parameters()
# -------------------------------
# Materials
# -------------------------------
Expand All @@ -46,7 +46,7 @@ function structure()
bc3 = FixedDof(:u, [3], bc3_label)
# Load
bc4_label = "tension"
bc4 = GlobalLoad(:u, t -> [tension(t), 0, 0], bc4_label)
bc4 = GlobalLoad(:u, t -> [tension_load(t), 0, 0], bc4_label)
# Get bc labels for the mesh
bc_labels = [bc1_label, bc2_label, bc3_label, bc4_label]
boundary_conditions = StructuralBoundaryCondition(bc1, bc2, bc3, bc4)
Expand All @@ -70,13 +70,9 @@ function structure()
end
gmsh_println(output)
msh_file = MshFile(mesh_path)
mesh = Mesh(msh_file, entities)
# -------------------------------
# Structure
# -------------------------------
apply!(materials, mesh)
apply!(boundary_conditions, mesh)

Structure(msh_file, materials, boundary_conditions, entities)
end;

Expand All @@ -91,7 +87,7 @@ function solve()
# -------------------------------
# Numerical solution
# -------------------------------
solve!(sa)
sol = solve!(sa)
end;

"Return random points to evaluate the solution"
Expand Down Expand Up @@ -143,14 +139,14 @@ end;
"Return analytical results for testing"
function analytic_solution(sol::AbstractSolution, p1::Point{dim}, p2::Point{dim},
e::AbstractElement{dim}) where {dim}
(; E, ν, λ, G, tension) = parameters()
(; 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},
x0::Real, y0::Real, z0::Real,
ν::Real=ν, E::Real=E)
C(t) = tension(t) * (1 - ν - 2ν^2) / (1 - ν)
C(t) = tension_load(t) * (1 - ν - 2ν^2) / (1 - ν)

ui(t) = C(t) / E * x0
uj(t) = 0.0
Expand All @@ -172,7 +168,7 @@ function analytic_solution(sol::AbstractSolution, p1::Point{dim}, p2::Point{dim}
"Compute strains νmeric solution ϵi, ϵj and ϵk for analytic validation."
function ϵ_ijk_analytic(λv::Vector{<:Real}, x0::Real, y0::Real, z0::Real, ν::Real=ν,
E::Real=E)
C(t) = tension(t) * (1 - ν - 2ν^2) / (1 - ν)
C(t) = tension_load(t) * (1 - ν - 2ν^2) / (1 - ν)

ϵi(t) = C(t) / E
ϵj(t) = 0.0
Expand All @@ -183,7 +179,7 @@ function analytic_solution(sol::AbstractSolution, p1::Point{dim}, p2::Point{dim}
## Stresses
"Compute strains νmeric solution ϵi, ϵj and ϵk for analytic validation."
function σ_ijk_analytic(λv::Vector{<:Real}, x0::Real, y0::Real, z0::Real, λ::Real, G::Real)
C(t) = tension(t) * (1 - ν - 2ν^2) / (1 - ν)
C(t) = tension_load(t) * (1 - ν - 2ν^2) / (1 - ν)

ϵi(t) = C(t) / E
ϵj(t) = 0.0
Expand Down Expand Up @@ -260,7 +256,7 @@ function test(sol::AbstractSolution)
@test σj_num ≈ σj_analy rtol = RTOL skip = true
@test σk_num ≈ σk_analy rtol = RTOL skip = true
end
end
end;

"Run the example"
function run()
Expand Down
6 changes: 3 additions & 3 deletions src/StructuralAnalyses/LinearStaticAnalyses.jl
Original file line number Diff line number Diff line change
Expand Up @@ -73,12 +73,12 @@ function Base.show(io::IO, sa::LinearStaticAnalysis)
end

"Solves a linear analysis problem mutating the state."
function _solve!(sa::LinearStaticAnalysis, alg::Nothing,
function _solve!(sa::LinearStaticAnalysis, ::Nothing,
linear_solver::SciMLBase.AbstractLinearAlgorithm)
s = structure(sa)

# Initialize solution.
solution = Solution(sa, alg)
solution = Solution(sa, linear_solver)

# Load factors iteration.
while !is_done(sa)
Expand All @@ -97,7 +97,7 @@ function _solve!(sa::LinearStaticAnalysis, alg::Nothing,
@debugtime "Step" step!(sa, linear_solver)

# Recompute σ and ε for the assembler
@debugtime "Update internal forces, stres and strains" assemble!(s, sa)
@debugtime "Update internal forces, stresses and strains" assemble!(s, sa)

# Save current state
store!(solution, current_state(sa), step)
Expand Down
19 changes: 19 additions & 0 deletions src/StructuralAnalyses/NonLinearStaticAnalyses.jl
Original file line number Diff line number Diff line change
Expand Up @@ -136,4 +136,23 @@ function step!(sa::NonLinearStaticAnalysis, ::NewtonRaphson,
update!(current_iteration(sa), norm_ΔU, rel_norm_ΔU, norm_r, rel_norm_r)
end

"Show the solution when solved with an in-house algorithm."
function Base.show(io::IO, ::MIME"text/plain",
solution::Solution{<:FullStaticState})
show(io, solution)

println("\nStats:")
println("----------")
# Check convergence
is_any_step_not_converged = any([criterion_step isa Union{NotConvergedYet,MaxIterCriterion}
for criterion_step in criterion(solution)])

num_iterations = reduce(+, iterations(solution))
avg_iterations = round(num_iterations / length(states(solution)); digits=1)
println("• Number of linear systems solved: $num_iterations")
println("• Average of iterations per step : $avg_iterations")
println("• Convergence success : $(!is_any_step_not_converged)")
_print_table(solution)
end

end # module
24 changes: 6 additions & 18 deletions src/StructuralSolvers/Solutions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ Of course the solution contains on the analysis and solver used to solve the pro
"""
module Solutions

using Reexport, PrettyTables, Dictionaries
using Reexport, PrettyTables, Dictionaries, LinearSolve

using ..Utils
using ..Entities
Expand Down Expand Up @@ -55,7 +55,7 @@ solver(sol::AbstractSolution) = sol.solver
Solution that stores all intermediate arrays during the analysis.
"""
struct Solution{ST<:AbstractStaticState,A<:AbstractStructuralAnalysis,
SS<:Union{AbstractSolver,Nothing}} <: AbstractSolution
SS<:Union{AbstractSolver,SciMLBase.AbstractLinearAlgorithm}} <: AbstractSolution
mvanzulli marked this conversation as resolved.
Show resolved Hide resolved
"Vector containing the converged structural states at each step."
states::Vector{ST}
"Analysis solved."
Expand All @@ -67,7 +67,7 @@ end
"Constructor with empty `AbstractStructuralState`s `Vector` and type `S`."
function Solution(analysis::A,
solver::SS) where {A<:AbstractStructuralAnalysis,
SS<:Union{AbstractSolver,Nothing}}
SS<:Union{AbstractSolver,SciMLBase.AbstractLinearAlgorithm}}

# TODO Use concrete types.
state = current_state(analysis)
Expand All @@ -88,28 +88,16 @@ function Solution(analysis::A,
Solution{StaticState,A,SS}(states, analysis, solver)
end

"Show the states solution."
function Base.show(io::IO, ::MIME"text/plain", solution::Solution)
"Generic minimal show method for a generic `solution`"
function Base.show(io::IO, ::MIME"text/plain",
solution::Solution)
println("Analysis solved:")
println("----------------\n")
show(io, analysis(solution))

println("\nSolver employed:")
println("----------------\n")
show(io, solver(solution))

println("\nStats:")
println("----------")
# Check convergence
is_any_step_not_converged = any([criterion_step isa Union{NotConvergedYet,MaxIterCriterion}
for criterion_step in criterion(solution)])

num_iterations = reduce(+, iterations(solution))
avg_iterations = round(num_iterations / length(states(solution)); digits=1)
println("• Number of linear systems solved: $num_iterations")
println("• Average of iterations per step : $avg_iterations")
println("• Convergence success : $(!is_any_step_not_converged)")
_print_table(solution)
end

"Print the solution table"
Expand Down
11 changes: 7 additions & 4 deletions src/StructuralSolvers/StructuralSolvers.jl
Original file line number Diff line number Diff line change
Expand Up @@ -63,8 +63,10 @@ struct ResidualForceCriterion <: AbstractConvergenceCriterion end
""" `ΔUCriterion` indicates displacements increment convergence criterion. """
struct ΔUCriterion <: AbstractConvergenceCriterion end

""" `ΔU_and_ResidualForce_Criteria` convergence criterion indicates that both ΔU and residual forces
converged . """
"""
`ΔU_and_ResidualForce_Criteria` convergence criterion indicates that both
ΔU and residual forces converged .
"""
struct ΔU_and_ResidualForce_Criteria <: AbstractConvergenceCriterion end

""" `MaxIterCriterion` criteria indicates that the maximum number of iterations has been reached. """
Expand Down Expand Up @@ -198,8 +200,9 @@ which holds the result and the algorithm used to obtain it.

This function mutates the state defined in the analysis problem.
Use [`solve`](@ref) to avoid mutation.
For linear analysis problems, the algorithm doesn't need to be provided.
Also a linar solver form the LinearSolve.jl package can by provided. By default the wrapper to IterativeSolvers Conjugatge Gradient is being used.
For linear analysis problems, the algorithm doesn't need to be provided since `linear_solve` is used.
Also a linear solver form the LinearSolve.jl package can by provided. By default the wrapper
to IterativeSolvers Conjugatge Gradient is being used.
"""
function solve!(problem::AbstractStructuralAnalysis, solver::Union{AbstractSolver,Nothing}=nothing,
linear_solve::SciMLBase.AbstractLinearAlgorithm=DEFAULT_LINEAR_SOLVER())
Expand Down
Loading