Cantilever Beam Test

General
1

Hi,
i tried to implement a cantilever beam test. So my problem would be the geometry which is fixed on the upper and lower left and a traction load downwards on the middle right. We then solve the unconstrained shape optimization problem

\underset{\Omega}{\mbox{min}} \;J(\Omega) + l Vol(\Omega) = \int_\Omega \sigma(u_\Omega):\epsilon(u_\Omega) \; dx + l \int_\Omega 1 \:dx

where u_\Omega is the solution to the linear elasticity pde. I then calculated the state and adjoint

def SolveStateEquation():
    u, v = fes.TnT()
    a = BilinearForm(fes)
    a += InnerProduct(sigma(u), eps(v)) * dx
    a.Assemble()
    ainv = a.mat.Inverse(fes.FreeDofs())
    
    g = ngs.CF((0,-1))
    f = LinearForm(fes)
    f+= g * v * ds("neumann")
    f.Assemble()
    
    gfu = GridFunction(fes)
    gfu.vec.data = ainv * f.vec

    p, w = fes.TnT()

    rhs_adjoint = LinearForm(fes)
    rhs_adjoint += (-1) * Cost(gfu, gfu).Diff(gfu, w)
    rhs_adjoint.Assemble()

    gfp = GridFunction(fes)
    gfp.vec.data = ainv * rhs_adjoint.vec
    
    return gfu, gfp

and shape derivative:

def SolveDeformationEquation():
    g = ngs.CF((0,-1))
    Lagrangian = Cost(gfu, gfu) + InnerProduct(sigma(gfu), eps(gfp))*dx - g*gfp*ds("neumann")
    
    VEC = H1(mesh, order=2, dim=2, dirichlet="dirichlet")
    (PHI, X) = VEC.TnT()
    
    dJ_Omega = LinearForm(VEC)
    dJ_Omega += Lagrangian.DiffShape(X)
    
    # Shape Diff glätten
    my_mesh_size = 0.3
    b = BilinearForm(VEC)
    reg_param = 3*my_mesh_size**2
    b += reg_param * InnerProduct(grad(X), grad(PHI)) * dx + InnerProduct(X, PHI) * dx
    # add Cauchy-Riemann equation for better angle-preserving deformations
    alpha = 1
    b += alpha * (PHI.Deriv()[0,0] - PHI.Deriv()[1,1])*(X.Deriv()[0,0] - X.Deriv()[1,1]) *dx
    b += alpha * (PHI.Deriv()[1,0] + PHI.Deriv()[0,1])*(X.Deriv()[1,0] + X.Deriv()[0,1]) *dx
    
    gfX = GridFunction(VEC)
    b.Assemble()
    dJ_Omega.Assemble()
    rhs = gfX.vec.CreateVector()
    rhs.data = dJ_Omega.vec
    gfX.vec.data = b.mat.Inverse(VEC.FreeDofs()) * rhs

    return gfX

Then i ran my algorithm which infers a descent direction, with and without line search, but the cost function is basically always increasing but i dont get why. I put the jupyter notebook also in here, could someone help or try to figure out with me what goes wrong?
Thank you in advance
Cantilever Beam Test.ipynb (18.5 KB)