boundary conditions by SubDomain vs. MeshFunction

Asked by Miro

Hi, I was trying to use two approaches to defining the domain boundary for a simple Laplace
equation with homog. Dirichlet boundary condition. The first approach uses a SubDomain class, while the other marks the boundary via MeshFunction. I am new to FEniCS, so I thought the two should yield identical results, but only the first one works. Any ideas about where the problem might be? Thank you.

=============================SUBDOMAIN APPROACH=======================
from dolfin import *
mesh = UnitSquare(10,10)
V = FunctionSpace(mesh, 'Lagrange', 1)

class Boundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary

boundary = Boundary()

u_0 = Expression('0')

bc = DirichletBC(V, u_0, Boundary())

u = TrialFunction(V)
v = TestFunction(V)
f = Expression('x[1]*sin(pi*x[0])')
a = inner(nabla_grad(u), nabla_grad(v))*dx
L = f*v*dx

A = assemble(a)
b = assemble(L)
bc.apply(A, b)

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

=========================MESH FUNCTION APPROACH=========================
from dolfin import *
mesh = UnitSquare(10,10)

V = FunctionSpace(mesh, 'Lagrange', 1)

boundary_parts = MeshFunction("uint", mesh, mesh.topology().dim()-1)

class Boundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary

boundary = Boundary()
boundary.mark(boundary_parts, 0)

u_0 = Expression('0')

bc = DirichletBC(V, u_0, boundary_parts, 0)

u = TrialFunction(V)
v = TestFunction(V)
f = Expression('x[1]*sin(pi*x[0])')
a = inner(nabla_grad(u), nabla_grad(v))*dx
L = f*v*dx

# Compute solution
A = assemble(a, exterior_facet_domains=boundary_parts)
b = assemble(L, exterior_facet_domains=boundary_parts)
bc.apply(A, b)

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

Question information

Language:
English Edit question
Status:
Solved
For:
Ubuntu fenics Edit question
Assignee:
No assignee Edit question
Solved by:
Miro
Solved:
Last query:
Last reply:
Revision history for this message
Miro (miroslav-kuchta) said :
#1

Sorry, this obviously belongs to DOLFIN.
anyways, the solution is

boundary = Boundary()
boundary_parts.set_all(4)#some number other than 0
boundary.mark(boundary_parts, 0)