From 496671e0d57b925fe451926bc3754605260e3e61 Mon Sep 17 00:00:00 2001 From: Z J Wegert <60646897+zjwegert@users.noreply.github.com> Date: Wed, 13 Nov 2024 21:48:39 +1100 Subject: [PATCH] BREAKING: Removed IntegrandWithMeasure type --- scripts/Benchmarks/benchmark.jl | 52 ++--- scripts/MPI/elastic_compliance_ALM.jl | 14 +- scripts/MPI/inverse_homenisation_ALM.jl | 18 +- .../MPI/nonlinear_thermal_compliance_ALM.jl | 10 +- scripts/MPI/thermal_compliance_ALM.jl | 14 +- scripts/MPI/thermal_compliance_ALM_PETSc.jl | 14 +- scripts/Serial/elastic_compliance_ALM.jl | 14 +- scripts/Serial/elastic_compliance_HPM.jl | 14 +- scripts/Serial/elastic_compliance_LM.jl | 10 +- scripts/Serial/hyperelastic_compliance_ALM.jl | 10 +- .../hyperelastic_compliance_neohook_ALM.jl | 10 +- scripts/Serial/inverse_homenisation_ALM.jl | 18 +- scripts/Serial/inverter_ALM.jl | 14 +- scripts/Serial/inverter_HPM.jl | 14 +- .../nonlinear_thermal_compliance_ALM.jl | 10 +- scripts/Serial/thermal_compliance_ALM.jl | 14 +- .../thermal_compliance_ALM_higherorderlsf.jl | 14 +- src/ChainRules.jl | 187 ++++-------------- src/GridapTopOpt.jl | 3 +- test/mpi/InverseHomogenisationALMTests.jl | 18 +- test/mpi/InverterHPMTests.jl | 14 +- test/mpi/NeohookAnalyticJacALMTests.jl | 16 +- .../mpi/NonlinearThermalComplianceALMTests.jl | 10 +- test/mpi/PZMultiFieldRepeatingStateTests.jl | 24 +-- test/mpi/ThermalComplianceALMTests.jl | 28 +-- test/mpi/ThermalComplianceALMTestsPETSc.jl | 28 +-- test/seq/InverseHomogenisationALMTests.jl | 18 +- test/seq/InverterHPMTests.jl | 14 +- test/seq/NeohookAnalyticJacALMTests.jl | 16 +- .../seq/NonlinearThermalComplianceALMTests.jl | 10 +- test/seq/PZMultiFieldRepeatingStateTests.jl | 24 +-- test/seq/ThermalComplianceALMTests.jl | 28 +-- test/seq/ThermalComplianceHPMTests.jl | 14 +- 33 files changed, 302 insertions(+), 414 deletions(-) diff --git a/scripts/Benchmarks/benchmark.jl b/scripts/Benchmarks/benchmark.jl index c31d395..13900d6 100644 --- a/scripts/Benchmarks/benchmark.jl +++ b/scripts/Benchmarks/benchmark.jl @@ -76,23 +76,23 @@ function nl_elast(mesh_partition,ranks,el_size,order,verbose) ## Volume change J(F) = sqrt(det(C(F))) ## Residual - res(u,v,φ,dΩ,dΓ_N) = ∫( (I ∘ φ)*((dE ∘ (∇(v),∇(u))) ⊙ (S ∘ ∇(u))) )*dΩ - ∫(g⋅v)dΓ_N + res(u,v,φ) = ∫( (I ∘ φ)*((dE ∘ (∇(v),∇(u))) ⊙ (S ∘ ∇(u))) )*dΩ - ∫(g⋅v)dΓ_N Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} Tv = Vector{PetscScalar} lin_solver = ElasticitySolver(V) nl_solver = NewtonSolver(lin_solver;maxiter=50,rtol=10^-8,verbose) state_map = NonlinearFEStateMap( - res,U,V,V_φ,U_reg,φh,dΩ,dΓ_N; + res,U,V,V_φ,U_reg,φh; assem_U = SparseMatrixAssembler(Tm,Tv,U,V), assem_adjoint = SparseMatrixAssembler(Tm,Tv,V,U), assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg), nls = nl_solver, adjoint_ls = lin_solver ) # Objective and constraints - J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*((dE ∘ (∇(u),∇(u))) ⊙ (S ∘ ∇(u))))dΩ + J(u,φ) = ∫((I ∘ φ)*((dE ∘ (∇(u),∇(u))) ⊙ (S ∘ ∇(u))))dΩ vol_D = sum(∫(1)dΩ) - C1(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ - dC1(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + C1(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ + dC1(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ pcfs = PDEConstrainedFunctionals(J,[C1],state_map,analytic_dC=[dC1]) # Velocity extension α = 4max_steps*γ*maximum(get_el_Δ(model)) @@ -152,24 +152,24 @@ function therm(mesh_partition,ranks,el_size,order,verbose) interp = SmoothErsatzMaterialInterpolation(η = 2*maximum(get_el_Δ(model))) I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ # Weak formulation - a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ - l(v,φ,dΩ,dΓ_N) = ∫(g*v)dΓ_N + a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ + l(v,φ) = ∫(g*v)dΓ_N Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} Tv = Vector{PetscScalar} solver = PETScLinearSolver() state_map = AffineFEStateMap( - a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N; + a,l,U,V,V_φ,U_reg,φh; assem_U = SparseMatrixAssembler(Tm,Tv,U,V), assem_adjoint = SparseMatrixAssembler(Tm,Tv,V,U), assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg), ls = solver,adjoint_ls = solver ) # Objective and constraints - J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ - dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ + dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ vol_D = sum(∫(1)dΩ) - C1(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ - dC1(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + C1(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ + dC1(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ pcfs = PDEConstrainedFunctionals(J,[C1],state_map, analytic_dJ=dJ,analytic_dC=[dC1]) # Velocity extension @@ -228,24 +228,24 @@ function elast(mesh_partition,ranks,el_size,order,verbose) interp = SmoothErsatzMaterialInterpolation(η = 2*maximum(get_el_Δ(model))) I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ # Weak formulation - a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ - l(v,φ,dΩ,dΓ_N) = ∫(v⋅g)dΓ_N + a(u,v,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + l(v,φ) = ∫(v⋅g)dΓ_N Tm = SparseMatrixCSR{0,PetscScalar,PetscInt} Tv = Vector{PetscScalar} solver = ElasticitySolver(V) state_map = AffineFEStateMap( - a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N; + a,l,U,V,V_φ,U_reg,φh; assem_U = SparseMatrixAssembler(Tm,Tv,U,V), assem_adjoint = SparseMatrixAssembler(Tm,Tv,V,U), assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg), ls = solver,adjoint_ls = solver ) # Objective and constraints - J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(u)))dΩ - dJ(q,u,φ,dΩ,dΓ_N) = ∫((C ⊙ ε(u) ⊙ ε(u))*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + J(u,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(u)))dΩ + dJ(q,u,φ) = ∫((C ⊙ ε(u) ⊙ ε(u))*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ vol_D = sum(∫(1)dΩ) - C1(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ - dC1(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + C1(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ + dC1(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ pcfs = PDEConstrainedFunctionals(J,[C1],state_map, analytic_dJ=dJ,analytic_dC=[dC1]) # Velocity extension @@ -320,15 +320,15 @@ function inverter_HPM(mesh_partition,ranks,el_size,order,verbose) 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 + a(u,v,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + ∫(ks*(u⋅v))dΓ_out + l(v,φ) = ∫(v⋅g)dΓ_in ## Optimisation functionals e₁ = VectorValue(1,0,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 + J(u,φ) = ∫((u⋅e₁)/vol_Γ_in)dΓ_in + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + UΓ_out(u,φ) = ∫((u⋅-e₁-δₓ)/vol_Γ_out)dΓ_out ## Finite difference solver ls_evo = HamiltonJacobiEvolution(FirstOrderStencil(3,Float64),model,V_φ,tol,max_steps) @@ -339,7 +339,7 @@ function inverter_HPM(mesh_partition,ranks,el_size,order,verbose) solver = ElasticitySolver(V) state_map = AffineFEStateMap( - a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_in,dΓ_out; + a,l,U,V,V_φ,U_reg,φh; assem_U = SparseMatrixAssembler(Tm,Tv,U,V), assem_adjoint = SparseMatrixAssembler(Tm,Tv,V,U), assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg), diff --git a/scripts/MPI/elastic_compliance_ALM.jl b/scripts/MPI/elastic_compliance_ALM.jl index 2cdafcc..ee18819 100644 --- a/scripts/MPI/elastic_compliance_ALM.jl +++ b/scripts/MPI/elastic_compliance_ALM.jl @@ -65,14 +65,14 @@ function main(mesh_partition,distribute,el_size,path) interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ)) I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ - a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ - l(v,φ,dΩ,dΓ_N) = ∫(v⋅g)dΓ_N + a(u,v,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + l(v,φ) = ∫(v⋅g)dΓ_N ## Optimisation functionals - J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(u)))dΩ - dJ(q,u,φ,dΩ,dΓ_N) = ∫((C ⊙ ε(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Ω + J(u,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(u)))dΩ + dJ(q,u,φ) = ∫((C ⊙ ε(u) ⊙ ε(u))*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ + dVol(q,u,φ) = ∫(-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) @@ -83,7 +83,7 @@ function main(mesh_partition,distribute,el_size,path) solver = ElasticitySolver(V) state_map = AffineFEStateMap( - a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N; + a,l,U,V,V_φ,U_reg,φh; assem_U = SparseMatrixAssembler(Tm,Tv,U,V), assem_adjoint = SparseMatrixAssembler(Tm,Tv,V,U), assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg), diff --git a/scripts/MPI/inverse_homenisation_ALM.jl b/scripts/MPI/inverse_homenisation_ALM.jl index 5946f58..0e36f2f 100644 --- a/scripts/MPI/inverse_homenisation_ALM.jl +++ b/scripts/MPI/inverse_homenisation_ALM.jl @@ -62,16 +62,16 @@ function main(mesh_partition,distribute,el_size,path) 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] + a(u,v,φ) = ∫((I ∘ φ) * C ⊙ ε(u) ⊙ ε(v))dΩ + l = [(v,φ) -> ∫(-(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Ω + Cᴴ(r,s,u,φ) = ∫((I ∘ φ)*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ εᴹ[s]))dΩ + dCᴴ(r,s,q,u,φ) = ∫(-q*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ (ε(u[s])+εᴹ[s]))*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + κ(u,φ) = -1/4*(Cᴴ(1,1,u,φ)+Cᴴ(2,2,u,φ)+2*Cᴴ(1,2,u,φ)) + dκ(q,u,φ) = -1/4*(dCᴴ(1,1,q,u,φ)+dCᴴ(2,2,q,u,φ)+2*dCᴴ(1,2,q,u,φ)) + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol(q,u,φ) = ∫(-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) @@ -82,7 +82,7 @@ function main(mesh_partition,distribute,el_size,path) solver = ElasticitySolver(V) state_map = RepeatingAffineFEStateMap( - 3,a,l,U,V,V_φ,U_reg,φh,dΩ; + 3,a,l,U,V,V_φ,U_reg,φh; assem_U = SparseMatrixAssembler(Tm,Tv,U,V), assem_adjoint = SparseMatrixAssembler(Tm,Tv,V,U), assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg), diff --git a/scripts/MPI/nonlinear_thermal_compliance_ALM.jl b/scripts/MPI/nonlinear_thermal_compliance_ALM.jl index beccd0e..0991d79 100644 --- a/scripts/MPI/nonlinear_thermal_compliance_ALM.jl +++ b/scripts/MPI/nonlinear_thermal_compliance_ALM.jl @@ -68,18 +68,18 @@ function main(mesh_partition,distribute,el_size,path) κ0 = 1 ξ = -1 κ(u) = κ0*(exp(ξ*u)) - res(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(v))dΩ - ∫(v)dΓ_N + res(u,v,φ) = ∫((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Ω + J(u,φ) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(u))dΩ + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol(q,u,φ) = ∫(-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) + state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh) pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dC=[dVol]) ## Hilbertian extension-regularisation problems diff --git a/scripts/MPI/thermal_compliance_ALM.jl b/scripts/MPI/thermal_compliance_ALM.jl index 9c78fbb..1e9ea3d 100644 --- a/scripts/MPI/thermal_compliance_ALM.jl +++ b/scripts/MPI/thermal_compliance_ALM.jl @@ -63,20 +63,20 @@ function main(mesh_partition,distribute,path) 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 + a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ + l(v,φ) = ∫(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Ω + J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ + dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ; + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol(q,u,φ) = ∫(-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 = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) + state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh) pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol]) ## Hilbertian extension-regularisation problems diff --git a/scripts/MPI/thermal_compliance_ALM_PETSc.jl b/scripts/MPI/thermal_compliance_ALM_PETSc.jl index 9248276..c5fb28e 100644 --- a/scripts/MPI/thermal_compliance_ALM_PETSc.jl +++ b/scripts/MPI/thermal_compliance_ALM_PETSc.jl @@ -63,14 +63,14 @@ function main(mesh_partition,distribute,el_size,path) 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 + a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ + l(v,φ) = ∫(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Ω + J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ + dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ; + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol(q,u,φ) = ∫(-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) @@ -81,7 +81,7 @@ function main(mesh_partition,distribute,el_size,path) solver = PETScLinearSolver() state_map = AffineFEStateMap( - a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N; + a,l,U,V,V_φ,U_reg,φh; assem_U = SparseMatrixAssembler(Tm,Tv,U,V), assem_adjoint = SparseMatrixAssembler(Tm,Tv,V,U), assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg), diff --git a/scripts/Serial/elastic_compliance_ALM.jl b/scripts/Serial/elastic_compliance_ALM.jl index 9f0b702..1e400ee 100644 --- a/scripts/Serial/elastic_compliance_ALM.jl +++ b/scripts/Serial/elastic_compliance_ALM.jl @@ -60,20 +60,20 @@ function main(path="./results/elastic_compliance_ALM/") interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ)) I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ - a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ - l(v,φ,dΩ,dΓ_N) = ∫(v⋅g)dΓ_N + a(u,v,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + l(v,φ) = ∫(v⋅g)dΓ_N ## Optimisation functionals - J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(u)))dΩ - dJ(q,u,φ,dΩ,dΓ_N) = ∫((C ⊙ ε(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Ω + J(u,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(u)))dΩ + dJ(q,u,φ) = ∫((C ⊙ ε(u) ⊙ ε(u))*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ; + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol(q,u,φ) = ∫(-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 = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) + state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh) pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol]) ## Hilbertian extension-regularisation problems diff --git a/scripts/Serial/elastic_compliance_HPM.jl b/scripts/Serial/elastic_compliance_HPM.jl index 9f787c3..eea2c88 100644 --- a/scripts/Serial/elastic_compliance_HPM.jl +++ b/scripts/Serial/elastic_compliance_HPM.jl @@ -61,20 +61,20 @@ function main(path="./results/elastic_compliance_HPM/") interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ)) I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ - a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ - l(v,φ,dΩ,dΓ_N) = ∫(v⋅g)dΓ_N + a(u,v,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + l(v,φ) = ∫(v⋅g)dΓ_N ## Optimisation functionals - J = (u,φ,dΩ,dΓ_N) -> ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(u)))dΩ - dJ = (q,u,φ,dΩ,dΓ_N) -> ∫((C ⊙ ε(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Ω + J = (u,φ) -> ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(u)))dΩ + dJ = (q,u,φ) -> ∫((C ⊙ ε(u) ⊙ ε(u))*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ; + Vol = (u,φ) -> ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol = (q,u,φ) -> ∫(-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 = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) + state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh) pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol]) ## Hilbertian extension-regularisation problems diff --git a/scripts/Serial/elastic_compliance_LM.jl b/scripts/Serial/elastic_compliance_LM.jl index b09e501..e4c090d 100644 --- a/scripts/Serial/elastic_compliance_LM.jl +++ b/scripts/Serial/elastic_compliance_LM.jl @@ -58,19 +58,19 @@ function main(path="./results/elastic_compliance_LM/") interp = SmoothErsatzMaterialInterpolation(η = η_coeff*maximum(el_Δ)) I,H,DH,ρ = interp.I,interp.H,interp.DH,interp.ρ - a(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ - l(v,φ,dΩ,dΓ_N) = ∫(v⋅g)dΓ_N + a(u,v,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + l(v,φ) = ∫(v⋅g)dΓ_N ## Optimisation functionals ξ = 2 - J = (u,φ,dΩ,dΓ_N) -> ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(u)) + ξ*(ρ ∘ φ))dΩ - dJ = (q,u,φ,dΩ,dΓ_N) -> ∫((- ξ + C ⊙ ε(u) ⊙ ε(u))*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ; + J = (u,φ) -> ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(u)) + ξ*(ρ ∘ φ))dΩ + dJ = (q,u,φ) -> ∫((- ξ + C ⊙ ε(u) ⊙ ε(u))*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 = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) + state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh) pcfs = PDEConstrainedFunctionals(J,state_map,analytic_dJ=dJ) ## Hilbertian extension-regularisation problems diff --git a/scripts/Serial/hyperelastic_compliance_ALM.jl b/scripts/Serial/hyperelastic_compliance_ALM.jl index dce96bd..42ef224 100644 --- a/scripts/Serial/hyperelastic_compliance_ALM.jl +++ b/scripts/Serial/hyperelastic_compliance_ALM.jl @@ -62,18 +62,18 @@ function main(path="./results/hyperelastic_compliance_ALM/") E(F) = 0.5*( F' ⋅ F - one(F) ) Σ(∇u) = λ*tr(E(F(∇u)))*one(∇u)+2*μ*E(F(∇u)) T(∇u) = F(∇u) ⋅ Σ(∇u) - res(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*((T ∘ ∇(u)) ⊙ ∇(v)))*dΩ - ∫(g⋅v)dΓ_N + res(u,v,φ) = ∫((I ∘ φ)*((T ∘ ∇(u)) ⊙ ∇(v)))*dΩ - ∫(g⋅v)dΓ_N ## Optimisation functionals - J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*((T ∘ ∇(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Ω + J(u,φ) = ∫((I ∘ φ)*((T ∘ ∇(u)) ⊙ ∇(u)))dΩ + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol(q,u,φ) = ∫(-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) + state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh) pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dC=[dVol]) ## Hilbertian extension-regularisation problems diff --git a/scripts/Serial/hyperelastic_compliance_neohook_ALM.jl b/scripts/Serial/hyperelastic_compliance_neohook_ALM.jl index eca9d1c..8927fcd 100644 --- a/scripts/Serial/hyperelastic_compliance_neohook_ALM.jl +++ b/scripts/Serial/hyperelastic_compliance_neohook_ALM.jl @@ -76,18 +76,18 @@ function main(path="./results/hyperelastic_compliance_neohook_ALM/") # Cauchy stress tensor and residual σ(∇u) = (1.0/J(F(∇u)))*F(∇u)⋅S(∇u)⋅(F(∇u))' - res(u,v,φ,dΩ,dΓ_N) = ∫( (I ∘ φ)*((dE∘(∇(v),∇(u))) ⊙ (S∘∇(u))) )*dΩ - ∫(g⋅v)dΓ_N + res(u,v,φ) = ∫( (I ∘ φ)*((dE∘(∇(v),∇(u))) ⊙ (S∘∇(u))) )*dΩ - ∫(g⋅v)dΓ_N ## Optimisation functionals - Obj(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*((dE∘(∇(u),∇(u))) ⊙ (S∘∇(u))))dΩ - Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; - dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + Obj(u,φ) = ∫((I ∘ φ)*((dE∘(∇(u),∇(u))) ⊙ (S∘∇(u))))dΩ + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol(q,u,φ) = ∫(-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) + state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh) pcfs = PDEConstrainedFunctionals(Obj,[Vol],state_map,analytic_dC=[dVol]) ## Hilbertian extension-regularisation problems diff --git a/scripts/Serial/inverse_homenisation_ALM.jl b/scripts/Serial/inverse_homenisation_ALM.jl index 972d432..0689e11 100644 --- a/scripts/Serial/inverse_homenisation_ALM.jl +++ b/scripts/Serial/inverse_homenisation_ALM.jl @@ -58,22 +58,22 @@ function main(path="./results/inverse_homenisation_ALM/") 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] + a(u,v,φ) = ∫((I ∘ φ) * C ⊙ ε(u) ⊙ ε(v) )dΩ + l = [(v,φ) -> ∫(-(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Ω + Cᴴ(r,s,u,φ) = ∫((I ∘ φ)*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ εᴹ[s]))dΩ + dCᴴ(r,s,q,u,φ) = ∫(-q*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ (ε(u[s])+εᴹ[s]))*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + κ(u,φ) = -1/4*(Cᴴ(1,1,u,φ)+Cᴴ(2,2,u,φ)+2*Cᴴ(1,2,u,φ)) + dκ(q,u,φ) = -1/4*(dCᴴ(1,1,q,u,φ)+dCᴴ(2,2,q,u,φ)+2*dCᴴ(1,2,q,u,φ)) + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ + dVol(q,u,φ) = ∫(-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Ω) + state_map = RepeatingAffineFEStateMap(3,a,l,U,V,V_φ,U_reg,φh) pcfs = PDEConstrainedFunctionals(κ,[Vol],state_map;analytic_dJ=dκ,analytic_dC=[dVol]) ## Hilbertian extension-regularisation problems diff --git a/scripts/Serial/inverter_ALM.jl b/scripts/Serial/inverter_ALM.jl index ad2d878..d1db290 100644 --- a/scripts/Serial/inverter_ALM.jl +++ b/scripts/Serial/inverter_ALM.jl @@ -74,21 +74,21 @@ function main(path="./results/inverter_ALM/") 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 + a(u,v,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + ∫(ks*(u⋅v))dΓ_out + l(v,φ) = ∫(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 + J(u,φ) = ∫((u⋅e₁)/vol_Γ_in)dΓ_in + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ + dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + UΓ_out(u,φ) = ∫((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) + state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh) pcfs = PDEConstrainedFunctionals(J,[Vol,UΓ_out],state_map,analytic_dC=[dVol,nothing]) ## Hilbertian extension-regularisation problems diff --git a/scripts/Serial/inverter_HPM.jl b/scripts/Serial/inverter_HPM.jl index fbe1312..cc20629 100644 --- a/scripts/Serial/inverter_HPM.jl +++ b/scripts/Serial/inverter_HPM.jl @@ -75,21 +75,21 @@ function main(path="./results/inverter_HPM/") 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 + a(u,v,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + ∫(ks*(u⋅v))dΓ_out + l(v,φ) = ∫(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 + J(u,φ) = ∫((u⋅e₁)/vol_Γ_in)dΓ_in + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ + dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + UΓ_out(u,φ) = ∫((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) + state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh) pcfs = PDEConstrainedFunctionals(J,[Vol,UΓ_out],state_map,analytic_dC=[dVol,nothing]) ## Hilbertian extension-regularisation problems diff --git a/scripts/Serial/nonlinear_thermal_compliance_ALM.jl b/scripts/Serial/nonlinear_thermal_compliance_ALM.jl index 2bf3248..b69a31b 100644 --- a/scripts/Serial/nonlinear_thermal_compliance_ALM.jl +++ b/scripts/Serial/nonlinear_thermal_compliance_ALM.jl @@ -65,18 +65,18 @@ function main(path="./results/nonlinear_thermal_compliance_ALM/") κ0 = 1 ξ = -1 κ(u) = κ0*(exp(ξ*u)) - res(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(v))dΩ - ∫(v)dΓ_N + res(u,v,φ) = ∫((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Ω + J(u,φ) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(u))dΩ + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol(q,u,φ) = ∫(-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) + state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh) pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dC=[dVol]) ## Hilbertian extension-regularisation problems diff --git a/scripts/Serial/thermal_compliance_ALM.jl b/scripts/Serial/thermal_compliance_ALM.jl index aaf45b2..a013cfa 100644 --- a/scripts/Serial/thermal_compliance_ALM.jl +++ b/scripts/Serial/thermal_compliance_ALM.jl @@ -61,20 +61,20 @@ function main(path="./results/thermal_compliance_ALM/") 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 + a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ + l(v,φ) = ∫(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Ω + J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ + dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ; + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol(q,u,φ) = ∫(-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 = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) + state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh) pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol]) ## Hilbertian extension-regularisation problems diff --git a/scripts/Serial/thermal_compliance_ALM_higherorderlsf.jl b/scripts/Serial/thermal_compliance_ALM_higherorderlsf.jl index 51f8e1a..4720998 100644 --- a/scripts/Serial/thermal_compliance_ALM_higherorderlsf.jl +++ b/scripts/Serial/thermal_compliance_ALM_higherorderlsf.jl @@ -67,20 +67,20 @@ function main(path="./results/thermal_compliance_ALM_higherorderlsf/") 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 + a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ + l(v,φ) = ∫(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Ω + J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ + dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ; + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol(q,u,φ) = ∫(-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 = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) + state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh) pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol]) ## Hilbertian extension-regularisation problems diff --git a/src/ChainRules.jl b/src/ChainRules.jl index 4661c6b..bcedbab 100644 --- a/src/ChainRules.jl +++ b/src/ChainRules.jl @@ -1,93 +1,38 @@ -abstract type AbstractIntegrandWithMeasure end - -""" - struct IntegrandWithMeasure{A,B<:Tuple} - -A wrapper to enable serial or parallel partial differentation of an -integral `F` using `Gridap.gradient`. This is required to allow automatic -differentation with `DistributedMeasure`. - -# Properties -- `F :: A`: A function that returns a `DomainContribution` or `DistributedDomainContribution`. -- `dΩ :: B`: A tuple of measures. -""" -struct IntegrandWithMeasure{A,B<:Tuple} <: AbstractIntegrandWithMeasure - F :: A - dΩ :: B -end - -""" - (F::IntegrandWithMeasure)(args...) - -Evaluate `F.F` given arguments `args`. """ -(F::IntegrandWithMeasure)(args...) = F.F(args...,F.dΩ...) + Gridap.gradient(F,uh::Vector,K::Int) -""" - Gridap.gradient(F::IntegrandWithMeasure,uh::Vector,K::Int) - -Given an an `IntegrandWithMeasure` `F` and a vector of `FEFunctions` `uh` (excluding measures) -evaluate the partial derivative of `F.F` with respect to `uh[K]`. +Given a function `F` that returns a DomainContribution when called, and a vector of +`FEFunctions` `uh`, evaluate the partial derivative of `F` with respect to `uh[K]`. # Example Suppose `uh` and `φh` are FEFunctions with measures `dΩ` and `dΓ_N`. Then the partial derivative of a function `J` wrt to `φh` is computed via ```` -J(u,φ,dΩ,dΓ_N) = ∫(f(u,φ))dΩ + ∫(g(u,φ))dΓ_N -J_iwm = IntegrandWithMeasure(J,(dΩ,dΓ_N)) -∂J∂φh = ∇(J_iwm,[uh,φh],2) +J(u,φ) = ∫(f(u,φ))dΩ + ∫(g(u,φ))dΓ_N +∂J∂φh = ∇(J,[uh,φh],2) ```` where `f` and `g` are user defined. """ -function Gridap.gradient(F::IntegrandWithMeasure,uh::Vector{<:FEFunction},K::Int) +function Gridap.gradient(F,uh::Vector{<:CellField},K::Int) @check 0 < K <= length(uh) - _f(uk) = F.F(uh[1:K-1]...,uk,uh[K+1:end]...,F.dΩ...) + _f(uk) = F(uh[1:K-1]...,uk,uh[K+1:end]...) return Gridap.gradient(_f,uh[K]) end -function Gridap.gradient(F::IntegrandWithMeasure,uh::Vector,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 - contribs = map(local_measures,local_fields) do dΩ,lf - # TODO: Remove second term below, this is a fix for the problem discussed in - # https://github.com/zjwegert/GridapTopOpt.jl/issues/46 - _f = u -> F.F(lf[1:K-1]...,u,lf[K+1:end]...,dΩ...) #+ ∑(∫(0)dΩ[i] for i = 1:length(dΩ)) - return Gridap.Fields.gradient(_f,lf[K]) - end - return DistributedDomainContribution(contribs) -end - -Gridap.gradient(F::IntegrandWithMeasure,uh) = Gridap.gradient(F,[uh],1) - """ - Gridap.jacobian(F::IntegrandWithMeasure,uh::Vector,K::Int) + Gridap.jacobian(F,uh::Vector,K::Int) -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]`. +Given a function `F` that returns a DomainContribution when called, and a +vector of `FEFunctions` or `CellField` `uh`, evaluate the Jacobian +`F` with respect to `uh[K]`. """ -function Gridap.jacobian(F::IntegrandWithMeasure,uh::Vector{<:Union{FEFunction,CellField}},K::Int) +function Gridap.jacobian(F,uh::Vector{<:CellField},K::Int) @check 0 < K <= length(uh) - _f(uk) = F.F(uh[1:K-1]...,uk,uh[K+1:end]...,F.dΩ...) + _f(uk) = F(uh[1:K-1]...,uk,uh[K+1:end]...) return Gridap.jacobian(_f,uh[K]) end -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 - contribs = map(local_measures,local_fields) do dΩ,lf - _f = u -> F.F(lf[1:K-1]...,u,lf[K+1:end]...,dΩ...) - return Gridap.jacobian(_f,lf[K]) - end - return DistributedDomainContribution(contribs) -end - -Gridap.jacobian(F::IntegrandWithMeasure,uh) = Gridap.jacobian(F,[uh],1) - function GridapDistributed.to_parray_of_arrays(a::NTuple{N,T}) where {N,T<:DebugArray} indices = linear_indices(first(a)) map(indices) do i @@ -107,21 +52,21 @@ function GridapDistributed.to_parray_of_arrays(a::NTuple{N,T}) where {N,T<:MPIAr end """ - struct StateParamIntegrandWithMeasure{A<:IntegrandWithMeasure,B,C,D} + struct StateParamIntegrandWithMeasure{A,B,C,D} -A wrapper to handle partial differentation of an [`IntegrandWithMeasure`](@ref) +A wrapper to handle partial differentation of a function F of a specific form (see below) in a `ChainRules.jl` compatible way with caching. # Assumptions -We assume that we have a `IntegrandWithMeasure` of the following form: +We assume that we have a function F of the following form: -`F(u,φ,dΩ₁,dΩ₂,...) = ∫(f(u,φ))dΩ₁ + ∫(g(u,φ))dΩ₂ + ...,`. +`F(u,φ) = ∫(f(u,φ))dΩ₁ + ∫(g(u,φ))dΩ₂ + ...,`. where `u` and `φ` are each expected to inherit from `Union{FEFunction,MultiFieldFEFunction}` or the GridapDistributed equivalent. """ -struct StateParamIntegrandWithMeasure{A<:AbstractIntegrandWithMeasure,B,C,D} +struct StateParamIntegrandWithMeasure{A,B,C,D} F :: A spaces :: B assems :: C @@ -129,17 +74,15 @@ struct StateParamIntegrandWithMeasure{A<:AbstractIntegrandWithMeasure,B,C,D} end """ - StateParamIntegrandWithMeasure(F::IntegrandWithMeasure,U::FESpace,V_φ::FESpace, + StateParamIntegrandWithMeasure(F,U::FESpace,V_φ::FESpace, U_reg::FESpace,assem_U::Assembler,assem_deriv::Assembler) Create an instance of `StateParamIntegrandWithMeasure`. """ function StateParamIntegrandWithMeasure( - F::AbstractIntegrandWithMeasure, - U::FESpace,V_φ::FESpace,U_reg::FESpace, + F,U::FESpace,V_φ::FESpace,U_reg::FESpace, assem_U::Assembler,assem_deriv::Assembler ) - # TODO: Find alternative to this later, is a fix for grad of embedded FEs φ₀, u₀ = interpolate(x->-sqrt((x[1]-1/2)^2+(x[2]-1/2)^2)+0.2,V_φ), zero(U) ∂j∂u_vecdata = collect_cell_vector(U,∇(F,[u₀,φ₀],1)) ∂j∂φ_vecdata = collect_cell_vector(U_reg,∇(F,[u₀,φ₀],2)) @@ -215,13 +158,6 @@ Return the solution/state `u` to the FE problem. """ get_state(::AbstractFEStateMap) = @abstractmethod -""" - get_measure(m::AbstractFEStateMap) - -Return the measures associated with the FE problem. -""" -get_measure(::AbstractFEStateMap) = @abstractmethod - """ get_spaces(m::AbstractFEStateMap) @@ -400,19 +336,7 @@ function ChainRulesCore.rrule(φ_to_u::AbstractFEStateMap,φ::AbstractVector) return ChainRulesCore.rrule(φ_to_u,φh) end -""" - StateParamIntegrandWithMeasure(f::Function,φ_to_u::AbstractFEStateMap) - -Create an instance of `StateParamIntegrandWithMeasure` given a `f` and -`φ_to_u`. -""" -function StateParamIntegrandWithMeasure(f::Function,φ_to_u::AbstractFEStateMap) - dΩ = get_measure(φ_to_u) - F = IntegrandWithMeasure(f,dΩ) - StateParamIntegrandWithMeasure(F,φ_to_u) -end - -function StateParamIntegrandWithMeasure(F::AbstractIntegrandWithMeasure,φ_to_u::AbstractFEStateMap) +function StateParamIntegrandWithMeasure(F::Function,φ_to_u::AbstractFEStateMap) U,V,V_φ,U_reg = φ_to_u.spaces assem_deriv = get_deriv_assembler(φ_to_u) assem_U = get_pde_assembler(φ_to_u) @@ -427,8 +351,8 @@ element operators `AffineFEOperator`. # Parameters -- `biform::A`: `IntegrandWithMeasure` defining the bilinear form. -- `liform::B`: `IntegrandWithMeasure` defining the linear form. +- `biform::A`: `Function` defining the bilinear form. +- `liform::B`: `Function` defining the linear form. - `spaces::C`: `Tuple` of finite element spaces. - `plb_caches::D`: A cache for the pullback operator. - `fwd_caches::E`: A cache for the forward problem. @@ -445,7 +369,7 @@ struct AffineFEStateMap{A,B,C,D,E,F} <: AbstractFEStateMap @doc """ AffineFEStateMap( a::Function,l::Function, - U,V,V_φ,U_reg,φh,dΩ...; + U,V,V_φ,U_reg,φh; assem_U = SparseMatrixAssembler(U,V), assem_adjoint = SparseMatrixAssembler(V,U), assem_deriv = SparseMatrixAssembler(U_reg,U_reg), @@ -460,8 +384,8 @@ struct AffineFEStateMap{A,B,C,D,E,F} <: AbstractFEStateMap Optional arguments enable specification of assemblers and linear solvers. """ function AffineFEStateMap( - biform::AbstractIntegrandWithMeasure,liform::AbstractIntegrandWithMeasure, - U,V,V_φ,U_reg,φh,dΩ...; + biform::Function,liform::Function, + U,V,V_φ,U_reg,φh; assem_U = SparseMatrixAssembler(U,V), assem_adjoint = SparseMatrixAssembler(V,U), assem_deriv = SparseMatrixAssembler(U_reg,U_reg), @@ -497,15 +421,8 @@ struct AffineFEStateMap{A,B,C,D,E,F} <: AbstractFEStateMap end end -function AffineFEStateMap(a::Function,l::Function,U,V,V_φ,U_reg,φh,dΩ...;kwargs...) - biform = IntegrandWithMeasure(a,dΩ) - liform = IntegrandWithMeasure(l,dΩ) - AffineFEStateMap(biform,liform,U,V,V_φ,U_reg,φh,dΩ;kwargs...) -end - # Getters get_state(m::AffineFEStateMap) = FEFunction(get_trial_space(m),m.fwd_caches[4]) -get_measure(m::AffineFEStateMap) = m.biform.dΩ get_spaces(m::AffineFEStateMap) = m.spaces get_assemblers(m::AffineFEStateMap) = (m.fwd_caches[6],m.plb_caches[2],m.adj_caches[4]) @@ -549,7 +466,7 @@ element operators. # Parameters -- `res::A`: an `IntegrandWithMeasure` defining the residual of the problem. +- `res::A`: a `Function` defining the residual of the problem. - `jac::B`: a `Function` defining Jacobian of the residual. - `spaces::C`: `Tuple` of finite element spaces. - `plb_caches::D`: A cache for the pullback operator. @@ -566,7 +483,7 @@ struct NonlinearFEStateMap{A,B,C,D,E,F} <: AbstractFEStateMap @doc """ NonlinearFEStateMap( - res::Function,U,V,V_φ,U_reg,φh,dΩ...; + res::Function,U,V,V_φ,U_reg,φh; assem_U = SparseMatrixAssembler(U,V), assem_adjoint = SparseMatrixAssembler(V,U), assem_deriv = SparseMatrixAssembler(U_reg,U_reg), @@ -581,7 +498,7 @@ struct NonlinearFEStateMap{A,B,C,D,E,F} <: AbstractFEStateMap Optional arguments enable specification of assemblers, nonlinear solver, and adjoint (linear) solver. """ function NonlinearFEStateMap( - res::AbstractIntegrandWithMeasure,jac,U,V,V_φ,U_reg,φh,dΩ...; + res::Function,jac::Function,U,V,V_φ,U_reg,φh; assem_U = SparseMatrixAssembler(U,V), assem_adjoint = SparseMatrixAssembler(V,U), assem_deriv = SparseMatrixAssembler(U_reg,U_reg), @@ -616,23 +533,14 @@ struct NonlinearFEStateMap{A,B,C,D,E,F} <: AbstractFEStateMap end end -function NonlinearFEStateMap(res::AbstractIntegrandWithMeasure,U,V,V_φ,U_reg,φh,dΩ...;jac::Nothing=nothing,kwargs...) - jacf = (u,du,v,φh) -> jacobian(res,[u,v,φh],1) - NonlinearFEStateMap(res,jacf,U,V,V_φ,U_reg,φh,dΩ;kwargs...) -end - -function NonlinearFEStateMap(res::Function,U,V,V_φ,U_reg,φh,dΩ...;jac=nothing,kwargs...) - _res = IntegrandWithMeasure(res,dΩ) +function NonlinearFEStateMap(res::Function,U,V,V_φ,U_reg,φh;jac=nothing,kwargs...) if isnothing(jac) - _jac = (u,du,v,φh) -> jacobian(_res,[u,v,φh],1) - else - _jac = IntegrandWithMeasure(jac,dΩ) + jac = (u,du,v,φh) -> jacobian(res,[u,v,φh],1) end - NonlinearFEStateMap(_res,_jac,U,V,V_φ,U_reg,φh,dΩ;kwargs...) + NonlinearFEStateMap(res,jac,U,V,V_φ,U_reg,φh;kwargs...) end get_state(m::NonlinearFEStateMap) = FEFunction(get_trial_space(m),m.fwd_caches[3]) -get_measure(m::NonlinearFEStateMap) = m.res.dΩ get_spaces(m::NonlinearFEStateMap) = m.spaces get_assemblers(m::NonlinearFEStateMap) = (m.fwd_caches[4],m.plb_caches[2],m.adj_caches[4]) @@ -676,8 +584,8 @@ a single bilinear form. # Parameters -- `biform`: `IntegrandWithMeasure` defining the bilinear form. -- `liform`: A vector of `IntegrandWithMeasure` defining the linear forms. +- `biform`: A `Function` defining the bilinear form. +- `liform`: A vector of `Function` defining the linear forms. - `spaces`: Repeated finite element spaces. - `spaces_0`: Original finite element spaces that are being repeated. - `plb_caches`: A cache for the pullback operator. @@ -696,7 +604,7 @@ struct RepeatingAffineFEStateMap{A,B,C,D,E,F,G} <: AbstractFEStateMap @doc """ RepeatingAffineFEStateMap( nblocks::Int,a::Function,l::Vector{<:Function}, - U0,V0,V_φ,U_reg,φh,dΩ...; + U0,V0,V_φ,U_reg,φh; assem_U = SparseMatrixAssembler(U0,V0), assem_adjoint = SparseMatrixAssembler(V0,U0), assem_deriv = SparseMatrixAssembler(U_reg,U_reg), @@ -717,8 +625,8 @@ struct RepeatingAffineFEStateMap{A,B,C,D,E,F,G} <: AbstractFEStateMap where each field corresponds to an entry in the vector of linear forms """ function RepeatingAffineFEStateMap( - nblocks::Int,biform::AbstractIntegrandWithMeasure,liforms::Vector{<:AbstractIntegrandWithMeasure}, - U0,V0,V_φ,U_reg,φh,dΩ...; + nblocks::Int,biform::Function,liforms::Vector{<:Function}, + U0,V0,V_φ,U_reg,φh; assem_U = SparseMatrixAssembler(U0,V0), assem_adjoint = SparseMatrixAssembler(V0,U0), assem_deriv = SparseMatrixAssembler(U_reg,U_reg), @@ -766,14 +674,6 @@ struct RepeatingAffineFEStateMap{A,B,C,D,E,F,G} <: AbstractFEStateMap end end -function RepeatingAffineFEStateMap( - nblocks::Int,a::Function,l::Vector{<:Function}, - U0,V0,V_φ,U_reg,φh,dΩ...;kwargs...) - biform = IntegrandWithMeasure(a,dΩ) - liforms = map(li -> IntegrandWithMeasure(li,dΩ),l) - RepeatingAffineFEStateMap(nblocks::Int,biform,liforms,U0,V0,V_φ,U_reg,φh,dΩ...;kwargs...) -end - const MultiFieldSpaceTypes = Union{<:MultiField.MultiFieldFESpace,<:GridapDistributed.DistributedMultiFieldFESpace} function repeat_spaces(nblocks::Integer,U0::FESpace,V0::FESpace) @@ -840,7 +740,6 @@ function repeated_allocate_in_domain(nblocks::Integer,M::AbstractBlockMatrix) end get_state(m::RepeatingAffineFEStateMap) = FEFunction(get_trial_space(m),m.fwd_caches[4]) -get_measure(m::RepeatingAffineFEStateMap) = m.biform.dΩ get_spaces(m::RepeatingAffineFEStateMap) = m.spaces get_assemblers(m::RepeatingAffineFEStateMap) = (m.fwd_caches[8],m.plb_caches[2],m.adj_caches[4]) @@ -1046,7 +945,6 @@ function Fields.evaluate!(pcf::PDEConstrainedFunctionals,φh) U_reg = get_deriv_space(pcf.state_map) deriv_assem = get_deriv_assembler(pcf.state_map) - dΩ = get_measure(pcf.state_map) ## Foward problem u, u_pullback = rrule(pcf.state_map,φh) @@ -1062,13 +960,6 @@ function Fields.evaluate!(pcf::PDEConstrainedFunctionals,φh) return j_val end function ∇!(F::StateParamIntegrandWithMeasure,dF,dF_analytic::Function) - # Analytic shape derivative - j_val = F(uh,φh) - _dF(q) = dF_analytic(q,uh,φh,dΩ...) - assemble_vector!(_dF,dF,deriv_assem,U_reg) - return j_val - end - function ∇!(F::StateParamIntegrandWithMeasure,dF,dF_analytic::AbstractIntegrandWithMeasure) # Analytic shape derivative j_val = F(uh,φh) _dF(q) = dF_analytic(q,uh,φh) @@ -1089,10 +980,6 @@ end # IO -function Base.show(io::IO,object::IntegrandWithMeasure) - print(io,"$(nameof(typeof(object)))") -end - function Base.show(io::IO,object::StateParamIntegrandWithMeasure) print(io,"$(nameof(typeof(object)))") end diff --git a/src/GridapTopOpt.jl b/src/GridapTopOpt.jl index 19eaa8d..ba410c0 100644 --- a/src/GridapTopOpt.jl +++ b/src/GridapTopOpt.jl @@ -23,7 +23,8 @@ using GridapDistributed using GridapDistributed: DistributedDiscreteModel, DistributedTriangulation, DistributedFESpace, DistributedDomainContribution, to_parray_of_arrays, allocate_in_domain, DistributedCellField, DistributedMultiFieldCellField, - DistributedMultiFieldFEBasis, BlockPMatrix, BlockPVector, change_ghost + DistributedMultiFieldFEBasis, BlockPMatrix, BlockPVector, change_ghost, + gradient, jacobian using PartitionedArrays using PartitionedArrays: getany, tuple_of_arrays, matching_ghost_indices diff --git a/test/mpi/InverseHomogenisationALMTests.jl b/test/mpi/InverseHomogenisationALMTests.jl index dcac064..f124c22 100644 --- a/test/mpi/InverseHomogenisationALMTests.jl +++ b/test/mpi/InverseHomogenisationALMTests.jl @@ -61,22 +61,22 @@ function main(distribute,mesh_partition;AD) 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] + a(u,v,φ) = ∫((I ∘ φ) * C ⊙ ε(u) ⊙ ε(v) )dΩ + l = [(v,φ) -> ∫(-(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Ω + Cᴴ(r,s,u,φ) = ∫((I ∘ φ)*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ εᴹ[s]))dΩ + dCᴴ(r,s,q,u,φ) = ∫(-q*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ (ε(u[s])+εᴹ[s]))*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + κ(u,φ) = -1/4*(Cᴴ(1,1,u,φ)+Cᴴ(2,2,u,φ)+2*Cᴴ(1,2,u,φ)) + dκ(q,u,φ) = -1/4*(dCᴴ(1,1,q,u,φ)+dCᴴ(2,2,q,u,φ)+2*dCᴴ(1,2,q,u,φ)) + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ + dVol(q,u,φ) = ∫(-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Ω) + state_map = RepeatingAffineFEStateMap(3,a,l,U,V,V_φ,U_reg,φh) pcfs = if AD PDEConstrainedFunctionals(κ,[Vol],state_map;analytic_dJ=dκ,analytic_dC=[dVol]) else diff --git a/test/mpi/InverterHPMTests.jl b/test/mpi/InverterHPMTests.jl index f4c614d..0450dfa 100644 --- a/test/mpi/InverterHPMTests.jl +++ b/test/mpi/InverterHPMTests.jl @@ -78,21 +78,21 @@ function main(distribute,mesh_partition) 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 + a(u,v,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + ∫(ks*(u⋅v))dΓ_out + l(v,φ) = ∫(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 + J(u,φ) = ∫((u⋅e₁)/vol_Γ_in)dΓ_in + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ + dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + UΓ_out(u,φ) = ∫((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) + state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh) pcfs = PDEConstrainedFunctionals(J,[Vol,UΓ_out],state_map,analytic_dC=[dVol,nothing]) ## Hilbertian extension-regularisation problems diff --git a/test/mpi/NeohookAnalyticJacALMTests.jl b/test/mpi/NeohookAnalyticJacALMTests.jl index df8041a..bf9936c 100644 --- a/test/mpi/NeohookAnalyticJacALMTests.jl +++ b/test/mpi/NeohookAnalyticJacALMTests.jl @@ -85,21 +85,21 @@ function main(distribute,mesh_partition) # Cauchy stress tensor and residual σ(∇u) = (1.0/J(F(∇u)))*F(∇u)⋅S(∇u)⋅(F(∇u))' - res(u,v,φ,dΩ,dΓ_N) = ∫( (I ∘ φ)*((dE∘(∇(v),∇(u))) ⊙ (S∘∇(u))) )*dΩ - ∫(g⋅v)dΓ_N - jac_mat(u,du,v,φ,dΩ,dΓ_N) = ∫( (I ∘ φ)*(dE∘(∇(v),∇(u))) ⊙ (dS∘(∇(du),∇(u))) )*dΩ - jac_geo(u,du,v,φ,dΩ,dΓ_N) = ∫( (I ∘ φ)*∇(v) ⊙ ( (S∘∇(u))⋅∇(du) ) )*dΩ - jac(u,du,v,φ,dΩ,dΓ_N) = jac_mat(u,du,v,φ,dΩ,dΓ_N) + jac_geo(u,du,v,φ,dΩ,dΓ_N) + res(u,v,φ) = ∫( (I ∘ φ)*((dE∘(∇(v),∇(u))) ⊙ (S∘∇(u))) )*dΩ - ∫(g⋅v)dΓ_N + jac_mat(u,du,v,φ) = ∫( (I ∘ φ)*(dE∘(∇(v),∇(u))) ⊙ (dS∘(∇(du),∇(u))) )*dΩ + jac_geo(u,du,v,φ) = ∫( (I ∘ φ)*∇(v) ⊙ ( (S∘∇(u))⋅∇(du) ) )*dΩ + jac(u,du,v,φ) = jac_mat(u,du,v,φ) + jac_geo(u,du,v,φ) ## Optimisation functionals - J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*((dE∘(∇(u),∇(u))) ⊙ (S∘∇(u))))dΩ - Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; - dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + J(u,φ) = ∫((I ∘ φ)*((dE∘(∇(u),∇(u))) ⊙ (S∘∇(u))))dΩ + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol(q,u,φ) = ∫(-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;jac) + state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh;jac) pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dC=[dVol]) ## Hilbertian extension-regularisation problems diff --git a/test/mpi/NonlinearThermalComplianceALMTests.jl b/test/mpi/NonlinearThermalComplianceALMTests.jl index 1b63a03..134a4cb 100644 --- a/test/mpi/NonlinearThermalComplianceALMTests.jl +++ b/test/mpi/NonlinearThermalComplianceALMTests.jl @@ -68,18 +68,18 @@ function main(distribute,mesh_partition) κ0 = 1 ξ = -1 κ(u) = κ0*(exp(ξ*u)) - res(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(v))dΩ - ∫(v)dΓ_N + res(u,v,φ) = ∫((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Ω + J(u,φ) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(u))dΩ + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol(q,u,φ) = ∫(-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) + state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh) pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dC=[dVol]) ## Hilbertian extension-regularisation problems diff --git a/test/mpi/PZMultiFieldRepeatingStateTests.jl b/test/mpi/PZMultiFieldRepeatingStateTests.jl index dbe6f25..c083b7f 100644 --- a/test/mpi/PZMultiFieldRepeatingStateTests.jl +++ b/test/mpi/PZMultiFieldRepeatingStateTests.jl @@ -67,21 +67,21 @@ function main(distribute,mesh_partition;AD) Eⁱ = (VectorValue(1.0,0.0,), VectorValue(0.0,1.0)) - a((u,ϕ),(v,q),φ,dΩ) = ∫((I ∘ φ) * (1/k0*((C ⊙ ε(u)) ⊙ ε(v)) - + a((u,ϕ),(v,q),φ) = ∫((I ∘ φ) * (1/k0*((C ⊙ ε(u)) ⊙ ε(v)) - 1/α0*((-∇(ϕ) ⋅ e) ⊙ ε(v)) + -1/α0*((e ⋅² ε(u)) ⋅ -∇(q)) + -γ0/β0*((κ ⋅ -∇(ϕ)) ⋅ -∇(q))) )dΩ; - l_ε = [((v,q),φ,dΩ) -> ∫(((I ∘ φ) * (-C ⊙ εᴹ[i] ⊙ ε(v) + k0/α0*(e ⋅² εᴹ[i]) ⋅ -∇(q))))dΩ for i = 1:3]; - l_E = [((v,q),φ,dΩ) -> ∫((I ∘ φ) * ((Eⁱ[i] ⋅ e ⊙ ε(v) + k0/α0*(κ ⋅ Eⁱ[i]) ⋅ -∇(q))))dΩ for i = 1:2]; + l_ε = [((v,q),φ) -> ∫(((I ∘ φ) * (-C ⊙ εᴹ[i] ⊙ ε(v) + k0/α0*(e ⋅² εᴹ[i]) ⋅ -∇(q))))dΩ for i = 1:3]; + l_E = [((v,q),φ) -> ∫((I ∘ φ) * ((Eⁱ[i] ⋅ e ⊙ ε(v) + k0/α0*(κ ⋅ Eⁱ[i]) ⋅ -∇(q))))dΩ for i = 1:2]; l = [l_ε; l_E] - function Cᴴ(r,s,uϕ,φ,dΩ) + function Cᴴ(r,s,uϕ,φ) u_s = uϕ[2s-1]; ϕ_s = uϕ[2s] ∫(1/k0 * (I ∘ φ) * (((C ⊙ (1/k0*ε(u_s) + εᴹ[s])) ⊙ εᴹ[r]) - ((-1/α0*∇(ϕ_s) ⋅ e) ⊙ εᴹ[r])))dΩ; end - function DCᴴ(r,s,q,uϕ,φ,dΩ) + function DCᴴ(r,s,q,uϕ,φ) u_r = uϕ[2r-1]; ϕ_r = uϕ[2r] u_s = uϕ[2s-1]; ϕ_s = uϕ[2s] ∫(- 1/k0 * q * ( @@ -93,19 +93,19 @@ function main(distribute,mesh_partition;AD) )dΩ; end - Bᴴ(uϕ,φ,dΩ) = 1/4*(Cᴴ(1,1,uϕ,φ,dΩ)+Cᴴ(2,2,uϕ,φ,dΩ)+2*Cᴴ(1,2,uϕ,φ,dΩ)) - DBᴴ(q,uϕ,φ,dΩ) = 1/4*(DCᴴ(1,1,q,uϕ,φ,dΩ)+DCᴴ(2,2,q,uϕ,φ,dΩ)+2*DCᴴ(1,2,q,uϕ,φ,dΩ)) + Bᴴ(uϕ,φ) = 1/4*(Cᴴ(1,1,uϕ,φ)+Cᴴ(2,2,uϕ,φ)+2*Cᴴ(1,2,uϕ,φ)) + DBᴴ(q,uϕ,φ) = 1/4*(DCᴴ(1,1,q,uϕ,φ)+DCᴴ(2,2,q,uϕ,φ)+2*DCᴴ(1,2,q,uϕ,φ)) - J(uϕ,φ,dΩ) = -1*Bᴴ(uϕ,φ,dΩ) - DJ(q,uϕ,φ,dΩ) = -1*DBᴴ(q,uϕ,φ,dΩ) - C1(uϕ,φ,dΩ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; - DC1(q,uϕ,φ,dΩ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + J(uϕ,φ) = -1*Bᴴ(uϕ,φ) + DJ(q,uϕ,φ) = -1*DBᴴ(q,uϕ,φ) + C1(uϕ,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + DC1(q,uϕ,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function stencil = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators - state_map = RepeatingAffineFEStateMap(5,a,l,UP,VQ,V_φ,U_reg,φh,dΩ) + state_map = RepeatingAffineFEStateMap(5,a,l,UP,VQ,V_φ,U_reg,φh) pcfs = if AD PDEConstrainedFunctionals(J,[C1],state_map;analytic_dJ=DJ,analytic_dC=[DC1]) else diff --git a/test/mpi/ThermalComplianceALMTests.jl b/test/mpi/ThermalComplianceALMTests.jl index 9b7e0a6..dee3d38 100644 --- a/test/mpi/ThermalComplianceALMTests.jl +++ b/test/mpi/ThermalComplianceALMTests.jl @@ -63,20 +63,20 @@ function main(distribute,mesh_partition;order,AD) 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 + a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ + l(v,φ) = ∫(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Ω + J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ + dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ; + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol(q,u,φ) = ∫(-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 = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) + state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh) pcfs = if AD PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol]) else @@ -157,20 +157,20 @@ function main_3d(distribute,mesh_partition,;order) 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 + a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ + l(v,φ) = ∫(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Ω + J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ + dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ; + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol(q,u,φ) = ∫(-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) + state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh) pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol]) ## Hilbertian extension-regularisation problems diff --git a/test/mpi/ThermalComplianceALMTestsPETSc.jl b/test/mpi/ThermalComplianceALMTestsPETSc.jl index db1b6d4..b3832bd 100644 --- a/test/mpi/ThermalComplianceALMTestsPETSc.jl +++ b/test/mpi/ThermalComplianceALMTestsPETSc.jl @@ -63,14 +63,14 @@ function main(distribute,mesh_partition;order,AD) 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 + a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ + l(v,φ) = ∫(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Ω + J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ + dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ; + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol(q,u,φ) = ∫(-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) @@ -81,7 +81,7 @@ function main(distribute,mesh_partition;order,AD) solver = PETScLinearSolver() state_map = AffineFEStateMap( - a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N; + a,l,U,V,V_φ,U_reg,φh; assem_U = SparseMatrixAssembler(Tm,Tv,U,V), assem_adjoint = SparseMatrixAssembler(Tm,Tv,V,U), assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg), @@ -171,14 +171,14 @@ function main_3d(distribute,mesh_partition,;order) 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 + a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ + l(v,φ) = ∫(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Ω + J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ + dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ; + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol(q,u,φ) = ∫(-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) @@ -189,7 +189,7 @@ function main_3d(distribute,mesh_partition,;order) solver = PETScLinearSolver() state_map = AffineFEStateMap( - a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N; + a,l,U,V,V_φ,U_reg,φh; assem_U = SparseMatrixAssembler(Tm,Tv,U,V), assem_adjoint = SparseMatrixAssembler(Tm,Tv,V,U), assem_deriv = SparseMatrixAssembler(Tm,Tv,U_reg,U_reg), diff --git a/test/seq/InverseHomogenisationALMTests.jl b/test/seq/InverseHomogenisationALMTests.jl index b0b79d4..3d1a264 100644 --- a/test/seq/InverseHomogenisationALMTests.jl +++ b/test/seq/InverseHomogenisationALMTests.jl @@ -59,22 +59,22 @@ function main(;AD) 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] + a(u,v,φ) = ∫((I ∘ φ) * C ⊙ ε(u) ⊙ ε(v) )dΩ + l = [(v,φ) -> ∫(-(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Ω + Cᴴ(r,s,u,φ) = ∫((I ∘ φ)*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ εᴹ[s]))dΩ + dCᴴ(r,s,q,u,φ) = ∫(-q*(C ⊙ (ε(u[r])+εᴹ[r]) ⊙ (ε(u[s])+εᴹ[s]))*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + κ(u,φ) = -1/4*(Cᴴ(1,1,u,φ)+Cᴴ(2,2,u,φ)+2*Cᴴ(1,2,u,φ)) + dκ(q,u,φ) = -1/4*(dCᴴ(1,1,q,u,φ)+dCᴴ(2,2,q,u,φ)+2*dCᴴ(1,2,q,u,φ)) + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ + dVol(q,u,φ) = ∫(-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Ω) + state_map = RepeatingAffineFEStateMap(3,a,l,U,V,V_φ,U_reg,φh) pcfs = if AD PDEConstrainedFunctionals(κ,[Vol],state_map;analytic_dJ=dκ,analytic_dC=[dVol]) else diff --git a/test/seq/InverterHPMTests.jl b/test/seq/InverterHPMTests.jl index 06fffac..1511cb8 100644 --- a/test/seq/InverterHPMTests.jl +++ b/test/seq/InverterHPMTests.jl @@ -76,21 +76,21 @@ function main() 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 + a(u,v,φ) = ∫((I ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ + ∫(ks*(u⋅v))dΓ_out + l(v,φ) = ∫(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 + J(u,φ) = ∫((u⋅e₁)/vol_Γ_in)dΓ_in + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ + dVol(q,u,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + UΓ_out(u,φ) = ∫((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) + state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh) pcfs = PDEConstrainedFunctionals(J,[Vol,UΓ_out],state_map,analytic_dC=[dVol,nothing]) ## Hilbertian extension-regularisation problems diff --git a/test/seq/NeohookAnalyticJacALMTests.jl b/test/seq/NeohookAnalyticJacALMTests.jl index 8ef5d70..e4b8bc0 100644 --- a/test/seq/NeohookAnalyticJacALMTests.jl +++ b/test/seq/NeohookAnalyticJacALMTests.jl @@ -83,21 +83,21 @@ function main() # Cauchy stress tensor and residual σ(∇u) = (1.0/J(F(∇u)))*F(∇u)⋅S(∇u)⋅(F(∇u))' - res(u,v,φ,dΩ,dΓ_N) = ∫( (I ∘ φ)*((dE∘(∇(v),∇(u))) ⊙ (S∘∇(u))) )*dΩ - ∫(g⋅v)dΓ_N - jac_mat(u,du,v,φ,dΩ,dΓ_N) = ∫( (I ∘ φ)*(dE∘(∇(v),∇(u))) ⊙ (dS∘(∇(du),∇(u))) )*dΩ - jac_geo(u,du,v,φ,dΩ,dΓ_N) = ∫( (I ∘ φ)*∇(v) ⊙ ( (S∘∇(u))⋅∇(du) ) )*dΩ - jac(u,du,v,φ,dΩ,dΓ_N) = jac_mat(u,du,v,φ,dΩ,dΓ_N) + jac_geo(u,du,v,φ,dΩ,dΓ_N) + res(u,v,φ) = ∫( (I ∘ φ)*((dE∘(∇(v),∇(u))) ⊙ (S∘∇(u))) )*dΩ - ∫(g⋅v)dΓ_N + jac_mat(u,du,v,φ) = ∫( (I ∘ φ)*(dE∘(∇(v),∇(u))) ⊙ (dS∘(∇(du),∇(u))) )*dΩ + jac_geo(u,du,v,φ) = ∫( (I ∘ φ)*∇(v) ⊙ ( (S∘∇(u))⋅∇(du) ) )*dΩ + jac(u,du,v,φ) = jac_mat(u,du,v,φ) + jac_geo(u,du,v,φ) ## Optimisation functionals - J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*((dE∘(∇(u),∇(u))) ⊙ (S∘∇(u))))dΩ - Vol(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; - dVol(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + J(u,φ) = ∫((I ∘ φ)*((dE∘(∇(u),∇(u))) ⊙ (S∘∇(u))))dΩ + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol(q,u,φ) = ∫(-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;jac) + state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh;jac) pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dC=[dVol]) ## Hilbertian extension-regularisation problems diff --git a/test/seq/NonlinearThermalComplianceALMTests.jl b/test/seq/NonlinearThermalComplianceALMTests.jl index c480870..c50f61c 100644 --- a/test/seq/NonlinearThermalComplianceALMTests.jl +++ b/test/seq/NonlinearThermalComplianceALMTests.jl @@ -66,18 +66,18 @@ function main() κ0 = 1 ξ = -1 κ(u) = κ0*(exp(ξ*u)) - res(u,v,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(v))dΩ - ∫(v)dΓ_N + res(u,v,φ) = ∫((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Ω + J(u,φ) = ∫((I ∘ φ)*(κ ∘ u)*∇(u)⋅∇(u))dΩ + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol(q,u,φ) = ∫(-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) + state_map = NonlinearFEStateMap(res,U,V,V_φ,U_reg,φh) pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dC=[dVol]) ## Hilbertian extension-regularisation problems diff --git a/test/seq/PZMultiFieldRepeatingStateTests.jl b/test/seq/PZMultiFieldRepeatingStateTests.jl index 5864d9f..8ed6c07 100644 --- a/test/seq/PZMultiFieldRepeatingStateTests.jl +++ b/test/seq/PZMultiFieldRepeatingStateTests.jl @@ -72,21 +72,21 @@ function main(;AD,use_mfs=false) Eⁱ = (VectorValue(1.0,0.0,), VectorValue(0.0,1.0)) - a((u,ϕ),(v,q),φ,dΩ) = ∫((I ∘ φ) * (1/k0*((C ⊙ ε(u)) ⊙ ε(v)) - + a((u,ϕ),(v,q),φ) = ∫((I ∘ φ) * (1/k0*((C ⊙ ε(u)) ⊙ ε(v)) - 1/α0*((-∇(ϕ) ⋅ e) ⊙ ε(v)) + -1/α0*((e ⋅² ε(u)) ⋅ -∇(q)) + -γ0/β0*((κ ⋅ -∇(ϕ)) ⋅ -∇(q))) )dΩ; - l_ε = [((v,q),φ,dΩ) -> ∫(((I ∘ φ) * (-C ⊙ εᴹ[i] ⊙ ε(v) + k0/α0*(e ⋅² εᴹ[i]) ⋅ -∇(q))))dΩ for i = 1:3]; - l_E = [((v,q),φ,dΩ) -> ∫((I ∘ φ) * ((Eⁱ[i] ⋅ e ⊙ ε(v) + k0/α0*(κ ⋅ Eⁱ[i]) ⋅ -∇(q))))dΩ for i = 1:2]; + l_ε = [((v,q),φ) -> ∫(((I ∘ φ) * (-C ⊙ εᴹ[i] ⊙ ε(v) + k0/α0*(e ⋅² εᴹ[i]) ⋅ -∇(q))))dΩ for i = 1:3]; + l_E = [((v,q),φ) -> ∫((I ∘ φ) * ((Eⁱ[i] ⋅ e ⊙ ε(v) + k0/α0*(κ ⋅ Eⁱ[i]) ⋅ -∇(q))))dΩ for i = 1:2]; l = [l_ε; l_E] - function Cᴴ(r,s,uϕ,φ,dΩ) + function Cᴴ(r,s,uϕ,φ) u_s = uϕ[2s-1]; ϕ_s = uϕ[2s] ∫(1/k0 * (I ∘ φ) * (((C ⊙ (1/k0*ε(u_s) + εᴹ[s])) ⊙ εᴹ[r]) - ((-1/α0*∇(ϕ_s) ⋅ e) ⊙ εᴹ[r])))dΩ; end - function DCᴴ(r,s,q,uϕ,φ,dΩ) + function DCᴴ(r,s,q,uϕ,φ) u_r = uϕ[2r-1]; ϕ_r = uϕ[2r] u_s = uϕ[2s-1]; ϕ_s = uϕ[2s] ∫(- 1/k0 * q * ( @@ -98,19 +98,19 @@ function main(;AD,use_mfs=false) )dΩ; end - Bᴴ(uϕ,φ,dΩ) = 1/4*(Cᴴ(1,1,uϕ,φ,dΩ)+Cᴴ(2,2,uϕ,φ,dΩ)+2*Cᴴ(1,2,uϕ,φ,dΩ)) - DBᴴ(q,uϕ,φ,dΩ) = 1/4*(DCᴴ(1,1,q,uϕ,φ,dΩ)+DCᴴ(2,2,q,uϕ,φ,dΩ)+2*DCᴴ(1,2,q,uϕ,φ,dΩ)) + Bᴴ(uϕ,φ) = 1/4*(Cᴴ(1,1,uϕ,φ)+Cᴴ(2,2,uϕ,φ)+2*Cᴴ(1,2,uϕ,φ)) + DBᴴ(q,uϕ,φ) = 1/4*(DCᴴ(1,1,q,uϕ,φ)+DCᴴ(2,2,q,uϕ,φ)+2*DCᴴ(1,2,q,uϕ,φ)) - J(uϕ,φ,dΩ) = -1*Bᴴ(uϕ,φ,dΩ) - DJ(q,uϕ,φ,dΩ) = -1*DBᴴ(q,uϕ,φ,dΩ) - C1(uϕ,φ,dΩ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; - DC1(q,uϕ,φ,dΩ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ + J(uϕ,φ) = -1*Bᴴ(uϕ,φ) + DJ(q,uϕ,φ) = -1*DBᴴ(q,uϕ,φ) + C1(uϕ,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + DC1(q,uϕ,φ) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ ## Finite difference solver and level set function stencil = HamiltonJacobiEvolution(FirstOrderStencil(2,Float64),model,V_φ,tol,max_steps) ## Setup solver and FE operators - state_map = RepeatingAffineFEStateMap(5,a,l,UP,VQ,V_φ,U_reg,φh,dΩ) + state_map = RepeatingAffineFEStateMap(5,a,l,UP,VQ,V_φ,U_reg,φh) pcfs = if AD PDEConstrainedFunctionals(J,[C1],state_map;analytic_dJ=DJ,analytic_dC=[DC1]) else diff --git a/test/seq/ThermalComplianceALMTests.jl b/test/seq/ThermalComplianceALMTests.jl index fcf3704..33591d4 100644 --- a/test/seq/ThermalComplianceALMTests.jl +++ b/test/seq/ThermalComplianceALMTests.jl @@ -61,20 +61,20 @@ function main(;order,AD) 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 + a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ + l(v,φ) = ∫(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Ω + J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ + dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ; + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol(q,u,φ) = ∫(-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 = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) + state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh) pcfs = if AD PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol]) else @@ -154,20 +154,20 @@ function main_3d(;order) 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 + a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ + l(v,φ) = ∫(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Ω + J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ + dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ; + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol(q,u,φ) = ∫(-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) + state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh) pcfs = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol]) ## Hilbertian extension-regularisation problems diff --git a/test/seq/ThermalComplianceHPMTests.jl b/test/seq/ThermalComplianceHPMTests.jl index 87fce14..2ebe2a5 100644 --- a/test/seq/ThermalComplianceHPMTests.jl +++ b/test/seq/ThermalComplianceHPMTests.jl @@ -61,20 +61,20 @@ function main(;order,AD_case) 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 + a(u,v,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(v))dΩ + l(v,φ) = ∫(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Ω + J(u,φ) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ + dJ(q,u,φ) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ; + Vol(u,φ) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ; + dVol(q,u,φ) = ∫(-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 = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N) + state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh) pcfs = if AD_case == :no_ad PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,analytic_dC=[dVol]) elseif AD_case == :with_ad