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)
# plt.show()

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

tol=1e-14
class Left_Boundary(SubDomain):
def inside(self, x, on_boundary):
if( near(x[0], 0, tol) and on_boundary ):
return True
else:
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)
meshFunc.set_all(0)
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):
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.

``````space_for_stress = TensorFunctionSpace(mesh, "DG", 0)