# assemble_system on subdomains

I am looking at the example demo_subdomains

If we just replace

solve(a==L, u, bcs)

by the algebra formulation, that is,

A,b=assemble_

solve(A,

it gives different results. Can anyone verify that? Aren't these two formulation supposed to give the same result? Thanks.

The program is attached:

from dolfin import *

# Create classes for defining parts of the boundaries and the interior

# of the domain

class Left(SubDomain):

def inside(self, x, on_boundary):

return near(x[0], 0.0)

class Right(SubDomain):

def inside(self, x, on_boundary):

return near(x[0], 1.0)

class Bottom(SubDomain):

def inside(self, x, on_boundary):

return near(x[1], 0.0)

class Top(SubDomain):

def inside(self, x, on_boundary):

return near(x[1], 1.0)

class Obstacle(

def inside(self, x, on_boundary):

return (between(x[1], (0.5, 0.7)) and between(x[0], (0.2, 1.0)))

# Initialize sub-domain instances

left = Left()

top = Top()

right = Right()

bottom = Bottom()

obstacle = Obstacle()

# Define mesh

mesh = UnitSquareMesh(64, 64)

# Initialize mesh function for interior domains

domains = CellFunction(

domains.set_all(0)

obstacle.

# Initialize mesh function for boundary domains

boundaries = FacetFunction(

boundaries.

left.mark(

top.mark(

right.mark(

bottom.

# Define input data

a0 = Constant(1.0)

a1 = Constant(0.01)

g_L = Expression("- 10*exp(- pow(x[1] - 0.5, 2))")

g_R = Constant("1.0")

f = Constant(1.0)

# Define function space and basis functions

V = FunctionSpace(mesh, "CG", 2)

u = TrialFunction(V)

v = TestFunction(V)

# Define Dirichlet boundary conditions at top and bottom boundaries

bcs = [DirichletBC(V, 5.0, boundaries, 2),

# Define new measures associated with the interior domains and

# exterior boundaries

dx = Measure(

ds = Measure(

# Define variational form

F = (inner(a0*grad(u), grad(v))*dx(0) + inner(a1*grad(u), grad(v))*dx(1)

- g_L*v*ds(1) - g_R*v*ds(3)

- f*v*dx(0) - f*v*dx(1))

# Separate left and right hand sides of equation

a, L = lhs(F), rhs(F)

# Solve problem

u = Function(V)

#solve(a == L, u, bcs)

A,b=assemble_

solve(A,

# Evaluate integral of normal gradient over top boundary

n = FacetNormal(mesh)

m1 = dot(grad(u), n)*ds(2)

v1 = assemble(m1)

print "\int grad(u) * n ds(2) = ", v1

# Evaluate integral of u over the obstacle

m2 = u*dx(1)

v2 = assemble(m2)

print "\int u dx(1) = ", v2

# Plot solution and gradient

plot(u, title="u")

plot(grad(u), title="Projected grad(u)")

interactive()

## Question information

- Language:
- English Edit question

- Status:
- Answered

- For:
- DOLFIN Edit question

- Assignee:
- No assignee Edit question

- Last query:
- 2013-02-09

- Last reply:
- 2013-02-11

Johannes Ring (johannr) said : | #1 |

Yes, I can confirm this for DOLFIN 1.1.0. Please report a bug.

It works fine with current trunk.

Johannes Ring (johannr) said : | #2 |

I see you have already made a bug report.

Jan Blechta (blechta) said : | #3 |

Until the bug is fixed you can workaround this by explicitly passing cell_domains to assemble_system() method:

A,b=assemble_

## Can you help with this problem?

Provide an answer of your own, or ask Shiyuan for more information if necessary.