I want to obtain the Reaction Force of the plate (10x10) subjected to a load at the right edge. For this, I wrote the following code to compute the matrix stiffness A and then the reaction force. However, it does not seem that the method is correct as the stiffness and the reaction force are not accurate comparing to the result of Abaqus. Is there any other way to compute the stiffness matrix and the reaction force?
Moreover how can I plot the contour of the vector Force? I tried Force[0] but I have an error "don't know how to plot Force[0]"
Thanks
Variational form
a = inner(eps(v), Ceps(u))dx
l = inner(t,v[0])*ds(0)
stiffness A
A= assemble(a)
Solving linear variational problem
u= Function(U)
problem = LinearVariationalProblem(a, l, , bcs=bcs)
solver = LinearVariationalSolver(problem)
solver.parameters["linear_solver"] = "direct"
solver.solve()
Reaction force
Force = A * u.vector().array()
reaction_force = 0.0
coor = mesh.coordinates()
for i in range(mesh.num_vertices()):
if coor[i][0] == 10:
reaction_force = reaction_force + Force[i]
print (reaction_force)
