Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Documenter and refactor #40

Merged
merged 75 commits into from
Mar 12, 2024
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
75 commits
Select commit Hold shift + click to select a range
4b3f304
VelocityExtension.jl docs & clean
zjwegert Feb 3, 2024
c75ecf5
Utilities.jl refactor & docs
zjwegert Feb 3, 2024
7912016
Merge pull request #39 from zjwegert/main
zjwegert Feb 5, 2024
71ed3b1
Started ChainRules.jl documentation and added some TODOs
zjwegert Feb 6, 2024
0059957
ChainRules.jl docs draft
zjwegert Feb 6, 2024
603b5e8
Advection.jl docs
zjwegert Feb 6, 2024
0c6e80f
Benchmarks.jl docs and small rearrange
zjwegert Feb 6, 2024
ecf1e2a
Remove export TODOs
zjwegert Feb 7, 2024
f5a3beb
Merge branch 'main' into src_refactor_and_docs
zjwegert Feb 7, 2024
39b27c0
Clean up exports and general clean
zjwegert Feb 7, 2024
0cac404
Started optimiser docs + clean up
zjwegert Feb 7, 2024
91782a4
Remove TODO
zjwegert Feb 7, 2024
178e0bf
Added Optimiser.jl docs and more cleaning
zjwegert Feb 7, 2024
e64f9ef
Add note
zjwegert Feb 7, 2024
8215892
Name change!
zjwegert Feb 7, 2024
3fc85e6
Replace all 'LSTO_Distributed' -> 'LevelSetTopOpt'
zjwegert Feb 7, 2024
39ad12c
This shouldn't be here!
zjwegert Feb 7, 2024
4cb1354
Typos
zjwegert Feb 7, 2024
11963f6
Fix type/func headers
zjwegert Feb 7, 2024
9eba725
Started creating docs pages
zjwegert Feb 7, 2024
a46f0e8
Collapse toc
zjwegert Feb 7, 2024
25733fa
bug fix
zjwegert Feb 7, 2024
33ae174
Some more IO docs
zjwegert Feb 7, 2024
b26aac3
chainrules.md docs
zjwegert Feb 8, 2024
18f7962
Additional reference docs
zjwegert Feb 8, 2024
3534546
Update readme, added developer notes
zjwegert Feb 8, 2024
eb96588
Added skeleton for tutorials
zjwegert Feb 8, 2024
c346df6
Added more docs
zjwegert Feb 18, 2024
f2f6a7c
adjust tut 1
zjwegert Feb 18, 2024
08c1fa0
IMPORTANT: Added Gamma_D dirichlet tag to match docs
zjwegert Feb 18, 2024
4aae2ee
Docs update & remove Zygote
zjwegert Feb 19, 2024
2ed9d42
Additional docs in tut 1
zjwegert Feb 19, 2024
d40ec83
Working on tutorial, fixed project wide typo
zjwegert Feb 19, 2024
1c02c76
Still writing tut 1
zjwegert Feb 19, 2024
bc7db53
Refactoring and renaming + verified tut 1 core
zjwegert Feb 19, 2024
99f77c5
POSSIBLY BREAKING: Need to check all problems
zjwegert Feb 19, 2024
9729d46
Undid prev change, adjust ALM conv crit
zjwegert Feb 19, 2024
5cccfea
Finished a draft of tut1
zjwegert Feb 20, 2024
873242e
IMPORTANT: Removed `dF .*= -1` in ChainRules to correct TODO. This re…
zjwegert Feb 20, 2024
9453d6c
Spelling
zjwegert Feb 20, 2024
2722df2
replace with num
zjwegert Feb 21, 2024
51e5a68
Added more options to HPM, added inverter ALM for paper
zjwegert Feb 21, 2024
4ca6d11
typo fix
zjwegert Feb 22, 2024
a392481
QoL - added initial parameters function
zjwegert Feb 23, 2024
0d3f7df
Added 2d thermal MPI + PETSc for paper
zjwegert Feb 23, 2024
5d4bc2e
added for paper
zjwegert Feb 23, 2024
f650d1c
Updated 3d inverter scripts
zjwegert Feb 25, 2024
c4034ef
Note for self
zjwegert Feb 27, 2024
5016ece
fixed minor typos in ChainRules.jl docs
zjwegert Feb 28, 2024
1f7aeb8
several corrections
zjwegert Feb 29, 2024
7547a31
Reuse serial `get_el_Δ`
zjwegert Feb 29, 2024
b1aa365
Merge branch 'src_refactor_and_docs' of https://github.com/zjwegert/L…
zjwegert Feb 29, 2024
0ecdb7a
Resolve #35
zjwegert Mar 4, 2024
64d12c9
Fix
zjwegert Mar 4, 2024
db9682e
Fix Neumann BC setup
zjwegert Mar 4, 2024
0e95b49
Removed unnecessary `prod` call
zjwegert Mar 5, 2024
7a9bfa5
Found typo in advection, need to re-run sims arg
zjwegert Mar 5, 2024
d1e6d80
Fix
zjwegert Mar 5, 2024
6845663
Testing different sign approximation
zjwegert Mar 6, 2024
dc274d0
CHANGE TO ADVECTION - may need to undo
zjwegert Mar 7, 2024
471a2ee
Ready for 3D testing
zjwegert Mar 7, 2024
041792d
Added `reinit_mod` to control how often we reinitialise
zjwegert Mar 7, 2024
7626979
file naming for testing
zjwegert Mar 7, 2024
363e99a
Important: Changed storage in `SmoothErsatzMaterialInterpolation`, al…
zjwegert Mar 8, 2024
f2e93e6
fix 3d job
zjwegert Mar 11, 2024
054c3a1
bug fix
zjwegert Mar 11, 2024
1b8fb1a
minor adjustments to paper scripts
zjwegert Mar 11, 2024
592cf65
rename
zjwegert Mar 11, 2024
1628970
Added stuff from ODE refactor branch (temp code)
JordiManyer Mar 11, 2024
f97976e
Added MultiField support for RepeatingAffineFEStateMaps
JordiManyer Mar 12, 2024
7fddabe
Merge pull request #41 from zjwegert/jordi-dev
JordiManyer Mar 12, 2024
bcfec72
Added CI structure
JordiManyer Mar 12, 2024
999e86c
Merge pull request #42 from zjwegert/jordi-dev
JordiManyer Mar 12, 2024
aa3ca0d
Commented doc deployment
JordiManyer Mar 12, 2024
192d8ee
import estimate_period
zjwegert Mar 12, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
added for paper
  • Loading branch information
