From 357e02280a18f7c63ee14cbba5cf28f5c3c377b6 Mon Sep 17 00:00:00 2001 From: zjwegert <60646897+zjwegert@users.noreply.github.com> Date: Tue, 30 Apr 2024 08:12:09 +1000 Subject: [PATCH] Added some serial tests + bug fix in ChainRules --- scripts/MPI/3d_thermal_compliance_ALM.jl | 2 +- src/ChainRules.jl | 6 +- src/LevelSetTopOpt.jl | 4 +- test/seq/InverseHomogenisationALMTests.jl | 103 ++++++++++++++++ test/seq/InverterHPMTests.jl | 114 ++++++++++++++++++ .../seq/NonlinearThermalComplianceALMTests.jl | 101 ++++++++++++++++ test/seq/ThermalComplianceALMTests.jl | 112 +++++++++++++++-- test/seq/runtests.jl | 3 + 8 files changed, 433 insertions(+), 12 deletions(-) create mode 100644 test/seq/InverseHomogenisationALMTests.jl create mode 100644 test/seq/InverterHPMTests.jl create mode 100644 test/seq/NonlinearThermalComplianceALMTests.jl diff --git a/scripts/MPI/3d_thermal_compliance_ALM.jl b/scripts/MPI/3d_thermal_compliance_ALM.jl index 2c63d3e0..ce4849db 100644 --- a/scripts/MPI/3d_thermal_compliance_ALM.jl +++ b/scripts/MPI/3d_thermal_compliance_ALM.jl @@ -29,7 +29,7 @@ function main(mesh_partition,distribute,el_size) dom = (0,xmax,0,ymax,0,zmax) γ = 0.1 γ_reinit = 0.5 - max_steps = floor(Int,order*minimum(el_size)/step) + max_steps = floor(Int,order*minimum(el_size)/10) tol = 1/(coef*order^2)/minimum(el_size) κ = 1 η_coeff = 2 diff --git a/src/ChainRules.jl b/src/ChainRules.jl index 398f74db..1a2a5f49 100644 --- a/src/ChainRules.jl +++ b/src/ChainRules.jl @@ -65,13 +65,15 @@ Gridap.gradient(F::IntegrandWithMeasure,uh) = Gridap.gradient(F,[uh],1) Given an an `IntegrandWithMeasure` `F` and a vector of `FEFunctions` or `CellField` `uh` (excluding measures) evaluate the Jacobian `F.F` with respect to `uh[K]`. """ -function Gridap.jacobian(F::IntegrandWithMeasure,uh::Vector{<:FEFunction},K::Int) +function Gridap.jacobian(F::IntegrandWithMeasure,uh::Vector{<:Union{FEFunction,CellField}},K::Int) @check 0 < K <= length(uh) _f(uk) = F.F(uh[1:K-1]...,uk,uh[K+1:end]...,F.dΩ...) return Gridap.jacobian(_f,uh[K]) end -function Gridap.jacobian(F::IntegrandWithMeasure,uh::Vector,K::Int) +DistributedFields = Union{DistributedCellField,DistributedMultiFieldCellField} + +function Gridap.jacobian(F::IntegrandWithMeasure,uh::Vector{<:DistributedFields},K::Int) @check 0 < K <= length(uh) local_fields = map(local_views,uh) |> to_parray_of_arrays local_measures = map(local_views,F.dΩ) |> to_parray_of_arrays diff --git a/src/LevelSetTopOpt.jl b/src/LevelSetTopOpt.jl index 52d28514..43f948f7 100644 --- a/src/LevelSetTopOpt.jl +++ b/src/LevelSetTopOpt.jl @@ -21,8 +21,8 @@ using Gridap: writevtk using GridapDistributed using GridapDistributed: DistributedDiscreteModel, DistributedTriangulation, DistributedFESpace, DistributedDomainContribution, to_parray_of_arrays, - allocate_in_domain, DistributedCellField, DistributedMultiFieldFEBasis, - BlockPMatrix, BlockPVector, change_ghost + allocate_in_domain, DistributedCellField, DistributedMultiFieldCellField, + DistributedMultiFieldFEBasis, BlockPMatrix, BlockPVector, change_ghost using PartitionedArrays using PartitionedArrays: getany, tuple_of_arrays, matching_ghost_indices diff --git a/test/seq/InverseHomogenisationALMTests.jl b/test/seq/InverseHomogenisationALMTests.jl new file mode 100644 index 00000000..82929b7d --- /dev/null +++ b/test/seq/InverseHomogenisationALMTests.jl @@ -0,0 +1,103 @@ +module InverseHomogenisationALMTests +using Test + +using Gridap, LevelSetTopOpt + +""" + (Serial) Maximum bulk modulus inverse homogenisation with augmented Lagrangian method in 2D. + + Optimisation problem: + Min J(Ω) = -κ(Ω) + Ω + s.t., Vol(Ω) = vf, + ⎡For unique εᴹᵢ, find uᵢ∈V=H¹ₚₑᵣ(Ω)ᵈ, + ⎣∫ ∑ᵢ C ⊙ ε(uᵢ) ⊙ ε(vᵢ) dΩ = ∫ -∑ᵢ C ⊙ ε⁰ᵢ ⊙ ε(vᵢ) dΩ, ∀v∈V. +""" +function main(;AD) + ## Parameters + order = 1 + xmax,ymax=(1.0,1.0) + dom = (0,xmax,0,ymax) + el_size = (20,20) + γ = 0.05 + γ_reinit = 0.5 + max_steps = floor(Int,order*minimum(el_size)/10) + tol = 1/(5order^2)/minimum(el_size) + C = isotropic_elast_tensor(2,1.,0.3) + η_coeff = 2 + α_coeff = 4max_steps*γ + vf = 0.5 + + ## FE Setup + model = CartesianDiscreteModel(dom,el_size,isperiodic=(true,true)) + el_Δ = get_el_Δ(model) + f_Γ_D(x) = iszero(x) + update_labels!(1,model,f_Γ_D,"origin") + + ## Triangulations and measures + Ω = Triangulation(model) + dΩ = Measure(Ω,2*order) + vol_D = sum(∫(1)dΩ) + + ## Spaces + reffe = ReferenceFE(lagrangian,VectorValue{2,Float64},order) + reffe_scalar = ReferenceFE(lagrangian,Float64,order) + V = TestFESpace(model,reffe;dirichlet_tags=["origin"]) + U = TrialFESpace(V,VectorValue(0.0,0.0)) + V_reg = V_φ = TestFESpace(model,reffe_scalar) + U_reg = TrialFESpace(V_reg) + + ## Create FE functions + lsf_fn = x->max(initial_lsf(2,0.4)(x),initial_lsf(2,0.4;b=VectorValue(0,0.5))(x)) + φh = interpolate(lsf_fn,V_φ) + + ## Interpolation and weak form + interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ)) + I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ + + εᴹ = (TensorValue(1.,0.,0.,0.), # ϵᵢⱼ⁽¹¹⁾≡ϵᵢⱼ⁽¹⁾ + TensorValue(0.,0.,0.,1.), # ϵᵢⱼ⁽²²⁾≡ϵᵢⱼ⁽²⁾ + TensorValue(0.,1/2,1/2,0.)) # ϵᵢⱼ⁽¹²⁾≡ϵᵢⱼ⁽³⁾ + + a(u,v,φ,dΩ) = ∫((I ∘ φ) * C ⊙ ε(u) ⊙ ε(v) )dΩ + l = [(v,φ,dΩ) -> ∫(-(I ∘ φ)* C ⊙ εᴹ[i] ⊙ ε(v))dΩ for i in 1:3] + + ## Optimisation functionals + Cᴴ(r,s,u,φ,dΩ) = ∫((I ∘ φ)*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ εᴹ[s]))dΩ + dCᴴ(r,s,q,u,φ,dΩ) = ∫(-q*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ (ε(u[s])+εᴹ[s]))*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + κ(u,φ,dΩ) = -1/4*(Cᴴ(1,1,u,φ,dΩ)+Cᴴ(2,2,u,φ,dΩ)+2*Cᴴ(1,2,u,φ,dΩ)) + dκ(q,u,φ,dΩ) = -1/4*(dCᴴ(1,1,q,u,φ,dΩ)+dCᴴ(2,2,q,u,φ,dΩ)+2*dCᴴ(1,2,q,u,φ,dΩ)) + Vol(u,φ,dΩ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ + dVol(q,u,φ,dΩ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + + ## Finite difference solver and level set function + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + + ## Setup solver and FE operators + state_map = RepeatingAffineFEStateMap(3,a,l,U,V,V_φ,U_reg,φh,dΩ) + pcfs = if AD + PDEConstrainedFunctionals(κ,[Vol],state_map;analytic_dJ=dκ,analytic_dC=[dVol]) + else + PDEConstrainedFunctionals(κ,[Vol],state_map) + end + + ## Hilbertian extension-regularisation problems + α = α_coeff*maximum(el_Δ) + a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ + vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) + + ## Optimiser + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; + γ,γ_reinit,verbose=true,constraint_names=[:Vol]) + + # Do a few iterations + vars, state = iterate(optimiser) + vars, state = iterate(optimiser,state) + true +end + +# Test that these run successfully +@test main(;AD=true) +@test main(;AD=false) + +end # module \ No newline at end of file diff --git a/test/seq/InverterHPMTests.jl b/test/seq/InverterHPMTests.jl new file mode 100644 index 00000000..ce94a389 --- /dev/null +++ b/test/seq/InverterHPMTests.jl @@ -0,0 +1,114 @@ +module InverterHPMTests +using Test + +using Gridap, LevelSetTopOpt + +""" + (Serial) Inverter mechanism with Hilbertian projection method in 2D. + + Optimisation problem: + Min J(Ω) = ηᵢₙ*∫ u⋅e₁ dΓᵢₙ/Vol(Γᵢₙ) + Ω + s.t., Vol(Ω) = vf, + C(Ω) = 0, + ⎡u∈V=H¹(Ω;u(Γ_D)=0)ᵈ, + ⎣∫ C ⊙ ε(u) ⊙ ε(v) dΩ + ∫ kₛv⋅u dΓₒᵤₜ = ∫ v⋅g dΓᵢₙ , ∀v∈V. + + where C(Ω) = ∫ -u⋅e₁-δₓ dΓₒᵤₜ/Vol(Γₒᵤₜ). We assume symmetry in the problem to aid + convergence. +""" +function main() + ## Parameters + order = 1 + dom = (0,1,0,0.5) + el_size = (20,20) + γ = 0.1 + γ_reinit = 0.5 + max_steps = floor(Int,order*minimum(el_size)/10) + tol = 1/(5order^2)/minimum(el_size) + C = isotropic_elast_tensor(2,1.0,0.3) + η_coeff = 2 + α_coeff = 4max_steps*γ + vf = 0.4 + δₓ = 0.2 + ks = 0.1 + g = VectorValue(0.5,0) + + ## FE Setup + model = CartesianDiscreteModel(dom,el_size) + el_Δ = get_el_Δ(model) + f_Γ_in(x) = (x[1] ≈ 0.0) && (x[2] <= 0.03 + eps()) + f_Γ_out(x) = (x[1] ≈ 1.0) && (x[2] <= 0.07 + eps()) + f_Γ_D(x) = (x[1] ≈ 0.0) && (x[2] >= 0.4) + f_Γ_D2(x) = (x[2] ≈ 0.0) + update_labels!(1,model,f_Γ_in,"Gamma_in") + update_labels!(2,model,f_Γ_out,"Gamma_out") + update_labels!(3,model,f_Γ_D,"Gamma_D") + update_labels!(4,model,f_Γ_D2,"SymLine") + + ## Triangulations and measures + Ω = Triangulation(model) + Γ_in = BoundaryTriangulation(model,tags="Gamma_in") + Γ_out = BoundaryTriangulation(model,tags="Gamma_out") + dΩ = Measure(Ω,2order) + dΓ_in = Measure(Γ_in,2order) + dΓ_out = Measure(Γ_out,2order) + vol_D = sum(∫(1)dΩ) + vol_Γ_in = sum(∫(1)dΓ_in) + vol_Γ_out = sum(∫(1)dΓ_out) + + ## Spaces + reffe = ReferenceFE(lagrangian,VectorValue{2,Float64},order) + reffe_scalar = ReferenceFE(lagrangian,Float64,order) + V = TestFESpace(model,reffe;dirichlet_tags=["Gamma_D","SymLine"], + dirichlet_masks=[(true,true),(false,true)]) + U = TrialFESpace(V,[VectorValue(0.0,0.0),VectorValue(0.0,0.0)]) + V_φ = TestFESpace(model,reffe_scalar) + V_reg = TestFESpace(model,reffe_scalar;dirichlet_tags=["Gamma_in","Gamma_out"]) + U_reg = TrialFESpace(V_reg,[0,0]) + + ## Create FE functions + lsf_fn(x) = min(max(initial_lsf(6,0.2)(x),-sqrt((x[1]-1)^2+(x[2]-0.5)^2)+0.2), + sqrt((x[1])^2+(x[2]-0.5)^2)-0.1) + φh = interpolate(lsf_fn,V_φ) + + ## Interpolation and weak form + interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ)) + I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ + + a(u,v,φ,dΩ,dΓ_in,dΓ_out) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + ∫(ks*(u⋅v))dΓ_out + l(v,φ,dΩ,dΓ_in,dΓ_out) = ∫(v⋅g)dΓ_in + + ## Optimisation functionals + e₁ = VectorValue(1,0) + J(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅e₁)/vol_Γ_in)dΓ_in + Vol(u,φ,dΩ,dΓ_in,dΓ_out) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ + dVol(q,u,φ,dΩ,dΓ_in,dΓ_out) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + UΓ_out(u,φ,dΩ,dΓ_in,dΓ_out) = ∫((u⋅-e₁-δₓ)/vol_Γ_out)dΓ_out + + ## Finite difference solver and level set function + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + + ## Setup solver and FE operators + state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_in,dΓ_out) + pcfs = PDEConstrainedFunctionals(J,[Vol,UΓ_out],state_map,analytic_dC=[dVol,nothing]) + + ## Hilbertian extension-regularisation problems + α = α_coeff*maximum(el_Δ) + a_hilb(p,q) = ∫(α^2*∇(p)⋅∇(q) + p*q)dΩ; + vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) + + ## Optimiser + optimiser = HilbertianProjection(pcfs,ls_evo,vel_ext,φh; + γ,γ_reinit,verbose=true,debug=true,constraint_names=[:Vol,:UΓ_out]) + + # Do a few iterations + vars, state = iterate(optimiser) + vars, state = iterate(optimiser,state) + true +end + +# Test that these run successfully +@test main() + +end # module \ No newline at end of file diff --git a/test/seq/NonlinearThermalComplianceALMTests.jl b/test/seq/NonlinearThermalComplianceALMTests.jl new file mode 100644 index 00000000..c624d997 --- /dev/null +++ b/test/seq/NonlinearThermalComplianceALMTests.jl @@ -0,0 +1,101 @@ +module NonlinearThermalComplianceALMTests +using Test + +using Gridap, LevelSetTopOpt + +""" + (Serial) Minimum thermal compliance with augmented Lagrangian method in 2D with nonlinear diffusivity. + + Optimisation problem: + Min J(Ω) = ∫ κ(u)*∇(u)⋅∇(u) dΩ + Ω + s.t., Vol(Ω) = vf, + ⎡u∈V=H¹(Ω;u(Γ_D)=0), + ⎣∫ κ(u)*∇(u)⋅∇(v) dΩ = ∫ v dΓ_N, ∀v∈V. + + In this example κ(u) = κ0*(exp(ξ*u)) +""" +function main() + ## Parameters + order = 1 + xmax=ymax=1.0 + prop_Γ_N = 0.2 + prop_Γ_D = 0.2 + dom = (0,xmax,0,ymax) + el_size = (20,20) + γ = 0.1 + γ_reinit = 0.5 + max_steps = floor(Int,order*minimum(el_size)/10) + tol = 1/(5order^2)/minimum(el_size) + η_coeff = 2 + α_coeff = 4max_steps*γ + vf = 0.4 + + ## FE Setup + model = CartesianDiscreteModel(dom,el_size); + el_Δ = get_el_Δ(model) + f_Γ_D(x) = (x[1] ≈ 0.0 && (x[2] <= ymax*prop_Γ_D + eps() || + x[2] >= ymax-ymax*prop_Γ_D - eps())); + f_Γ_N(x) = (x[1] ≈ xmax && ymax/2-ymax*prop_Γ_N/2 - eps() <= x[2] <= + ymax/2+ymax*prop_Γ_N/2 + eps()); + update_labels!(1,model,f_Γ_D,"Gamma_D") + update_labels!(2,model,f_Γ_N,"Gamma_N") + + ## Triangulations and measures + Ω = Triangulation(model) + Γ_N = BoundaryTriangulation(model,tags="Gamma_N") + dΩ = Measure(Ω,2*order) + dΓ_N = Measure(Γ_N,2*order) + vol_D = sum(∫(1)dΩ) + + ## Spaces + reffe_scalar = ReferenceFE(lagrangian,Float64,order) + V = TestFESpace(model,reffe_scalar;dirichlet_tags=["Gamma_D"]) + U = TrialFESpace(V,0.0) + V_φ = TestFESpace(model,reffe_scalar) + V_reg = TestFESpace(model,reffe_scalar;dirichlet_tags=["Gamma_N"]) + U_reg = TrialFESpace(V_reg,0) + + ## Create FE functions + φh = interpolate(initial_lsf(4,0.2),V_φ) + + ## Interpolation and weak form + interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ)) + I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ + + κ0 = 1 + ξ = -1 + κ(u) = κ0*(exp(ξ*u)) + res(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(v))dΩ - ∫(v)dΓ_N + + ## Optimisation functionals + J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(u))dΩ + Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + + ## Finite difference solver and level set function + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + + ## Setup solver and FE operators + state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) + pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dC=[dVol]) + + ## Hilbertian extension-regularisation problems + α = α_coeff*maximum(el_Δ) + a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ + vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) + + ## Optimiser + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; + γ,γ_reinit,verbose=true,constraint_names=[:Vol]) + + # Do a few iterations + vars, state = iterate(optimiser) + vars, state = iterate(optimiser,state) + true +end + +# Test that these run successfully +@test main() + +end # module \ No newline at end of file diff --git a/test/seq/ThermalComplianceALMTests.jl b/test/seq/ThermalComplianceALMTests.jl index a5891838..c96225f8 100644 --- a/test/seq/ThermalComplianceALMTests.jl +++ b/test/seq/ThermalComplianceALMTests.jl @@ -1,4 +1,5 @@ module ThermalComplianceALMTests +using Test using Gridap, LevelSetTopOpt @@ -12,22 +13,21 @@ using Gridap, LevelSetTopOpt ⎡u∈V=H¹(Ω;u(Γ_D)=0), ⎣∫ κ*∇(u)⋅∇(v) dΩ = ∫ v dΓ_N, ∀v∈V. """ -function main() - ## Parameters| - order = 1 - xmax=ymax=1.0 +function main(;order,AD) + ## Parameters + xmax = ymax = 1.0 prop_Γ_N = 0.2 prop_Γ_D = 0.2 dom = (0,xmax,0,ymax) el_size = (20,20) γ = 0.1 γ_reinit = 0.5 - max_steps = floor(Int,minimum(el_size)/10) + max_steps = floor(Int,order*minimum(el_size)/10) tol = 1/(5*order^2)/minimum(el_size) κ = 1 vf = 0.4 η_coeff = 2 - α_coeff = 4 + α_coeff = 4max_steps*γ ## FE Setup model = CartesianDiscreteModel(dom,el_size); @@ -73,6 +73,99 @@ function main() ## Finite difference solver and level set function ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) + ## Setup solver and FE operators + state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) + pcfs = if AD + PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol]) + else + PDEConstrainedFunctionals(J,[Vol],state_map) + end + + ## Hilbertian extension-regularisation problems + α = α_coeff*maximum(el_Δ) + a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ; + vel_ext = VelocityExtension(a_hilb,U_reg,V_reg) + + ## Optimiser + optimiser = AugmentedLagrangian(pcfs,ls_evo,vel_ext,φh; + γ,γ_reinit,verbose=true,constraint_names=[:Vol]) + + # Do a few iterations + vars, state = iterate(optimiser) + vars, state = iterate(optimiser,state) + true +end + +""" + (Serial) Minimum thermal compliance with augmented Lagrangian method in 3D. + + Optimisation problem: + Min J(Ω) = ∫ κ*∇(u)⋅∇(u) dΩ + Ω + s.t., Vol(Ω) = vf, + ⎡u∈V=H¹(Ω;u(Γ_D)=0), + ⎣∫ κ*∇(u)⋅∇(v) dΩ = ∫ v dΓ_N, ∀v∈V. +""" +function main_3d(;order) + ## Parameters + xmax = ymax = zmax = 1.0 + prop_Γ_N = 0.2 + prop_Γ_D = 0.2 + dom = (0,xmax,0,ymax,0,zmax) + el_size = (20,20,20) + γ = 0.1 + γ_reinit = 0.5 + max_steps = floor(Int,order*minimum(el_size)/10) + tol = 1/(5*order^2)/minimum(el_size) + κ = 1 + vf = 0.4 + η_coeff = 2 + α_coeff = 4max_steps*γ + + ## FE Setup + model = CartesianDiscreteModel(dom,el_size); + el_Δ = get_el_Δ(model) + f_Γ_D(x) = (x[1] ≈ 0.0) && (x[2] <= ymax*prop_Γ_D + eps() || x[2] >= ymax-ymax*prop_Γ_D - eps()) && + (x[3] <= zmax*prop_Γ_D + eps() || x[3] >= zmax-zmax*prop_Γ_D - eps()) + f_Γ_N(x) = (x[1] ≈ xmax) && (ymax/2-ymax*prop_Γ_N/2 - eps() <= x[2] <= ymax/2+ymax*prop_Γ_N/2 + eps()) && + (zmax/2-zmax*prop_Γ_N/2 - eps() <= x[3] <= zmax/2+zmax*prop_Γ_N/2 + eps()) + update_labels!(1,model,f_Γ_D,"Gamma_D") + update_labels!(2,model,f_Γ_N,"Gamma_N") + + ## Triangulations and measures + Ω = Triangulation(model) + Γ_N = BoundaryTriangulation(model,tags="Gamma_N") + dΩ = Measure(Ω,2*order) + dΓ_N = Measure(Γ_N,2*order) + vol_D = sum(∫(1)dΩ) + + ## Spaces + reffe_scalar = ReferenceFE(lagrangian,Float64,order) + V = TestFESpace(model,reffe_scalar;dirichlet_tags=["Gamma_D"]) + U = TrialFESpace(V,0.0) + V_φ = TestFESpace(model,reffe_scalar) + V_reg = TestFESpace(model,reffe_scalar;dirichlet_tags=["Gamma_N"]) + U_reg = TrialFESpace(V_reg,0) + + ## Create FE functions + φh = interpolate(initial_lsf(4,0.2),V_φ) + + ## Interpolation and weak form + interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ)) + I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ + + a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ + l(v,φ,dΩ,dΓ_N) = ∫(v)dΓ_N + + ## Optimisation functionals + J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ + dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ; + Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + + ## Finite difference solver and level set function + ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) + ## Setup solver and FE operators state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol]) @@ -89,8 +182,13 @@ function main() # Do a few iterations vars, state = iterate(optimiser) vars, state = iterate(optimiser,state) + true end -main() +# Test that these run successfully +@test main(;order=1,AD=true) +@test main(;order=2,AD=true) +@test main(;order=1,AD=false) +@test main_3d(;order=1) end # module \ No newline at end of file diff --git a/test/seq/runtests.jl b/test/seq/runtests.jl index 92ccacfa..20fd54e7 100644 --- a/test/seq/runtests.jl +++ b/test/seq/runtests.jl @@ -3,5 +3,8 @@ module LSTOSequentialTests using Test @time @testset "Thermal Compliance - ALM" begin include("ThermalComplianceALMTests.jl") end +@time @testset "Nonlinear Thermal Compliance - ALM" begin include("NonLinearThermalComplianceALMTests.jl") end +@time @testset "Inverse Homogenisation - ALM" begin include("InverseHomogenisationALMTests.jl") end +@time @testset "Inverter - HPM" begin include("InverterHPMTests.jl") end end # module \ No newline at end of file