The following code represents a linear elasticity problem with a uniform traction force applied in the normal direction on boundary 1. How should I modify this code if I want to set the traction force as a function of the z-coordinate (e.g., T = z + 1)?
T = Constant(mesh, 10)
n = FacetNormal(mesh)
dU = TrialFunction(VU)
U_ = TestFunction(VU)
ds = Measure('ds', domain=mesh, subdomain_data=facet_markers)
Wint = inner(sigma(dU, Delta_T), eps(U_)) * dx
aM = lhs(Wint)
LM = rhs(Wint) + dot(T * n, U_) * ds(1)