Solving topology optimization with external top-center force not working

variational formulation
1

Sorry in advance for my poor knowledge in the matter.

The topology I’m seeking as shown in the figure:

image
image844×230 2.41 KB

I’m following the guide written here for cantilever and did the tweaks but whenever I move the force vector(the load) inside the object or at a point on its boundary from the top, the solver rests at 0. Can someone tell me why and guide me through this?
Here is the code:

# Problem parameters
thetamoy = 0.4 # target average material density
E = Constant(1)
nu = Constant(0.3)
lamda = E*nu/(1+nu)/(1-2*nu)
mu = E/(2*(1+nu))
f = Constant((0, -1)) # vertical downwards force
lx,ly = [4.0,1.0]    
Lag,Nx,Ny = [130.0,150,50]
# Mesh
mesh  = RectangleMesh(Point(0.0,0.0),Point(lx,ly),Nx,Ny,'crossed') 
# Boundaries
tol = 1E-14
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0) and abs(x[1]) > tol
class Right(SubDomain):
    def inside(self, x, on_boundary):
        return abs(x[0]-lx) < tol and abs(x[1]) < tol
class Load(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 2.0) and near(x[1], 1.05, 0.05)

left, right = [Left(), Right()]
load = Load()
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)

load.mark(boundaries, 1)

ds = Measure("ds", subdomain_data=boundaries)
# Function space for density field
V0 = FunctionSpace(mesh, "DG", 0)
# Function space for displacement
V2 = VectorFunctionSpace(mesh, "CG", 2)
# Fixed boundary condtions
bc = [DirichletBC(V2, Constant((0, 0)), left, method='pointwise'), DirichletBC(V2, Constant((0, 0)), right, method='pointwise')]

p = Constant(1) # SIMP penalty exponent

exponent_counter = 0 # exponent update counter
lagrange = Constant(1) # Lagrange multiplier for volume constraint

thetaold = Function(V0, name="Density")
thetaold.interpolate(Constant(thetamoy))
coeff = thetaold**p
theta = Function(V0)

volume = assemble(Constant(1.)*dx(domain=mesh))
avg_density_0 = assemble(thetaold*dx)/volume # initial average density
avg_density = 0.

def eps(v):
    return sym(grad(v))
def sigma(v):
    return coeff*(lamda*div(v)*Identity(2)+2*mu*eps(v))
def energy_density(u, v):
    return inner(sigma(u), eps(v))

# Inhomogeneous elastic variational problem
u_ = TestFunction(V2)
du = TrialFunction(V2)
a = inner(sigma(u_), eps(du))*dx
L = dot(f, u_)*ds(1)

Thanks in advance.