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

General
1

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):
    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 ?
2

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.

3

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.

4

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).

5

Got it. Thank you very much.