From d16be9f3d74b807a075d6d867b4f70d4f5115c7a Mon Sep 17 00:00:00 2001 From: Tolulope Oladele <59896301+tolulade2023@users.noreply.github.com> Date: Fri, 15 Mar 2024 14:45:32 -0400 Subject: [PATCH 01/14] Goodness intialization Adding functions relevant to the optimzaiton of mechanical test data richness and experimental certainty. --- src/Functions/Functions.jl | 17 +- src/Functions/defgrad.jl | 436 ++++++++++++++++++++++++++++ src/Functions/goodness.jl | 277 ++++++++++++++++++ src/TopOptProblems/problem_types.jl | 98 ++++++- 4 files changed, 826 insertions(+), 2 deletions(-) create mode 100644 src/Functions/defgrad.jl create mode 100644 src/Functions/goodness.jl diff --git a/src/Functions/Functions.jl b/src/Functions/Functions.jl index 38a706e0..84456461 100644 --- a/src/Functions/Functions.jl +++ b/src/Functions/Functions.jl @@ -13,6 +13,8 @@ using SparseArrays, Statistics, ChainRulesCore, Zygote using Nonconvex: Nonconvex using Flux using AbstractDifferentiation: AbstractDifferentiation +using StatsBase, Statistics +using Plots:heatmap const AD = AbstractDifferentiation export Volume, @@ -45,7 +47,16 @@ export Volume, MaterialInterpolation, MultiMaterialVariables, element_densities, - tounit + tounit, + DefGradTensor, + GoodnessTensor, + ElementDefGradTensor, + ElementGoodnessTensor, + FToK2AndK3, + Entropy_Calc, + Entropy, + SMu_gen, + SAlpha_gen const to = TimerOutput() @@ -68,6 +79,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..df3a34e8 --- /dev/null +++ b/src/Functions/defgrad.jl @@ -0,0 +1,436 @@ +@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 (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) + return sum( + map(1:n_basefuncs, 1:n_quad) do a, q_point + _u = cellu[dim * (a - 1) .+ (1:dim)] + return tensor_kernel(f, q_point, a)(DisplacementResult(_u)) + end, + ) + I(3) +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] +# @unpack E, ν, q_point, a, cellvalues = f +# ∇ϕ = Vector(shape_gradient(cellvalues, q_point, a)) +# ϵ = (u.u .* ∇ϕ' .+ ∇ϕ .* u.u') ./ 2 +# c1 = E * ν / (1 - ν^2) * sum(diag(ϵ)) +# c2 = E * ν * (1 + ν) +# return c1 * I + c2 * ϵ +# end +function (f::ElementDefGradTensorKernel)(u::DisplacementResult) # ----------------------------------------------------------------------------------------------------------------[C8] ---- nifty + @unpack E, ν, q_point, a, cellvalues = f + ∇ϕ = Vector(shape_gradient(cellvalues, q_point, a)) + # ϵ = (u.u .* ∇ϕ' .+ ∇ϕ .* u.u') ./ 2 + # c1 = E * ν / (1 - ν^2) * sum(diag(ϵ)) + # c2 = E * ν * (1 + ν) + F_quad = u.u * ∇ϕ' # should this I(3) go somewhere else + return F_quad +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(Δ))) +end + +# function tensor_kernel(f::DefGradTensor, quad, basef) # -------------------------------------------------------------------------------------------------------------------------[C6] +# return ElementDefGradTensorKernel( +# f.problem.E, +# f.problem.ν, +# quad, +# basef, +# f.cellvalues, +# TopOptProblems.getdim(f.problem), +# ) +# 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 + + + +# function von_mises(σ::AbstractMatrix) +# if size(σ, 1) == 2 +# t1 = σ[1, 1]^2 - σ[1, 1] * σ[2, 2] + σ[2, 2]^2 +# t2 = 3 * σ[1, 2]^2 +# elseif size(σ, 1) == 3 +# t1 = ((σ[1, 1] - σ[2, 2])^2 + (σ[2, 2] - σ[3, 3])^2 + (σ[3, 3] - σ[1, 1])^2) / 2 +# t2 = 3 * (σ[1, 2]^2 + σ[2, 3]^2 + σ[3, 1]^2) +# else +# throw(ArgumentError("Unsupported stress tensor type.")) +# end +# return sqrt(t1 + t2) +# end + +# function von_mises_stress_function(solver::AbstractFEASolver) +# st = StressTensor(solver) +# dp = Displacement(solver) +# return x -> von_mises.(st(dp(x))) +# end + + + + + + + + + + + + + + + + + + + + + + +# function FToK2AndK3(F::Vector{Matrix{Float64}}) +# k1_list=[] +# k2_list=[] +# k3_list=[] +# lambda_list=[] +# #Catches potential error of an empty F value +# if length(F)==0 +# throw("F is empty") +# end +# #Catches the error of F being the wrong size +# x=[0 0 0; 0 0 0; 0 0 0] +# for l in 1:length(F) +# if size(F[l])!=size(x) +# msg = "Error at index $l: Deformation gradient shape is not (3,3)" +# throw(ArgumentError(msg)) +# end + +# #Here begins calculations to get k1, k2, and k3 +# C=transpose(F[l])*F[l] +# #creates object R which holds eigenvalues and eigenvectors, which are then extracted +# R=eigen(C) +# lam2 = R.values +# n = R.vectors + +# Q=n +# g=[0.0,0.0,0.0] +# lam=sqrt.(lam2) +# #Tries to sort lam vector, and attempts to catch errors in trying to do so +# try +# lam=sort([lam[1],lam[2],lam[3]],rev=true) +# catch +# has_complex=false +# for s in 1:3 +# if lam2[s]<0 +# has_complex=true +# end +# end +# if has_complex==true +# throw("Lambda values include complex numbers") +# else +# throw("Lambda values cannot be sorted") +# end +# end + + +# kproduct=lam[1]*lam[2]*lam[3] +# k1=log(lam[1]*lam[2]*lam[3]) + +# for r in 1:3 +# g[r] = float(log(lam[r])-(k1/3)) +# end + +# k2=sqrt(g[1]^2+g[2]^2+g[3]^2) +# k3=3*sqrt(6)*g[1]*g[2]*g[3]/(k2^3); + +# #Adds k1, k2, and k3 valus for element to the lists of values of all elements +# push!(k1_list,k1) +# push!(k2_list,k2) +# push!(k3_list,k3) + +# map(Float64,k1_list) +# map(Float64,k2_list) +# map(Float64,k3_list) + +# end +# return(map(Float64,k1_list),map(Float64,k2_list),map(Float64,k3_list)) +# end + +# 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 + +# # function Sensitivity(F::Vector{Matrix{Float64}},alpha::Float64=2.0,mu::Float64=1.0) diff --git a/src/Functions/goodness.jl b/src/Functions/goodness.jl new file mode 100644 index 00000000..56c6d457 --- /dev/null +++ b/src/Functions/goodness.jl @@ -0,0 +1,277 @@ +function FToK2AndK3(F::Vector{Matrix{Float64}}) + k1_list=[] + k2_list=[] + k3_list=[] + lambda_list=[] + #Catches potential error of an empty F value + if length(F)==0 + throw("F is empty") + end + #Catches the error of F being the wrong size + x=[0 0 0; 0 0 0; 0 0 0] + for l in 1:length(F) + if size(F[l])!=size(x) + msg = "Error at index $l: Deformation gradient shape is not (3,3)" + throw(ArgumentError(msg)) + end + + #Here begins calculations to get k1, k2, and k3 + C=transpose(F[l])*F[l] + #creates object R which holds eigenvalues and eigenvectors, which are then extracted + R=eigen(C) + lam2 = R.values + n = R.vectors + + Q=n + g=[0.0,0.0,0.0] + lam=sqrt.(lam2) + + #Tries to sort lam vector, and attempts to catch errors in trying to do so + #Not necessary + # try + # lam=sort([lam[1],lam[2],lam[3]],rev=true) + # catch + # has_complex=false + # for s in 1:3 + # if lam2[s]<0 + # has_complex=true + # end + # end + # if has_complex==true + # throw("Lambda values include complex numbers") + # else + # throw("Lambda values cannot be sorted") + # end + # end + + + kproduct=lam[1]*lam[2]*lam[3] + k1=log(lam[1]*lam[2]*lam[3]) + + + + + # for r in 1:3 + # g[r] = float(log(lam[r])-(k1/3)) + # end + + g1=float(log(lam[1])-(k1/3)) + g2=float(log(lam[2])-(k1/3)) + g3=float(log(lam[3])-(k1/3)) + g=[g1,g2,g3] + + + + k2=sqrt(g[1]^2+g[2]^2+g[3]^2) + k3=3*sqrt(6)*g[1]*g[2]*g[3]/(k2^3); + + + + + + + + #Adds k1, k2, and k3 valus for element to the lists of values of all elements + k1_list=vcat(k1_list,k1) + k2_list=vcat(k2_list,k2) + k3_list= vcat(k3_list,k3) + + end + return(map(Float64,k1_list),map(Float64,k2_list),map(Float64,k3_list)) +end + +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 + +# function Sensitivity(F::Vector{Matrix{Float64}},alpha::Float64=2.0,mu::Float64=1.0) + +function Sensitivity(x,dispfun,defgrad) + sim1_u= dispfun(x) + F = defgrad(sim1_u); + k1,k2,k3=FToK2AndK3(F) + S=SAlpha_gen(k2,k3); + Sens_Corrected=ceil(x).*S[1] + return Sens_Corrected +end \ No newline at end of file diff --git a/src/TopOptProblems/problem_types.jl b/src/TopOptProblems/problem_types.jl index 79ebcf6e..9d8454e3 100644 --- a/src/TopOptProblems/problem_types.jl +++ b/src/TopOptProblems/problem_types.jl @@ -341,7 +341,7 @@ 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) +#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] fnode = Tuple(getnodeset(p.rect_grid.grid, "down_force"))[1] @@ -718,3 +718,99 @@ function RayProblem( end nnodespercell(p::RayProblem) = nnodespercell(p.rect_grid) getcloaddict(p::RayProblem) = p.loads + +@params struct JosephLoadCantilever{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")}, ::JosephLoadCantilever) + return println("TopOpt joseph load cantilever beam problem") +end + +function JosephLoadCantilever( + ::Type{Val{CellType}}, + nels::NTuple{dim,Int}, + sizes::NTuple{dim}, + E=1.0, + ν=0.3, + disp=1.0, +) where {dim,CellType} + iseven(nels[2]) && (length(nels) < 3 || iseven(nels[3])) || + throw("Grid does not have an even number of elements along the y and/or z axes.") + + _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 + #addfaceset!(rect_grid.grid, "fixed_all", x -> left(rect_grid, x)); + 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 = Dirichlet(:u, getfaceset(rect_grid.grid, "fixed_all"), (x,t) -> zeros(T, dim), collect(1:dim)) + 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) -> fill(T(disp), dim), collect(1:dim) + ) + add!(ch, dbc_right) + close!(ch) + t = T(0) + update!(ch, t) + + metadata = Metadata(dh) + + #fnode = Tuple(getnodeset(rect_grid.grid, "down_force"))[1] + #node_dofs = metadata.node_dofs + #force_dof = node_dofs[2, fnode] + + N = nnodespercell(rect_grid) + M = nfacespercell(rect_grid) + + black, white = find_black_and_white(dh) + varind = find_varind(black, white) + + return JosephLoadCantilever( + rect_grid, E, ν, ch, disp, black, white, varind, metadata + ) +end + +nnodespercell(p::Union{PointLoadCantilever,HalfMBB,JosephLoadCantilever}) = nnodespercell(p.rect_grid) \ No newline at end of file From da3cf0d9d04c21437d854bdc5ac46e0838d0e12e Mon Sep 17 00:00:00 2001 From: jbeckett Date: Tue, 2 Apr 2024 17:05:35 -0400 Subject: [PATCH 02/14] Add DBC Problem Type and 2D DefGrad/Goodness support --- Project.toml | 4 + src/Functions/defgrad.jl | 347 ++------------------------- src/Functions/goodness.jl | 50 ++-- src/TopOptProblems/TopOptProblems.jl | 3 +- 4 files changed, 45 insertions(+), 359 deletions(-) diff --git a/Project.toml b/Project.toml index 7a801710..ecff3bc3 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" @@ -37,6 +40,7 @@ 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" Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/src/Functions/defgrad.jl b/src/Functions/defgrad.jl index df3a34e8..06d7745c 100644 --- a/src/Functions/defgrad.jl +++ b/src/Functions/defgrad.jl @@ -67,7 +67,7 @@ function (f::ElementDefGradTensor)(u::DisplacementResult; element_dofs=false) #- _u = cellu[dim * (a - 1) .+ (1:dim)] return tensor_kernel(f, q_point, a)(DisplacementResult(_u)) end, - ) + I(3) + ) + I(dim) end @params struct ElementDefGradTensorKernel{T} <: AbstractFunction{T} # ------------------------------------------------------------------------------------------------------------[C7] @@ -78,43 +78,23 @@ end cellvalues::Any dim::Int end -# function (f::ElementDefGradTensorKernel)(u::DisplacementResult) # ----------------------------------------------------------------------------------------------------------------[C8] -# @unpack E, ν, q_point, a, cellvalues = f -# ∇ϕ = Vector(shape_gradient(cellvalues, q_point, a)) -# ϵ = (u.u .* ∇ϕ' .+ ∇ϕ .* u.u') ./ 2 -# c1 = E * ν / (1 - ν^2) * sum(diag(ϵ)) -# c2 = E * ν * (1 + ν) -# return c1 * I + c2 * ϵ -# end function (f::ElementDefGradTensorKernel)(u::DisplacementResult) # ----------------------------------------------------------------------------------------------------------------[C8] ---- nifty @unpack E, ν, q_point, a, cellvalues = f ∇ϕ = Vector(shape_gradient(cellvalues, q_point, a)) # ϵ = (u.u .* ∇ϕ' .+ ∇ϕ .* u.u') ./ 2 # c1 = E * ν / (1 - ν^2) * sum(diag(ϵ)) # c2 = E * ν * (1 + ν) - F_quad = u.u * ∇ϕ' # should this I(3) go somewhere else - return F_quad + 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(Δ))) + 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] -# return ElementDefGradTensorKernel( -# f.problem.E, -# f.problem.ν, -# quad, -# basef, -# f.cellvalues, -# TopOptProblems.getdim(f.problem), -# ) -# end - function tensor_kernel(f::DefGradTensor, quad, basef) # -------------------------------------------------------------------------------------------------------------------------[C6* Altered ] - if string(typeof(f.problem))[1:12]!="InpStiffness" + #if string(typeof(f.problem))[1:12]!="InpStiffness" return ElementDefGradTensorKernel( f.problem.E, f.problem.ν, @@ -123,314 +103,17 @@ function tensor_kernel(f::DefGradTensor, 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 + #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 - - - -# function von_mises(σ::AbstractMatrix) -# if size(σ, 1) == 2 -# t1 = σ[1, 1]^2 - σ[1, 1] * σ[2, 2] + σ[2, 2]^2 -# t2 = 3 * σ[1, 2]^2 -# elseif size(σ, 1) == 3 -# t1 = ((σ[1, 1] - σ[2, 2])^2 + (σ[2, 2] - σ[3, 3])^2 + (σ[3, 3] - σ[1, 1])^2) / 2 -# t2 = 3 * (σ[1, 2]^2 + σ[2, 3]^2 + σ[3, 1]^2) -# else -# throw(ArgumentError("Unsupported stress tensor type.")) -# end -# return sqrt(t1 + t2) -# end - -# function von_mises_stress_function(solver::AbstractFEASolver) -# st = StressTensor(solver) -# dp = Displacement(solver) -# return x -> von_mises.(st(dp(x))) -# end - - - - - - - - - - - - - - - - - - - - - - -# function FToK2AndK3(F::Vector{Matrix{Float64}}) -# k1_list=[] -# k2_list=[] -# k3_list=[] -# lambda_list=[] -# #Catches potential error of an empty F value -# if length(F)==0 -# throw("F is empty") -# end -# #Catches the error of F being the wrong size -# x=[0 0 0; 0 0 0; 0 0 0] -# for l in 1:length(F) -# if size(F[l])!=size(x) -# msg = "Error at index $l: Deformation gradient shape is not (3,3)" -# throw(ArgumentError(msg)) -# end - -# #Here begins calculations to get k1, k2, and k3 -# C=transpose(F[l])*F[l] -# #creates object R which holds eigenvalues and eigenvectors, which are then extracted -# R=eigen(C) -# lam2 = R.values -# n = R.vectors - -# Q=n -# g=[0.0,0.0,0.0] -# lam=sqrt.(lam2) -# #Tries to sort lam vector, and attempts to catch errors in trying to do so -# try -# lam=sort([lam[1],lam[2],lam[3]],rev=true) -# catch -# has_complex=false -# for s in 1:3 -# if lam2[s]<0 -# has_complex=true -# end -# end -# if has_complex==true -# throw("Lambda values include complex numbers") -# else -# throw("Lambda values cannot be sorted") -# end -# end - - -# kproduct=lam[1]*lam[2]*lam[3] -# k1=log(lam[1]*lam[2]*lam[3]) - -# for r in 1:3 -# g[r] = float(log(lam[r])-(k1/3)) -# end - -# k2=sqrt(g[1]^2+g[2]^2+g[3]^2) -# k3=3*sqrt(6)*g[1]*g[2]*g[3]/(k2^3); - -# #Adds k1, k2, and k3 valus for element to the lists of values of all elements -# push!(k1_list,k1) -# push!(k2_list,k2) -# push!(k3_list,k3) - -# map(Float64,k1_list) -# map(Float64,k2_list) -# map(Float64,k3_list) - -# end -# return(map(Float64,k1_list),map(Float64,k2_list),map(Float64,k3_list)) -# end - -# 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 - -# # function Sensitivity(F::Vector{Matrix{Float64}},alpha::Float64=2.0,mu::Float64=1.0) +end \ No newline at end of file diff --git a/src/Functions/goodness.jl b/src/Functions/goodness.jl index 56c6d457..98b56035 100644 --- a/src/Functions/goodness.jl +++ b/src/Functions/goodness.jl @@ -1,22 +1,31 @@ function FToK2AndK3(F::Vector{Matrix{Float64}}) - k1_list=[] - k2_list=[] - k3_list=[] - lambda_list=[] + k1_list=Float64[] + k2_list=Float64[] + k3_list=Float64[] + lambda_list=Float64[] #Catches potential error of an empty F value - if length(F)==0 - throw("F is empty") + #if length(F)==0 + # throw("F is empty") + #end + + function Fto3by3(F) + return Vector{Matrix{Float64}}(map(eachindex(F)) do i + F_tmp = F[i] + if size(F_tmp) == (2, 2) + # Build an immutable 3x3 matrix from the 2x2 F_item + [F_tmp zeros(2, 1); 0 0 1] + elseif size(F_tmp) == (3, 3) + F_tmp + else + error("Unexpected deformation gradient size at index $l") + end + end) end - #Catches the error of F being the wrong size - x=[0 0 0; 0 0 0; 0 0 0] - for l in 1:length(F) - if size(F[l])!=size(x) - msg = "Error at index $l: Deformation gradient shape is not (3,3)" - throw(ArgumentError(msg)) - end + F3by3 = Fto3by3(F) # 3 by 3 + for l in 1:length(F3by3) #Here begins calculations to get k1, k2, and k3 - C=transpose(F[l])*F[l] + C=transpose(F3by3[l])*F3by3[l] #creates object R which holds eigenvalues and eigenvectors, which are then extracted R=eigen(C) lam2 = R.values @@ -48,9 +57,6 @@ function FToK2AndK3(F::Vector{Matrix{Float64}}) kproduct=lam[1]*lam[2]*lam[3] k1=log(lam[1]*lam[2]*lam[3]) - - - # for r in 1:3 # g[r] = float(log(lam[r])-(k1/3)) # end @@ -60,16 +66,8 @@ function FToK2AndK3(F::Vector{Matrix{Float64}}) g3=float(log(lam[3])-(k1/3)) g=[g1,g2,g3] - - k2=sqrt(g[1]^2+g[2]^2+g[3]^2) - k3=3*sqrt(6)*g[1]*g[2]*g[3]/(k2^3); - - - - - - + k3=3*sqrt(6)*g[1]*g[2]*g[3]/(k2^3) #Adds k1, k2, and k3 valus for element to the lists of values of all elements k1_list=vcat(k1_list,k1) diff --git a/src/TopOptProblems/TopOptProblems.jl b/src/TopOptProblems/TopOptProblems.jl index 6fc7ee7b..3edab198 100644 --- a/src/TopOptProblems/TopOptProblems.jl +++ b/src/TopOptProblems/TopOptProblems.jl @@ -49,6 +49,7 @@ export RayProblem, bcmatrix, save_mesh, RandomMagnitude, - MultiLoad + MultiLoad, + JosephLoadCantilever end # module From baa5295ad33c066db1bbe4e05f1b0228bd5b2101 Mon Sep 17 00:00:00 2001 From: jbeckett Date: Thu, 25 Apr 2024 15:56:46 -0400 Subject: [PATCH 03/14] Clean up TensionBar problem type Converted the name from JosephLoadCantilever to TensionBar, added a schematic for the documentation, and deleted a requirements related to enforcing an even number of elements along the y/z directions that was unnecessary. --- src/TopOptProblems/TopOptProblems.jl | 2 +- src/TopOptProblems/problem_types.jl | 70 +++++++++++++++++++--------- 2 files changed, 50 insertions(+), 22 deletions(-) diff --git a/src/TopOptProblems/TopOptProblems.jl b/src/TopOptProblems/TopOptProblems.jl index 3edab198..f22a1d37 100644 --- a/src/TopOptProblems/TopOptProblems.jl +++ b/src/TopOptProblems/TopOptProblems.jl @@ -50,6 +50,6 @@ export RayProblem, save_mesh, RandomMagnitude, MultiLoad, - JosephLoadCantilever + TensionBar end # module diff --git a/src/TopOptProblems/problem_types.jl b/src/TopOptProblems/problem_types.jl index 9d8454e3..a054e9fc 100644 --- a/src/TopOptProblems/problem_types.jl +++ b/src/TopOptProblems/problem_types.jl @@ -719,7 +719,43 @@ end nnodespercell(p::RayProblem) = nnodespercell(p.rect_grid) getcloaddict(p::RayProblem) = p.loads -@params struct JosephLoadCantilever{dim,T,N,M} <: StiffnessTopOptProblem{dim,T} +""" +``` +///**********************************-> +///* *-> +///* *-> 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 @@ -731,21 +767,18 @@ getcloaddict(p::RayProblem) = p.loads metadata::Metadata end -function Base.show(::IO, ::MIME{Symbol("text/plain")}, ::JosephLoadCantilever) - return println("TopOpt joseph load cantilever beam problem") +function Base.show(::IO, ::MIME{Symbol("text/plain")}, ::TensionBar) + return println("TopOpt tension bar problem") end -function JosephLoadCantilever( +function TensionBar( ::Type{Val{CellType}}, nels::NTuple{dim,Int}, sizes::NTuple{dim}, - E=1.0, - ν=0.3, - disp=1.0, + E, + ν, + disp, ) where {dim,CellType} - iseven(nels[2]) && (length(nels) < 3 || iseven(nels[3])) || - throw("Grid does not have an even number of elements along the y and/or z axes.") - _T = promote_type(eltype(sizes), typeof(E), typeof(ν), typeof(disp)) if _T <: Integer T = Float64 @@ -767,9 +800,7 @@ function JosephLoadCantilever( 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) - ) + addnodeset!(rect_grid.grid, "disp_right", x -> right(rect_grid, x)) # Create displacement field u dh = DofHandler(rect_grid.grid) @@ -784,13 +815,10 @@ function JosephLoadCantilever( ch = ConstraintHandler(dh) #dbc = Dirichlet(:u, getfaceset(rect_grid.grid, "fixed_all"), (x,t) -> zeros(T, dim), collect(1:dim)) - dbc_left = Dirichlet( - :u, getnodeset(rect_grid.grid, "fixed_left"), (x, t) -> zeros(T, dim), collect(1:dim) - ) + 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) -> fill(T(disp), dim), collect(1:dim) - ) + #dbc_right = Dirichlet(:u, getnodeset(rect_grid.grid, "disp_right"), (x, t) -> fill(T(disp), dim), collect(1:dim)) + dbc_right = Dirichlet(:u, getnodeset(rect_grid.grid, "disp_right"), (x, t) -> T[disp; zeros(dim-1)], collect(1:dim)) add!(ch, dbc_right) close!(ch) t = T(0) @@ -808,9 +836,9 @@ function JosephLoadCantilever( black, white = find_black_and_white(dh) varind = find_varind(black, white) - return JosephLoadCantilever( + return TensionBar( rect_grid, E, ν, ch, disp, black, white, varind, metadata ) end -nnodespercell(p::Union{PointLoadCantilever,HalfMBB,JosephLoadCantilever}) = nnodespercell(p.rect_grid) \ No newline at end of file +nnodespercell(p::Union{PointLoadCantilever,HalfMBB,TensionBar}) = nnodespercell(p.rect_grid) \ No newline at end of file From 32ba84e71c3dfb5b13985cb3854336cf3c103598 Mon Sep 17 00:00:00 2001 From: jbeckett Date: Thu, 25 Apr 2024 16:33:11 -0400 Subject: [PATCH 04/14] Fixed quadrature formulation in defgrad Fixes bug where all of the quadrature points were not iterated over all of the basis functions. Also adds missing quadrature weight and determinant of the Jacobian to the formulation --- src/Functions/defgrad.jl | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/src/Functions/defgrad.jl b/src/Functions/defgrad.jl index 06d7745c..d6400df1 100644 --- a/src/Functions/defgrad.jl +++ b/src/Functions/defgrad.jl @@ -62,12 +62,13 @@ function (f::ElementDefGradTensor)(u::DisplacementResult; element_dofs=false) #- n_basefuncs = getnbasefunctions(st.cellvalues) n_quad = getnquadpoints(st.cellvalues) dim = TopOptProblems.getdim(st.problem) - return sum( - map(1:n_basefuncs, 1:n_quad) do a, q_point + return 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, - ) + I(dim) + end) * dΩ + end) + I(dim) end @params struct ElementDefGradTensorKernel{T} <: AbstractFunction{T} # ------------------------------------------------------------------------------------------------------------[C7] @@ -81,9 +82,6 @@ end function (f::ElementDefGradTensorKernel)(u::DisplacementResult) # ----------------------------------------------------------------------------------------------------------------[C8] ---- nifty @unpack E, ν, q_point, a, cellvalues = f ∇ϕ = Vector(shape_gradient(cellvalues, q_point, a)) - # ϵ = (u.u .* ∇ϕ' .+ ∇ϕ .* u.u') ./ 2 - # c1 = E * ν / (1 - ν^2) * sum(diag(ϵ)) - # c2 = E * ν * (1 + ν) return u.u * ∇ϕ' end function ChainRulesCore.rrule(f::ElementDefGradTensorKernel, u::DisplacementResult) From 2cdeae5523c76c7b4247c1267b0f307c8e5b202e Mon Sep 17 00:00:00 2001 From: jbeckett Date: Fri, 26 Apr 2024 11:44:48 -0400 Subject: [PATCH 05/14] Update .gitignore for dev work Ignore .ipynb, .jl, and .vtu files frequently developed during testing. --- .gitignore | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index 281f6a83..49832ec0 100644 --- a/.gitignore +++ b/.gitignore @@ -14,4 +14,8 @@ Manifest.toml */__pycache__ -examples/benchmark/iterations \ No newline at end of file +examples/benchmark/iterations + +*.ipynb +*.jl +*.vtu \ No newline at end of file From ee1ea1a113848151f963500e54b1128362581545 Mon Sep 17 00:00:00 2001 From: jbeckett Date: Fri, 10 May 2024 13:23:40 -0400 Subject: [PATCH 06/14] Add volume correction in defgrad --- src/Functions/defgrad.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Functions/defgrad.jl b/src/Functions/defgrad.jl index d6400df1..467f11a0 100644 --- a/src/Functions/defgrad.jl +++ b/src/Functions/defgrad.jl @@ -62,13 +62,14 @@ function (f::ElementDefGradTensor)(u::DisplacementResult; element_dofs=false) #- n_basefuncs = getnbasefunctions(st.cellvalues) n_quad = getnquadpoints(st.cellvalues) dim = TopOptProblems.getdim(st.problem) + V = sum(st.cellvalues.detJdV) return 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) + I(dim) + end) ./ V + I(dim) end @params struct ElementDefGradTensorKernel{T} <: AbstractFunction{T} # ------------------------------------------------------------------------------------------------------------[C7] From 5a93b469f8bc1c0102e9fe558018b4b705d0959f Mon Sep 17 00:00:00 2001 From: jbeckett Date: Mon, 17 Jun 2024 13:46:43 -0400 Subject: [PATCH 07/14] Add TensionRoller problem type Added a tension bar problem type where x displacement are the only enforced condition on the bar. Two additional points have contraints to prevent translation of the sample and to prevent rigid body rotation about the x-axis. --- src/TopOptProblems/TopOptProblems.jl | 3 +- src/TopOptProblems/problem_types.jl | 145 ++++++++++++++++++++++++--- 2 files changed, 135 insertions(+), 13 deletions(-) diff --git a/src/TopOptProblems/TopOptProblems.jl b/src/TopOptProblems/TopOptProblems.jl index f22a1d37..a8b1ec8b 100644 --- a/src/TopOptProblems/TopOptProblems.jl +++ b/src/TopOptProblems/TopOptProblems.jl @@ -50,6 +50,7 @@ export RayProblem, save_mesh, RandomMagnitude, MultiLoad, - TensionBar + TensionBar, + TensionRoller end # module diff --git a/src/TopOptProblems/problem_types.jl b/src/TopOptProblems/problem_types.jl index a054e9fc..ec0cb740 100644 --- a/src/TopOptProblems/problem_types.jl +++ b/src/TopOptProblems/problem_types.jl @@ -341,7 +341,6 @@ 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] fnode = Tuple(getnodeset(p.rect_grid.grid, "down_force"))[1] @@ -716,7 +715,6 @@ 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 """ @@ -794,7 +792,6 @@ function TensionBar( if haskey(rect_grid.grid.facesets, "fixed_left") pop!(rect_grid.grid.facesets, "fixed_left") end - #addfaceset!(rect_grid.grid, "fixed_all", x -> left(rect_grid, x)); addnodeset!(rect_grid.grid, "fixed_left", x -> left(rect_grid, x)) if haskey(rect_grid.grid.nodesets, "disp_right") @@ -814,10 +811,8 @@ function TensionBar( ch = ConstraintHandler(dh) - #dbc = Dirichlet(:u, getfaceset(rect_grid.grid, "fixed_all"), (x,t) -> zeros(T, dim), collect(1:dim)) 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) -> fill(T(disp), dim), collect(1:dim)) dbc_right = Dirichlet(:u, getnodeset(rect_grid.grid, "disp_right"), (x, t) -> T[disp; zeros(dim-1)], collect(1:dim)) add!(ch, dbc_right) close!(ch) @@ -826,9 +821,137 @@ function TensionBar( metadata = Metadata(dh) - #fnode = Tuple(getnodeset(rect_grid.grid, "down_force"))[1] - #node_dofs = metadata.node_dofs - #force_dof = node_dofs[2, fnode] + 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[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(0) + update!(ch, t) + + metadata = Metadata(dh) N = nnodespercell(rect_grid) M = nfacespercell(rect_grid) @@ -836,9 +959,7 @@ function TensionBar( black, white = find_black_and_white(dh) varind = find_varind(black, white) - return TensionBar( - rect_grid, E, ν, ch, disp, black, white, varind, metadata - ) + return TensionRoller(rect_grid, E, ν, ch, disp, black, white, varind, metadata) end -nnodespercell(p::Union{PointLoadCantilever,HalfMBB,TensionBar}) = nnodespercell(p.rect_grid) \ No newline at end of file +nnodespercell(p::Union{PointLoadCantilever,HalfMBB,RayProblem,TensionBar,TensionRoller}) = nnodespercell(p.rect_grid) \ No newline at end of file From 4e408a4fa69a4c3f7448d51e7c537968e24d720f Mon Sep 17 00:00:00 2001 From: jbeckett Date: Fri, 5 Jul 2024 13:56:09 -0400 Subject: [PATCH 08/14] Add plane strain components 2D defgrad tensors All FEA currently done in TopOpt assumes plane strain. This small fix just clears up the ambiguous 3D terms of defgrad tensors of 2D problem types (i.e., F33=0 and F31=F32=F13=23=0) --- src/Functions/defgrad.jl | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/src/Functions/defgrad.jl b/src/Functions/defgrad.jl index 467f11a0..90a47b81 100644 --- a/src/Functions/defgrad.jl +++ b/src/Functions/defgrad.jl @@ -51,6 +51,16 @@ function ChainRulesCore.rrule(::typeof(reinit!), st::ElementDefGradTensor, celli 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 @@ -63,13 +73,18 @@ function (f::ElementDefGradTensor)(u::DisplacementResult; element_dofs=false) #- n_quad = getnquadpoints(st.cellvalues) dim = TopOptProblems.getdim(st.problem) V = sum(st.cellvalues.detJdV) - return sum(map(1:n_quad) do q_point + 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] From e5e1ed9c8051c88c52b83cb2804c9d712df5f6b0 Mon Sep 17 00:00:00 2001 From: jbeckett Date: Fri, 5 Jul 2024 14:36:46 -0400 Subject: [PATCH 09/14] Clean and restructure goodness.jl Add detailed comments and make code more concise. The primary focus of the edits were to the sensitivity related functions (i.e., not the entropy related functions). --- Project.toml | 1 + src/Functions/Functions.jl | 10 +- src/Functions/goodness.jl | 213 +++++++++++++++++++++++-------------- 3 files changed, 139 insertions(+), 85 deletions(-) diff --git a/Project.toml b/Project.toml index ecff3bc3..08f7a8dc 100644 --- a/Project.toml +++ b/Project.toml @@ -42,6 +42,7 @@ 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/Functions/Functions.jl b/src/Functions/Functions.jl index 84456461..cebf8dca 100644 --- a/src/Functions/Functions.jl +++ b/src/Functions/Functions.jl @@ -14,6 +14,7 @@ using Nonconvex: Nonconvex using Flux using AbstractDifferentiation: AbstractDifferentiation using StatsBase, Statistics +using Symbolics using Plots:heatmap const AD = AbstractDifferentiation @@ -49,14 +50,15 @@ export Volume, element_densities, tounit, DefGradTensor, - GoodnessTensor, ElementDefGradTensor, - ElementGoodnessTensor, - FToK2AndK3, + FToK123, + orth_decomp, + sensitivityFieldFncs, + sensitivityPVW, Entropy_Calc, Entropy, SMu_gen, - SAlpha_gen + SAlpha_gen const to = TimerOutput() diff --git a/src/Functions/goodness.jl b/src/Functions/goodness.jl index 98b56035..04198e51 100644 --- a/src/Functions/goodness.jl +++ b/src/Functions/goodness.jl @@ -1,83 +1,149 @@ -function FToK2AndK3(F::Vector{Matrix{Float64}}) - k1_list=Float64[] - k2_list=Float64[] - k3_list=Float64[] - lambda_list=Float64[] - #Catches potential error of an empty F value - #if length(F)==0 - # throw("F is empty") - #end +#= 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 Fto3by3(F) - return Vector{Matrix{Float64}}(map(eachindex(F)) do i - F_tmp = F[i] - if size(F_tmp) == (2, 2) - # Build an immutable 3x3 matrix from the 2x2 F_item - [F_tmp zeros(2, 1); 0 0 1] - elseif size(F_tmp) == (3, 3) - F_tmp - else - error("Unexpected deformation gradient size at index $l") - end - end) - end - F3by3 = Fto3by3(F) # 3 by 3 +function FToK123(F::AbstractMatrix{T}) where {T} + @assert size(F) == (3, 3) + K = Vector{T}(undef,3) + g = Vector{T}(undef,3) + + C = transpose(F)*F + λ = sqrt.(eigen(C).values) + + K[1] = log(λ[1]*λ[2]*λ[3]) + g = log.(λ) .- (K[1]/3) + K[2] = sqrt(sum(g.^2)) + K[3] = 3*sqrt(6)*g[1]*g[2]*g[3]/(K[2]^3) + return 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) - for l in 1:length(F3by3) - #Here begins calculations to get k1, k2, and k3 - C=transpose(F3by3[l])*F3by3[l] - #creates object R which holds eigenvalues and eigenvectors, which are then extracted - R=eigen(C) - lam2 = R.values - n = R.vectors + 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 - Q=n - g=[0.0,0.0,0.0] - lam=sqrt.(lam2) +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 - #Tries to sort lam vector, and attempts to catch errors in trying to do so - #Not necessary - # try - # lam=sort([lam[1],lam[2],lam[3]],rev=true) - # catch - # has_complex=false - # for s in 1:3 - # if lam2[s]<0 - # has_complex=true - # end - # end - # if has_complex==true - # throw("Lambda values include complex numbers") - # else - # throw("Lambda values cannot be sorted") - # end - # 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)) - kproduct=lam[1]*lam[2]*lam[3] - k1=log(lam[1]*lam[2]*lam[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 - # for r in 1:3 - # g[r] = float(log(lam[r])-(k1/3)) - # 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 - g1=float(log(lam[1])-(k1/3)) - g2=float(log(lam[2])-(k1/3)) - g3=float(log(lam[3])-(k1/3)) - g=[g1,g2,g3] +# 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! - k2=sqrt(g[1]^2+g[2]^2+g[3]^2) - k3=3*sqrt(6)*g[1]*g[2]*g[3]/(k2^3) +# Outputs +# sens_IVW (sensitivity metric across all elements) - #Adds k1, k2, and k3 valus for element to the lists of values of all elements - k1_list=vcat(k1_list,k1) - k2_list=vcat(k2_list,k2) - k3_list= vcat(k3_list,k3) +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(map(Float64,k1_list),map(Float64,k2_list),map(Float64,k3_list)) + 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)) @@ -153,9 +219,6 @@ function Entropy(F::Vector{Vector{Matrix{Float64}}},n_bins::Int64=100,offset::Fl 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] @@ -219,7 +282,6 @@ function Entropy(F::Vector{Vector{Matrix{Float64}}},n_bins::Int64=100,offset::Fl 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=[] @@ -261,15 +323,4 @@ function SAlpha_gen(k2_list::Vector{Float64},k3_list::Vector{Float64},alpha::Flo end Sum=sum(SAlpha_List) return (map(Float64,SAlpha_List)), Sum -end - -# function Sensitivity(F::Vector{Matrix{Float64}},alpha::Float64=2.0,mu::Float64=1.0) - -function Sensitivity(x,dispfun,defgrad) - sim1_u= dispfun(x) - F = defgrad(sim1_u); - k1,k2,k3=FToK2AndK3(F) - S=SAlpha_gen(k2,k3); - Sens_Corrected=ceil(x).*S[1] - return Sens_Corrected end \ No newline at end of file From a4511d4240d30deb3799e46a056d37c24e090ab2 Mon Sep 17 00:00:00 2001 From: jbeckett Date: Fri, 5 Jul 2024 15:07:51 -0400 Subject: [PATCH 10/14] Patch for autodiff compatibility of FToK123() Making a simple indexing fix to fix a mutation error from Zygote when calling autodiff. --- src/Functions/goodness.jl | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/Functions/goodness.jl b/src/Functions/goodness.jl index 04198e51..b8cf8c58 100644 --- a/src/Functions/goodness.jl +++ b/src/Functions/goodness.jl @@ -7,17 +7,16 @@ 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) - K = Vector{T}(undef,3) g = Vector{T}(undef,3) C = transpose(F)*F λ = sqrt.(eigen(C).values) - K[1] = log(λ[1]*λ[2]*λ[3]) - g = log.(λ) .- (K[1]/3) - K[2] = sqrt(sum(g.^2)) - K[3] = 3*sqrt(6)*g[1]*g[2]*g[3]/(K[2]^3) - return K + 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} From b04ad393440f11a6657107c3d115ba4887e54a88 Mon Sep 17 00:00:00 2001 From: jbeckett Date: Mon, 5 Aug 2024 17:30:32 -0400 Subject: [PATCH 11/14] Hyperelastic upgrade draft Adding basic fundamentals of a hyperelastic solver. Still needs plenty of cleanings and comments. The main issue currently is autodiff still needs to be upgraded to use ImplicitDiff.jl to handle the iterative nature of the solver. The code still needs additional validation, but the bones are there. --- src/FEA/FEA.jl | 3 + src/FEA/iterative_hyperelastic_solver.jl | 404 +++++++++++++++++++++ src/FEA/solvers_api.jl | 5 + src/TopOpt.jl | 1 + src/TopOptProblems/TopOptProblems.jl | 5 +- src/TopOptProblems/assemble.jl | 80 +++- src/TopOptProblems/elementinfo.jl | 59 +++ src/TopOptProblems/elementmatrix.jl | 5 +- src/TopOptProblems/matrices_and_vectors.jl | 306 +++++++++++++++- 9 files changed, 863 insertions(+), 5 deletions(-) create mode 100644 src/FEA/iterative_hyperelastic_solver.jl diff --git a/src/FEA/FEA.jl b/src/FEA/FEA.jl index 0c001d4e..5996dfef 100644 --- a/src/FEA/FEA.jl +++ b/src/FEA/FEA.jl @@ -12,10 +12,12 @@ export AbstractFEASolver, DirectDisplacementSolver, PCGDisplacementSolver, StaticMatrixFreeDisplacementSolver, + HyperelasticDisplacementSolver, Displacement, Direct, CG, MatrixFree, + Hyperelastic, FEASolver, Assembly, DefaultCriteria, @@ -34,6 +36,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..04ee7576 --- /dev/null +++ b/src/FEA/iterative_hyperelastic_solver.jl @@ -0,0 +1,404 @@ +abstract type AbstractHyperelasticSolver <: AbstractDisplacementSolver end + +mutable struct HyperelasticDisplacementSolver{ + T, + dim, + TP1<:AbstractPenalty{T}, + TP2<:StiffnessTopOptProblem{dim,T}, + TG<:GlobalFEAInfo_hyperelastic{T}, + TE<:ElementFEAInfo_hyperelastic{dim,T}, + Tu<:AbstractVector{T}, +} <: AbstractHyperelasticSolver + mp + problem::TP2 + globalinfo::TG + elementinfo::TE + u::Tu # JGB: u --> u0 + lhs::Tu + rhs::Tu + vars::Tu + penalty::TP1 + prev_penalty::TP1 + xmin::T + qr::Bool +end +function Base.show(::IO, ::MIME{Symbol("text/plain")}, x::HyperelasticDisplacementSolver) + return println("TopOpt 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), + qr=false, +) where {dim,T} + u = zeros(T, ndofs(sp.ch.dh)) + apply!(u,sp.ch) # apply dbc for initial guess + elementinfo = ElementFEAInfo_hyperelastic(mp, sp, u, quad_order, Val{:Static}) # JGB: add u + globalinfo = GlobalFEAInfo_hyperelastic(sp) # JGB: small issue this leads to symmetric K initialization + #u = zeros(T, ndofs(sp.ch.dh)) # JGB + lhs = similar(u) + rhs = similar(u) + vars = fill(one(T), getncells(sp.ch.dh.grid) - sum(sp.black) - sum(sp.white)) + varind = sp.varind + return HyperelasticDisplacementSolver( + mp, sp, globalinfo, elementinfo, u, lhs, rhs, vars, penalty, prev_penalty, xmin, qr # u to u0 + ) +end +function (s::HyperelasticDisplacementSolver{T})( + ::Type{Val{safe}}=Val{false}, + ::Type{newT}=T; + assemble_f=true, + reuse_fact=false, + #rhs=assemble_f ? s.globalinfo.f : s.rhs, + rhs=assemble_f ? s.globalinfo.g : s.rhs, + lhs=assemble_f ? s.u : s.lhs, + kwargs..., +) where {T,safe,newT} + #=globalinfo = s.globalinfo + assemble!( + globalinfo, + s.problem, + s.elementinfo, + s.vars, + getpenalty(s), + s.xmin; + assemble_f=assemble_f, + ) + K = globalinfo.K + if safe + m = meandiag(K) + for i in 1:size(K, 1) + if K[i, i] ≈ zero(T) + K[i, i] = m + end + end + end + nans = false + if !reuse_fact + newK = T === newT ? K : newT.(K) + if s.qr + globalinfo.qrK = qr(newK.data) + else + cholK = cholesky(Symmetric(K); check=false) + if issuccess(cholK) + globalinfo.cholK = cholK + else + @warn "The global stiffness matrix is not positive definite. Please check your boundary conditions." + lhs .= T(NaN) + nans = true + end + end + end + nans && return nothing + new_rhs = T === newT ? rhs : newT.(rhs) + fact = s.qr ? globalinfo.qrK : globalinfo.cholK + lhs .= fact \ new_rhs + =# + # NEW STUFF + globalinfo = s.globalinfo + #assemble_hyperelastic!( + # globalinfo, + # s.problem, + # s.elementinfo, + # s.vars, + # getpenalty(s), + # s.xmin; + # assemble_f=assemble_f, + #) + #K = globalinfo.K + dh = s.problem.ch.dh + ch = s.problem.ch + + # Pre-allocation of vectors for the solution and Newton increments + _ndofs = ndofs(dh) + un = s.u + #un = zeros(_ndofs) # previous solution vector + u = zeros(_ndofs) + #u = s.u + Δu = zeros(_ndofs) + ΔΔu = zeros(_ndofs) + #apply!(un, ch) + + # Create sparse matrix and residual vector + #K = create_sparsity_pattern(dh) + #K = globalinfo.K + #g = zeros(_ndofs) + #g = globalinfo.g + + # Perform Newton iterations + newton_itr = -1 + NEWTON_TOL = 1e-8 + NEWTON_MAXITER = 1000 # OG is 30 + #prog = ProgressMeter.ProgressThresh(NEWTON_TOL, "Solving:") + + while true; newton_itr += 1 + # Construct the current guess + u .= un .+ Δu + # Compute residual and tangent for current guess + s.elementinfo = ElementFEAInfo_hyperelastic(s.mp, s.problem, u, default_quad_order(s.problem), Val{:Static}) # JGB: add u + #s.globalinfo = GlobalFEAInfo_hyperelastic(s.problem) # JGB: small issue this leads to symmetric K initialization + #assemble_global!(K, g, dh, cv, fv, mp, u, ΓN) + assemble_hyperelastic!(s.globalinfo,s.problem,s.elementinfo,s.vars,getpenalty(s),s.xmin,assemble_f=assemble_f) + #K = s.globalinfo.K + #g = s.globalinfo.g + # Apply boundary conditions + #apply_zero!(s.globalinfo.K, s.globalinfo.g, ch) + # Compute the residual norm and compare with tolerance + normg = norm(s.globalinfo.g) + if normg < NEWTON_TOL + break + elseif newton_itr > NEWTON_MAXITER + error("Reached maximum Newton iterations, aborting") + end + + # Compute increment using conjugate gradients + IterativeSolvers.cg!(ΔΔu, s.globalinfo.K, s.globalinfo.g; maxiter=1000) + + apply_zero!(ΔΔu, ch) + Δu .-= ΔΔu #Δu = Δu - Δ(Δu) + end + + s.u .= u #, F_storage, F_avg + return nothing +end + + + + + + + + + + + + + + + + + + +# ALL OF THE NEW STUFF TO BE PUT PLACES +# ATTRIBUTION: the following code was primarily sourced from sample code provided in Ferrite.jl documentation +#= +struct NeoHooke # JGB: add to .ipynb + μ::Float64 + λ::Float64 +end + +function Ψ(C, mp::NeoHooke) # JGB: add to .ipynb + μ = mp.μ + λ = mp.λ + Ic = tr(C) + J = sqrt(det(C)) + return μ / 2 * (Ic - 3) - μ * log(J) + λ / 2 * log(J)^2 +end + +function constitutive_driver(C, mp::NeoHooke) + # 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; + +function assemble_element!(ke, ge, cell, cv, fv, mp, ue, ΓN, F_storage, F_avg) + # Reinitialize cell values, and reset output arrays + reinit!(cv, cell) + fill!(ke, 0.0) + fill!(ge, 0.0) + + b = Vec{3}((0.0, -0.5, 0.0)) # Body force + tn = 0.1 # Traction (to be scaled with surface normal) + ndofs = getnbasefunctions(cv) + cell_id=cell.current_cellid.x + + for qp in 1:getnquadpoints(cv) + dΩ = getdetJdV(cv, qp) + # Compute deformation gradient F and right Cauchy-Green tensor C + ∇u = function_gradient(cv, qp, ue) + F = one(∇u) + ∇u #*note that one() is the identity tensor times whatever you put in + C = tdot(F) # F' ⋅ F + # Compute stress and tangent + S, ∂S∂C = constitutive_driver(C, mp) + P = F ⋅ S + I = one(S) + ∂P∂F = otimesu(I, S) + 2 * otimesu(F, I) ⊡ ∂S∂C ⊡ otimesu(F', I) + + # Store F in our data structure + #F_storage[(cell.current_cellid, qp)] = F + F_storage[cell_id][qp] = F + + F_avg[cell_id] = F_avg[cell_id] + cv.qr_weights[1]*F + + # Loop over test functions + for i in 1:ndofs + # Test function and gradient + δui = shape_value(cv, qp, i) + ∇δui = shape_gradient(cv, qp, i) + # Add contribution to the residual from this test function + ge[i] += ( ∇δui ⊡ P - δui ⋅ b ) * dΩ + + ∇δui∂P∂F = ∇δui ⊡ ∂P∂F # Hoisted computation + for j in 1:ndofs + ∇δuj = shape_gradient(cv, qp, j) + # Add contribution to the tangent + ke[i, j] += ( ∇δui∂P∂F ⊡ ∇δuj ) * dΩ + end + end + end + + # Surface integral for the traction + for face in 1:nfaces(cell) # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! dloads start + if (cellid(cell), face) in ΓN + reinit!(fv, cell, face) + for q_point in 1:getnquadpoints(fv) + t = tn * getnormal(fv, q_point) + dΓ = getdetJdV(fv, q_point) + for i in 1:ndofs + δui = shape_value(fv, q_point, i) + ge[i] -= (δui ⋅ t) * dΓ + end + end + end + end +end; + +function assemble_global!(K, g, dh, cv, fv, mp, u, ΓN, F_storage, F_avg) + n = ndofs_per_cell(dh) + ke = zeros(n, n) + ge = zeros(n) + + # start_assemble resets K and g + assembler = start_assemble(K, g) + + # Loop over all cells in the grid + for cell in CellIterator(dh) + global_dofs = celldofs(cell) + ue = u[global_dofs] # element dofs + assemble_element!(ke, ge, cell, cv, fv, mp, ue, ΓN, F_storage, F_avg) + assemble!(assembler, global_dofs, ge, ke) + end +end; + +function solve() + # Generate a grid + N = 10 + L = 1.0 + left = zero(Vec{3}) + right = L * ones(Vec{3}) + grid = generate_grid(Tetrahedron, (N, N, N), left, right) + + # Material parameters + E = 10.0 + ν = 0.3 + μ = E / (2(1 + ν)) + λ = (E * ν) / ((1 + ν) * (1 - 2ν)) + mp = NeoHooke(μ, λ) + + # Finite element base + ip = Lagrange{3, RefTetrahedron, 1}() + qr = QuadratureRule{3, RefTetrahedron}(1) + qr_face = QuadratureRule{2, RefTetrahedron}(1) + cv = CellVectorValues(qr, ip) + fv = FaceVectorValues(qr_face, ip) + + # DofHandler + dh = DofHandler(grid) + push!(dh, :u, 3) # Add a displacement field + close!(dh) + #F_storage = Dict{Tuple{Int, Int}, Tensor{2,3,Float64}}() + ncells = Ferrite.getncells(dh.grid) + nqpoints_per_cell = getnquadpoints(cv) + #F_storage = [Vector{Matrix{Float64}(undef,3,3)}(undef, nqpoints_per_cell) for _ in 1:ncells] + #F_storage = [Vector{Tensor{2,3,Float64}}(undef, nqpoints_per_cell) for _ in 1:ncells] + + # Creating a nested array where each cell will hold an array of Tensors + F_storage = Vector{Vector{Tensor{2, 3, Float64}}}(undef, ncells) + # Initialize the inner arrays in the nested array for each cell + for i in 1:ncells + F_storage[i] = Vector{Tensor{2, 3, Float64}}(undef, nqpoints_per_cell) + end + + F_avg = Vector{Tensor{2, 3, Float64}}(undef, ncells) # Array for averaged deformation gradients. + for i in 1:ncells + F_avg[i] = zero(Tensor{2, 3}) + end + + function rotation(X, t) + θ = pi / 3 # 60° + x, y, z = X + return t * Vec{3}(( + 0.0, + L/2 - y + (y-L/2)*cos(θ) - (z-L/2)*sin(θ), + L/2 - z + (y-L/2)*sin(θ) + (z-L/2)*cos(θ) + )) + end + + dbcs = ConstraintHandler(dh) + # Add a homogeneous boundary condition on the "clamped" edge + dbc = Dirichlet(:u, getfaceset(grid, "right"), (x,t) -> [0.0, 0.0, 0.0], [1, 2, 3]) + add!(dbcs, dbc) + dbc = Dirichlet(:u, getfaceset(grid, "left"), (x,t) -> rotation(x, t), [1, 2, 3]) + add!(dbcs, dbc) + close!(dbcs) + t = 0.5 + Ferrite.update!(dbcs, t) + + # Neumann part of the boundary + ΓN = union( + getfaceset(grid, "top"), + getfaceset(grid, "bottom"), + getfaceset(grid, "front"), + getfaceset(grid, "back"), + ) + + # Pre-allocation of vectors for the solution and Newton increments + _ndofs = ndofs(dh) + un = zeros(_ndofs) # previous solution vector + u = zeros(_ndofs) + Δu = zeros(_ndofs) + ΔΔu = zeros(_ndofs) + apply!(un, dbcs) + + # Create sparse matrix and residual vector + K = create_sparsity_pattern(dh) + g = zeros(_ndofs) + + # Perform Newton iterations + newton_itr = -1 + NEWTON_TOL = 1e-8 + NEWTON_MAXITER = 30 + prog = ProgressMeter.ProgressThresh(NEWTON_TOL, "Solving:") + + while true; newton_itr += 1 + # Construct the current guess + u .= un .+ Δu + # Compute residual and tangent for current guess + assemble_global!(K, g, dh, cv, fv, mp, u, ΓN, F_storage, F_avg) + # Apply boundary conditions + apply_zero!(K, g, dbcs) + # Compute the residual norm and compare with tolerance + normg = norm(g) + ProgressMeter.update!(prog, normg; showvalues = [(:iter, newton_itr)]) + if normg < NEWTON_TOL + break + elseif newton_itr > NEWTON_MAXITER + error("Reached maximum Newton iterations, aborting") + end + + # Compute increment using conjugate gradients + IterativeSolvers.cg!(ΔΔu, K, g; maxiter=1000) + + apply_zero!(ΔΔu, dbcs) + Δu .-= ΔΔu #Δu = Δu - Δ(Δu) + end + + return u, F_storage, F_avg +end + +u, F_storage, F_avg, LogEH = solve()=# \ 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/TopOpt.jl b/src/TopOpt.jl index 9ce76f22..c44fbe8d 100644 --- a/src/TopOpt.jl +++ b/src/TopOpt.jl @@ -97,6 +97,7 @@ export TopOpt, Compliance, CG, Direct, + Hyperelastic, Assembly, MatrixFree, FEASolver, diff --git a/src/TopOptProblems/TopOptProblems.jl b/src/TopOptProblems/TopOptProblems.jl index a8b1ec8b..bc77f02c 100644 --- a/src/TopOptProblems/TopOptProblems.jl +++ b/src/TopOptProblems/TopOptProblems.jl @@ -39,7 +39,9 @@ export RayProblem, StiffnessTopOptProblem, AbstractTopOptProblem, GlobalFEAInfo, + GlobalFEAInfo_hyperelastic, ElementFEAInfo, + ElementFEAInfo_hyperelastic, YoungsModulus, assemble, assemble_f!, @@ -51,6 +53,7 @@ export RayProblem, RandomMagnitude, MultiLoad, TensionBar, - TensionRoller + 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 619aef72..70194820 100644 --- a/src/TopOptProblems/elementinfo.jl +++ b/src/TopOptProblems/elementinfo.jl @@ -38,6 +38,21 @@ An instance of the `ElementFEAInfo` type stores element information such as: cells::Any end +@params struct ElementFEAInfo_hyperelastic{dim,T} + Kes::AbstractVector{<:AbstractMatrix{T}} + fes::AbstractVector{<:AbstractVector{T}} + ges::AbstractVector{<:AbstractVector{T}} + fixedload::AbstractVector{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, @@ -86,6 +101,39 @@ function ElementFEAInfo( ) end +function ElementFEAInfo_hyperelastic( + mp, sp, u, quad_order=2, ::Type{Val{mat_type}}=Val{:Static} +) where {mat_type} + Kes, weights, dloads, ges, cellvalues, facevalues = make_Kes_and_fes_hyperelastic( + mp, sp, u, quad_order, Val{mat_type} + ) + 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)) + 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!!!!!!!! + 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 @@ -146,6 +194,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 f0020167..088ca508 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..5918254d 100644 --- a/src/TopOptProblems/matrices_and_vectors.jl +++ b/src/TopOptProblems/matrices_and_vectors.jl @@ -48,7 +48,7 @@ function gettypes( return Matrix{T}, Vector{T} 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 +103,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}}) 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 = gettypes(T, Val{mat_type}, Val{Kesize}) + Kes, weights, ges = _make_Kes_and_weights_hyperelastic( + dh, + Tuple{MatrixType,VectorType}, + Val{n_basefuncs}, + Val{dim * n_basefuncs}, + #C, + mp, + u, + ρ, + quadrature_rule, + cellvalues, + cellvaluesV, + ) + dloads, ges2 = _make_dloads_hyperelastic(weights, problem, facevalues, facevaluesV, ges) + ges3 = ges + ges2 + + return Kes, weights, dloads, ges3, 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 +262,173 @@ function _make_Kes_and_weights( return Kes, weights end +function Ψ(C, mp) # JGB: add to .ipynb + μ = mp.μ + λ = mp.λ + Ic = tr(C) + J = sqrt(det(C)) + return μ / 2 * (Ic - 3) - μ * log(J) + λ / 2 * log(J)^2 +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}}, + ::Type{Val{n_basefuncs}}, + ::Type{Val{Kesize}}, + mp, + u, + ρ, + quadrature_rule, + cellvalues, + cellvaluesV, +) where {dim,N,T,MatrixType<:StaticArray,VectorType,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 = MatrixType{T}[] + #Kes = Vector{MatrixType{T}}() + #Kes = Vector{MatrixType}[] + #Kes = Vector{T,MatrixType}() + Kes = Vector{MatrixType}() + sizehint!(Kes, nel) + weights = [zeros(VectorType) for i in 1:nel] + ges = [zeros(VectorType) for i in 1:nel] + Ke_e = zeros(T, dim, dim) + ge = zeros(T, Kesize) + ge2 = zeros(T, Kesize) + fe = zeros(T, Kesize) + fe2 = zeros(T, Kesize) + Ke_0 = Matrix{T}(undef, Kesize, Kesize) + Ke_02 = Matrix{T}(undef, Kesize, Kesize) # JGB: just for checking + 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] + #fe2 = weights[k] + ge = ges[k] + #ge2 = ges[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 + 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 + #=for b in 1:n_basefuncs + ∇ϕb = shape_gradient(cellvalues, q_point, b) # JGB: like ∇δui 3x3 + ϕb = shape_value(cellvalues, q_point, b) # JGB: like δui + for d2 in 1:dim + #ge[b] += ( ∇ϕb ⊡ P - ϕb ⋅ body_force ) * dΩ # Add contribution to the residual from this test function + #∇ϕb∂P∂F = ∇ϕb ⊡ ∂P∂F # Hoisted computation # JGB add + fe2 = @set fe2[(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) + ge2 = @set ge2[(b - 1) * dim + d2] += ( ∇ϕb ⋅ P[d2,:] - ϕb * body_force[d2]) * dΩ + for a in 1:n_basefuncs + ∇ϕa = shape_gradient(cellvalues, q_point, a) # JGB: like ∇δuj + Ke_e .= dotdot(∇ϕa, SymmetricTensor{4,dim}(∂P∂F), ∇ϕb) * dΩ#(∇ϕb∂P∂F ⊡ ∇ϕa) * dΩ -------error: doesn't understand method + for d1 in 1:dim + # C = SymmetricTensor{4,dim}(g) + #Ke_0[a,b] += (∇ϕb∂P∂F ⊡ ∇ϕa) * dΩ + #Ke_02[dim * (a - 1) + d1, dim * (b - 1) + d2] += Ke_e[d1, d2] + Ke_02[dim * (a - 1) + d1, dim * (b - 1) + d2] += Ke_e[d1,d2] + end + end + end + end=# + end + # ge and fe still technically need checked with mean, norm, etc. to verify they are working correctly because they are trivial the first time through + weights[k] = fe + ges[k] = ge + 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 # 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 +482,69 @@ function _make_dloads(fes, problem, facevalues) return dloads end +function _make_dloads_hyperelastic(fes, problem, facevalues, facevaluesV, ges) + 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) + 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 +570,24 @@ function make_cload(problem) end return sparsevec(inds, vals, ndofs(dh)) end + +function make_cload_hyperelastic(problem) + T = floattype(problem) + dim = getdim(problem) + cloads = getcloaddict(problem) + 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 From 02a72fd6737c31fa5658fa14ff2bbb540d959a45 Mon Sep 17 00:00:00 2001 From: jbeckett Date: Mon, 5 Aug 2024 17:40:56 -0400 Subject: [PATCH 12/14] Update .gitignore Only ignore .jl files in top level of directory. Also ignore .png and .zip files. --- .gitignore | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/.gitignore b/.gitignore index 49832ec0..63395e93 100644 --- a/.gitignore +++ b/.gitignore @@ -17,5 +17,7 @@ Manifest.toml examples/benchmark/iterations *.ipynb -*.jl -*.vtu \ No newline at end of file +/*.jl +*.vtu +*.png +*.zip \ No newline at end of file From 70824729abcc03efedc76f1a5c3bd11d84027cc4 Mon Sep 17 00:00:00 2001 From: Beckett Date: Thu, 10 Oct 2024 09:34:50 -0400 Subject: [PATCH 13/14] Add time stepping to hyperelastic solver Additional validation is pending regarding the time stepping code, but confidence is fairly high that it is implemented correctly. Stability seems to be excellent so long as sufficiently small time steps are used to achieve the solution. This also includes some of the foundation for implenting a nearly-incompressible hyperelastic solver that is not implemented yet. --- src/FEA/FEA.jl | 1 + src/FEA/iterative_hyperelastic_solver.jl | 193 +++++++++++++-------- src/Functions/displacement.jl | 38 ++++ src/TopOptProblems/elementinfo.jl | 6 +- src/TopOptProblems/matrices_and_vectors.jl | 25 +-- src/TopOptProblems/problem_types.jl | 42 ++--- 6 files changed, 195 insertions(+), 110 deletions(-) diff --git a/src/FEA/FEA.jl b/src/FEA/FEA.jl index 5996dfef..d13a4176 100644 --- a/src/FEA/FEA.jl +++ b/src/FEA/FEA.jl @@ -9,6 +9,7 @@ using Parameters: @unpack export AbstractFEASolver, AbstractDisplacementSolver, + AbstractHyperelasticSolver, DirectDisplacementSolver, PCGDisplacementSolver, StaticMatrixFreeDisplacementSolver, diff --git a/src/FEA/iterative_hyperelastic_solver.jl b/src/FEA/iterative_hyperelastic_solver.jl index 04ee7576..90fb4962 100644 --- a/src/FEA/iterative_hyperelastic_solver.jl +++ b/src/FEA/iterative_hyperelastic_solver.jl @@ -1,6 +1,6 @@ -abstract type AbstractHyperelasticSolver <: AbstractDisplacementSolver end +abstract type AbstractHyperelasticSolver <: AbstractFEASolver end -mutable struct HyperelasticDisplacementSolver{ +mutable struct HyperelasticCompressibleDisplacementSolver{ T, dim, TP1<:AbstractPenalty{T}, @@ -14,119 +14,158 @@ mutable struct HyperelasticDisplacementSolver{ globalinfo::TG elementinfo::TE u::Tu # JGB: u --> u0 - lhs::Tu - rhs::Tu vars::Tu penalty::TP1 prev_penalty::TP1 xmin::T - qr::Bool + tsteps::Int + ntsteps::Int end -function Base.show(::IO, ::MIME{Symbol("text/plain")}, x::HyperelasticDisplacementSolver) - return println("TopOpt hyperelastic solver") +mutable struct HyperelasticNearlyIncompressibleDisplacementSolver{ + T, + dim, + TP1<:AbstractPenalty{T}, + TP2<:StiffnessTopOptProblem{dim,T}, + TG<:GlobalFEAInfo_hyperelastic{T}, + TE<:ElementFEAInfo_hyperelastic{dim,T}, + Tu<:AbstractVector{T}, +} <: AbstractHyperelasticSolver + mp + problem::TP2 + globalinfo::TG + elementinfo::TE + u::Tu # JGB: u --> u0 + 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, + xmin=T(1)/1000, penalty=PowerPenalty{T}(1), prev_penalty=deepcopy(penalty), quad_order=default_quad_order(sp), - qr=false, + 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}) # JGB: add u + elementinfo = ElementFEAInfo_hyperelastic(mp, sp, u, quad_order, Val{:Static}, nearlyincompressible; ts = ts0) # JGB: add u globalinfo = GlobalFEAInfo_hyperelastic(sp) # JGB: small issue this leads to symmetric K initialization #u = zeros(T, ndofs(sp.ch.dh)) # JGB - lhs = similar(u) - rhs = similar(u) vars = fill(one(T), getncells(sp.ch.dh.grid) - sum(sp.black) - sum(sp.white)) varind = sp.varind - return HyperelasticDisplacementSolver( - mp, sp, globalinfo, elementinfo, u, lhs, rhs, vars, penalty, prev_penalty, xmin, qr # u to u0 - ) + if nearlyincompressible + return HyperelasticNearlyIncompressibleDisplacementSolver(mp, sp, globalinfo, elementinfo, u, vars, penalty, prev_penalty, xmin, tstep, ntsteps) + else + return HyperelasticCompressibleDisplacementSolver(mp, sp, globalinfo, elementinfo, u, vars, penalty, prev_penalty, xmin, tstep, ntsteps) + end end -function (s::HyperelasticDisplacementSolver{T})( +function (s::HyperelasticCompressibleDisplacementSolver{T})( ::Type{Val{safe}}=Val{false}, ::Type{newT}=T; assemble_f=true, - reuse_fact=false, - #rhs=assemble_f ? s.globalinfo.f : s.rhs, - rhs=assemble_f ? s.globalinfo.g : s.rhs, - lhs=assemble_f ? s.u : s.lhs, kwargs..., ) where {T,safe,newT} - #=globalinfo = s.globalinfo - assemble!( - globalinfo, - s.problem, - s.elementinfo, - s.vars, - getpenalty(s), - s.xmin; - assemble_f=assemble_f, - ) - K = globalinfo.K - if safe - m = meandiag(K) - for i in 1:size(K, 1) - if K[i, i] ≈ zero(T) - K[i, i] = m - end - end - end - nans = false - if !reuse_fact - newK = T === newT ? K : newT.(K) - if s.qr - globalinfo.qrK = qr(newK.data) - else - cholK = cholesky(Symmetric(K); check=false) - if issuccess(cholK) - globalinfo.cholK = cholK - else - @warn "The global stiffness matrix is not positive definite. Please check your boundary conditions." - lhs .= T(NaN) - nans = true + globalinfo = s.globalinfo + dh = s.problem.ch.dh + ch = s.problem.ch + + # Pre-allocation of vectors for the solution and Newton increments + _ndofs = ndofs(dh) + #un = s.u + un = zeros(_ndofs) # previous solution vector + #apply!(un, ch) + + # Perform Newton iterations + NEWTON_TOL = 1e-8 + NEWTON_MAXITER =30 # OG is 30 + #prog = ProgressMeter.ProgressThresh(NEWTON_TOL, "Solving:") + + 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 = -1 + while true; newton_itr += 1 + # Construct the current guess + u .= un .+ Δu + #u += Δu + # Compute residual and tangent for current guess + s.elementinfo = ElementFEAInfo_hyperelastic(s.mp, s.problem, u, default_quad_order(s.problem), Val{:Static}; ts=ts) # JGB: add u + #s.globalinfo = GlobalFEAInfo_hyperelastic(s.problem) # JGB: small issue this leads to symmetric K initialization + #assemble_global!(K, g, dh, cv, fv, mp, u, ΓN) + assemble_hyperelastic!(s.globalinfo,s.problem,s.elementinfo,s.vars,getpenalty(s),s.xmin,assemble_f=assemble_f) + #K = s.globalinfo.K + #g = s.globalinfo.g + # Apply boundary conditions + #if newton_itr == 1 + # BSON.bson("C:\\Users\\jbecktt\\.julia\\juliaup\\julia-1.10.5+0.x64.w64.mingw32\\dev\\TopOpt\\test2.bson",K=s.globalinfo.K,g=s.globalinfo.g) + #end + #apply_zero!(s.globalinfo.K, s.globalinfo.g, ch) # why is this not active currently!! + #if newton_itr == 1 + # BSON.bson("C:\\Users\\jbecktt\\.julia\\juliaup\\julia-1.10.5+0.x64.w64.mingw32\\dev\\TopOpt\\test2.bson",K1=s.globalinfo.K,g1=s.globalinfo.g) + #end + # Compute the residual norm and compare with tolerance + normg = norm(s.globalinfo.g) + println("Tstep: $tstep / $ntsteps. Iteration: $newton_itr. normg is equal to " * string(normg)) + if normg < NEWTON_TOL + break + elseif newton_itr > NEWTON_MAXITER + error("Reached maximum Newton iterations, aborting") end + + # Compute increment using conjugate gradients + IterativeSolvers.cg!(ΔΔu, s.globalinfo.K, s.globalinfo.g; maxiter=1000) + + apply_zero!(ΔΔu, ch) + Δu .-= ΔΔu #Δu = Δu - Δ(Δu) + #u = un + Δu + #println("Finished iteration $newton_itr of tstep $tstep of $ntsteps") end + un = u end - nans && return nothing - new_rhs = T === newT ? rhs : newT.(rhs) - fact = s.qr ? globalinfo.qrK : globalinfo.cholK - lhs .= fact \ new_rhs - =# - # NEW STUFF + # worst case scenario I can save F here for now + s.u .= un #, F_storage, F_avg + 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} globalinfo = s.globalinfo - #assemble_hyperelastic!( - # globalinfo, - # s.problem, - # s.elementinfo, - # s.vars, - # getpenalty(s), - # s.xmin; - # assemble_f=assemble_f, - #) - #K = globalinfo.K dh = s.problem.ch.dh ch = s.problem.ch # Pre-allocation of vectors for the solution and Newton increments _ndofs = ndofs(dh) un = s.u - #un = zeros(_ndofs) # previous solution vector u = zeros(_ndofs) - #u = s.u Δu = zeros(_ndofs) ΔΔu = zeros(_ndofs) - #apply!(un, ch) - - # Create sparse matrix and residual vector - #K = create_sparsity_pattern(dh) - #K = globalinfo.K - #g = zeros(_ndofs) - #g = globalinfo.g # Perform Newton iterations newton_itr = -1 diff --git a/src/Functions/displacement.jl b/src/Functions/displacement.jl index 4ecd2d45..3eafc7de 100644 --- a/src/Functions/displacement.jl +++ b/src/Functions/displacement.jl @@ -7,10 +7,24 @@ 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 @@ -37,6 +51,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(3, 3) 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 @@ -54,6 +82,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) # I need to add F support +end + """ rrule for autodiff. diff --git a/src/TopOptProblems/elementinfo.jl b/src/TopOptProblems/elementinfo.jl index 70194820..3e751ac3 100644 --- a/src/TopOptProblems/elementinfo.jl +++ b/src/TopOptProblems/elementinfo.jl @@ -102,10 +102,10 @@ function ElementFEAInfo( end function ElementFEAInfo_hyperelastic( - mp, sp, u, quad_order=2, ::Type{Val{mat_type}}=Val{:Static} + mp, sp, u, quad_order=2, ::Type{Val{mat_type}}=Val{:Static}, nearlyincompressible=false; ts = 1.0, ) where {mat_type} Kes, weights, dloads, ges, cellvalues, facevalues = make_Kes_and_fes_hyperelastic( - mp, sp, u, quad_order, Val{mat_type} + mp, sp, u, quad_order, Val{mat_type}, ts ) element_Kes = convert( # make sure this isn't going to symmetric Vector{<:ElementMatrix}, @@ -113,7 +113,7 @@ function ElementFEAInfo_hyperelastic( bc_dofs=sp.ch.prescribed_dofs, dof_cells=sp.metadata.dof_cells, ) - fixedload = Vector(make_cload_hyperelastic(sp)) + 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) diff --git a/src/TopOptProblems/matrices_and_vectors.jl b/src/TopOptProblems/matrices_and_vectors.jl index 5918254d..0e7f35a2 100644 --- a/src/TopOptProblems/matrices_and_vectors.jl +++ b/src/TopOptProblems/matrices_and_vectors.jl @@ -103,7 +103,7 @@ 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}}) where {mat_type} +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) @@ -150,10 +150,10 @@ function make_Kes_and_fes_hyperelastic(mp, problem, u, quad_order, ::Type{Val{ma cellvalues, cellvaluesV, ) - dloads, ges2 = _make_dloads_hyperelastic(weights, problem, facevalues, facevaluesV, ges) - ges3 = ges + ges2 + dloads, ges2 = _make_dloads_hyperelastic(weights, problem, facevalues, facevaluesV, ges, ts) + ges += ges2 - return Kes, weights, dloads, ges3, cellvalues, facevalues # switched to ges2 to solve error where 2x ges2 with no contribution from dload + return Kes, weights, dloads, ges, 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 @@ -265,9 +265,14 @@ end function Ψ(C, mp) # JGB: add to .ipynb μ = mp.μ λ = mp.λ - Ic = tr(C) + I1 = tr(C) J = sqrt(det(C)) - return μ / 2 * (Ic - 3) - μ * log(J) + λ / 2 * log(J)^2 + 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 @@ -482,7 +487,7 @@ function _make_dloads(fes, problem, facevalues) return dloads end -function _make_dloads_hyperelastic(fes, problem, facevalues, facevaluesV, ges) +function _make_dloads_hyperelastic(fes, problem, facevalues, facevaluesV, ges, ts) dim = getdim(problem) N = nnodespercell(problem) T = floattype(problem) @@ -502,7 +507,7 @@ function _make_dloads_hyperelastic(fes, problem, facevalues, facevaluesV, ges) ges2[i] .= 0 end end - pressuredict = getpressuredict(problem) + pressuredict = getpressuredict(problem; ts=ts) dh = getdh(problem) grid = dh.grid boundary_matrix = grid.boundary_matrix @@ -571,10 +576,10 @@ function make_cload(problem) return sparsevec(inds, vals, ndofs(dh)) end -function make_cload_hyperelastic(problem) +function make_cload_hyperelastic(problem,ts) T = floattype(problem) dim = getdim(problem) - cloads = getcloaddict(problem) + cloads = getcloaddict(problem; ts=ts) dh = getdh(problem) metadata = getmetadata(problem) node_dofs = metadata.node_dofs diff --git a/src/TopOptProblems/problem_types.jl b/src/TopOptProblems/problem_types.jl index ec0cb740..41474cbf 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) """ @@ -162,7 +162,8 @@ function PointLoadCantilever( ) add!(ch, dbc) close!(ch) - t = T(0) + + t = T(1) update!(ch, t) metadata = Metadata(dh) @@ -323,7 +324,7 @@ function HalfMBB( add!(ch, dbc2) close!(ch) - t = T(0) + t = T(1) update!(ch, t) metadata = Metadata(dh) @@ -341,8 +342,8 @@ function HalfMBB( return HalfMBB(rect_grid, E, ν, ch, force, force_dof, black, white, varind, metadata) end -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 @@ -461,7 +462,7 @@ function LBeam( add!(ch, dbc) close!(ch) - t = T(0) + t = T(1) update!(ch, t) metadata = Metadata(dh) @@ -521,8 +522,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 @@ -612,7 +613,7 @@ function TieBeam( add!(ch, dbc2) close!(ch) - t = T(0) + t = T(1) update!(ch, t) metadata = Metadata(dh) @@ -625,8 +626,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 @@ -699,7 +700,7 @@ function RayProblem( add!(ch, dbc) end close!(ch) - t = T(0) + t = T(1) update!(ch, t) metadata = Metadata(dh) @@ -715,8 +716,9 @@ function RayProblem( return RayProblem(rect_grid, 1.0, 0.3, ch, loadsdict, black, white, varind, metadata) end -getcloaddict(p::RayProblem) = p.loads - +function getcloaddict(p::RayProblem; ts=1.0) + return Dict(k => v .* ts for (k, v) in p.loads) +end """ ``` ///**********************************-> @@ -813,10 +815,10 @@ function TensionBar( 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[disp; zeros(dim-1)], collect(1:dim)) + 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(0) + t = T(1) update!(ch, t) metadata = Metadata(dh) @@ -939,7 +941,7 @@ function TensionRoller( 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[disp], [1]) # set u1 to 'disp' for right face + 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) @@ -948,7 +950,7 @@ function TensionRoller( add!(ch, dbc_fixed_midpt2) end close!(ch) - t = T(0) + t = T(1) update!(ch, t) metadata = Metadata(dh) From edef8ae6679260c9d031a7736dd33ccea8d65640 Mon Sep 17 00:00:00 2001 From: Beckett Date: Thu, 10 Oct 2024 21:05:35 -0400 Subject: [PATCH 14/14] Add def grad export and smart time stepping Nearly incompressible HE is still only partially implemented. The smart time stepping is not adjustable yet via kwargs. --- src/FEA/iterative_hyperelastic_solver.jl | 439 ++++++--------------- src/Functions/displacement.jl | 4 +- src/TopOptProblems/elementinfo.jl | 4 +- src/TopOptProblems/matrices_and_vectors.jl | 112 +++--- 4 files changed, 182 insertions(+), 377 deletions(-) diff --git a/src/FEA/iterative_hyperelastic_solver.jl b/src/FEA/iterative_hyperelastic_solver.jl index 90fb4962..9fa43ea8 100644 --- a/src/FEA/iterative_hyperelastic_solver.jl +++ b/src/FEA/iterative_hyperelastic_solver.jl @@ -8,12 +8,14 @@ mutable struct HyperelasticCompressibleDisplacementSolver{ 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 @@ -29,12 +31,14 @@ mutable struct HyperelasticNearlyIncompressibleDisplacementSolver{ 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 @@ -64,15 +68,19 @@ function HyperelasticDisplacementSolver( 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 - if nearlyincompressible - return HyperelasticNearlyIncompressibleDisplacementSolver(mp, sp, globalinfo, elementinfo, u, vars, penalty, prev_penalty, xmin, tstep, ntsteps) - else - return HyperelasticCompressibleDisplacementSolver(mp, sp, globalinfo, elementinfo, u, vars, penalty, prev_penalty, xmin, tstep, ntsteps) - end + 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}, @@ -80,72 +88,99 @@ function (s::HyperelasticCompressibleDisplacementSolver{T})( assemble_f=true, kwargs..., ) where {T,safe,newT} + elementinfo = s.elementinfo globalinfo = s.globalinfo dh = s.problem.ch.dh ch = s.problem.ch - # Pre-allocation of vectors for the solution and Newton increments _ndofs = ndofs(dh) - #un = s.u un = zeros(_ndofs) # previous solution vector - #apply!(un, ch) - # Perform Newton iterations NEWTON_TOL = 1e-8 - NEWTON_MAXITER =30 # OG is 30 - #prog = ProgressMeter.ProgressThresh(NEWTON_TOL, "Solving:") + NEWTON_MAXITER = 30 + CG_MAXITER = 1000 + TS_MAXITER_ABS = 200 + TS_MAXITER_REL = 10 - ntsteps = s.ntsteps - for tstep ∈ 1:ntsteps - ts = tstep/ntsteps + function HyperelasticSolverCore(ts) update!(ch,ts) - apply!(un,ch) - println(maximum(un)) + apply!(un,ch) + #println(maximum(un)) u = zeros(_ndofs) Δu = zeros(_ndofs) ΔΔu = zeros(_ndofs) - newton_itr = -1 + newton_itr = 0 + normg = zeros(NEWTON_MAXITER) while true; newton_itr += 1 - # Construct the current guess - u .= un .+ Δu - #u += Δu - # Compute residual and tangent for current guess - s.elementinfo = ElementFEAInfo_hyperelastic(s.mp, s.problem, u, default_quad_order(s.problem), Val{:Static}; ts=ts) # JGB: add u - #s.globalinfo = GlobalFEAInfo_hyperelastic(s.problem) # JGB: small issue this leads to symmetric K initialization - #assemble_global!(K, g, dh, cv, fv, mp, u, ΓN) - assemble_hyperelastic!(s.globalinfo,s.problem,s.elementinfo,s.vars,getpenalty(s),s.xmin,assemble_f=assemble_f) - #K = s.globalinfo.K - #g = s.globalinfo.g - # Apply boundary conditions - #if newton_itr == 1 - # BSON.bson("C:\\Users\\jbecktt\\.julia\\juliaup\\julia-1.10.5+0.x64.w64.mingw32\\dev\\TopOpt\\test2.bson",K=s.globalinfo.K,g=s.globalinfo.g) - #end - #apply_zero!(s.globalinfo.K, s.globalinfo.g, ch) # why is this not active currently!! - #if newton_itr == 1 - # BSON.bson("C:\\Users\\jbecktt\\.julia\\juliaup\\julia-1.10.5+0.x64.w64.mingw32\\dev\\TopOpt\\test2.bson",K1=s.globalinfo.K,g1=s.globalinfo.g) - #end - # Compute the residual norm and compare with tolerance - normg = norm(s.globalinfo.g) - println("Tstep: $tstep / $ntsteps. Iteration: $newton_itr. normg is equal to " * string(normg)) - if normg < NEWTON_TOL + 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, s.globalinfo.K, s.globalinfo.g; maxiter=1000) - + IterativeSolvers.cg!(ΔΔu, globalinfo.K, globalinfo.g; maxiter=CG_MAXITER) apply_zero!(ΔΔu, ch) - Δu .-= ΔΔu #Δu = Δu - Δ(Δu) - #u = un + Δu - #println("Finished iteration $newton_itr of tstep $tstep of $ntsteps") + Δu .-= ΔΔu end un = u end - # worst case scenario I can save F here for now - s.u .= un #, F_storage, F_avg + + 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 @@ -156,288 +191,52 @@ function (s::HyperelasticNearlyIncompressibleDisplacementSolver{T})( reuse_fact=false, kwargs..., ) where {T,safe,newT} + elementinfo = s.elementinfo globalinfo = s.globalinfo dh = s.problem.ch.dh ch = s.problem.ch - # Pre-allocation of vectors for the solution and Newton increments - _ndofs = ndofs(dh) - un = s.u - u = zeros(_ndofs) - Δu = zeros(_ndofs) - ΔΔu = zeros(_ndofs) - - # Perform Newton iterations - newton_itr = -1 - NEWTON_TOL = 1e-8 - NEWTON_MAXITER = 1000 # OG is 30 - #prog = ProgressMeter.ProgressThresh(NEWTON_TOL, "Solving:") - - while true; newton_itr += 1 - # Construct the current guess - u .= un .+ Δu - # Compute residual and tangent for current guess - s.elementinfo = ElementFEAInfo_hyperelastic(s.mp, s.problem, u, default_quad_order(s.problem), Val{:Static}) # JGB: add u - #s.globalinfo = GlobalFEAInfo_hyperelastic(s.problem) # JGB: small issue this leads to symmetric K initialization - #assemble_global!(K, g, dh, cv, fv, mp, u, ΓN) - assemble_hyperelastic!(s.globalinfo,s.problem,s.elementinfo,s.vars,getpenalty(s),s.xmin,assemble_f=assemble_f) - #K = s.globalinfo.K - #g = s.globalinfo.g - # Apply boundary conditions - #apply_zero!(s.globalinfo.K, s.globalinfo.g, ch) - # Compute the residual norm and compare with tolerance - normg = norm(s.globalinfo.g) - if normg < NEWTON_TOL - break - elseif newton_itr > NEWTON_MAXITER - error("Reached maximum Newton iterations, aborting") - end - - # Compute increment using conjugate gradients - IterativeSolvers.cg!(ΔΔu, s.globalinfo.K, s.globalinfo.g; maxiter=1000) - - apply_zero!(ΔΔu, ch) - Δu .-= ΔΔu #Δu = Δu - Δ(Δu) - end - - s.u .= u #, F_storage, F_avg - return nothing -end - - - - - - - - - - - - - - - - - - -# ALL OF THE NEW STUFF TO BE PUT PLACES -# ATTRIBUTION: the following code was primarily sourced from sample code provided in Ferrite.jl documentation -#= -struct NeoHooke # JGB: add to .ipynb - μ::Float64 - λ::Float64 -end - -function Ψ(C, mp::NeoHooke) # JGB: add to .ipynb - μ = mp.μ - λ = mp.λ - Ic = tr(C) - J = sqrt(det(C)) - return μ / 2 * (Ic - 3) - μ * log(J) + λ / 2 * log(J)^2 -end - -function constitutive_driver(C, mp::NeoHooke) - # 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; - -function assemble_element!(ke, ge, cell, cv, fv, mp, ue, ΓN, F_storage, F_avg) - # Reinitialize cell values, and reset output arrays - reinit!(cv, cell) - fill!(ke, 0.0) - fill!(ge, 0.0) - - b = Vec{3}((0.0, -0.5, 0.0)) # Body force - tn = 0.1 # Traction (to be scaled with surface normal) - ndofs = getnbasefunctions(cv) - cell_id=cell.current_cellid.x - - for qp in 1:getnquadpoints(cv) - dΩ = getdetJdV(cv, qp) - # Compute deformation gradient F and right Cauchy-Green tensor C - ∇u = function_gradient(cv, qp, ue) - F = one(∇u) + ∇u #*note that one() is the identity tensor times whatever you put in - C = tdot(F) # F' ⋅ F - # Compute stress and tangent - S, ∂S∂C = constitutive_driver(C, mp) - P = F ⋅ S - I = one(S) - ∂P∂F = otimesu(I, S) + 2 * otimesu(F, I) ⊡ ∂S∂C ⊡ otimesu(F', I) - - # Store F in our data structure - #F_storage[(cell.current_cellid, qp)] = F - F_storage[cell_id][qp] = F - - F_avg[cell_id] = F_avg[cell_id] + cv.qr_weights[1]*F - - # Loop over test functions - for i in 1:ndofs - # Test function and gradient - δui = shape_value(cv, qp, i) - ∇δui = shape_gradient(cv, qp, i) - # Add contribution to the residual from this test function - ge[i] += ( ∇δui ⊡ P - δui ⋅ b ) * dΩ - - ∇δui∂P∂F = ∇δui ⊡ ∂P∂F # Hoisted computation - for j in 1:ndofs - ∇δuj = shape_gradient(cv, qp, j) - # Add contribution to the tangent - ke[i, j] += ( ∇δui∂P∂F ⊡ ∇δuj ) * dΩ - end - end - end - - # Surface integral for the traction - for face in 1:nfaces(cell) # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! dloads start - if (cellid(cell), face) in ΓN - reinit!(fv, cell, face) - for q_point in 1:getnquadpoints(fv) - t = tn * getnormal(fv, q_point) - dΓ = getdetJdV(fv, q_point) - for i in 1:ndofs - δui = shape_value(fv, q_point, i) - ge[i] -= (δui ⋅ t) * dΓ - end - end - end - end -end; - -function assemble_global!(K, g, dh, cv, fv, mp, u, ΓN, F_storage, F_avg) - n = ndofs_per_cell(dh) - ke = zeros(n, n) - ge = zeros(n) - - # start_assemble resets K and g - assembler = start_assemble(K, g) - - # Loop over all cells in the grid - for cell in CellIterator(dh) - global_dofs = celldofs(cell) - ue = u[global_dofs] # element dofs - assemble_element!(ke, ge, cell, cv, fv, mp, ue, ΓN, F_storage, F_avg) - assemble!(assembler, global_dofs, ge, ke) - end -end; - -function solve() - # Generate a grid - N = 10 - L = 1.0 - left = zero(Vec{3}) - right = L * ones(Vec{3}) - grid = generate_grid(Tetrahedron, (N, N, N), left, right) - - # Material parameters - E = 10.0 - ν = 0.3 - μ = E / (2(1 + ν)) - λ = (E * ν) / ((1 + ν) * (1 - 2ν)) - mp = NeoHooke(μ, λ) - - # Finite element base - ip = Lagrange{3, RefTetrahedron, 1}() - qr = QuadratureRule{3, RefTetrahedron}(1) - qr_face = QuadratureRule{2, RefTetrahedron}(1) - cv = CellVectorValues(qr, ip) - fv = FaceVectorValues(qr_face, ip) - - # DofHandler - dh = DofHandler(grid) - push!(dh, :u, 3) # Add a displacement field - close!(dh) - #F_storage = Dict{Tuple{Int, Int}, Tensor{2,3,Float64}}() - ncells = Ferrite.getncells(dh.grid) - nqpoints_per_cell = getnquadpoints(cv) - #F_storage = [Vector{Matrix{Float64}(undef,3,3)}(undef, nqpoints_per_cell) for _ in 1:ncells] - #F_storage = [Vector{Tensor{2,3,Float64}}(undef, nqpoints_per_cell) for _ in 1:ncells] - - # Creating a nested array where each cell will hold an array of Tensors - F_storage = Vector{Vector{Tensor{2, 3, Float64}}}(undef, ncells) - # Initialize the inner arrays in the nested array for each cell - for i in 1:ncells - F_storage[i] = Vector{Tensor{2, 3, Float64}}(undef, nqpoints_per_cell) - end - - F_avg = Vector{Tensor{2, 3, Float64}}(undef, ncells) # Array for averaged deformation gradients. - for i in 1:ncells - F_avg[i] = zero(Tensor{2, 3}) - end - - function rotation(X, t) - θ = pi / 3 # 60° - x, y, z = X - return t * Vec{3}(( - 0.0, - L/2 - y + (y-L/2)*cos(θ) - (z-L/2)*sin(θ), - L/2 - z + (y-L/2)*sin(θ) + (z-L/2)*cos(θ) - )) - end - - dbcs = ConstraintHandler(dh) - # Add a homogeneous boundary condition on the "clamped" edge - dbc = Dirichlet(:u, getfaceset(grid, "right"), (x,t) -> [0.0, 0.0, 0.0], [1, 2, 3]) - add!(dbcs, dbc) - dbc = Dirichlet(:u, getfaceset(grid, "left"), (x,t) -> rotation(x, t), [1, 2, 3]) - add!(dbcs, dbc) - close!(dbcs) - t = 0.5 - Ferrite.update!(dbcs, t) - - # Neumann part of the boundary - ΓN = union( - getfaceset(grid, "top"), - getfaceset(grid, "bottom"), - getfaceset(grid, "front"), - getfaceset(grid, "back"), - ) - - # Pre-allocation of vectors for the solution and Newton increments _ndofs = ndofs(dh) un = zeros(_ndofs) # previous solution vector - u = zeros(_ndofs) - Δu = zeros(_ndofs) - ΔΔu = zeros(_ndofs) - apply!(un, dbcs) - - # Create sparse matrix and residual vector - K = create_sparsity_pattern(dh) - g = zeros(_ndofs) - # Perform Newton iterations - newton_itr = -1 NEWTON_TOL = 1e-8 NEWTON_MAXITER = 30 - prog = ProgressMeter.ProgressThresh(NEWTON_TOL, "Solving:") + CG_MAXITER = 1000 - while true; newton_itr += 1 - # Construct the current guess - u .= un .+ Δu - # Compute residual and tangent for current guess - assemble_global!(K, g, dh, cv, fv, mp, u, ΓN, F_storage, F_avg) - # Apply boundary conditions - apply_zero!(K, g, dbcs) - # Compute the residual norm and compare with tolerance - normg = norm(g) - ProgressMeter.update!(prog, normg; showvalues = [(:iter, newton_itr)]) - if normg < NEWTON_TOL - break - elseif newton_itr > NEWTON_MAXITER - error("Reached maximum Newton iterations, aborting") + 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 - - # Compute increment using conjugate gradients - IterativeSolvers.cg!(ΔΔu, K, g; maxiter=1000) - - apply_zero!(ΔΔu, dbcs) - Δu .-= ΔΔu #Δu = Δu - Δ(Δu) + un = u end - - return u, F_storage, F_avg -end - -u, F_storage, F_avg, LogEH = solve()=# \ No newline at end of file + s.u .= un + s.F .= elementinfo.Fes + return nothing +end \ No newline at end of file diff --git a/src/Functions/displacement.jl b/src/Functions/displacement.jl index 3eafc7de..fc5ce9e2 100644 --- a/src/Functions/displacement.jl +++ b/src/Functions/displacement.jl @@ -60,7 +60,7 @@ function Displacement(solver::AbstractHyperelasticSolver; maxfevals=10^8) global_dofs = zeros(Int, k) total_ndof = ndofs(dh) u = zeros(T, total_ndof) - F = [zeros(3, 3) for _ in 1:total_ndof/dim] + 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 @@ -89,7 +89,7 @@ function (dp::HyperelasticDisplacement{T})(x::PseudoDensities) where {T} @assert length(global_dofs) == ndofs_per_cell(solver.problem.ch.dh) solver.vars .= x.x solver() - return DisplacementResult(copy(solver.u)) #, copy(solver.F) # I need to add F support + return DisplacementResult(copy(solver.u)), copy(solver.F) #, copy(solver.F) # I need to add F support end """ diff --git a/src/TopOptProblems/elementinfo.jl b/src/TopOptProblems/elementinfo.jl index 3e751ac3..d6ced066 100644 --- a/src/TopOptProblems/elementinfo.jl +++ b/src/TopOptProblems/elementinfo.jl @@ -43,6 +43,7 @@ end 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} @@ -104,7 +105,7 @@ 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, cellvalues, facevalues = make_Kes_and_fes_hyperelastic( + 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 @@ -123,6 +124,7 @@ function ElementFEAInfo_hyperelastic( weights, ges, fixedload, # this is g version now!!!!!!!! + Fes, cellvolumes, cellvalues, facevalues, diff --git a/src/TopOptProblems/matrices_and_vectors.jl b/src/TopOptProblems/matrices_and_vectors.jl index 0e7f35a2..d3eea8da 100644 --- a/src/TopOptProblems/matrices_and_vectors.jl +++ b/src/TopOptProblems/matrices_and_vectors.jl @@ -1,51 +1,84 @@ 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::Bool=true) = symmetric ? Symmetric(create_sparsity_pattern(sp.ch.dh)) : create_sparsity_pattern(sp.ch.dh) @@ -136,10 +169,10 @@ function make_Kes_and_fes_hyperelastic(mp, problem, u, quad_order, ::Type{Val{ma n_basefuncs = getnbasefunctions(cellvalues) Kesize = dim * n_basefuncs - MatrixType, VectorType = gettypes(T, Val{mat_type}, Val{Kesize}) - Kes, weights, ges = _make_Kes_and_weights_hyperelastic( + 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}, + Tuple{MatrixType,VectorType,MatrixTypeF}, Val{n_basefuncs}, Val{dim * n_basefuncs}, #C, @@ -153,7 +186,7 @@ function make_Kes_and_fes_hyperelastic(mp, problem, u, quad_order, ::Type{Val{ma dloads, ges2 = _make_dloads_hyperelastic(weights, problem, facevalues, facevaluesV, ges, ts) ges += ges2 - return Kes, weights, dloads, ges, cellvalues, facevalues # switched to ges2 to solve error where 2x ges2 with no contribution from dload + 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 @@ -287,7 +320,7 @@ end # `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}}, + ::Type{Tuple{MatrixType,VectorType,MatrixTypeF}}, ::Type{Val{n_basefuncs}}, ::Type{Val{Kesize}}, mp, @@ -296,26 +329,17 @@ function _make_Kes_and_weights_hyperelastic( quadrature_rule, cellvalues, cellvaluesV, -) where {dim,N,T,MatrixType<:StaticArray,VectorType,n_basefuncs,Kesize} +) 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 = MatrixType{T}[] - #Kes = Vector{MatrixType{T}}() - #Kes = Vector{MatrixType}[] - #Kes = Vector{T,MatrixType}() Kes = Vector{MatrixType}() sizehint!(Kes, nel) weights = [zeros(VectorType) for i in 1:nel] ges = [zeros(VectorType) for i in 1:nel] - Ke_e = zeros(T, dim, dim) - ge = zeros(T, Kesize) - ge2 = zeros(T, Kesize) - fe = zeros(T, Kesize) - fe2 = zeros(T, Kesize) + Fes = [zeros(MatrixTypeF) for _ in 1:nel] Ke_0 = Matrix{T}(undef, Kesize, Kesize) - Ke_02 = Matrix{T}(undef, Kesize, Kesize) # JGB: just for checking celliterator = CellIterator(dh) for (k, cell) in enumerate(celliterator) Ke_0 .= 0 @@ -325,13 +349,13 @@ function _make_Kes_and_weights_hyperelastic( reinit!(cellvalues, cell) reinit!(cellvaluesV, cell) fe = weights[k] - #fe2 = weights[k] ge = ges[k] - #ge2 = 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 @@ -349,30 +373,10 @@ function _make_Kes_and_weights_hyperelastic( Ke_0[a,b] += (∇ϕb∂P∂F ⊡ ∇ϕa) * dΩ end end - #=for b in 1:n_basefuncs - ∇ϕb = shape_gradient(cellvalues, q_point, b) # JGB: like ∇δui 3x3 - ϕb = shape_value(cellvalues, q_point, b) # JGB: like δui - for d2 in 1:dim - #ge[b] += ( ∇ϕb ⊡ P - ϕb ⋅ body_force ) * dΩ # Add contribution to the residual from this test function - #∇ϕb∂P∂F = ∇ϕb ⊡ ∂P∂F # Hoisted computation # JGB add - fe2 = @set fe2[(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) - ge2 = @set ge2[(b - 1) * dim + d2] += ( ∇ϕb ⋅ P[d2,:] - ϕb * body_force[d2]) * dΩ - for a in 1:n_basefuncs - ∇ϕa = shape_gradient(cellvalues, q_point, a) # JGB: like ∇δuj - Ke_e .= dotdot(∇ϕa, SymmetricTensor{4,dim}(∂P∂F), ∇ϕb) * dΩ#(∇ϕb∂P∂F ⊡ ∇ϕa) * dΩ -------error: doesn't understand method - for d1 in 1:dim - # C = SymmetricTensor{4,dim}(g) - #Ke_0[a,b] += (∇ϕb∂P∂F ⊡ ∇ϕa) * dΩ - #Ke_02[dim * (a - 1) + d1, dim * (b - 1) + d2] += Ke_e[d1, d2] - Ke_02[dim * (a - 1) + d1, dim * (b - 1) + d2] += Ke_e[d1,d2] - end - end - end - end=# end - # ge and fe still technically need checked with mean, norm, etc. to verify they are working correctly because they are trivial the first time through 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)) @@ -381,7 +385,7 @@ function _make_Kes_and_weights_hyperelastic( push!(Kes, MatrixType(Ke_0)) end end - return Kes, weights, ges # weights is fes + return Kes, weights, ges, Fes # weights is fes end # Fallback #= bring this up to speed at a later time to match with