How to apply a force in a point of a 2D mesh

variational formulation
1

hello, I am trying to solve a 2D variational problem in dolfin legacy with a mesh that has stiffening buster on top (treated as a timoshenko beam). At this point I should apply a force at the top right point, but the integral of the vertical stresses of the constrained side does not return the total force applied.
The energy formulation works; in fact, the same displacement-controlled model returns the expected results.
I have also tried the Pointwise function, but got no results, thanks for the help in advance

n_mesh=10
Ea=1000
nia=.3
Eb=70000
nib=.3

from dolfin import*
mesh=UnitSquareMesh(n_mesh,n_mesh)
U_ = VectorElement("CG", mesh.ufl_cell(), 2, dim=2)
T_ = FiniteElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, U_*T_)
U=FunctionSpace(mesh,U_)
D=FunctionSpace(mesh,"CG",1)
Mat=TensorFunctionSpace(mesh,"CG",1,(2,2))

v_Test = TestFunction(V)
v_Trial = TrialFunction(V)
v=Function(V)

x3=Expression("x[1]",degree=1)
def boundary_left(x, on_boundary): return near(x[0],0)and on_boundary
def boundary_right(x, on_boundary): return near(x[0],1) and on_boundary
def boundary_bottom(x, on_boundary): return near(x[1],0) and on_boundary
def boundary_top(x, on_boundary): return near(x[1],1) and on_boundary
def boundary_point(x): return near(x[1],1) and near(x[0],1) 

boundary_parts = MeshFunction('size_t', mesh, mesh.topology().dim()-1,0)
ds = ds(domain=mesh, subdomain_data=boundary_parts)
class Right(SubDomain):
    def inside(self, x, on_boundary):
        return boundary_right(x, on_boundary)
Right().mark(boundary_parts, 1)
class Top(SubDomain):
    def inside(self, x, on_boundary):
        return boundary_top(x, on_boundary)
Top().mark(boundary_parts, 2)
class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return boundary_bottom(x, on_boundary)
Bottom().mark(boundary_parts, 3)

class Boundary_point(SubDomain):
    def inside(self, x,on_boundary):
        return boundary_point(x)
Boundary_point().mark(boundary_parts, 4)

mua=.5*Ea/(1+nia)
mub=.5*Eb/(1+nib)
betaa= nia*Ea/((1+nia)*(1-(2*nia)))


def epsilon(u): return sym(grad(u))
def sigma(u): return 2*mua*epsilon(u)+betaa*tr(epsilon(u))*Identity(2)
def elastic_energy(u): return 0.5*inner(sigma(u),epsilon(u))

def line_energy(v):
  u,theta= split(v)
  axial_energy=.5*Eb*A_b*(grad(u[0])[0]*grad(u[0])[0])
  combined_energy=-.5*Eb*A_b*H*(grad(u[0])[0]*grad(theta)[0])
  momentum_enery=.5*Eb*J_base*grad(theta)[0]*grad(theta)[0]
  shear_energy=.5*mub*kappa_cal*A_b*((grad(u[1])[0]-theta)*(grad(u[1])[0]-theta))
  return axial_energy+momentum_enery+shear_energy+combined_energy

def strain_energy(v):
  u,theta=split(v)
  return elastic_energy(u)+(.00000001*inner(grad(theta),grad(theta)))

force_points=interpolate(Constant((0,0,0)),V)
force_pointsnarray=force_points.vector().get_local()
for i in range(len(V.tabulate_dof_coordinates()-1)):
  if boundary_point(V.tabulate_dof_coordinates()[i])==True:
    force_pointsnarray[i]=1
force_points.vector().set_local(force_pointsnarray)

force_intensity=interpolate(Constant((0,1,0)),V)
forcedir=interpolate(Constant((0,0,0)),V)
forcedirarray=forcedir.vector().get_local()

for i in range(len(forcedirarray)-1):
  if boundary_point(V.tabulate_dof_coordinates()[i])==True:
    forcedirarray[i]=force_points.vector().get_local()[i]*force_intensity.vector().get_local()[i]
forcedir.vector().set_local(forcedirarray)

#forcedir=interpolate(Constant((0,0,0)),V)
bc_bottom=DirichletBC(V.sub(0),Constant((0,0)),boundary_bottom)
#bc_point=DirichletBC(V.sub(0).sub(1),Constant(1),boundary_point,method='pointwise')
ene1=strain_energy(v)*dx
ene2=line_energy(v)*ds(2)
ene=ene1+ene2
var1=derivative(ene,v,v_Test)
var2=derivative(var1,v,v_Trial)
#solve(var2==dot(forcedir,v_Test)*dx,v,[bc_bottom])#,bc_point])    # way 1
solve(var2==dot(force_intensity,v_Test)*ds(4),v,[bc_bottom])#,bc_point])    # way 2
u,theta=split(v)
T=project(sigma(u),Mat)
us=project(u,U)
tetas=project(theta,D)
print(assemble(T[1,0]*ds(3)))