Clamped Beam under the End Load Problem

Question

Hello
I have solved a linear elasticity problem : A clamped beam under the end load. The left side is clamped and I want to apply a vertical end load on the right edge. A part of my code:

# Create mesh and define function space
mesh = RectangleMesh(Point(0, 0), Point(L, W), 10, 3)
V = VectorFunctionSpace(mesh, 'P', 1)

class LeftEdge(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1e-6
        return on_boundary and abs(x[0]) < tol

class RightEdge(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1e-6
       return on_boundary and abs(x[0] - L) < tol

left_edge = LeftEdge()
right_edge = RightEdge()
sub_domains = FacetFunction("size_t", mesh, mesh.topology().dim() - 1)
sub_domains.set_all(0)
right_edge.mark(sub_domains, 1)
left_edge.mark(sub_domains, 2)
ds = Measure("ds")[sub_domains]

#Zero Dirichlet B.C on the left side
bc = DirichletBC(V, Constant((0, 0)), left_edge)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

#End Load
T = Constant((0, -10))

#Solve
a = inner(sigma(u), eps(v))*dx
L =dot(T, v)*ds(1)

# Compute solution
u = Function(V)
solve(a == L, u, bc)

The question is, when I plot the end load:

p=interpolate(T,V)
plot(p, title='Load')                          

This is what comes up:

As you can see this is not a load applied on only one edge. I mean this is like a distributed load started from the one end to the other end. In addition when I import my results to the Paraview I get the below pattern for the end load:

Again, as you can see it looks like a distributed load (some arrows started from the beginning to the end of the beam), while I just expect one arrow on the right edge.
Could you please tell me why it looks like this?

Answer 1

You've interpolated the constant T over the space V, which results in a Function defined over all of V, which you are then plotting over the whole mesh. Your solution should be ok (at least wrt. how the load is applied) since you are integrating this boundary condition properly with ds(1).

If you want your p to truly be a point load, you could:

1. Create a zero constant over V, then edit its value for the dofs at the vertex you are interested in (see e.g. this question)

2. (Better) Subclass Expression:

   class PointLoad(Expression):
       def __init__(self, **kwargs):
           self.point = kwargs['point']
           self.value = kwargs['value']
       def eval(self, value, x):
           if near (x[0], self.point[0]) and near(x[1], self.point[1]):
               value[0] = self.value[0]
               value[1] = self.value[1]
           else:  # Reset the input variable because it's likely to be reused
               values[0] = 0
               values[1] = 0
       def value_shape(self):
           return (2,)
   

   Then use this with

    T = PointLoad(point=(1.0,.2), value=(0,-10), degree=0)
   

3. Or use dolfin.PointSource(). Note that you can apply it to one component of a vector:

   u = Function(V)
   ps = PointSource(V.sub(1), L)
   ps.apply(u.vector()) 
   

Comments

But of course! (slap myself in the forehead) We are overwriting value with self.value, when we should be assigning to its components! As it was the code was just discarding the input variable. I'm fixing my answer. I should've run this before posting, sorry

---

Plot T with

plot(T, mesh=V.mesh())

---

Ok, that one was subtle: for efficiency reasons the parameter values is owned by the caller to eval(), so that one needs not allocate memory on each call. But then it is most likely going to be reused! After PointLoad.eval() sets it to some value, it will keep it, so we need to reset it:

if near...:
    # stuff
else:
    values[0] = 0
    values[1] = 0

This is pretty wasteful, but I see no other workaround.

---

Thank you so much for your help. It really works this time! It was not as easy as I thought it would be!

---

Glad it helped! I did learn stuff too :)