diff --git a/.gitignore b/.gitignore index 281f6a83..63395e93 100644 --- a/.gitignore +++ b/.gitignore @@ -14,4 +14,10 @@ Manifest.toml */__pycache__ -examples/benchmark/iterations \ No newline at end of file +examples/benchmark/iterations + +*.ipynb +/*.jl +*.vtu +*.png +*.zip \ No newline at end of file diff --git a/Project.toml b/Project.toml index 9fad67a5..3a362040 100644 --- a/Project.toml +++ b/Project.toml @@ -9,6 +9,7 @@ ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Einsum = "b7d42ee7-0b51-5a75-98ca-779d3107e4c0" +FFMPEG = "c87230d0-a227-11e9-1b43-d7ebe4e7570a" Ferrite = "c061ca5d-56c9-439f-9c0e-210fe06d3992" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" @@ -28,6 +29,8 @@ NonconvexPercival = "4296f080-e499-44ba-a02c-ae40015c44d5" NonconvexSemidefinite = "9e21ff56-eb25-42d3-b86f-5b0612f555e7" Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Preconditioners = "af69fa37-3177-5a40-98ee-561f696e4fcd" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" @@ -36,7 +39,9 @@ Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c" +Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" diff --git a/src/FEA/FEA.jl b/src/FEA/FEA.jl index 36544f9c..6a8172e9 100644 --- a/src/FEA/FEA.jl +++ b/src/FEA/FEA.jl @@ -9,13 +9,16 @@ using Parameters: @unpack export AbstractFEASolver, AbstractDisplacementSolver, + AbstractHyperelasticSolver, DirectDisplacementSolver, PCGDisplacementSolver, StaticMatrixFreeDisplacementSolver, + HyperelasticDisplacementSolver, Displacement, Direct, CG, MatrixFree, + Hyperelastic, FEASolver, Assembly, DefaultCriteria, @@ -34,6 +37,7 @@ include("direct_displacement_solver.jl") include("assembly_cg_displacement_solvers.jl") include("matrix_free_cg_displacement_solvers.jl") include("matrix_free_apply_bcs.jl") +include("iterative_hyperelastic_solver.jl") include("simulate.jl") include("solvers_api.jl") diff --git a/src/FEA/iterative_hyperelastic_solver.jl b/src/FEA/iterative_hyperelastic_solver.jl new file mode 100644 index 00000000..9fa43ea8 --- /dev/null +++ b/src/FEA/iterative_hyperelastic_solver.jl @@ -0,0 +1,242 @@ +abstract type AbstractHyperelasticSolver <: AbstractFEASolver end + +mutable struct HyperelasticCompressibleDisplacementSolver{ + T, + dim, + TP1<:AbstractPenalty{T}, + TP2<:StiffnessTopOptProblem{dim,T}, + TG<:GlobalFEAInfo_hyperelastic{T}, + TE<:ElementFEAInfo_hyperelastic{dim,T}, + Tu<:AbstractVector{T}, + TF<:AbstractVector{<:AbstractMatrix{T}}, +} <: AbstractHyperelasticSolver + mp + problem::TP2 + globalinfo::TG + elementinfo::TE + u::Tu # JGB: u --> u0 + F::TF + vars::Tu + penalty::TP1 + prev_penalty::TP1 + xmin::T + tsteps::Int + ntsteps::Int +end +mutable struct HyperelasticNearlyIncompressibleDisplacementSolver{ + T, + dim, + TP1<:AbstractPenalty{T}, + TP2<:StiffnessTopOptProblem{dim,T}, + TG<:GlobalFEAInfo_hyperelastic{T}, + TE<:ElementFEAInfo_hyperelastic{dim,T}, + Tu<:AbstractVector{T}, + TF<:AbstractVector{<:AbstractMatrix{T}}, +} <: AbstractHyperelasticSolver + mp + problem::TP2 + globalinfo::TG + elementinfo::TE + u::Tu # JGB: u --> u0 + F::TF + vars::Tu + penalty::TP1 + prev_penalty::TP1 + xmin::T + tsteps::Int + ntsteps::Int +end +function Base.show(::IO, ::MIME{Symbol("text/plain")}, x::HyperelasticCompressibleDisplacementSolver) + return println("TopOpt compressible hyperelastic solver") +end +function Base.show(::IO, ::MIME{Symbol("text/plain")}, x::HyperelasticNearlyIncompressibleDisplacementSolver) + return println("TopOpt nearly-incompressible hyperelastic solver") +end +function HyperelasticDisplacementSolver( + mp, # JGB: add type later + sp::StiffnessTopOptProblem{dim,T}; # JGB: eventually add ::HyperelaticParam type + xmin=T(1)/1000, + penalty=PowerPenalty{T}(1), + prev_penalty=deepcopy(penalty), + quad_order=default_quad_order(sp), + tstep = 1, + ntsteps = 20, + nearlyincompressible=false +) where {dim,T} + u = zeros(T, ndofs(sp.ch.dh)) + ts0 = tstep/ntsteps + update!(sp.ch,ts0) # set initial time-step (adjusts dirichlet bcs) + apply!(u,sp.ch) # apply dbc for initial guess + elementinfo = ElementFEAInfo_hyperelastic(mp, sp, u, quad_order, Val{:Static}, nearlyincompressible; ts = ts0) # JGB: add u + F = [zeros(typeof(elementinfo.Fes[1])) for _ in 1:getncells(sp.ch.dh.grid)] + globalinfo = GlobalFEAInfo_hyperelastic(sp) # JGB: small issue this leads to symmetric K initialization + #u = zeros(T, ndofs(sp.ch.dh)) # JGB + vars = fill(one(T), getncells(sp.ch.dh.grid) - sum(sp.black) - sum(sp.white)) + varind = sp.varind + return nearlyincompressible ? + HyperelasticNearlyIncompressibleDisplacementSolver(mp, sp, globalinfo, elementinfo, u, F, vars, penalty, prev_penalty, xmin, tstep, ntsteps) : + HyperelasticCompressibleDisplacementSolver(mp, sp, globalinfo, elementinfo, u, F, vars, penalty, prev_penalty, xmin, tstep, ntsteps) + #if nearlyincompressible + # return HyperelasticNearlyIncompressibleDisplacementSolver(mp, sp, globalinfo, elementinfo, u, F, vars, penalty, prev_penalty, xmin, tstep, ntsteps) + #else + # return HyperelasticCompressibleDisplacementSolver(mp, sp, globalinfo, elementinfo, u, F, vars, penalty, prev_penalty, xmin, tstep, ntsteps) + #end +end +function (s::HyperelasticCompressibleDisplacementSolver{T})( + ::Type{Val{safe}}=Val{false}, + ::Type{newT}=T; + assemble_f=true, + kwargs..., +) where {T,safe,newT} + elementinfo = s.elementinfo + globalinfo = s.globalinfo + dh = s.problem.ch.dh + ch = s.problem.ch + + _ndofs = ndofs(dh) + un = zeros(_ndofs) # previous solution vector + + NEWTON_TOL = 1e-8 + NEWTON_MAXITER = 30 + CG_MAXITER = 1000 + TS_MAXITER_ABS = 200 + TS_MAXITER_REL = 10 + + function HyperelasticSolverCore(ts) + update!(ch,ts) + apply!(un,ch) + #println(maximum(un)) + u = zeros(_ndofs) + Δu = zeros(_ndofs) + ΔΔu = zeros(_ndofs) + + newton_itr = 0 + normg = zeros(NEWTON_MAXITER) + while true; newton_itr += 1 + u .= un .+ Δu # current trial solution + # Compute residual norm for current guess + elementinfo = ElementFEAInfo_hyperelastic(s.mp, s.problem, u, default_quad_order(s.problem), Val{:Static}; ts=ts) + assemble_hyperelastic!(globalinfo,s.problem,elementinfo,s.vars,getpenalty(s),s.xmin,assemble_f=assemble_f) + apply_zero!(globalinfo.K,globalinfo.g,ch) + normg[newton_itr] = norm(globalinfo.g) + println("Tstep: $ts / 1]. Iteration: $newton_itr. normg is equal to " * string(normg[newton_itr])) + # Check for convergence + if normg[newton_itr] < NEWTON_TOL + break + elseif newton_itr > 1 && normg[newton_itr] > normg[newton_itr-1] + error("Newton iteration resulted in an increase in normg, aborting") + elseif newton_itr > NEWTON_MAXITER + error("Reached maximum Newton iterations, aborting") + end + # Compute increment using conjugate gradients + IterativeSolvers.cg!(ΔΔu, globalinfo.K, globalinfo.g; maxiter=CG_MAXITER) + apply_zero!(ΔΔu, ch) + Δu .-= ΔΔu + end + un = u + end + + inc_up = 1.1 + inc_down = 0.5 + delay_up = 5 + delay_down = 5 + inc_delay = 10 # integer + ntsteps = s.ntsteps + iter_ts = 0 + ts = 0 + Δts0 = 1/ntsteps + Δts = Δts0 + conv = zeros(TS_MAXITER_ABS) + #for tstep ∈ 1:ntsteps + # ts = tstep/ntsteps + #end + while ts < 1 + iter_ts += 1 + ts += Δts + if ts > 1 + Δts = 1 - ts + ts = 1 + elseif iter_ts > TS_MAXITER_REL && sum(conv[iter_ts-(TS_MAXITER_REL):iter_ts-1]) == 0 + error("Reached maximum number of successive failed ts iterations ($TS_MAXITER_REL), aborting") + elseif iter_ts > TS_MAXITER_ABS + error("Reached maximum number of allowed ts iterations ($TS_MAXITER_ABS), aborting") + end + try + HyperelasticSolverCore(ts) + conv[iter_ts] = 1 + if sum(conv[1:iter_ts]) == iter_ts || (iter_ts - findlast(x -> x == 0, conv[1:iter_ts-1])) > delay_up + 1 # increase Δts if it's never failed or if it's been 'inc_delay' or more successful iterations since last divergence + Δts *= inc_up + end + catch + conv[iter_ts] = 0 + ts -= Δts # revert to previous successful ts + println("REMINDER TO CHECK THIS!!") + if any(x -> x == 1, conv) && (iter_ts - findlast(x -> x == 1, conv[1:iter_ts-1])) < delay_down # decrease Δts a little if it's been 'down_delay' or less failures in a row + Δts *= 1/inc_up + else + Δts *= inc_down # otherwise make a big decrease in Δts + end + end + println("ts = $ts. Δts/Δts0 = $(Δts/Δts0)") + end + s.u .= un + s.F .= elementinfo.Fes + return nothing +end + +function (s::HyperelasticNearlyIncompressibleDisplacementSolver{T})( + ::Type{Val{safe}}=Val{false}, + ::Type{newT}=T; + assemble_f=true, + reuse_fact=false, + kwargs..., +) where {T,safe,newT} + elementinfo = s.elementinfo + globalinfo = s.globalinfo + dh = s.problem.ch.dh + ch = s.problem.ch + + _ndofs = ndofs(dh) + un = zeros(_ndofs) # previous solution vector + + NEWTON_TOL = 1e-8 + NEWTON_MAXITER = 30 + CG_MAXITER = 1000 + + ntsteps = s.ntsteps + for tstep ∈ 1:ntsteps + ts = tstep/ntsteps + update!(ch,ts) + apply!(un,ch) + println(maximum(un)) + u = zeros(_ndofs) + Δu = zeros(_ndofs) + ΔΔu = zeros(_ndofs) + + newton_itr = 0 + normg = zeros(NEWTON_MAXITER) + while true; newton_itr += 1 + u .= un .+ Δu # current trial solution + # Compute residual norm for current guess + elementinfo = ElementFEAInfo_hyperelastic(s.mp, s.problem, u, default_quad_order(s.problem), Val{:Static}; ts=ts) + assemble_hyperelastic!(globalinfo,s.problem,elementinfo,s.vars,getpenalty(s),s.xmin,assemble_f=assemble_f) + apply_zero!(globalinfo.K,globalinfo.g,ch) + normg[newton_itr] = norm(globalinfo.g) + println("Tstep: $tstep / $ntsteps. Iteration: $newton_itr. normg is equal to " * string(normg[newton_itr])) + # Check for convergence + if normg[newton_itr] < NEWTON_TOL + break + elseif newton_itr > NEWTON_MAXITER + error("Reached maximum Newton iterations, aborting") + end + # Compute increment using conjugate gradients + IterativeSolvers.cg!(ΔΔu, globalinfo.K, globalinfo.g; maxiter=CG_MAXITER) + apply_zero!(ΔΔu, ch) + Δu .-= ΔΔu + end + un = u + end + s.u .= un + s.F .= elementinfo.Fes + return nothing +end \ No newline at end of file diff --git a/src/FEA/solvers_api.jl b/src/FEA/solvers_api.jl index fe1b081e..4002030e 100644 --- a/src/FEA/solvers_api.jl +++ b/src/FEA/solvers_api.jl @@ -4,6 +4,7 @@ abstract type SolverSubtype end struct Direct <: SolverType end struct CG <: SolverType end +struct Hyperelastic <: SolverType end struct MatrixFree <: SolverSubtype end struct Assembly <: SolverSubtype end @@ -38,6 +39,10 @@ function FEASolver(::Type{CG}, ::Type{Assembly}, problem; kwargs...) return PCGDisplacementSolver(problem; kwargs...) end +function FEASolver(::Type{Hyperelastic}, mp, problem; kwargs...) # JGB: add place to add NeoHookean <: :HyperelasticParam + return HyperelasticDisplacementSolver(mp, problem; kwargs...) +end + function default_quad_order(problem) if TopOptProblems.getdim(problem) == 2 && TopOptProblems.nnodespercell(problem) in (3, 6) || diff --git a/src/Functions/Functions.jl b/src/Functions/Functions.jl index 6e27b817..f493c0fa 100644 --- a/src/Functions/Functions.jl +++ b/src/Functions/Functions.jl @@ -13,6 +13,9 @@ using SparseArrays, Statistics, ChainRulesCore, Zygote using Nonconvex: Nonconvex using Flux using AbstractDifferentiation: AbstractDifferentiation +using StatsBase, Statistics +using Symbolics +using Plots:heatmap const AD = AbstractDifferentiation export Volume, @@ -45,7 +48,17 @@ export Volume, MaterialInterpolation, MultiMaterialVariables, element_densities, - tounit + tounit, + DefGradTensor, + ElementDefGradTensor, + FToK123, + orth_decomp, + sensitivityFieldFncs, + sensitivityPVW, + Entropy_Calc, + Entropy, + SMu_gen, + SAlpha_gen const to = TimerOutput() @@ -68,6 +81,10 @@ include("assemble_K.jl") include("element_ksigma.jl") include("element_k.jl") +# Goodness related +include("defgrad.jl") +include("goodness.jl") + # TODO no rrules yet include("truss_stress.jl") diff --git a/src/Functions/defgrad.jl b/src/Functions/defgrad.jl new file mode 100644 index 00000000..90a47b81 --- /dev/null +++ b/src/Functions/defgrad.jl @@ -0,0 +1,133 @@ +@params struct DefGradTensor{T} <: AbstractFunction{T} + problem::Any + solver::Any + global_dofs::Vector{Int} + cellvalues::Any + cells::Any + _::T +end +function DefGradTensor(solver) + problem = solver.problem + dh = problem.ch.dh + n = ndofs_per_cell(dh) + global_dofs = zeros(Int, n) + cellvalues = solver.elementinfo.cellvalues + return DefGradTensor( + problem, solver, global_dofs, cellvalues, collect(CellIterator(dh)), 0.0 + ) +end + +function Ferrite.reinit!(s::DefGradTensor, cellidx) + reinit!(s.cellvalues, s.cells[cellidx]) + celldofs!(s.global_dofs, s.problem.ch.dh, cellidx) + return s +end +function ChainRulesCore.rrule(::typeof(reinit!), st::DefGradTensor, cellidx) + return reinit!(st, cellidx), _ -> (NoTangent(), NoTangent(), NoTangent()) +end + +function (f::DefGradTensor)(dofs::DisplacementResult) # This is called in the PkgTest.ipynb by good(sim1_u) ----------------------------------------------------[C1] root for interation being studied + return map(1:length(f.cells)) do cellidx + cf = f[cellidx] # This runs C2 in order to makea ElementDefGradTensor struct for this cell with indiex cellidx + return cf(dofs) # cf is a ElementDefGradTensor struct that will run C4 to yield a __________ + end +end + +@params struct ElementDefGradTensor{T} <: AbstractFunction{T} # C2 creates one of these objects ----------------------------------------------------------------[C3] + defgrad_tensor::DefGradTensor{T} + cell::Any + cellidx::Any +end +function Base.getindex(f::DefGradTensor{T}, cellidx) where {T} # This is encounter in C1 -----------------------------------------------------------------------[C2] + reinit!(f, cellidx) + return ElementDefGradTensor(f, f.cells[cellidx], cellidx) +end + +function Ferrite.reinit!(s::ElementDefGradTensor, cellidx) + reinit!(s.defgrad_tensor, cellidx) + return s +end +function ChainRulesCore.rrule(::typeof(reinit!), st::ElementDefGradTensor, cellidx) + return reinit!(st, cellidx), _ -> (NoTangent(), NoTangent(), NoTangent()) +end + +function Fto3by3(F::Matrix{T}) where {T} + if size(F) == (2,2) + F = [F zeros(T,2); zeros(T,1,2) T(1)] # plane strain assumption dictates that F₃₃ = 1 + elseif size(F) == (3,3) + F + else + error("Invalid dimensions of deformation gradient tensor") + end +end + +function (f::ElementDefGradTensor)(u::DisplacementResult; element_dofs=false) #---------------------------------------------------------------------------------[C4] summing from C8 + st = f.defgrad_tensor + reinit!(f, f.cellidx) # refreshing f + if element_dofs # i think this is just choosing between local and global dofs + cellu = u.u + else + cellu = u.u[copy(st.global_dofs)] + end + n_basefuncs = getnbasefunctions(st.cellvalues) + n_quad = getnquadpoints(st.cellvalues) + dim = TopOptProblems.getdim(st.problem) + V = sum(st.cellvalues.detJdV) + F = sum(map(1:n_quad) do q_point + dΩ = getdetJdV(st.cellvalues, q_point) + sum(map(1:n_basefuncs) do a + _u = cellu[dim * (a - 1) .+ (1:dim)] + return tensor_kernel(f, q_point, a)(DisplacementResult(_u)) + end) * dΩ + end) ./ V + I(dim) + + if dim == 2 + F = Fto3by3(F) + end + return F +end + +@params struct ElementDefGradTensorKernel{T} <: AbstractFunction{T} # ------------------------------------------------------------------------------------------------------------[C7] + E::T + ν::T + q_point::Int + a::Int + cellvalues::Any + dim::Int +end +function (f::ElementDefGradTensorKernel)(u::DisplacementResult) # ----------------------------------------------------------------------------------------------------------------[C8] ---- nifty + @unpack E, ν, q_point, a, cellvalues = f + ∇ϕ = Vector(shape_gradient(cellvalues, q_point, a)) + return u.u * ∇ϕ' +end +function ChainRulesCore.rrule(f::ElementDefGradTensorKernel, u::DisplacementResult) + v, (∇,) = AD.value_and_jacobian( + AD.ForwardDiffBackend(), u -> vec(f(DisplacementResult(u))), u.u + ) + return reshape(v, f.dim, f.dim), Δ -> (NoTangent(), Tangent{typeof(u)}(; u=∇' * vec(Δ))) # Need to verify that this is true +end + +function tensor_kernel(f::DefGradTensor, quad, basef) # -------------------------------------------------------------------------------------------------------------------------[C6* Altered ] + #if string(typeof(f.problem))[1:12]!="InpStiffness" + return ElementDefGradTensorKernel( + f.problem.E, + f.problem.ν, + quad, + basef, + f.cellvalues, + TopOptProblems.getdim(f.problem), + ) + #else + # return ElementDefGradTensorKernel( + # f.problem.inp_content.E, + # f.problem.inp_content.ν, + # quad, + # basef, + # f.cellvalues, + # TopOptProblems.getdim(f.problem), + #) + #end +end +function tensor_kernel(f::ElementDefGradTensor, quad, basef) # --------------------------------------------------------------------------------------------------------------------[C5] + return tensor_kernel(f.defgrad_tensor, quad, basef) +end \ No newline at end of file diff --git a/src/Functions/displacement.jl b/src/Functions/displacement.jl index 2ea56115..ff23cddf 100644 --- a/src/Functions/displacement.jl +++ b/src/Functions/displacement.jl @@ -13,10 +13,24 @@ mutable struct Displacement{ maxfevals::Int end +@params mutable struct HyperelasticDisplacement{T} <: AbstractFunction{T} + u::AbstractVector{T} # displacement vector + F::AbstractVector # deformation gradient tensor + dudx_tmp::AbstractVector # directional derivative + solver::AbstractHyperelasticSolver + global_dofs::AbstractVector{<:Integer} + fevals::Int + maxfevals::Int +end + function Base.show(::IO, ::MIME{Symbol("text/plain")}, ::Displacement) return println("TopOpt displacement function") end +function Base.show(::IO, ::MIME{Symbol("text/plain")}, ::HyperelasticDisplacement) + return println("TopOpt displacement function for hyperelastic strain regimes") +end + struct DisplacementResult{T,N,A<:AbstractArray{T,N}} <: AbstractArray{T,N} u::A end @@ -43,6 +57,20 @@ function Displacement(solver::AbstractFEASolver; maxfevals=10^8) return Displacement(u, dudx_tmp, solver, global_dofs, 0, maxfevals) end +function Displacement(solver::AbstractHyperelasticSolver; maxfevals=10^8) + dim = TopOptProblems.getdim(solver.problem) + dim == 3 || throw("2D hyperelastic FEA is not supported yet.") + T = eltype(solver.u) + dh = solver.problem.ch.dh + k = ndofs_per_cell(dh) + global_dofs = zeros(Int, k) + total_ndof = ndofs(dh) + u = zeros(T, total_ndof) + F = [zeros(typeof(solver.elementinfo.Fes[1])) for _ in 1:total_ndof/dim] + dudx_tmp = zeros(T, length(solver.vars)) + return HyperelasticDisplacement(u, F, dudx_tmp, solver, global_dofs, 0, maxfevals) +end + """ # Arguments `x` = design variables @@ -60,6 +88,16 @@ function (dp::Displacement{T})(x::PseudoDensities) where {T} return DisplacementResult(copy(solver.u)) end +function (dp::HyperelasticDisplacement{T})(x::PseudoDensities) where {T} + @unpack solver, global_dofs = dp + @unpack penalty, problem, xmin = solver + dp.fevals += 1 + @assert length(global_dofs) == ndofs_per_cell(solver.problem.ch.dh) + solver.vars .= x.x + solver() + return DisplacementResult(copy(solver.u)), copy(solver.F) #, copy(solver.F) # I need to add F support +end + """ rrule for autodiff. diff --git a/src/Functions/goodness.jl b/src/Functions/goodness.jl new file mode 100644 index 00000000..b8cf8c58 --- /dev/null +++ b/src/Functions/goodness.jl @@ -0,0 +1,325 @@ +#= Orthogonal decomposition function to convert F to an invariant basis (K₁, K₂, K₃) for isotropic hyperelastic strain (Criscione JC et al. JMPS, 2000) +K = FToK123(F) and K = orth_decomp(F) +In the current formulation orth_decomp is not autodiff compatible, but FToK123(F) ≈ orth_decomp(F) +Benchmarking suggests that FToK123 is also about twice as fast as orth_decomp +K is of dims a Vector{Vector{T}} of dim [1:3][1:nels] +Outputs: K = [K₁ K₂ K₃] where each term is a vector of length nels =# + +function FToK123(F::AbstractMatrix{T}) where {T} + @assert size(F) == (3, 3) + g = Vector{T}(undef,3) + + C = transpose(F)*F + λ = sqrt.(eigen(C).values) + + K₁ = log(λ[1]*λ[2]*λ[3]) + g = log.(λ) .- (K₁/3) + K₂ = sqrt(sum(g.^2)) + K₃ = 3*sqrt(6)*g[1]*g[2]*g[3]/(K₂^3) + return [K₁ K₂ K₃] +end + +function FToK123(F::AbstractVector{<:AbstractMatrix{T}}) where {T} + nels = length(F) + Kel = map(i -> FToK123(F[i]), 1:nels) + K = [[Kel[i][j] for i in 1:nels] for j in 1:3] + return K +end + +function orth_decomp(F::AbstractMatrix{T}) where {T} + @assert size(F) == (3, 3) + K = Vector{T}(undef,3) + + J = det(F) + B = F*transpose(F) + V = sqrt(B) + η = log(V) + devη = η - (1/3)*(tr(η))*Matrix(I,3,3) + + K[1] = log(J) + K[2] = sqrt(tr(devη^2)) + Φ = devη/K[2] + K[3] = 3*sqrt(6)*det(Φ) + return K # warning: final K values may need to be modified to ensure they are within physical limits (e.g., K₃ ∈ [-1 1]) +end + +function orth_decomp(F::AbstractVector{<:AbstractMatrix{T}}) where {T} + nels = length(F) + Kel = map(i -> orth_decomp(F[i]), 1:nels) + K = [[Kel[i][j] for i in 1:nels] for j in 1:3] # transposing the vector of vectors to yield K[1:3][el] + return K +end + +#= Function that produces stress and kinematic terms used to construct the virtual fields +stressTerms, kinTerms = sensitivityFieldFunctions(:ConstitutiveModel) +stressTerms = fns[1] +kinTerms = fns[2] +Outputs: stressTerms[i,j] (i.e. ∂²W(K₁,K₂,K₃)/∂Kᵢ∂ξⱼ) and kinTerms[i,j,k] (this is equivalent to d^3W/dK_i/dXi_j/dK_kk) =# + +function sensitivityFieldFncs(matlModel::Symbol) + @variables α K[1:3] # K[1:3] are the orthogonal strain invariants + λᵅ = Array{Any}(undef,3) + λᵅ[1] = exp(α*K[2]*sqrt(2/3)*sin((-asin(K[3])+2π)/3)) + λᵅ[2] = exp(α*K[2]*sqrt(2/3)*sin((-asin(K[3]))/3)) + λᵅ[3] = exp(α*K[2]*sqrt(2/3)*sin((-asin(K[3])-2π)/3)) + + # construct the work function for called constitutive model, W + Ī₁ = sum(substitute(λᵅ[i], Dict(α => 2)) for i in 1:3) + bulk = (exp(K[1])-1)^2 # equivalent to (J-1)² + if matlModel in (:MooneyRivlin, :MR) # note: bulk modulus (κ) must be the final parameter in ξ for all matlModel + ξ = @variables C₁₀ C₀₁ κ + Ī₂ = sum(substitute(λᵅ[i], Dict(α => -2)) for i in 1:3) + W = C₁₀*(Ī₁-3) + C₀₁*(Ī₂-3) + (κ/2)*bulk + elseif matlModel in (:NeoHookean, :NH) + ξ = @variables µ κ + W = (µ/2)*(Ī₁-3) + (κ/2)*bulk + elseif matlModel in (:Yeoh2, :Y2) + ξ = @variables C₁₀ C₂₀ κ + W = C₁₀*(Ī₁-3) + C₂₀*(Ī₁-3)^2 + (κ/2)*bulk + elseif matlModel in (:Yeoh3, :Y3) + ξ = @variables C₁₀ C₂₀ C₃₀ κ + W = C₁₀*(Ī₁-3) + C₂₀*(Ī₁-3)^2 + C₃₀*(Ī₁-3)^3 + (κ/2)*bulk + end + + # solve for terms used to in stress and kinematic field calculations + stressTerms = Array{Any}(undef,length(K),length(ξ)) + for i = 1:length(K) + for j = 1:length(ξ) + ∂Kᵢ∂ξⱼ = Differential(K[i])*Differential(ξ[j]) # ∂²/∂Kᵢ∂ξⱼ + if i == 1 && j == length(ξ) + stressTerms[i,j] = eval(build_function(expand_derivatives(∂Kᵢ∂ξⱼ(W)),K[1])) # ∂²W(K₁)/∂K₁∂κ + elseif i != 1 && j != length(ξ) + stressTerms[i,j] = eval(build_function(expand_derivatives(∂Kᵢ∂ξⱼ(W)),K[2],K[3])) # ∂²W(K₂,K₃)/∂Kᵢ∂ξⱼ + end + end + end + kinTerms = Array{Any}(undef,length(K),length(ξ),length(K)) + for i = 1:length(K) + for j = 1:length(ξ) + for k = 1:length(K) + ∂Kᵢ∂ξⱼ∂Kₖ = Differential(K[i])*Differential(ξ[j])*Differential(K[k]) # ∂³/∂Kᵢ∂ξⱼ∂Kₖ + if i == 1 && j == length(ξ) && k == 1 + kinTerms[i,j,k] = eval(build_function(expand_derivatives(∂Kᵢ∂ξⱼ∂Kₖ(W)),K[1])) # ∂³W(K₁)/∂²K₁∂κ + elseif i != 1 && j != length(ξ) && k != 1 + kinTerms[i,j,k] = eval(build_function(expand_derivatives(∂Kᵢ∂ξⱼ∂Kₖ(W)),K[2],K[3])) # ∂³W(K₂,K₃)/∂Kᵢ∂ξⱼ∂Kₖ + end + end + end + end + return stressTerms, kinTerms +end + +# Function that produces a volume-weighted element-wise sensitivity metric +# sens_IVW = sensitivity_PVW(x,K,V0,stressTerms,kinTerms,ξ) + # x is the + # K₁, K₂, and K₃ are all vectors of size [nels] (acquired from FtoK123 or orth_decomp) + # Where matl_props is a vector, e.g. for Mooney-Rivlin matl_props could be [0.04 0.001 0.5] + # Vol is a vector of size [nels] that contains information of the volume for each element + # stressTerms and kinTerms are acquired from sensitivityFieldFunctions + # This function only needs to be run once! + +# Outputs +# sens_IVW (sensitivity metric across all elements) + +function sensitivityPVW(x,K,V0,stressTerms,kinTerms,ξ) + # Stress term evaluation at K + ξᵢ∂²W_∂K₁∂ξᵢ = ξ[end]*stressTerms[1,end].(K[1]) + ξᵢ∂²W_∂K₂∂ξᵢ = sum(map(i -> ξ[i]*stressTerms[2,i].(K[2],K[3]), 1:length(ξ)-1)) + ξᵢ∂²W_∂K₃∂ξᵢ = sum(map(i -> ξ[i]*stressTerms[3,i].(K[2],K[3]), 1:length(ξ)-1)) + + # Internal virtual work (IVW) evaluation + IVW_bar = zeros(length(x)) + for i = 1:length(ξ) # ξᵢ + for j = 1:length(K) # Kⱼ + if i == length(ξ) && j == 1 + IVW_bar += 3*V0.*(x.^2).*ξᵢ∂²W_∂K₁∂ξᵢ.*kinTerms[1,i,j].(K[1]) + elseif i != length(ξ) && j != 1 + IVW_bar += V0.*(x.^2).*(ξᵢ∂²W_∂K₂∂ξᵢ.*kinTerms[2,i,j].(K[2],K[3]) + (9*(ones(length(x))-(K[3].^2))./(K[2].^2)).*ξᵢ∂²W_∂K₃∂ξᵢ.*kinTerms[3,i,j].(K[2],K[3])) + end + end + end + return IVW_bar/sum(V0.*x) # Volume-weighted element-wise sensitivity metric of all IVW fields +end + +# NOTE: Functions below are legacy material. They should work properly, however, they are less well vetted than the sensitivity related functions. + +function Entropy_Calc(h::Matrix{Int64}) + #Calcuates entropy based on h matrix given + p=transpose(h)./(sum(h)) + ElementWise_Entropy=(p.*log.(p)) + temp_p=deepcopy(ElementWise_Entropy) + replace!(temp_p, -NaN=>0) + H_pwise=-sum(temp_p) + #Entropy of Gaussian was not considered + ElementWise_Entropy*=-1 + return (H_pwise,ElementWise_Entropy) +end + +function Entropy(F::Vector{Vector{Matrix{Float64}}},n_bins::Int64=100,offset::Float64=0.005,make_plot::Bool=false,save_plot::Bool=false,saveplot_name::String="",saveplot_path::String="") + #Creates a vector of entropy values, in the case multiple time signatures are given. If not and F is of length one, returns a vector of length one + #Elemetnt-Wise value returned by function retains NaN values and is therefore not able to be differenciated without further tinkering + Entropy_Value_List=[] + ElementWise_Entropy_List=[] + for t in 1:length(F) + #Pre-processing of F vector + F_copy=deepcopy(F[t]) + for i in 1:length(F_copy) + no_nan=true + for q in 1:3 + for w=1:3 + if isnan(F_copy[i][q,w])==true + no_nan=false + end + end + end + F_copy[i]=reshape(F_copy[i],(3,3)) + F_copy[i]=F_copy[i]/(det(F_copy[i])^(1/3)) + end + + #Gets vectors of k1, k2, and k3 + k1,k2,k3=FToK2AndK3(F_copy) + K3_Vector=k3 + K2_Vector=k2 + + #Generates 2 histogram, and the h matrix which represents the counts in matrix form + K2_edges = (0.0-offset, 1.0-offset) + K3_edges=(-1-offset, 1.0-offset) + Hist_K2_Edges=K2_edges[1]:((K2_edges[2]-K2_edges[1])/n_bins):K2_edges[2] + Hist_K3_Edges=K3_edges[1]:((K3_edges[2]-K3_edges[1])/n_bins):K3_edges[2] + Hist_Matrix=fit(Histogram, (map(Float64,K3_Vector), map(Float64,K2_Vector)), (Hist_K2_Edges, Hist_K3_Edges)) + #Below line defines bins differently, in a way which doesn't quite work but it worth keep around + #Hist_Matrix=fit(Histogram,(map(Float64,K2_Vector),map(Float64,K3_Vector)),nbins=(n_bins,n_bins)) + h=Hist_Matrix.weights + #Makes plots, and subsequently saves plots if inputs are set ot do so + if make_plot==true + #Below line creates figure slightly differently, using a method that is not the same as what was used to get h + # hist = histogram2d(K2_Vector, K3_Vector, nbins=n_bins, color=:viridis, xlims=K2_edges, ylims=K3_edges, background="black") + p=(heatmap(h)) + title!(p,"K2 v. K3"); + xlabel!(p,"K2"); + ylabel!(p,"K3"); + # xlims!(p,(K2_edges[1], K2_edges[2])); + # ylims!(p,(K3_edges[1], K3_edges[2])); + display(p) + if save_plot==true + name=string(saveplot_name,".png") + path=(string(saveplot_path,name)) + savefig(p,path) + end + end + (Entropy_Value,ElementWise_Entropy)=Entropy_Calc(h) + push!(Entropy_Value_List,Entropy_Value) + push!(ElementWise_Entropy_List,ElementWise_Entropy) + end + if length(Entropy_Value_List)==1 + Entropy_Value_List=Entropy_Value_List[1] + ElementWise_Entropy_List=ElementWise_Entropy_List[1] + end + return (map(Float64,Entropy_Value_List),map(Float64,ElementWise_Entropy_List)) + end + + function Entropy(F::Vector{Matrix{Float64}},n_bins::Int64=100,offset::Float64=0.005,make_plot::Bool=false,save_plot::Bool=false,saveplot_name::String="",saveplot_path::String="") + #Creates a vector of entropy values, in the case multiple time signatures are given. If not and F is of length one, returns a vector of length one + F=[F] + Entropy_Value_List=[] + ElementWise_Entropy_List=[] + for t in 1:length(F) + #Pre-processing of F vector + F_copy=deepcopy(F[t]) + for i in 1:length(F_copy) + no_nan=true + for q in 1:3 + for w=1:3 + if isnan(F_copy[i][q,w])==true + no_nan=false + end + end + end + F_copy[i]=reshape(F_copy[i],(3,3)) + F_copy[i]=F_copy[i]/(det(F_copy[i])^(1/3)) + end + + #Gets vectors of k1, k2, and k3 + k1,k2,k3=FToK2AndK3(F_copy) + K3_Vector=k3 + K2_Vector=k2 + + #Generates 2 histogram, and the h matrix which represents the counts in matrix form + K2_edges = (0.0-offset, 1.0-offset) + K3_edges=(-1-offset, 1.0-offset) + Hist_K2_Edges=K2_edges[1]:((K2_edges[2]-K2_edges[1])/n_bins):K2_edges[2] + Hist_K3_Edges=K3_edges[1]:((K3_edges[2]-K3_edges[1])/n_bins):K3_edges[2] + Hist_Matrix=fit(Histogram, (map(Float64,K3_Vector), map(Float64,K2_Vector)), (Hist_K2_Edges, Hist_K3_Edges)) + #Below line defines bins differently, in a way which doesn't quite work but it worth keep around + #Hist_Matrix=fit(Histogram,(map(Float64,K2_Vector),map(Float64,K3_Vector)),nbins=(n_bins,n_bins)) + h=Hist_Matrix.weights + #Makes plots, and subsequently saves plots if inputs are set ot do so + if make_plot==true + #Below line creates figure slightly differently, using a method that is not the same as what was used to get h + # hist = histogram2d(K2_Vector, K3_Vector, nbins=n_bins, color=:viridis, xlims=K2_edges, ylims=K3_edges, background="black") + p=(heatmap(h)) + title!(p,"K2 v. K3"); + xlabel!(p,"K2"); + ylabel!(p,"K3"); + # xlims!(p,(K2_edges[1], K2_edges[2])); + # ylims!(p,(K3_edges[1], K3_edges[2])); #Error on this line?? + display(p) + if save_plot==true + name=string(saveplot_name,".png") + path=(string(saveplot_path,name)) + savefig(p,path) + end + end + (Entropy_Value,ElementWise_Entropy)=Entropy_Calc(h) + push!(Entropy_Value_List,Entropy_Value) + push!(ElementWise_Entropy_List,ElementWise_Entropy) + end + if length(Entropy_Value_List)==1 + Entropy_Value_List=Entropy_Value_List[1] + ElementWise_Entropy_List=ElementWise_Entropy_List[1] + end + return (map(Float64,Entropy_Value_List),map.(Float64,ElementWise_Entropy_List)) +end + +function SMu_gen(k2_list::Vector{Float64},k3_list::Vector{Float64},alpha::Float64=2.0,mu::Float64=1.0) + #Takes in a list of k2 and k3 values, and returns both total and element-wise sensitivity in terms of mu + SMu_List=[] + for i in 1:length(k3_list) + k3=k3_list[i] + k2=k2_list[i] + g1 = sqrt(2/3)*sin(-asin(k3)/3+(2*pi/3)) + g2 = sqrt(2/3)*sin(-asin(k3)/3) + g3 = sqrt(2/3)*sin(-asin(k3)/3-(2*pi/3)) + + val1=g1*exp(g1*alpha*k2) + val2=g2*exp((g2*alpha*k2)) + val3=g3*exp((g3*alpha*k2)) + SMu=val1+val2+val3 + + push!(SMu_List,SMu) + end + Sum=sum(SMu_List) + return map(Float64,SMu_List), Sum +end + +function SAlpha_gen(k2_list::Vector{Float64},k3_list::Vector{Float64},alpha::Float64=2.0,mu::Float64=1.0) + #Takes in a list of k2 and k3 values, and returns both total and element-wise sensitivity in terms of Alpha + SAlpha_List=[] + for i in 1:length(k3_list) + k3=k3_list[i] + k2=k2_list[i] + + g1 = sqrt(2/3)*sin(-asin(k3)/3+(2*pi/3)) + g2 = sqrt(2/3)*sin(-asin(k3)/3) + g3 = sqrt(2/3)*sin(-asin(k3)/3-(2*pi/3)) + + val1=(g1^2)*k2*(exp(g1*alpha*k2)) + val2=(g2^2)*k2*(exp(g2*alpha*k2)) + val3=(g3^2)*k2*(exp(g3*alpha*k2)) + + SAlpha=mu*(val1+val2+val3) + SAlpha_List=vcat(SAlpha_List,[SAlpha]) + end + Sum=sum(SAlpha_List) + return (map(Float64,SAlpha_List)), Sum +end \ No newline at end of file diff --git a/src/TopOpt.jl b/src/TopOpt.jl index c747e407..856521d3 100644 --- a/src/TopOpt.jl +++ b/src/TopOpt.jl @@ -103,6 +103,7 @@ export TopOpt, Compliance, CG, Direct, + Hyperelastic, Assembly, MatrixFree, FEASolver, diff --git a/src/TopOptProblems/TopOptProblems.jl b/src/TopOptProblems/TopOptProblems.jl index e8bbea4f..4f89e70a 100644 --- a/src/TopOptProblems/TopOptProblems.jl +++ b/src/TopOptProblems/TopOptProblems.jl @@ -36,7 +36,9 @@ export RayProblem, StiffnessTopOptProblem, AbstractTopOptProblem, GlobalFEAInfo, + GlobalFEAInfo_hyperelastic, ElementFEAInfo, + ElementFEAInfo_hyperelastic, YoungsModulus, assemble, assemble_f!, @@ -46,6 +48,9 @@ export RayProblem, bcmatrix, save_mesh, RandomMagnitude, - MultiLoad + MultiLoad, + TensionBar, + TensionRoller, + assemble_hyperelastic! end # module diff --git a/src/TopOptProblems/assemble.jl b/src/TopOptProblems/assemble.jl index f2c5f58e..1e9f0407 100644 --- a/src/TopOptProblems/assemble.jl +++ b/src/TopOptProblems/assemble.jl @@ -32,7 +32,7 @@ function assemble!( _K = K isa Symmetric ? K.data : K _K.nzval .= 0 - assembler = Ferrite.AssemblerSparsityPattern(_K, f, Int[], Int[]) + assembler = Ferrite.AssemblerSparsityPattern(_K, f, Int[], Int[]) # For TopOpt folks https://github.com/Ferrite-FEM/Ferrite.jl/blob/98e2b674a8c6d740a779c022de31f73051cb353c/src/assembler.jl#L131 global_dofs = zeros(Int, ndofs_per_cell(dh)) fe = zeros(typeof(fes[1])) @@ -87,6 +87,84 @@ function assemble!( return nothing end +function assemble_hyperelastic!( + globalinfo::GlobalFEAInfo_hyperelastic{T}, + problem::StiffnessTopOptProblem{dim,T}, + elementinfo::ElementFEAInfo_hyperelastic{dim,T,TK}, + vars=ones(T, getncells(getdh(problem).grid)), + penalty=PowerPenalty(T(1)), + xmin=T(0.001); + assemble_f=true, +) where {dim,T,TK} + ch = problem.ch + dh = ch.dh + K, g = globalinfo.K, globalinfo.g + if assemble_f # temporarily ignored + g .= elementinfo.fixedload + end + Kes, ges = elementinfo.Kes, elementinfo.ges + black = problem.black + white = problem.white + varind = problem.varind + + _K = K isa Symmetric ? K.data : K + _K.nzval .= 0 + assembler = Ferrite.AssemblerSparsityPattern(_K, g, Int[], Int[]) # For TopOpt folks https://github.com/Ferrite-FEM/Ferrite.jl/blob/98e2b674a8c6d740a779c022de31f73051cb353c/src/assembler.jl#L131 + + global_dofs = zeros(Int, ndofs_per_cell(dh)) + ge = zeros(typeof(ges[1])) + Ke = zeros(T, size(rawmatrix(Kes[1]))) + + celliterator = CellIterator(dh) + for (i, cell) in enumerate(celliterator) + # get global_dofs for cell#i + celldofs!(global_dofs, dh, i) + ge = ges[i] + _Ke = rawmatrix(Kes[i]) + Ke = _Ke isa Symmetric ? _Ke.data : _Ke + if black[i] + if assemble_f + Ferrite.assemble!(assembler, global_dofs, Ke, ge) + else + Ferrite.assemble!(assembler, global_dofs, Ke) + end + elseif white[i] + if PENALTY_BEFORE_INTERPOLATION + px = xmin + else + px = penalty(xmin) + end + Ke = px * Ke + if assemble_f + ge = px * ge + Ferrite.assemble!(assembler, global_dofs, Ke, ge) + else + Ferrite.assemble!(assembler, global_dofs, Ke) + end + else + if PENALTY_BEFORE_INTERPOLATION + px = density(penalty(vars[varind[i]]), xmin) + else + px = penalty(density(vars[varind[i]], xmin)) + end + Ke = px * Ke + if assemble_f + ge = px * ge + Ferrite.assemble!(assembler, global_dofs, Ke, ge) + else + Ferrite.assemble!(assembler, global_dofs, Ke) + end + end + end + + #* apply boundary condition + _K = TK <: Symmetric ? K.data : K + #apply!(_K, g, ch) + apply_zero!(_K, g, ch) + + return nothing +end + function assemble_f( problem::StiffnessTopOptProblem{dim,T}, elementinfo::ElementFEAInfo{dim,T}, diff --git a/src/TopOptProblems/elementinfo.jl b/src/TopOptProblems/elementinfo.jl index 39ef1521..42453003 100644 --- a/src/TopOptProblems/elementinfo.jl +++ b/src/TopOptProblems/elementinfo.jl @@ -52,6 +52,22 @@ struct ElementFEAInfo{ cells::Tc3 end +@params struct ElementFEAInfo_hyperelastic{dim,T} + Kes::AbstractVector{<:AbstractMatrix{T}} + fes::AbstractVector{<:AbstractVector{T}} + ges::AbstractVector{<:AbstractVector{T}} + fixedload::AbstractVector{T} + Fes::AbstractVector{<:AbstractMatrix{T}} + cellvolumes::AbstractVector{T} + cellvalues::CellValues{dim,T,<:Any} + facevalues::FaceValues{<:Any,T,<:Any} + metadata::Metadata + black::AbstractVector + white::AbstractVector + varind::AbstractVector{Int} + cells::Any +end + function Base.show(io::Base.IO, ::MIME"text/plain", efeainfo::ElementFEAInfo) return print( io, @@ -100,6 +116,40 @@ function ElementFEAInfo( ) end +function ElementFEAInfo_hyperelastic( + mp, sp, u, quad_order=2, ::Type{Val{mat_type}}=Val{:Static}, nearlyincompressible=false; ts = 1.0, +) where {mat_type} + Kes, weights, dloads, ges, Fes, cellvalues, facevalues = make_Kes_and_fes_hyperelastic( + mp, sp, u, quad_order, Val{mat_type}, ts + ) + element_Kes = convert( # make sure this isn't going to symmetric + Vector{<:ElementMatrix}, + Kes; + bc_dofs=sp.ch.prescribed_dofs, + dof_cells=sp.metadata.dof_cells, + ) + fixedload = Vector(make_cload_hyperelastic(sp,ts)) + assemble_f!(fixedload, sp, dloads) # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! GUT FEELING + #assemble_f!(fixedload, sp, ges) # it would seem that this was an issue + cellvolumes = get_cell_volumes(sp, cellvalues) + cells = sp.ch.dh.grid.cells + return ElementFEAInfo_hyperelastic( + element_Kes, + weights, + ges, + fixedload, # this is g version now!!!!!!!! + Fes, + cellvolumes, + cellvalues, + facevalues, + sp.metadata, + sp.black, + sp.white, + sp.varind, + cells, + ) +end + """ struct GlobalFEAInfo{T, TK<:AbstractMatrix{T}, Tf<:AbstractVector{T}, Tchol} K::TK @@ -160,6 +210,17 @@ function GlobalFEAInfo(K, f) return GlobalFEAInfo(K, f, chol, qrfact) end +@params mutable struct GlobalFEAInfo_hyperelastic{T} + K::AbstractMatrix{T} + g::AbstractVector{T} +end + +function GlobalFEAInfo_hyperelastic(sp::StiffnessTopOptProblem) + K = initialize_K(sp; symmetric=false) + g = initialize_f(sp) + return GlobalFEAInfo_hyperelastic(K, g) +end + """ get_cell_volumes(sp::StiffnessTopOptProblem{dim, T}, cellvalues) diff --git a/src/TopOptProblems/elementmatrix.jl b/src/TopOptProblems/elementmatrix.jl index 26211026..69847e83 100644 --- a/src/TopOptProblems/elementmatrix.jl +++ b/src/TopOptProblems/elementmatrix.jl @@ -65,14 +65,15 @@ function Base.convert( ) where {N,T,TM<:StaticMatrix{N,N,T}} fill_matrix = zero(TM) fill_mask = ones(SVector{N,Bool}) - element_Kes = fill(ElementMatrix(fill_matrix, fill_mask), length(Kes)) + element_Kes = fill(ElementMatrix(fill_matrix, fill_mask), length(Kes)) #lil different for i in bc_dofs d_cells = dof_cells[i] for c in d_cells (cellid, localdof) = c Ke = element_Kes[cellid] new_Ke = @set Ke.mask[localdof] = false - element_Kes[cellid] = Symmetric(new_Ke) + #element_Kes[cellid] = Symmetric(new_Ke) ######################################### DOOOOOOOOOOOOO NOOTTTTTTTTTTTTTTTTTTT FORGETTTTTTTT + element_Kes[cellid] = new_Ke end end for e in 1:length(element_Kes) diff --git a/src/TopOptProblems/matrices_and_vectors.jl b/src/TopOptProblems/matrices_and_vectors.jl index 0d87f74b..d3eea8da 100644 --- a/src/TopOptProblems/matrices_and_vectors.jl +++ b/src/TopOptProblems/matrices_and_vectors.jl @@ -1,54 +1,87 @@ function gettypes( ::Type{T}, # number type ::Type{Val{:Static}}, # matrix type - ::Type{Val{Kesize}}, # matrix size + ::Type{Val{Kesize}}; # matrix size + hyperelastic = false, ) where {T,Kesize} - return SMatrix{Kesize,Kesize,T,Kesize^2}, SVector{Kesize,T} + return hyperelastic ? + (SMatrix{Kesize,Kesize,T,Kesize^2}, SVector{Kesize,T}, SMatrix{3,3,T,3^2}) : + (SMatrix{Kesize,Kesize,T,Kesize^2}, SVector{Kesize,T}) end function gettypes( ::Type{T}, # number type ::Type{Val{:SMatrix}}, # matrix type - ::Type{Val{Kesize}}, # matrix size + ::Type{Val{Kesize}}; # matrix size + hyperelastic = false, ) where {T,Kesize} - return SMatrix{Kesize,Kesize,T,Kesize^2}, SVector{Kesize,T} + if hyperelastic + return SMatrix{Kesize,Kesize,T,Kesize^2}, SVector{Kesize,T}, SMatrix{3,3,T,3^2} + else + return SMatrix{Kesize,Kesize,T,Kesize^2}, SVector{Kesize,T} + end end function gettypes( ::Type{T}, # number type ::Type{Val{:MMatrix}}, # matrix type - ::Type{Val{Kesize}}, # matrix size + ::Type{Val{Kesize}}; # matrix size + hyperelastic = false, ) where {T,Kesize} - return MMatrix{Kesize,Kesize,T,Kesize^2}, MVector{Kesize,T} + if hyperelastic + return MMatrix{Kesize,Kesize,T,Kesize^2}, MVector{Kesize,T}, MMatrix{3,3,T,3^2} + else + return MMatrix{Kesize,Kesize,T,Kesize^2}, MVector{Kesize,T} + end end function gettypes( ::Type{BigFloat}, # number type ::Type{Val{:Static}}, # matrix type - ::Type{Val{Kesize}}, # matrix size + ::Type{Val{Kesize}}; # matrix size + hyperelastic = false, ) where {Kesize} - return SizedMatrix{Kesize,Kesize,BigFloat,Kesize^2}, SizedVector{Kesize,BigFloat} + if hyperelastic + return SizedMatrix{Kesize,Kesize,BigFloat,Kesize^2}, SizedVector{Kesize,BigFloat}, SizedMatrix{3,3,BigFloat,3^2} + else + return SizedMatrix{Kesize,Kesize,BigFloat,Kesize^2}, SizedVector{Kesize,BigFloat} + end end function gettypes( ::Type{BigFloat}, # number type ::Type{Val{:SMatrix}}, # matrix type - ::Type{Val{Kesize}}, # matrix size + ::Type{Val{Kesize}}; # matrix size + hyperelastic = false, ) where {Kesize} - return SizedMatrix{Kesize,Kesize,BigFloat,Kesize^2}, SizedVector{Kesize,BigFloat} + if hyperelastic + return SizedMatrix{Kesize,Kesize,BigFloat,Kesize^2}, SizedVector{Kesize,BigFloat}, SizedMatrix{3,3,BigFloat,3^2} + else + return SizedMatrix{Kesize,Kesize,BigFloat,Kesize^2}, SizedVector{Kesize,BigFloat} + end end function gettypes( ::Type{BigFloat}, # number type ::Type{Val{:MMatrix}}, # matrix type - ::Type{Val{Kesize}}, # matrix size + ::Type{Val{Kesize}}; # matrix size + hyperelastic = false, ) where {Kesize} - return SizedMatrix{Kesize,Kesize,BigFloat,Kesize^2}, SizedVector{Kesize,BigFloat} + if hyperelastic + return SizedMatrix{Kesize,Kesize,BigFloat,Kesize^2}, SizedVector{Kesize,BigFloat}, SizedMatrix{3,3,BigFloat,3^2} + else + return SizedMatrix{Kesize,Kesize,BigFloat,Kesize^2}, SizedVector{Kesize,BigFloat} + end end function gettypes( ::Type{T}, # number type ::Any, # matrix type - ::Any, # matrix size + ::Any; # matrix size + hyperelastic = false, ) where {T} - return Matrix{T}, Vector{T} + if hyperelastic + return Matrix{T}, Vector{T}, Matrix{T} + else + return Matrix{T}, Vector{T} + end end -initialize_K(sp::StiffnessTopOptProblem) = Symmetric(create_sparsity_pattern(sp.ch.dh)) +initialize_K(sp::StiffnessTopOptProblem;symmetric::Bool=true) = symmetric ? Symmetric(create_sparsity_pattern(sp.ch.dh)) : create_sparsity_pattern(sp.ch.dh) initialize_f(sp::StiffnessTopOptProblem{dim,T}) where {dim,T} = zeros(T, ndofs(sp.ch.dh)) @@ -103,6 +136,59 @@ function make_Kes_and_fes(problem, quad_order, ::Type{Val{mat_type}}) where {mat return Kes, weights, dloads, cellvalues, facevalues end +function make_Kes_and_fes_hyperelastic(mp, problem, u, quad_order, ::Type{Val{mat_type}}, ts) where {mat_type} + T = floattype(problem) + dim = getdim(problem) + geom_order = getgeomorder(problem) + dh = getdh(problem) + E = getE(problem) + ν = getν(problem) + ρ = getdensity(problem) + + refshape = Ferrite.getrefshape(dh.field_interpolations[1]) + + #λ = E * ν / ((1 + ν) * (1 - 2 * ν)) + #μ = E / (2 * (1 + ν)) + #δ(i, j) = i == j ? T(1) : T(0) + #g(i, j, k, l) = λ * δ(i, j) * δ(k, l) + μ * (δ(i, k) * δ(j, l) + δ(i, l) * δ(j, k)) + #C = SymmetricTensor{4,dim}(g) + + # Shape functions and quadrature rule + interpolation_space = Lagrange{dim,refshape,geom_order}() + quadrature_rule = QuadratureRule{dim,refshape}(quad_order) + cellvalues = CellScalarValues(quadrature_rule, interpolation_space) # JGB: change from CellScalarValues + cellvaluesV = CellVectorValues(quadrature_rule, interpolation_space) # JGB: change from CellScalarValues + facevalues = FaceScalarValues( + QuadratureRule{dim - 1,refshape}(quad_order), interpolation_space + ) # JGB: change from CellScalarValues + facevaluesV = FaceVectorValues( + QuadratureRule{dim - 1,refshape}(quad_order), interpolation_space + ) + + # Calculate element stiffness matrices + n_basefuncs = getnbasefunctions(cellvalues) + + Kesize = dim * n_basefuncs + MatrixType, VectorType, MatrixTypeF = gettypes(T, Val{mat_type}, Val{Kesize}; hyperelastic = true) + Kes, weights, ges, Fes = _make_Kes_and_weights_hyperelastic( + dh, + Tuple{MatrixType,VectorType,MatrixTypeF}, + Val{n_basefuncs}, + Val{dim * n_basefuncs}, + #C, + mp, + u, + ρ, + quadrature_rule, + cellvalues, + cellvaluesV, + ) + dloads, ges2 = _make_dloads_hyperelastic(weights, problem, facevalues, facevaluesV, ges, ts) + ges += ges2 + + return Kes, weights, dloads, ges, Fes, cellvalues, facevalues # switched to ges2 to solve error where 2x ges2 with no contribution from dload +end + const g = [0.0, 9.81, 0.0] # N/kg or m/s^2 # Element stiffness matrices are StaticArrays @@ -209,6 +295,149 @@ function _make_Kes_and_weights( return Kes, weights end +function Ψ(C, mp) # JGB: add to .ipynb + μ = mp.μ + λ = mp.λ + I1 = tr(C) + J = sqrt(det(C)) + return μ / 2 * (I1 - 3) - μ * log(J) + λ / 2 * (J - 1)^2 # Ferrite.jl version + #return μ / 2 * (Ic - 3 - 2 * log(J)) + λ / 2 * (J-1)^2 # Bower version + #Cnew = @MArray C + #I1bar = Ic*J^-2/3 + #I1bar = Ic*det(C)^-1/3 + #return μ / 2 * (I1bar - 3) + 0.5*(λ + 2μ/3)*(J-1)^2 # ABAQUS/Bower version +end + +function constitutive_driver(C, mp) # JGB removed type ::NeoHook from mp + # Compute all derivatives in one function call + ∂²Ψ∂C², ∂Ψ∂C = Tensors.hessian(y -> Ψ(y, mp), C, :all) + S = 2.0 * ∂Ψ∂C + ∂S∂C = 2.0 * ∂²Ψ∂C² + return S, ∂S∂C +end + +# Element stiffness matrices are StaticArrays +# `weights` : a vector of `xdim` vectors, element_id => self-weight load vector +function _make_Kes_and_weights_hyperelastic( + dh::DofHandler{dim,N,T}, + ::Type{Tuple{MatrixType,VectorType,MatrixTypeF}}, + ::Type{Val{n_basefuncs}}, + ::Type{Val{Kesize}}, + mp, + u, + ρ, + quadrature_rule, + cellvalues, + cellvaluesV, +) where {dim,N,T,MatrixType<:StaticArray,VectorType,MatrixTypeF<:StaticArray,n_basefuncs,Kesize} + # Calculate element stiffness matrices + nel = getncells(dh.grid) + body_force = ρ .* g # Force per unit volume + #Kes = Symmetric{T,MatrixType}[] # JGB: scary (is this sucker symmetric?) + Kes = Vector{MatrixType}() + sizehint!(Kes, nel) + weights = [zeros(VectorType) for i in 1:nel] + ges = [zeros(VectorType) for i in 1:nel] + Fes = [zeros(MatrixTypeF) for _ in 1:nel] + Ke_0 = Matrix{T}(undef, Kesize, Kesize) + celliterator = CellIterator(dh) + for (k, cell) in enumerate(celliterator) + Ke_0 .= 0 + #Ke_02 .= 0 + global_dofs = celldofs(cell) # new + ue = u[global_dofs] # new + reinit!(cellvalues, cell) + reinit!(cellvaluesV, cell) + fe = weights[k] + ge = ges[k] + Fe = Fes[k] + for q_point in 1:getnquadpoints(cellvalues) + dΩ = getdetJdV(cellvalues, q_point) + ∇u = function_gradient(cellvaluesV, q_point, ue) # JGB add (NEEDS TO BE CHECKED!!) + F = one(∇u) + ∇u # JGB add + Fe += F*dΩ + C = tdot(F) # JGB add + S, ∂S∂C = constitutive_driver(C, mp) # JGB add + P = F ⋅ S # JGB add + I = one(S) # JGB add + ∂P∂F = otimesu(I, S) + 2 * otimesu(F, I) ⊡ ∂S∂C ⊡ otimesu(F', I) # JGB add (neither P or F is symmetric so ∂P∂F will not either) + for b in 1:Kesize + ∇ϕb = shape_gradient(cellvaluesV, q_point, b) # JGB: like ∇δui + ϕb = shape_value(cellvaluesV, q_point, b) # JGB: like δui + ∇ϕb∂P∂F = ∇ϕb ⊡ ∂P∂F # Hoisted computation # JGB add + #fe = @set fe[(b - 1) * dim + d2] += ϕb * body_force[d2] * dΩ # weird but probably fine... just not in ferrite.jl example code (leaving in case of zygote issues) + fe = @set fe[b] += ϕb ⋅ body_force * dΩ + ge = @set ge[b] += ( ∇ϕb ⊡ P - ϕb ⋅ body_force ) * dΩ # Add contribution to the residual from this test function + for a in 1:Kesize + ∇ϕa = shape_gradient(cellvaluesV, q_point, a) # JGB: like ∇δuj + Ke_0[a,b] += (∇ϕb∂P∂F ⊡ ∇ϕa) * dΩ + end + end + end + weights[k] = fe + ges[k] = ge + Fes[k] = Fe + if MatrixType <: SizedMatrix # Work around because full constructor errors + #push!(Kes, Symmetric(SizedMatrix{Kesize,Kesize,T}(Ke_0))) + push!(Kes, SizedMatrix{Kesize,Kesize,T}(Ke_0)) + else + #push!(Kes, Symmetric(MatrixType(Ke_0))) + push!(Kes, MatrixType(Ke_0)) + end + end + return Kes, weights, ges, Fes # weights is fes +end +# Fallback +#= bring this up to speed at a later time to match with +function _make_Kes_and_weights( + dh::DofHandler{dim,N,T}, + ::Type{Tuple{MatrixType,VectorType}}, + ::Type{Val{n_basefuncs}}, + ::Type{Val{Kesize}}, + C, + ρ, + quadrature_rule, + cellvalues, +) where {dim,N,T,MatrixType,VectorType,n_basefuncs,Kesize} + # Calculate element stiffness matrices + nel = getncells(dh.grid) + body_force = ρ .* g # Force per unit volume + Kes = let Kesize = Kesize, nel = nel + [Symmetric(zeros(T, Kesize, Kesize), :U) for i in 1:nel] + end + weights = let Kesize = Kesize, nel = nel + [zeros(T, Kesize) for i in 1:nel] + end + Ke_e = zeros(T, dim, dim) + celliterator = CellIterator(dh) + for (k, cell) in enumerate(celliterator) + reinit!(cellvalues, cell) + fe = weights[k] + for q_point in 1:getnquadpoints(cellvalues) + dΩ = getdetJdV(cellvalues, q_point) + for b in 1:n_basefuncs + ∇ϕb = shape_gradient(cellvalues, q_point, b) + ϕb = shape_value(cellvalues, q_point, b) + for d2 in 1:dim + fe[(b - 1) * dim + d2] += ϕb * body_force[d2] * dΩ + for a in 1:n_basefuncs + ∇ϕa = shape_gradient(cellvalues, q_point, a) + Ke_e .= dotdot(∇ϕa, C, ∇ϕb) * dΩ + for d1 in 1:dim + #if dim*(b-1) + d2 >= dim*(a-1) + d1 + Kes[k].data[dim * (a - 1) + d1, dim * (b - 1) + d2] += Ke_e[ + d1, d2 + ] + #end + end + end + end + end + end + end + return Kes, weights +end=# + """ _make_dload(problem) @@ -262,6 +491,69 @@ function _make_dloads(fes, problem, facevalues) return dloads end +function _make_dloads_hyperelastic(fes, problem, facevalues, facevaluesV, ges, ts) + dim = getdim(problem) + N = nnodespercell(problem) + T = floattype(problem) + dloads = deepcopy(fes) + ges2 = deepcopy(ges) + for i in 1:length(dloads) + if eltype(dloads) <: SArray + dloads[i] = zero(eltype(dloads)) + else + dloads[i] .= 0 + end + end + for i in 1:length(ges2) + if eltype(ges2) <: SArray + ges2[i] = zero(eltype(ges2)) + else + ges2[i] .= 0 + end + end + pressuredict = getpressuredict(problem; ts=ts) + dh = getdh(problem) + grid = dh.grid + boundary_matrix = grid.boundary_matrix + cell_coords = zeros(Ferrite.Vec{dim,T}, N) + n_basefuncs = getnbasefunctions(facevalues) + for k in keys(pressuredict) + t = -pressuredict[k] # traction = negative the pressure + faceset = getfacesets(problem)[k] + for (cellid, faceid) in faceset + boundary_matrix[faceid, cellid] || + throw("Face $((cellid, faceid)) not on boundary.") + fe = dloads[cellid] + ge = ges2[cellid] + getcoordinates!(cell_coords, grid, cellid) + reinit!(facevalues, cell_coords, faceid) + for q_point in 1:getnquadpoints(facevalues) + dΓ = getdetJdV(facevalues, q_point) # Face area + normal = getnormal(facevalues, q_point) # Nomral vector at quad point + for i in 1:n_basefuncs + ϕ = shape_value(facevalues, q_point, i) # Shape function value + for d in 1:dim + if fe isa SArray + fe = @set fe[(i - 1) * dim + d] += ϕ * t * normal[d] * dΓ + else + fe[(i - 1) * dim + d] += ϕ * t * normal[d] * dΓ + end + if ge isa SArray + ge = @set ge[(i - 1) * dim + d] -= ϕ * t * normal[d] * dΓ + else + ge[(i - 1) * dim + d] -= ϕ * t * normal[d] * dΓ + end + end + end + end + dloads[cellid] = fe + ges2[cellid] = ge + end + end + + return dloads, ges2 +end + """ make_cload(problem) @@ -287,3 +579,24 @@ function make_cload(problem) end return sparsevec(inds, vals, ndofs(dh)) end + +function make_cload_hyperelastic(problem,ts) + T = floattype(problem) + dim = getdim(problem) + cloads = getcloaddict(problem; ts=ts) + dh = getdh(problem) + metadata = getmetadata(problem) + node_dofs = metadata.node_dofs + inds = Int[] + vals = T[] + for nodeidx in keys(cloads) + for (dofidx, force) in enumerate(cloads[nodeidx]) + if force != 0 + dof = node_dofs[(nodeidx - 1) * dim + dofidx] + push!(inds, dof) + push!(vals, force) + end + end + end + return sparsevec(inds, vals, ndofs(dh)) +end diff --git a/src/TopOptProblems/problem_types.jl b/src/TopOptProblems/problem_types.jl index 917d07ad..d5958ed2 100644 --- a/src/TopOptProblems/problem_types.jl +++ b/src/TopOptProblems/problem_types.jl @@ -21,9 +21,9 @@ getgeomorder(p::StiffnessTopOptProblem) = nnodespercell(p) == 9 ? 2 : 1 getdensity(::StiffnessTopOptProblem{dim,T}) where {dim,T} = T(0) getmetadata(p::StiffnessTopOptProblem) = p.metadata getdh(p::StiffnessTopOptProblem) = p.ch.dh -getcloaddict(p::StiffnessTopOptProblem{dim,T}) where {dim,T} = Dict{String,Vector{T}}() -getpressuredict(p::StiffnessTopOptProblem{dim,T}) where {dim,T} = Dict{String,T}() -getfacesets(p::StiffnessTopOptProblem{dim,T}) where {dim,T} = Dict{String,Tuple{Int,T}}() +getcloaddict(p::StiffnessTopOptProblem{dim,T}; kwargs...) where {dim,T} = Dict{String,Vector{T}}() +getpressuredict(p::StiffnessTopOptProblem{dim,T}; kwargs...) where {dim,T} = Dict{String,T}() +getfacesets(p::StiffnessTopOptProblem{dim,T}; kwargs...) where {dim,T} = Dict{String,Tuple{Int,T}}() Ferrite.getncells(problem::StiffnessTopOptProblem) = Ferrite.getncells(getdh(problem).grid) """ @@ -174,7 +174,8 @@ function PointLoadCantilever( ) add!(ch, dbc) close!(ch) - t = T(0) + + t = T(1) update!(ch, t) metadata = Metadata(dh) @@ -347,7 +348,7 @@ function HalfMBB( add!(ch, dbc2) close!(ch) - t = T(0) + t = T(1) update!(ch, t) metadata = Metadata(dh) @@ -365,9 +366,8 @@ function HalfMBB( return HalfMBB(rect_grid, E, ν, ch, force, force_dof, black, white, varind, metadata) end -nnodespercell(p::Union{PointLoadCantilever,HalfMBB}) = nnodespercell(p.rect_grid) -function getcloaddict(p::Union{PointLoadCantilever{dim,T},HalfMBB{dim,T}}) where {dim,T} - f = T[0, -p.force, 0] +function getcloaddict(p::Union{PointLoadCantilever{dim,T},HalfMBB{dim,T}}; ts=1.0) where {dim,T} + f = T[0, -p.force*ts, 0] fnode = Tuple(getnodeset(p.rect_grid.grid, "down_force"))[1] return Dict{Int,Vector{T}}(fnode => f) end @@ -496,7 +496,7 @@ function LBeam( add!(ch, dbc) close!(ch) - t = T(0) + t = T(1) update!(ch, t) metadata = Metadata(dh) @@ -556,8 +556,8 @@ end nnodespercell(p::LBeam{T,N}) where {T,N} = N getdim(::LBeam) = 2 -function getcloaddict(p::LBeam{T}) where {T} - f = T[0, -p.force] +function getcloaddict(p::LBeam{T};ts=1.0) where {T} + f = T[0, -p.force*ts] fnode = Tuple(getnodeset(getdh(p).grid, "load"))[1] return Dict{Int,Vector{T}}(fnode => f) end @@ -656,7 +656,7 @@ function TieBeam( add!(ch, dbc2) close!(ch) - t = T(0) + t = T(1) update!(ch, t) metadata = Metadata(dh) @@ -669,8 +669,8 @@ end getdim(::TieBeam) = 2 nnodespercell(::TieBeam{T,N}) where {T,N} = N -function getpressuredict(p::TieBeam{T}) where {T} - return Dict{String,T}("rightload" => 2 * p.force, "bottomload" => -p.force) +function getpressuredict(p::TieBeam{T}; ts=1.0) where {T} + return Dict{String,T}("rightload" => 2 * p.force * ts, "bottomload" => -p.force * ts) end getfacesets(p::TieBeam) = getdh(p).grid.facesets @@ -754,7 +754,7 @@ function RayProblem( add!(ch, dbc) end close!(ch) - t = T(0) + t = T(1) update!(ch, t) metadata = Metadata(dh) @@ -770,5 +770,252 @@ function RayProblem( return RayProblem(rect_grid, 1.0, 0.3, ch, loadsdict, black, white, varind, metadata) end -nnodespercell(p::RayProblem) = nnodespercell(p.rect_grid) -getcloaddict(p::RayProblem) = p.loads +function getcloaddict(p::RayProblem; ts=1.0) + return Dict(k => v .* ts for (k, v) in p.loads) +end +""" +``` +///**********************************-> +///* *-> +///* *-> disp +///* *-> +///**********************************-> + + +@params struct TensionBar{dim, T, N, M} <: StiffnessTopOptProblem{dim, T} + rect_grid::RectilinearGrid{dim, T, N, M} + E::T + ν::T + ch::ConstraintHandler{<:DofHandler{dim, <:Cell{dim,N,M}, T}, T} + disp::T + black::AbstractVector + white::AbstractVector + varind::AbstractVector{Int} + metadata::Metadata +end +``` + +- `dim`: dimension of the problem +- `T`: number type for computations and coordinates +- `N`: number of nodes in a cell of the grid +- `M`: number of faces in a cell of the grid +- `rect_grid`: a RectilinearGrid struct +- `E`: Young's modulus +- `ν`: Poisson's ration +- `disp`: displacement over the right face of the beam (positive is right) +- `ch`: a `Ferrite.ConstraintHandler` struct +- `metadata`: Metadata having various cell-node-dof relationships +- `black`: a `BitVector` of length equal to the number of elements where `black[e]` is 1 iff the `e`^th element must be part of the final design +- `white`: a `BitVector` of length equal to the number of elements where `white[e]` is 1 iff the `e`^th element must not be part of the final design +- `varind`: an `AbstractVector{Int}` of length equal to the number of elements where `varind[e]` gives the index of the decision variable corresponding to element `e`. Because some elements can be fixed to be black or white, not every element has a decision variable associated. +""" +@params struct TensionBar{dim,T,N,M} <: StiffnessTopOptProblem{dim,T} + rect_grid::RectilinearGrid{dim,T,N,M} + E::T + ν::T + ch::ConstraintHandler{<:DofHandler{dim,<:Cell{dim,N,M},T},T} + disp::T + black::AbstractVector + white::AbstractVector + varind::AbstractVector{Int} + metadata::Metadata +end + +function Base.show(::IO, ::MIME{Symbol("text/plain")}, ::TensionBar) + return println("TopOpt tension bar problem") +end + +function TensionBar( + ::Type{Val{CellType}}, + nels::NTuple{dim,Int}, + sizes::NTuple{dim}, + E, + ν, + disp, +) where {dim,CellType} + _T = promote_type(eltype(sizes), typeof(E), typeof(ν), typeof(disp)) + if _T <: Integer + T = Float64 + else + T = _T + end + if CellType === :Linear || dim === 3 + rect_grid = RectilinearGrid(Val{:Linear}, nels, T.(sizes)) + else + rect_grid = RectilinearGrid(Val{:Quadratic}, nels, T.(sizes)) + end + + if haskey(rect_grid.grid.facesets, "fixed_left") + pop!(rect_grid.grid.facesets, "fixed_left") + end + addnodeset!(rect_grid.grid, "fixed_left", x -> left(rect_grid, x)) + + if haskey(rect_grid.grid.nodesets, "disp_right") + pop!(rect_grid.grid.nodesets, "disp_right") + end + addnodeset!(rect_grid.grid, "disp_right", x -> right(rect_grid, x)) + + # Create displacement field u + dh = DofHandler(rect_grid.grid) + if CellType === :Linear || dim === 3 + push!(dh, :u, dim) # Add a displacement field + else + ip = Lagrange{2,RefCube,2}() + push!(dh, :u, dim, ip) # Add a displacement field + end + close!(dh) + + ch = ConstraintHandler(dh) + + dbc_left = Dirichlet(:u, getnodeset(rect_grid.grid, "fixed_left"), (x, t) -> zeros(T, dim), collect(1:dim)) + add!(ch, dbc_left) + dbc_right = Dirichlet(:u, getnodeset(rect_grid.grid, "disp_right"), (x, t) -> t*T[disp; zeros(dim-1)], collect(1:dim)) + add!(ch, dbc_right) + close!(ch) + t = T(1) + update!(ch, t) + + metadata = Metadata(dh) + + N = nnodespercell(rect_grid) + M = nfacespercell(rect_grid) + + black, white = find_black_and_white(dh) + varind = find_varind(black, white) + + return TensionBar(rect_grid, E, ν, ch, disp, black, white, varind, metadata) +end + +""" +``` +O**********************************-> +O* *-> +O* *-> disp x +O* *-> +O**********************************-> +// + +@params struct TensionRoller{dim, T, N, M} <: StiffnessTopOptProblem{dim, T} + rect_grid::RectilinearGrid{dim, T, N, M} + E::T + ν::T + ch::ConstraintHandler{<:DofHandler{dim, <:Cell{dim,N,M}, T}, T} + disp::T + black::AbstractVector + white::AbstractVector + varind::AbstractVector{Int} + metadata::Metadata +end +``` + +- `dim`: dimension of the problem +- `T`: number type for computations and coordinates +- `N`: number of nodes in a cell of the grid +- `M`: number of faces in a cell of the grid +- `rect_grid`: a RectilinearGrid struct +- `E`: Young's modulus +- `ν`: Poisson's ration +- `disp`: displacement over the right face of the beam (positive is right) +- `ch`: a `Ferrite.ConstraintHandler` struct +- `metadata`: Metadata having various cell-node-dof relationships +- `black`: a `BitVector` of length equal to the number of elements where `black[e]` is 1 iff the `e`^th element must be part of the final design +- `white`: a `BitVector` of length equal to the number of elements where `white[e]` is 1 iff the `e`^th element must not be part of the final design +- `varind`: an `AbstractVector{Int}` of length equal to the number of elements where `varind[e]` gives the index of the decision variable corresponding to element `e`. Because some elements can be fixed to be black or white, not every element has a decision variable associated. +""" +@params struct TensionRoller{dim,T,N,M} <: StiffnessTopOptProblem{dim,T} + rect_grid::RectilinearGrid{dim,T,N,M} + E::T + ν::T + ch::ConstraintHandler{<:DofHandler{dim,<:Cell{dim,N,M},T},T} + disp::T + black::AbstractVector + white::AbstractVector + varind::AbstractVector{Int} + metadata::Metadata +end + +function Base.show(::IO, ::MIME{Symbol("text/plain")}, ::TensionRoller) + return println("TopOpt tension bar problem with unconstrained y and z displacement at the grips") +end + +function TensionRoller( + ::Type{Val{CellType}}, + nels::NTuple{dim,Int}, + sizes::NTuple{dim}, + E, + ν, + disp, +) where {dim,CellType} + _T = promote_type(eltype(sizes), typeof(E), typeof(ν), typeof(disp)) + if _T <: Integer + T = Float64 + else + T = _T + end + if CellType === :Linear || dim === 3 + rect_grid = RectilinearGrid(Val{:Linear}, nels, T.(sizes)) + else + rect_grid = RectilinearGrid(Val{:Quadratic}, nels, T.(sizes)) + end + + if haskey(rect_grid.grid.facesets, "fixed_left") + pop!(rect_grid.grid.facesets, "fixed_left") + end + addnodeset!(rect_grid.grid, "fixed_left", x -> left(rect_grid, x)) + + if haskey(rect_grid.grid.nodesets, "disp_right") + pop!(rect_grid.grid.nodesets, "disp_right") + end + addnodeset!(rect_grid.grid, "disp_right", x -> right(rect_grid, x)) + + if haskey(rect_grid.grid.nodesets, "fixed_pt1") + pop!(rect_grid.grid.nodesets, "fixed_pt1") + end + if haskey(rect_grid.grid.nodesets, "fixed_pt2") + pop!(rect_grid.grid.nodesets, "fixed_pt2") + end + if dim == 2 + addnodeset!(rect_grid.grid, "fixed_pt1", x -> left(rect_grid, x) && bottom(rect_grid, x)) + elseif dim == 3 + addnodeset!(rect_grid.grid, "fixed_pt1", x -> left(rect_grid, x) && bottom(rect_grid, x) && front(rect_grid, x)) + addnodeset!(rect_grid.grid, "fixed_pt2", x -> left(rect_grid, x) && bottom(rect_grid, x) && back(rect_grid, x)) + end + + # Create displacement field u + dh = DofHandler(rect_grid.grid) + if CellType === :Linear || dim === 3 + push!(dh, :u, dim) # Add a displacement field + else + ip = Lagrange{2,RefCube,2}() + push!(dh, :u, dim, ip) # Add a displacement field + end + close!(dh) + + ch = ConstraintHandler(dh) + + dbc_left = Dirichlet(:u, getnodeset(rect_grid.grid, "fixed_left"), (x, t) -> zeros(T, 1), [1]) # set u1 to 0 for left face + add!(ch, dbc_left) + dbc_right = Dirichlet(:u, getnodeset(rect_grid.grid, "disp_right"), (x, t) -> t*T[disp], [1]) # set u1 to 'disp' for right face + add!(ch, dbc_right) + dbc_fixed_midpt1 = Dirichlet(:u, getnodeset(rect_grid.grid, "fixed_pt1"), (x, t) -> zeros(T, dim), collect(1:dim)) # fix left bottom (front) point to prevent translation in y + add!(ch, dbc_fixed_midpt1) + if dim == 3 + dbc_fixed_midpt2 = Dirichlet(:u, getnodeset(rect_grid.grid, "fixed_pt2"), (x, t) -> zeros(T, 1), [2]) # set u2 = 0 for left bottom back point to prevent moment about x-axis + add!(ch, dbc_fixed_midpt2) + end + close!(ch) + t = T(1) + update!(ch, t) + + metadata = Metadata(dh) + + N = nnodespercell(rect_grid) + M = nfacespercell(rect_grid) + + black, white = find_black_and_white(dh) + varind = find_varind(black, white) + + return TensionRoller(rect_grid, E, ν, ch, disp, black, white, varind, metadata) +end + +nnodespercell(p::Union{PointLoadCantilever,HalfMBB,RayProblem,TensionBar,TensionRoller}) = nnodespercell(p.rect_grid) \ No newline at end of file