Newbie: Mixed Neumann on top of a box

Question

Dear all,

I am studying Fenics, and I find it fascinating. I am playing right now with a simple elastic problem. However, I am finding not quite easy to impose a mixed Neumann-Dirichlet and I hope you can guide me.

The problem is quite easy: a box, with clamping on the bottom side, and a force on the top side.

My attempt right now is below. I've tried to put a force on the boundary, but I don't know if (or how) I could define the force as a dirac delta on the top side, but I think this approach will cause more harm than good... am I right?

My question is then: do I have to decompose the box into five domains, with the 2D domain, and the four sides, or just two? I am reading this documentation, and I think I should mark the top, defining a class, and just use a function for the Dirichlet condition on the bottom. I've tried this, but nothing happens on the deformation side.

EDIT: I've changed my code to apply a subdomain condition. It compiles and runs, but I cannot see any deformation: u is constant zero on the domain.

EDIT: Now the on-top condition handles the top boundary actually. However, nothing changes, and I don't know how to proceed...

Can you help me in figuring out what I am doing wrong?

from dolfin import *

# Constants
Young   = 30e6
Poisson = 0.25
Force   = 1.0
Width   = 7.0
Height  = 14.0
Mu      = Young / (2.0 * (1.0 + Poisson))
Lambda  = Young * Poisson / ((1.0 + Poisson) * (1.0 - 2.0 * Poisson))

# Floating point tolerance
Epsilon = 10e-7

# Create mesh and define function space (x0, y0)-(x1, y1), nx, ny
mesh = RectangleMesh(0, 0, Width, Height, 10, 20)

# 1st order Lagrange shape functions
V = VectorFunctionSpace(mesh, 'CG', 1)

# Displacement function
u = TrialFunction(V)

# Test function
v = TestFunction(V)

# Strain
def eps(v):
    return sym(grad(v))

# Stress
def sigma(v):
    return 2.0 * Mu * eps(v) + Lambda * tr(eps(v)) * Identity(v.geometric_dimension())

# Bulk Force
f = Constant((0, 0))

# Set up boundary conditions
u0 = Constant((0.0, 0.0))
f1 = Constant((0.0, Force * 1000))

def on_bottom(x, on_boundary):
    return on_boundary and abs(x[1]) < Epsilon

def on_top(x, on_boundary):
    return on_boundary and abs(x[1] - Height) < Epsilon

force_vector = Expression(("0.0", "1000"))

# Now let's move to the boundary
boundary_parts = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)

# The top boundary class
class TopBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0] - Height) < Epsilon

# Mark the top boundary with 0
NeumannBoundary = TopBoundary()
NeumannBoundary.mark(boundary_parts, 0)

bcBottom = DirichletBC(V, u0, on_bottom)
bcs      = [bcBottom]

# Variational form
a = inner(sigma(u), eps(v)) * dx
L = inner(f, v) * dx + 10000 * inner(force_vector, v) * ds(0)

Comments

I've changed the class to

class TopBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and abs(x[0] - Height) < Epsilon

Nothing changes... What am I doing wrong here? Thanks for any help! :)

Answer 1

Solved! Following this example, I've managed to solve my problem. The main differences are: don't solve, assemble & solve:

# Variational form
a = inner(sigma(u), eps(v)) * dx
L = inner(f, v) * dx + inner(force_vector, v) * ds(1)

# Compute solution
A = assemble(a, exterior_facet_domains = boundary_parts)
b = assemble(L, exterior_facet_domains = boundary_parts)
for condition in bcs: condition.apply(A, b)

# SOLVE!
u = Function(V)
solve(A, u.vector(), b, 'lu')

Damn, that wasn't quite easy.

Comments

*_domains of assemble has been deprecated in DOLFIN version 1.4.0.