From b36da58d90b6a1f18db888cf91d599cdd090aef5 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Sat, 9 May 2026 17:23:24 +0300 Subject: [PATCH] chore(legacy): remove orphaned elasticity Problem module Delete the Dict-based continuum elasticity `Problem` implementation that was never wired into `Legacy.jl`. - Drop `src/legacy/problems_elasticity.jl` (assembly over legacy elements). --- src/legacy/problems_elasticity.jl | 594 ------------------------------ 1 file changed, 594 deletions(-) delete mode 100644 src/legacy/problems_elasticity.jl diff --git a/src/legacy/problems_elasticity.jl b/src/legacy/problems_elasticity.jl deleted file mode 100644 index 74ab676..0000000 --- a/src/legacy/problems_elasticity.jl +++ /dev/null @@ -1,594 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - -""" Elasticity equations. - -Field equation is: - - m∂²u/∂t² = ∇⋅σ - b - -Weak form is: find u∈U such that ∀v in V - - δW := ∫ρ₀∂²u/∂t²⋅δu dV₀ + ∫S:δE dV₀ - ∫b₀⋅δu dV₀ - ∫t₀⋅δu dA₀ = 0 - -where - - ρ₀ = density - b₀ = displacement load - t₀ = displacement traction - -Formulations ------------- -plane stress, plane strain, 3D - -References ----------- - -https://en.wikipedia.org/wiki/Linear_elasticity -https://en.wikipedia.org/wiki/Finite_strain_theory -https://en.wikipedia.org/wiki/Stress_measures -https://en.wikipedia.org/wiki/Mooney%E2%80%93Rivlin_solid -https://en.wikipedia.org/wiki/Strain_energy_density_function -https://en.wikipedia.org/wiki/Plane_stress -https://en.wikipedia.org/wiki/Hooke's_law - -""" -mutable struct Elasticity <: FieldProblem - # these are found from problem.properties for type Problem{Elasticity} - formulation::Symbol - finite_strain::Bool - geometric_stiffness::Bool - store_fields::Vector{Symbol} -end -function Elasticity() - # formulations: plane_stress, plane_strain, continuum - return Elasticity(:continuum, false, false, []) -end - -function get_unknown_field_name(problem::Problem{Elasticity}) - return "displacement" -end - -function get_formulation_type(problem::Problem{Elasticity}) - return :incremental -end - -""" - assemble!(assembly:Assembly, problem::Problem{Elasticity}, elements, time) - -Start finite element assembly procedure for Elasticity problem. - -Function groups elements to arrays by their type and assembles one element type -at time. This makes it possible to pre-allocate matrices common to same type -of elements. -""" -function assemble!(assembly::Assembly, problem::Problem{Elasticity}, - elements::Vector{Element}, time) - formulation = Val{problem.properties.formulation} - for (element_type, elements_subset) in group_by_element_type(elements) - assemble!(assembly, problem, elements_subset, time, formulation) - end -end - -function assemble!(assembly::Assembly, problem::Problem{Elasticity}, - elements::Vector{T}, time, formulation) where {T<:Element} - - # Normal assembly (parallel assembly disabled for now - needs property fields) - local_buffer = allocate_buffer(problem, elements) - for i in 1:length(elements) - assemble_element!(assembly, problem, elements[i], local_buffer, time, formulation) - end -end - -include("problems_elasticity_2d.jl") - -const Elasticity3DSurfaceElements = Union{Poi1,Tri3,Tri6,Quad4,Quad8,Quad9} -const Elasticity3DVolumeElements = Union{Tet4,Pyr5,Wedge6,Wedge15,Hex8,Tet10,Hex20,Hex27} - -function initialize_internal_params!(params, ip, type_) #::Type{Val{:type_2d}}) - param_keys = keys(params) - all_keys = ip.fields.keys - ip_fields = filter(x -> isassigned(all_keys, x), collect(1:length(all_keys))) - - if !("params_initialized" in ip_fields) - for key in param_keys - update!(ip, key, 0.0 => params[key]) - end - if type_ == Val{:type_2d} - update!(ip, "stress", 0.0 => [0.0, 0.0, 0.0]) - update!(ip, "strain", 0.0 => [0.0, 0.0, 0.0]) - elseif type_ == Val{:type_3d} - update!(ip, "stress", 0.0 => [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]) - update!(ip, "strain", 0.0 => [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]) - else - error("daa") - end - update!(ip, "prev_time", 0.0 => 0.0) - update!(ip, "params_initialized", 0.0 => true) - end -end - -# Buffer for elasticity assembly (keyword constructor without Parameters.jl) -struct Elasticity3DLocalBuffers{B,T} - ndofs::Int - dim::Int - bi::BasisInfo{B,T} - BL::Matrix{T} - BNL::Matrix{T} - Km::Matrix{T} - Kg::Matrix{T} - f_int::Vector{T} - f_ext::Vector{T} - f_buffer::Vector{T} - f_buffer_dim::Vector{T} - gdofs::Vector{Int} - gradu::Matrix{T} - strain::Matrix{T} - strain_vec::Vector{T} - stress_vec::Vector{T} - F::Matrix{T} - D::Matrix{T} - Dtan::Matrix{T} - Bt_mul_D::Matrix{T} - Bt_mul_D_mul_B::Matrix{T} - Bt_mul_S::Vector{T} - - # Keyword constructor - function Elasticity3DLocalBuffers(; ndofs::Int, dim::Int, bi::BasisInfo{B,T}) where {B,T} - new{B,T}( - ndofs, dim, bi, - zeros(T, 6, ndofs), # BL - zeros(T, 9, ndofs), # BNL - zeros(T, ndofs, ndofs), # Km - zeros(T, ndofs, ndofs), # Kg - zeros(T, ndofs), # f_int - zeros(T, ndofs), # f_ext - zeros(T, ndofs), # f_buffer - zeros(T, div(ndofs, dim)), # f_buffer_dim - zeros(Int, ndofs), # gdofs - zeros(T, dim, dim), # gradu - zeros(T, dim, dim), # strain - zeros(T, 6), # strain_vec - zeros(T, 6), # stress_vec - zeros(T, dim, dim), # F - zeros(T, 6, 6), # D - zeros(T, 6, 6), # Dtan - zeros(T, ndofs, 6), # Bt_mul_D - zeros(T, ndofs, ndofs), # Bt_mul_D_mul_B - zeros(T, ndofs) # Bt_mul_S - ) - end -end - -function allocate_buffer(problem::Problem{Elasticity}, ::Vector{Element{El}}) where El<:Elasticity3DVolumeElements - dim = get_unknown_field_dimension(problem) - - nnodes = length(El) - ndofs = dim * nnodes - - return Elasticity3DLocalBuffers(ndofs=ndofs, dim=dim, bi=BasisInfo(El)) -end - -function reset_element!(buf::Elasticity3DLocalBuffers) - fill!(buf.Km, 0.0) - fill!(buf.Kg, 0.0) - fill!(buf.f_int, 0.0) - fill!(buf.f_ext, 0.0) - return -end - -function reset_integration_point!(buf::Elasticity3DLocalBuffers) - fill!(buf.F, 0.0) - fill!(buf.strain, 0.0) - fill!(buf.D, 0.0) - fill!(buf.BL, 0.0) - fill!(buf.BNL, 0.0) - return -end - -function to_voigt!(strain_vec, strain) - strain_vec[1] = strain[1, 1] - strain_vec[2] = strain[2, 2] - strain_vec[3] = strain[3, 3] - strain_vec[4] = 2.0 * strain[1, 2] - strain_vec[5] = 2.0 * strain[2, 3] - strain_vec[6] = 2.0 * strain[1, 3] - return -end - -const u = ([0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]) -const X = ([-93.7197, -93.7197, 150.883], [-91.657, -85.8251, 157.885], [-100.523, -88.8309, 157.883], [-91.6593, -88.8309, 157.883], [-92.6883, -89.7724, 154.384], [-96.0902, -87.328, 157.883], [-97.1216, -91.2753, 154.383], [-92.6895, -91.2753, 154.383], [-91.6581, -87.328, 157.883], [-96.0914, -88.8309, 157.883]) -const displacement_load_string = [string("displacement load ", i) for i in 1:3] -""" Assemble 3d continuum elements in general solid mechanics problem. """ -function assemble_element!(assembly::Assembly, - assembler::AssemblerSparsityPattern, - problem::Problem{Elasticity}, - element::Element{El}, - local_buffer::Elasticity3DLocalBuffers, - time, ::Type{Val{:continuum}}, - use_csc=false) where El<:Elasticity3DVolumeElements - cheating = false - props = problem.properties - dim = get_unknown_field_dimension(problem) - - nnodes = length(El) - ndofs = dim * nnodes - - # Unpack buffer fields (replaced Parameters.@unpack) - bi = local_buffer.bi - BL = local_buffer.BL - BNL = local_buffer.BNL - Km = local_buffer.Km - Kg = local_buffer.Kg - f_int = local_buffer.f_int - f_ext = local_buffer.f_ext - f_buffer = local_buffer.f_buffer - f_buffer_dim = local_buffer.f_buffer_dim - gdofs = local_buffer.gdofs - gradu = local_buffer.gradu - strain = local_buffer.strain - strain_vec = local_buffer.strain_vec - stress_vec = local_buffer.stress_vec - F = local_buffer.F - D = local_buffer.D - Dtan = local_buffer.Dtan - Bt_mul_D = local_buffer.Bt_mul_D - Bt_mul_D_mul_B = local_buffer.Bt_mul_D_mul_B - Bt_mul_S = local_buffer.Bt_mul_S - - if !cheating - u = element("displacement", time) - X = element("geometry", time) - end - reset_element!(local_buffer) - - for ip in get_integration_points(element) - reset_integration_point!(local_buffer) - eval_basis!(bi, X, ip) - w = ip.weight * bi.detJ - N = bi.N - dN = bi.grad # deriatives of basis functions w.r.t. X, i.e. ∂N/∂X - grad!(bi, gradu, u) # displacement gradient ∇u - - # calculate strain tensor and deformation gradient - #F[:,:] += I - for i in 1:dim - F[i, i] += 1.0 - end - if props.finite_strain - strain[:, :] = 1 / 2 * (gradu + gradu' + gradu' * gradu) - F[:, :] += gradu - else - strain[:, :] .= 1 / 2 .* (gradu .+ gradu') - end - - to_voigt!(strain_vec, strain) - - # material stiffness start - - if props.finite_strain - for i = 1:nnodes - BL[1, 3*(i-1)+1] = F[1, 1] * dN[1, i] - BL[1, 3*(i-1)+2] = F[2, 1] * dN[1, i] - BL[1, 3*(i-1)+3] = F[3, 1] * dN[1, i] - BL[2, 3*(i-1)+1] = F[1, 2] * dN[2, i] - BL[2, 3*(i-1)+2] = F[2, 2] * dN[2, i] - BL[2, 3*(i-1)+3] = F[3, 2] * dN[2, i] - BL[3, 3*(i-1)+1] = F[1, 3] * dN[3, i] - BL[3, 3*(i-1)+2] = F[2, 3] * dN[3, i] - BL[3, 3*(i-1)+3] = F[3, 3] * dN[3, i] - BL[4, 3*(i-1)+1] = F[1, 1] * dN[2, i] + F[1, 2] * dN[1, i] - BL[4, 3*(i-1)+2] = F[2, 1] * dN[2, i] + F[2, 2] * dN[1, i] - BL[4, 3*(i-1)+3] = F[3, 1] * dN[2, i] + F[3, 2] * dN[1, i] - BL[5, 3*(i-1)+1] = F[1, 2] * dN[3, i] + F[1, 3] * dN[2, i] - BL[5, 3*(i-1)+2] = F[2, 2] * dN[3, i] + F[2, 3] * dN[2, i] - BL[5, 3*(i-1)+3] = F[3, 2] * dN[3, i] + F[3, 3] * dN[2, i] - BL[6, 3*(i-1)+1] = F[1, 3] * dN[1, i] + F[1, 1] * dN[3, i] - BL[6, 3*(i-1)+2] = F[2, 3] * dN[1, i] + F[2, 1] * dN[3, i] - BL[6, 3*(i-1)+3] = F[3, 3] * dN[1, i] + F[3, 1] * dN[3, i] - end - else - for i = 1:nnodes - BL[1, 3*(i-1)+1] = dN[1, i] - BL[2, 3*(i-1)+2] = dN[2, i] - BL[3, 3*(i-1)+3] = dN[3, i] - BL[4, 3*(i-1)+1] = dN[2, i] - BL[4, 3*(i-1)+2] = dN[1, i] - BL[5, 3*(i-1)+2] = dN[3, i] - BL[5, 3*(i-1)+3] = dN[2, i] - BL[6, 3*(i-1)+1] = dN[3, i] - BL[6, 3*(i-1)+3] = dN[1, i] - end - end - - # calculate stress - - if cheating - E = 200e3 - nu = 0.3 - else - E = element("youngs modulus", ip, time)::Float64 - nu = element("poissons ratio", ip, time)::Float64 - end - la = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)) - mu = E / (2.0 * (1.0 + nu)) - D[1, 1] = D[2, 2] = D[3, 3] = 2 * mu + la - D[4, 4] = D[5, 5] = D[6, 6] = mu - D[1, 2] = D[2, 1] = D[2, 3] = D[3, 2] = D[1, 3] = D[3, 1] = la - - # determine material model - - - material_model = :linear_elasticity - if haskey(element, "plasticity") - material_model = :ideal_plasticity - end - - # calculate stress vector based on material model - - if material_model == :linear_elasticity - copyto!(Dtan, D) - mul!(stress_vec, Dtan, strain_vec) - end - - if material_model == :ideal_plasticity - plastic_def = element("plasticity")[ip.id] - - calculate_stress! = plastic_def["type"] - yield_surface_ = plastic_def["yield_surface"] - params = plastic_def["params"] - initialize_internal_params!(params, ip, Val{:type_3d}) - - if time == 0.0 - error("Given step time = $(time). Please select time > 0.0") - end - - t_last = ip("prev_time", time) - update!(ip, "prev_time", time => t_last) - dt = time - t_last - stress_last = ip("stress", t_last) - strain_last = ip("strain", t_last) - dstrain_vec = strain_vec - strain_last - fill!(stress_vec, 0.0) - fill!(Dtan, 0.0) - plastic_strain = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - calculate_stress!(stress_vec, stress_last, dstrain_vec, plastic_strain, D, params, Dtan, yield_surface_, time, dt, Val{:type_3d}) - - end - - :strain in props.store_fields && update!(ip, "strain", time => strain_vec) - :stress in props.store_fields && update!(ip, "stress", time => stress_vec) - :stress11 in props.store_fields && update!(ip, "stress11", time => stress_vec[1]) - :stress22 in props.store_fields && update!(ip, "stress22", time => stress_vec[2]) - :stress33 in props.store_fields && update!(ip, "stress33", time => stress_vec[3]) - :stress12 in props.store_fields && update!(ip, "stress12", time => stress_vec[4]) - :stress23 in props.store_fields && update!(ip, "stress23", time => stress_vec[5]) - :stress13 in props.store_fields && update!(ip, "stress13", time => stress_vec[6]) - :plastic_strain in props.store_fields && update!(ip, "plastic_strain", time => plastic_strain) - - #Km += w*BL'*Dtan*BL - mul!(Bt_mul_D, transpose(BL), Dtan) - mul!(Bt_mul_D_mul_B, Bt_mul_D, BL) - rmul!(Bt_mul_D_mul_B, w) - for i = 1:ndofs^2 - @inbounds Km[i] += Bt_mul_D_mul_B[i] - end - - # material stiffness end - if props.geometric_stiffness - # take geometric stiffness into account - - for i = 1:size(dN, 2) - BNL[1, 3*(i-1)+1] = dN[1, i] - BNL[2, 3*(i-1)+1] = dN[2, i] - BNL[3, 3*(i-1)+1] = dN[3, i] - BNL[4, 3*(i-1)+2] = dN[1, i] - BNL[5, 3*(i-1)+2] = dN[2, i] - BNL[6, 3*(i-1)+2] = dN[3, i] - BNL[7, 3*(i-1)+3] = dN[1, i] - BNL[8, 3*(i-1)+3] = dN[2, i] - BNL[9, 3*(i-1)+3] = dN[3, i] - end - - S3 = zeros(3 * dim, 3 * dim) - S3[1, 1] = stress_vec[1] - S3[2, 2] = stress_vec[2] - S3[3, 3] = stress_vec[3] - S3[1, 2] = S3[2, 1] = stress_vec[4] - S3[2, 3] = S3[3, 2] = stress_vec[5] - S3[1, 3] = S3[3, 1] = stress_vec[6] - S3[4:6, 4:6] = S3[7:9, 7:9] = S3[1:3, 1:3] - - Kg += w * BNL' * S3 * BNL - - end - - # internal load - mul!(Bt_mul_S, transpose(BL), stress_vec) - rmul!(Bt_mul_S, w) - f_int .+= Bt_mul_S - - # external load start - if haskey(element, "displacement load") - T = element("displacement load", ip, time)::Vector{Float64} - mul!(f_buffer, w, vec(T * N)) - f_ext .+= f_buffer - end - - for i = 1:dim - if haskey(element, displacement_load_string[i]) - b = element(displacement_load_string[i], ip, time)::Float64 - mul!(f_buffer_dim, w, N) - for (i, j) in enumerate(1:dim:length(f_ext)) - f_ext[j] = b * f_buffer_dim[i] - end - end - end - # external load end - end - - FEMBase.get_gdofs!(gdofs, problem, element) - - # Update f_ext in place to be f_ext - f_int - f_ext .-= f_int - - if use_csc - # add contributions to K, Kg, f - @inbounds assemble_local!(assembler, gdofs, Km, f_ext) - - if props.geometric_stiffness - @inbounds assemble_local_matrix!(assembler, gdofs, Kg) - end - else - add!(assembly.f, gdofs, f_ext) - add!(assembly.K, gdofs, gdofs, Km) - if props.geometric_stiffness - add!(assembly.Kg, gdofs, gdofs, Kg) - end - end - - return nothing -end - -""" Elasticity equations, surface traction for continuum formulation. """ -function assemble!(assembly::Assembly, - problem::Problem{Elasticity}, - elements::Vector{Element{El}}, - time, ::Type{Val{:continuum}}) where El<:Elasticity3DSurfaceElements - - props = problem.properties - dim = get_unknown_field_dimension(problem) - - for element in elements - nnodes = size(element, 2) - f = zeros(dim * nnodes) - - has_concentrated_forces = false - for ip in get_integration_points(element) - detJ = element(ip, time, Val{:detJ}) - w = ip.weight * detJ - N = element(ip, time) - if haskey(element, "displacement traction force") - T = element("displacement traction force", ip, time) - f += w * vec(T * N) - end - for i in 1:dim - if haskey(element, "displacement traction force $i") - T = element("displacement traction force $i", ip, time) - f[i:dim:end] += w * vec(T * N) - end - if haskey(element, "concentrated force $i") - has_concentrated_forces = true - T = element("concentrated force $i", ip, time) - f[i:dim:end] += w * vec(T * N) - end - end - if haskey(element, "surface pressure") - J = element(ip, time, Val{:Jacobian})' - n = cross(J[:, 1], J[:, 2]) - n /= norm(n) - # sign convention, positive pressure is towards surface - p = -element("surface pressure", ip, time) - f += w * p * vec(n * N) - end - end - if has_concentrated_forces - update!(element, "concentrated force", time => Any[f]) - end - - gdofs = get_gdofs(problem, element) - add!(assembly.f, gdofs, f) - end -end - -""" - assemble!(assembly, problem, elements, time, ::Type{Val{:continuum}}) - -Assemble all other elements for continuum elasticity problems. Basically, throw -an exception telling to filter invalid elements out from the element set. -""" -function assemble!(assembly::Assembly, - problem::Problem{Elasticity}, - elements::Vector{Element{El}}, - time, ::Type{Val{:continuum}}) where El - @info("It looks that you are trying to assemble elements of type $El to 3d continuum " * - "problem. However, they are not supported yet. To filter out elements from a " * - "element set, try `filter(element->!isa(element, Element{$El}), elements)`") - error("Tried to assemble unsupported elements of type $El to 3d continuum problem.") -end - -""" Return strain tensor. """ -function get_strain_tensor(problem, element, ip, time) - gradu = element("displacement", ip, time, Val{:Grad}) - eps = 0.5 * (gradu' + gradu) - return eps -end - -""" Return stress tensor. """ -function get_stress_tensor(problem, element, ip, time) - eps = get_strain_tensor(problem, element, ip, time) - E = element("youngs modulus", ip, time) - nu = element("poissons ratio", ip, time) - mu = E / (2.0 * (1.0 + nu)) - la = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)) - S = la * tr(eps) * I + 2.0 * mu * eps - return S -end - -""" Return stain vector in "ABAQUS" order 11, 22, 33, 12, 23, 13. """ -function get_strain_vector(problem, element, ip, time) - eps = get_strain_tensor(problem, element, ip, time) - return [eps[1, 1], eps[2, 2], eps[3, 3], eps[1, 2], eps[2, 3], eps[1, 3]] -end - -""" Return stress vector in "ABAQUS" order 11, 22, 33, 12, 23, 13. """ -function get_stress_vector(problem, element, ip, time) - S = get_stress_tensor(problem, element, ip, time) - return [S[1, 1], S[2, 2], S[3, 3], S[1, 2], S[2, 3], S[1, 3]] -end - -""" Make least squares fit for some field to nodes. """ -function lsq_fit(problem, elements, field, time) - A = SparseMatrixCOO() - b = SparseMatrixCOO() - volume = 0.0 - for element in elements - gdofs = get_connectivity(element) - for ip in get_integration_points(element) - detJ = element(ip, time, Val{:detJ}) - w = ip.weight * detJ - N = element(ip, time) - f = field(problem, element, ip, time) - add!(A, gdofs, gdofs, w * kron(N', N)) - for i = 1:length(f) - add!(b, gdofs, w * f[i] * N, i) - end - volume += w - end - end - A = sparse(A) - b = sparse(b) - A = 1 / 2 * (A + A') - - nz = get_nonzero_rows(A) - F = ldlt(A[nz, nz]) - - x = F \ b[nz, :] - - nodal_values = Dict(node_id => Vector(x[idx, :]) for (idx, node_id) in enumerate(nz)) - return nodal_values -end - - -""" Postprocessing, extrapolate strain to nodes using least-squares fit. """ -function postprocess!(problem::Problem{Elasticity}, time::Float64, ::Type{Val{:strain}}) - elements = get_elements(problem) - strain = lsq_fit(problem, elements, get_strain_vector, time) - update!(elements, "strain", time => strain) -end - -function postprocess!(problem::Problem{Elasticity}, time::Float64, ::Type{Val{:stress}}) - elements = get_elements(problem) - stress = lsq_fit(problem, elements, get_stress_vector, time) - update!(elements, "stress", time => stress) -end