Mark elements on boundary - minimal code included

Question

Dear all,

this is a very straightforward code which is not working properly as it marks only one boundary element. You can see everything just by running the code below. I suppose I made a foolish mistake somewhere as I am a novice in both python and fenics...

from dolfin import *

mesh = UnitSquare(12,12)

# Create FunctionSpaces
V =  FunctionSpace(mesh, 'CG', 1)

# Boundaries case
def dOmega(x, on_boundary): return x[0] > (1.0 - DOLFIN_EPS) or x[0] < DOLFIN_EPS or x[1] > (1.0 - DOLFIN_EPS) or x[1] < DOLFIN_EPS
g0 = Constant(0.0)

# No-slip boundary condition for velocity
bc0 = DirichletBC(V, g0, dOmega)

# Parameters and functions
f = Constant(1.)
uh = TrialFunction(V)
vh = TestFunction(V)
a = dot(grad(uh), grad(vh))*dx
L = f*vh*dx

# Compute solution
uh = Function(V)
solve(a == L, uh, bc0)

# Mark boundary adjacent cells
boundary_adjacent_cells = [myCell for myCell in cells(mesh)
                                  if any([myFacet.exterior() for myFacet in facets(myCell)])]
cell_domains = CellFunction('int', mesh)
cell_domains.set_all(1)
for myCell in boundary_adjacent_cells:
    cell_domains[myCell] = 0

# Plot cell_domains
plot(cell_domains, interactive=True)

# Try to integrate over the whole domain and over interior cells
integral = assemble(uh*dx)
print 'Integral = %e\r'%integral
integral2 = assemble(uh*dx(1))
print 'Integral2 = %e\r'%integral2

# Plot solution and mesh
plot(uh)
plot(mesh)

# Hold plot
interactive()

Other problem is, that the second integral is zero. Because if we integrate over the only marked boundary cell, the integral should be a positive number...

Comments

I ran your code as it appeared in the window above, and I got more than one boundary element marked - in fact all elements that touched the boundary were marked as "0" in the first plot. What version of FEniCS are you running? To find out, from a Python shell, you could type:

from dolfin import *
dolfin_version()
'1.3.0+'

The output of your code on my machine was:
Integral = 3.436498e-02
Integral2 = 0.000000e+00

---

Thank you very much for your reply, Timm.

The result of

from dolfin import *
dolfin_version()

is

Traceback (most recent call last):
  File "python.py", line 2, in <module>
    dolfin_version()
NameError: name 'dolfin_version' is not defined

I will now install a newer version, since according to Ubuntu software center my current version is fenics 1:1.0.0-1.
The current output is the same on my machine. The current result from the first plot is in my case:

---

So. My version of dolphin now is 1.3.0.
There was a problem with update, which is described here:
http://fenicsproject.org/download/ubuntu_details.html
And what about my example?
Now I got the boundary elements properly:

Now, it left to solve the issue, that integral2 is still 0. I thought that assemble(uh*dx(1)) means, that we are integrating over the "orange" domain from the image...

Answer 1

Hi, you need to define this new measure associated with the cell function

dx = Measure('dx')[cell_domains]
# Try to integrate over the whole domain
integral = assemble(uh*(dx(0)+dx(1)))
print 'Integral = %e\r'%integral
# and over interior cells
integral2 = assemble(uh*dx(1))
print 'Integral2 = %e\r'%integral2

Comments

Great, now the whole problem is solved.

Just to repeat, basically the problem was in 2 different things:
- old version of fenics from main Ubuntu repository
- incorrect usage of measures

Thank you both very much! (timm and MiroK)

Fully working python script:

from dolfin import *

mesh = UnitSquareMesh(12,12)

# Create FunctionSpaces
V =  FunctionSpace(mesh, 'CG', 1)

# Boundaries case
def dOmega(x, on_boundary): return x[0] > (1.0 - DOLFIN_EPS) or x[0] < DOLFIN_EPS or x[1] > (1.0 - DOLFIN_EPS) or x[1] < DOLFIN_EPS
g0 = Constant(0.0)

# No-slip boundary condition for velocity
bc0 = DirichletBC(V, g0, dOmega)

# Parameters and functions
f = Constant(1.)
uh = TrialFunction(V)
vh = TestFunction(V)
a = dot(grad(uh), grad(vh))*dx
L = f*vh*dx

# Compute solution
uh = Function(V)
solve(a == L, uh, bc0)

# Mark boundary adjacent cells
boundary_adjacent_cells = [myCell for myCell in cells(mesh)
                                  if any([myFacet.exterior() for myFacet in facets(myCell)])]
cell_domains = CellFunction('size_t', mesh)
cell_domains.set_all(1)
for myCell in boundary_adjacent_cells:
    cell_domains[myCell] = 0
dx = Measure('dx')[cell_domains]

# Plot cell_domains
plot(cell_domains, interactive=True)

# Try to integrate over the whole domain and over interior cells
integral = assemble(uh*(dx(0)+dx(1)))
print 'Integral = %e\r'%integral
integral2 = assemble(uh*dx(1))
print 'Integral2 = %e\r'%integral2

# Plot solution and mesh
plot(uh)
plot(mesh)

# Hold plot
interactive()