Computing the self-weight of a truss

Hello community,
I have a truss like structure. I am trying to apply the self-weight (-g * rho) to the back cells only. My current implementation only applies force on the hollow areas (white cells) of the section (shown in Figure 1).
Could someone help me with properly applying the self-weight given the density?

image

Fig. 1 The black areas is the actual truss (where the material is and forces should be applied here) and the yellow areas are the hollow part of the truss where the forces are being applied (where there is actually no material)

Here is a working example:

from dolfin import *
from fenics_adjoint import *
import torch_fenics
from ufl import nabla_div


def left(x, on_boundary):
    return near(x[0], -2) and on_boundary

def eps(v):
    return sym(grad(v))

def sigma(v):
    return coeff*(lamda*div(v)*Identity(2)+2*mu*eps(v))


# Problem parameters
p = Constant(5) 
E = Constant(203 * 1e8) # kg/m2
nu = Constant(0.265)

g = Constant(9.8) # m/s2
rho = Constant(7985) # kg/m3

lamda = E*nu/(1+nu)/(1-2*nu)
mu = E/(2*(1+nu))

# Mesh
mesh = RectangleMesh(Point(-2, 0), Point(2, 1), 20, 10, "crossed")
facets = MeshFunction("size_t", mesh, 1)


V0 = FunctionSpace(mesh, "DG", 0)        # Function space for density field
V2 = VectorFunctionSpace(mesh, "CG", 2)  # Function space for displacement


# Assigning the density
# Peta is a numpy array containing the values of the truss; here i used a random value to have a working example
peta = np.ones((800)) * 0.5
theta = Function(V0, name = 'Density')
theta.vector().set_local(peta)

# Inhomogeneous elastic variational problem
u_ = TestFunction(V2)
du = TrialFunction(V2)

u = Function(V2, name="Displacement")
coeff = theta ** p

f = theta * Constant((0, -rho*g))

a = inner(sigma(u_), eps(du))*dx
L = inner(f, u_)* dx

# Applying clamped end
bc = DirichletBC(V2, Constant((0, 0)), left)

# solve
solve(a == L, u, bc)

# Plotting the displacement
plot(u/1e9, mode="displacement")

Feel free to message me so I can also send you the mesh file if needed.

Thank you for your help.

I was curious if anyone had a chance to take a look at this problem, I’d really appreciate any pointers to see how can I deal with this situation.

P.S. I am using an exponent p in the variational formulation (see variable coeff). My initial guess is that as I increase the value for the exponent p, both the elastic modulus and the effect of the density of my material increase. Any ideas on this sentence as well?

Thank you for your time in advance.