Skip to content

Commit

Permalink
BREAKING: Removed IntegrandWithMeasure type
Browse files Browse the repository at this point in the history
  • Loading branch information
@zjwegert
zjwegert committed Nov 13, 2024
1 parent 88440c8 commit 496671e
Show file tree
Hide file tree
Showing 33 changed files with 302 additions and 414 deletions.
52 changes: 26 additions & 26 deletions scripts/Benchmarks/benchmark.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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),
Expand Down
14 changes: 7 additions & 7 deletions scripts/MPI/elastic_compliance_ALM.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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),
Expand Down
18 changes: 9 additions & 9 deletions scripts/MPI/inverse_homenisation_ALM.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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),
Expand Down
10 changes: 5 additions & 5 deletions scripts/MPI/nonlinear_thermal_compliance_ALM.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
14 changes: 7 additions & 7 deletions scripts/MPI/thermal_compliance_ALM.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
14 changes: 7 additions & 7 deletions scripts/MPI/thermal_compliance_ALM_PETSc.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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),
Expand Down
14 changes: 7 additions & 7 deletions scripts/Serial/elastic_compliance_ALM.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
14 changes: 7 additions & 7 deletions scripts/Serial/elastic_compliance_HPM.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
10 changes: 5 additions & 5 deletions scripts/Serial/elastic_compliance_LM.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
10 changes: 5 additions & 5 deletions scripts/Serial/hyperelastic_compliance_ALM.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
10 changes: 5 additions & 5 deletions scripts/Serial/hyperelastic_compliance_neohook_ALM.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading

0 comments on commit 496671e

Please sign in to comment.