zjwegert committed Feb 23, 2024
commit 5d4bc2eb1263a2ea29ad827e42afd69df3559a9f
Original file line number Diff line number Diff line change
@@ -0,0 +1,193 @@
using Gridap, LevelSetTopOpt

"""
(Serial) Minimum hyperelastic (neohookean) compliance with augmented Lagrangian method in 2D.
"""
function main()
## Parameters
order = 1
xmax,ymax=2.0,1.0
prop_Γ_N = 0.4
dom = (0,xmax,0,ymax)
el_size = (400,200)
γ = 0.1
γ_reinit = 0.5
max_steps = floor(Int,minimum(el_size)/10)
tol = 1/(10order^2)*prod(inv,minimum(el_size))
η_coeff = 2
α_coeff = 4
vf = 0.4
path = dirname(dirname(@__DIR__))*"/results/highres_hyperelastic_compliance_neohook_ALM"
mkdir(path)

## FE Setup
model = CartesianDiscreteModel(dom,el_size);
el_Δ = get_el_Δ(model)
f_Γ_D(x) = (x[1] ≈ 0.0)
f_Γ_N(x) = (x[1] ≈ xmax && ymax/2-ymax*prop_Γ_N/4 - eps() <= x[2] <= ymax/2+ymax*prop_Γ_N/4 + 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 = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
reffe_scalar = ReferenceFE(lagrangian,Float64,order)
V = TestFESpace(model,reffe;dirichlet_tags=["Gamma_D"])
U = TrialFESpace(V,VectorValue(0.0,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.ρ

## Material properties
_E = 1000
ν = 0.3
μ, λ = _E/(2*(1 + ν)), _E*ν/((1 + ν)*(1 - 2*ν))
g = VectorValue(0,-10)

## Neohookean hyperelastic material
# Deformation gradient
F(∇u) = one(∇u) + ∇u'
J(F) = sqrt(det(C(F)))

# Derivative of green Strain
dE(∇du,∇u) = 0.5*( ∇du⋅F(∇u) + (∇du⋅F(∇u))' )

# Right Caughy-green deformation tensor
C(F) = (F')⋅F

# Constitutive law (Neo hookean)
function S(∇u)
Cinv = inv(C(F(∇u)))
μ*(one(∇u)-Cinv) + λ*log(J(F(∇u)))*Cinv
end

# 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

## 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Ω

## Finite difference solver and level set function
stencil = AdvectionStencil(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(Obj,[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,stencil,vel_ext,φh;
γ,γ_reinit,verbose=true,constraint_names=[:Vol])
for (it,uh,φh) in optimiser
write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
write_history(path*"/history.txt",optimiser.history)
end
it = get_history(optimiser).niter; uh = get_state(pcfs)
write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1)
end

main()

function main_LE_comparison()
## Parameters
order = 1
xmax,ymax=(2.0,1.0)
prop_Γ_N = 0.4
dom = (0,xmax,0,ymax)
el_size = (400,200)
γ = 0.1
γ_reinit = 0.5
max_steps = floor(Int,minimum(el_size)/10)
tol = 1/(10order^2)*prod(inv,minimum(el_size))
C = isotropic_elast_tensor(2,1000,0.3)
η_coeff = 2
α_coeff = 4
vf = 0.4
g = VectorValue(0,-10)
path = dirname(dirname(@__DIR__))*"/results/highres_LE_comparison"
mkdir(path)

## FE Setup
model = CartesianDiscreteModel(dom,el_size)
el_Δ = get_el_Δ(model)
f_Γ_D(x) = (x[1] ≈ 0.0)
f_Γ_N(x) = (x[1] ≈ xmax && ymax/2-ymax*prop_Γ_N/4 - eps() <= x[2] <= ymax/2+ymax*prop_Γ_N/4 + 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 = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
reffe_scalar = ReferenceFE(lagrangian,Float64,order)
V = TestFESpace(model,reffe;dirichlet_tags=["Gamma_D"])
U = TrialFESpace(V,VectorValue(0.0,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 ∘ φ)*(C ⊙ ε(u) ⊙ ε(v)))dΩ
l(v,φ,dΩ,dΓ_N) = ∫(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Ω

## Finite difference solver and level set function
stencil = AdvectionStencil(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 = PDEConstrainedFunctionals(J,[Vol],state_map,analytic_dJ=dJ,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,stencil,vel_ext,φh;
γ,γ_reinit,verbose=true,constraint_names=[:Vol])
for (it,uh,φh) in optimiser
write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
write_history(path*"/history.txt",optimiser.history)
end
it = get_history(optimiser).niter; uh = get_state(pcfs)
write_vtk(Ω,path*"/struc_$it",it,["phi"=>φh,"H(phi)"=>(H ∘ φh),"|nabla(phi)|"=>(norm ∘ ∇(φh)),"uh"=>uh];iter_mod=1)
end

main_LE_comparison()
104 changes: 104 additions & 0 deletions scripts/_dev/Added_for_paper/thermal_MPI_PETSc.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,104 @@
using LevelSetTopOpt, Gridap, GridapDistributed, GridapPETSc, PartitionedArrays,
SparseMatricesCSR

function main(mesh_partition,distribute)
ranks = distribute(LinearIndices((prod(mesh_partition),)))
# FE parameters
order = 1 # Finite element order
xmax = ymax = 1.0 # Domain size
dom = (0,xmax,0,ymax) # Bounding domain
el_size = (400,400) # Mesh partition size
prop_Γ_N = 0.4 # Γ_N size parameter
prop_Γ_D = 0.2 # Γ_D size parameter
f_Γ_N(x) = (x[1] ≈ xmax && # Γ_N indicator function
ymax/2-ymax*prop_Γ_N/4 - eps() <= x[2] <= ymax/2+ymax*prop_Γ_N/4 + eps())
f_Γ_D(x) = (x[1] ≈ 0.0 && # Γ_D indicator function
(x[2] <= ymax*prop_Γ_D + eps() || x[2] >= ymax-ymax*prop_Γ_D - eps())) # @$\lvert\lvert$@
# FD parameters
γ = 0.1 # HJ eqn time step coeff
γ_reinit = 0.5 # Reinit. eqn time step coeff
max_steps = floor(Int,minimum(el_size)/10) # Max steps for advection
tol = 1/(10order^2)*prod(inv,minimum(el_size)) # Advection tolerance
# Problem parameters
κ = 1 # Diffusivity
g = 1 # Heat flow in
vf = 0.4 # Volume fraction constraint
lsf_func = initial_lsf(4,0.2) # Initial level set function
iter_mod = 10 # VTK Output modulo
path = "./results/tut1_MPI_PETSc_rtol12/" # Output path
i_am_main(ranks) && mkpath(path) # Create path
# Model
model = CartesianDiscreteModel(ranks,mesh_partition,dom,el_size);
update_labels!(1,model,f_Γ_D,"Gamma_D")
update_labels!(2,model,f_Γ_N,"Gamma_N")
# Triangulation and measures
Ω = Triangulation(model)
Γ_N = BoundaryTriangulation(model,tags="Gamma_N")
dΩ = Measure(Ω,2*order)
dΓ_N = Measure(Γ_N,2*order)
# Spaces
reffe = ReferenceFE(lagrangian,Float64,order)
V = TestFESpace(model,reffe;dirichlet_tags=["Gamma_D"])
U = TrialFESpace(V,0.0)
V_φ = TestFESpace(model,reffe)
V_reg = TestFESpace(model,reffe;dirichlet_tags=["Gamma_N"])
U_reg = TrialFESpace(V_reg,0)
# Level set and interpolator
φh = interpolate(lsf_func,V_φ)
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
Tm = SparseMatrixCSR{0,PetscScalar,PetscInt}
Tv = Vector{PetscScalar}
solver = PETScLinearSolver()
state_map = AffineFEStateMap(
a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N;
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Ω
vol_D = sum(∫(1)dΩ)
C(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
dC(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
pcfs = PDEConstrainedFunctionals(J,[C],state_map,
analytic_dJ=dJ,analytic_dC=[dC])
# Velocity extension
α = 4*maximum(get_el_Δ(model))
a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ
vel_ext = VelocityExtension(
a_hilb, U_reg, V_reg;
assem = SparseMatrixAssembler(Tm,Tv,U_reg,V_reg),
ls = solver
)
# Finite difference scheme
scheme = FirstOrderStencil(length(el_size),Float64)
stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps)
# Optimiser
optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
γ,γ_reinit,verbose=i_am_main(ranks),constraint_names=[:Vol])
# Solve
for (it,uh,φh) in optimiser
data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|nabla(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
iszero(it % iter_mod) && (writevtk(Ω,path*"out$it",cellfields=data);GC.gc())
write_history(path*"/history.txt",get_history(optimiser);ranks)
end
# Final structure
it = get_history(optimiser).niter; uh = get_state(pcfs)
writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,
"H(φ)"=>(H ∘ φh),"|nabla(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
end

with_mpi() do distribute
mesh_partition = (2,2)
solver_options = "-pc_type gamg -ksp_type cg -ksp_error_if_not_converged true
-ksp_converged_reason -ksp_rtol 1.0e-12"
GridapPETSc.with(args=split(solver_options)) do
main(mesh_partition,distribute)
end
end
82 changes: 82 additions & 0 deletions scripts/_dev/Added_for_paper/thermal_serial.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,82 @@
using LevelSetTopOpt, Gridap

function main()
# FE parameters
order = 1 # Finite element order
xmax = ymax = 1.0 # Domain size
dom = (0,xmax,0,ymax) # Bounding domain
el_size = (400,400) # Mesh partition size
prop_Γ_N = 0.4 # Γ_N size parameter
prop_Γ_D = 0.2 # Γ_D size parameter
f_Γ_N(x) = (x[1] ≈ xmax && # Γ_N indicator function
ymax/2-ymax*prop_Γ_N/4 - eps() <= x[2] <= ymax/2+ymax*prop_Γ_N/4 + eps())
f_Γ_D(x) = (x[1] ≈ 0.0 && # Γ_D indicator function
(x[2] <= ymax*prop_Γ_D + eps() || x[2] >= ymax-ymax*prop_Γ_D - eps())) # @$\lvert\lvert$@
# FD parameters
γ = 0.1 # HJ eqn time step coeff
γ_reinit = 0.5 # Reinit. eqn time step coeff
max_steps = floor(Int,minimum(el_size)/10) # Max steps for advection
tol = 1/(10order^2)*prod(inv,minimum(el_size)) # Advection tolerance
# Problem parameters
κ = 1 # Diffusivity
g = 1 # Heat flow in
vf = 0.4 # Volume fraction constraint
lsf_func = initial_lsf(4,0.2) # Initial level set function
iter_mod = 10 # VTK Output modulo
path = "./results/tut1/" # Output path
mkpath(path) # Create path
# Model
model = CartesianDiscreteModel(dom,el_size);
update_labels!(1,model,f_Γ_D,"Gamma_D")
update_labels!(2,model,f_Γ_N,"Gamma_N")
# Triangulation and measures
Ω = Triangulation(model)
Γ_N = BoundaryTriangulation(model,tags="Gamma_N")
dΩ = Measure(Ω,2*order)
dΓ_N = Measure(Γ_N,2*order)
# Spaces
reffe = ReferenceFE(lagrangian,Float64,order)
V = TestFESpace(model,reffe;dirichlet_tags=["Gamma_D"])
U = TrialFESpace(V,0.0)
V_φ = TestFESpace(model,reffe)
V_reg = TestFESpace(model,reffe;dirichlet_tags=["Gamma_N"])
U_reg = TrialFESpace(V_reg,0)
# Level set and interpolator
φh = interpolate(lsf_func,V_φ)
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
state_map = AffineFEStateMap(a,l,U,V,V_φ,U_reg,φh,dΩ,dΓ_N)
# Objective and constraints
J(u,φ,dΩ,dΓ_N) = ∫((I ∘ φ)*κ*∇(u)⋅∇(u))dΩ
dJ(q,u,φ,dΩ,dΓ_N) = ∫(κ*∇(u)⋅∇(u)*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
vol_D = sum(∫(1)dΩ)
C(u,φ,dΩ,dΓ_N) = ∫(((ρ ∘ φ) - vf)/vol_D)dΩ
dC(q,u,φ,dΩ,dΓ_N) = ∫(-1/vol_D*q*(DH ∘ φ)*(norm ∘ ∇(φ)))dΩ
pcfs = PDEConstrainedFunctionals(J,[C],state_map,
analytic_dJ=dJ,analytic_dC=[dC])
# Velocity extension
α = 4*maximum(get_el_Δ(model))
a_hilb(p,q) =∫(α^2*∇(p)⋅∇(q) + p*q)dΩ
vel_ext = VelocityExtension(a_hilb,U_reg,V_reg)
# Finite difference scheme
scheme = FirstOrderStencil(length(el_size),Float64)
stencil = AdvectionStencil(scheme,model,V_φ,tol,max_steps)
# Optimiser
optimiser = AugmentedLagrangian(pcfs,stencil,vel_ext,φh;
γ,γ_reinit,verbose=true,constraint_names=[:Vol])
# Solve
for (it,uh,φh) in optimiser
data = ["φ"=>φh,"H(φ)"=>(H ∘ φh),"|nabla(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh]
iszero(it % iter_mod) && (writevtk(Ω,path*"out$it",cellfields=data);GC.gc())
write_history(path*"/history.txt",get_history(optimiser))
end
# Final structure
it = get_history(optimiser).niter; uh = get_state(pcfs)
writevtk(Ω,path*"out$it",cellfields=["φ"=>φh,
"H(φ)"=>(H ∘ φh),"|nabla(φ)|"=>(norm ∘ ∇(φh)),"uh"=>uh])
end

main()