How to calculate the force on a "point Dirichlet boundary"?

Dear all,

I did a simple 2D elastic simulation: a 2D elastic body whose left boundary is fixed and an upward displacement is imposed on its right top corner. I know how to calculate the force on a Dirichlet boundary. but when the Dirichlet boundary is imposed on a point, I don’t know how to calculate the force on it. Below is my code. Do you guys have any idea how to do this?

I appreciate any help. Thank you very much in advance.

from dolfin import *
import numpy as np
import matplotlib.pyplot as plt

#length and width of the 2D elastic body
L = 10
H = 2

mesh = RectangleMesh(Point(0, 0), Point(L, H), 20, 4, "crossed")
# plot(mesh)

#space for displacement
V = VectorFunctionSpace(mesh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)

class Left_Boundary(SubDomain):
    def inside(self, x, on_boundary):
        if( near(x[0], 0, tol) and on_boundary ):
            return True
            return False
class Right_Top_Point(SubDomain):
    def inside(self, x, on_boundary):
        return np.isclose(x[0], L) & np.isclose(x[1], H)

left = Left_Boundary()
right_top_point = Right_Top_Point()

meshFunc = MeshFunction("size_t", mesh, 1)
left.mark(meshFunc, 1) 
ds = Measure("ds", domain=mesh, subdomain_data=meshFunc)  

#material properties
E = Constant(1e5)
nu = Constant(0.3)
lamda_ = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)

def eps(v):
    return sym(grad(v))
def sigma(v):
    return lamda_*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)

f = Constant((0, 0))    #ignore body force
a = inner(sigma(u), eps(v)) *dx
l = inner(f, v)*dx

dbc_left_boundary = DirichletBC(V, Constant((0, 0)), left)
dbc_right_top_point = DirichletBC(V.sub(1), Constant(0.2), right_top_point, method = "pointwise")
dbcs = [dbc_left_boundary, dbc_right_top_point]

u_function = Function(V, name="displacement")
solve(a == l, u_function, dbcs)
File("displacement.pvd") << u_function

#output the force on the left Dirichlet boundary:
force_left_boundary = assemble( (-sigma(u_function)[0,0])*ds(1) )
print("force on the left Dirichlet boundary:", force_left_boundary)

#question: how to calculate the force on the "right_top_point" Dirichlet boundary ?

Project -sigma(u_function)[0,0] into a suitable function space and use point evaluation of that function to get the force at a point.

Thank you for your reply. Do you mean like this?

space_for_stress = TensorFunctionSpace(mesh, "DG", 0)
sigma_function = project( sigma(u_function), space_for_stress )
force_right_top_point = sigma_function(L, H)[3]
print("force on the right_top_point Dirichlet boundary:", force_right_top_point)

But I am not sure if it is correct.

Yes, something along those lines. Please note that if the node of interest is shared between cells, the result is not well defined (unique stress in each cell having the point).

Got it. Thank you very much.