Zero results with Neumann and Dirichlet mixed boundary conditions

Question

Dear all,

I am trying to solve the 2-d static elasticity problem in a rectangle domain. left and bottom are fixed. When I applied displacement at the top in y-direction, I obtain the correct results. But when I applied a traction at the top in y-direction, I always got the results as zero. I think the problem is related with the form of integral on the boundaries, but can't figure out the bug. Anyone's help will be greatly appreciated.

from dolfin import *

# Create mesh
x0, y0, x1, y1, nx, ny = 0.0, 0.0, 1.0, 2.0, 50, 100
mesh = RectangleMesh(x0, y0, x1, y1, nx, ny, 'crossed')

# Create function space
V = VectorFunctionSpace(mesh, "Lagrange", 2)

# Create test and trial functions, and source term
u, w = TrialFunction(V), TestFunction(V)
b = Constant((0.0, 0.0))

# Elasticity parameters
E, nu = 100.0e9, 0.3
mu, lmbda = E/(2.0*(1.0 + nu)), E*nu/((1.0 + nu)*(1.0 - 2.0*nu))

# Stress
sigma = 2*mu*sym(grad(u)) + lmbda*tr(grad(u))*Identity(w.cell().d)

#define boundary left side
class left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0.0)

#define boundary right side
class right(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 1.0)

#define boundary top side
class top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 2.0)

#define boundary bottom side
class bottom(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0.0)

#initialize subdomains
left = left()
right = right()
top = top()
bottom = bottom()

#initialize mesh function for boundary functions and tag boundaries
boundaries = FacetFunction('size_t', mesh)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)

#define new measure associated with boundaries
ds = Measure('ds')[boundaries]

# Define boundary condiition 
# **Apply displacement**                                 
zero = Constant(0.0)
bc1 = DirichletBC(V.sub(0), zero, left)     # u.n = 0
bc2 = DirichletBC(V.sub(1), zero, bottom)   # u.n = 0
bc3 = DirichletBC(V.sub(0), zero, right)    # u.n = 0
bc4 = DirichletBC(V.sub(1), 2.0*0.05, top)     # u.n = 0
bcs = [bc1, bc2, bc4]

# **Apply traction**
#t = Constant((0.0, 1.0))
#bcs = [bc1, bc2]

# Governing balance equation
a = inner(sigma, grad(w))*dx 
L = dot(b, w)*dx #+ dot(t, w)*ds(2)

# Set up PDE and solve
u = Function(V)
problem = LinearVariationalProblem(a, L, u, bcs)
solver = LinearVariationalSolver(problem)
solver.parameters["symmetric"] = True
solver.solve()

sigma_w = project(2*mu*sym(grad(u)) + lmbda*tr(grad(u))*Identity(w.cell().d))

# Plot Results
plot(mesh, axes=True, title = "Meshed Geometry" )
plot(u, mode="displacement", axes=True,title = "Displacement Vectors" )
plot(sigma_w[1, 1], axes=True, title = "sigma22", interactive=True) 

Comments

Did you get this message?

*** Warning: Bilinear and linear forms do not have same exterior facet subdomains in SystemAssembler. Taking subdomains from bilinear form

I think you need to have a term in ds on the LHS (even if it is zero). Maybe this is a "feature" of ffc...

---

Maybe try adding a term like

Constant(0.0)*dot(u, w)*ds(2)

to a.

---

Dear Chris,

Thank you. It works now. Can you tell me why should I add that term? That's very strange because in a there is no surface integral actually. Does it mean that a and L both should have surface integral, if it exists?

---

I think the subdomains need to be defined for both forms, 'a' and 'L'. If not, it looks at 'a' by default. There may be another way to do this, but adding a zero valued term seems to work.

---

Dear Chris_richardson:

Would you like to post your answer again? I will then select it as the answer.

---

I see. Thank you so much.