Problem about extract nodal force

1

hi, I want to extract nodal force and i write the following code but do not know is it correct for extracting nodal force

L = 1.; H = 1.
Nx = 10; Ny = 10
mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny)
nodal_coordinates = mesh.coordinates()

Vm = FunctionSpace(mesh,'DG',1)
W = VectorFunctionSpace(mesh,'P',1)
F = TensorFunctionSpace(mesh, "Lagrange", 2)
u , v = TrialFunction(W), TestFunction(W)

E = 10e3
nu = 0.3
# plane problem type
model = "plane_stress"
# lame's parameters for plane strain
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

# lame's parameters for plane stress
if model == "plane_stress":
    lmbda = 2*mu*lmbda/(lmbda+2*mu)

bottom = CompiledSubDomain("near(x[1], 0) && on_boundary")
left = CompiledSubDomain("near(x[0], 0) && on_boundary")

bc_u = [DirichletBC(W.sub(1), Constant(0.0), bottom),
        DirichletBC(W.sub(0), Constant(0.0), left),]

# traction boundary condition
top = CompiledSubDomain("near(x[1], 1) && on_boundary")
load = Constant((0, 1e-3))

# index boundary no.
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
top.mark(boundaries,1)
ds = Measure("ds")(subdomain_data=boundaries)
n = FacetNormal(mesh)

def epsilon(u):
    return 0.5*(nabla_grad(u)+nabla_grad(u).T)
def sigma(u):
    return 2.0*mu*epsilon(u)+lmbda*tr(epsilon(u))*Identity(len(u))

LHS = inner(sigma(u), epsilon(v))*dx
RHS = dot(load, v)*ds(1)
u = Function(W, name="Displacement")

K = assemble(LHS)
f = assemble(RHS)
[bc.apply(K, f) for bc in bc_u]
solve(K,u.vector(),f)

v2d = vertex_to_dof_map(W)
nodal_solution__matrix = np.zeros((mesh.num_vertices(),7))

nodal_force = K*u.vector()
nodal_force = np.array(nodal_force)
nodal_force = np.array(nodal_force[v2d])
for node_j in range(mesh.num_vertices()):
    nodal_j_x_force = nodal_force[2 * node_j]
    nodal_j_y_force = nodal_force[2 * node_j + 1]

    nodal_solution__matrix[node_j, 5] = nodal_j_x_force
    nodal_solution__matrix[node_j, 6] = nodal_j_y_force

thanks