Integrals over (sub)domain

Question

Dear all,

I haven't found the answer in the book, so I am resorting to the forum. I have a domain, and, apart from solving a problem on it, I'd like to perform some integrations over it, and over its subdomains.

So, If I have dx and ds, like in the code below, how can I perform an integration over ds(1)?

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

The bigger picture is that I will implement an iterative solver, and at each iteration, I need to compute some properties on the domain and on its parts.

Thanks!

Answer 1

Hi, subdomain integration is demonstrated in this demo.

Comments

Thanks MiroK, I see how I can do it. But my code won't work:

mesh = RectangleMesh(0, 0, Width, Height, 10, 20)

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

# Force to be applied to the top
force_vector = Expression(("0", "1000000"))

# 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[1] - Height) < Epsilon

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

# Dirichlet condition on the bottom boundary
bcBottom = DirichletBC(V, u0, on_bottom)
bcs      = [bcBottom]

# Variational form
a = inner(sigma(u), eps(v)) * dx
L = inner(f, v) * dx + 1*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(A, u.vector(), b, 'gmres')

# Integrals
markers = FacetFunctionSizet(mesh, 1)
TopBoundary().mark(markers, 1);

ds = ds[markers]
lx = 1.0 * ds(1)

totlen = assemble(lx)

print 'Length of top boundary by integration: ', totlen

The error I'm getting is this:

Traceback (most recent call last):
  File "linelastic.py", line 251, in <module>
    totlen = assemble(lx)
  File "/Applications/FEniCS.app/Contents/Resources/lib/python2.7/site-packages/dolfin/fem/assembling.py", line 230, in assemble
    cell_domains, exterior_facet_domains, interior_facet_domains)
  File "/Applications/FEniCS.app/Contents/Resources/lib/python2.7/site-packages/dolfin/fem/assembling.py", line 80, in _create_dolfin_form
    mesh=mesh)
  File "/Applications/FEniCS.app/Contents/Resources/lib/python2.7/site-packages/dolfin/fem/form.py", line 95, in __init__
    assert form.domains(), "Expecting a completed form with domains at this point."
AssertionError: Expecting a completed form with domains at this point.

What am I doing wrong?

---

Could you provide your the dolfin version? Btw, you declare the facet function markers such that all facets are marked as 1 and then you mark the facets by TopBoundary as 1 again, i.e. nothing changes and the assembled form returns circumference of the entire mesh.

---

I just downloaded the app for MacOS X, 1.4.0.

Since I am quite a newbie, do you mean that my code regarding marking boundaries can be reduced to something simpler? I've tried the following but I always have the same error...

ds = ds[1]
lx = 1.0 * ds(1)

totlen = assemble(lx)

---

This shoud fix your problem and also illustrates the error with facet function initialization I mentioned above

from dolfin import *

Width, Height = 5, 10
mesh = RectangleMesh(0, 0, Width, Height, 10, 20)

class TopBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], Height)

# Your old way
markers = FacetFunctionSizet(mesh, 1)  # Note initialization to 1
ds = ds[markers]
TopBoundary().mark(markers, 1) # marking to 1, but all facets are 1 already
form = 1*ds(1, domain=mesh)
ans = assemble(form)   # circumnference
assert abs(ans - 2*(Width + Height)) < 1E-13
print ans

# Correct way
markers = FacetFunctionSizet(mesh, 0) # Note initialization to 0
ds = ds[markers]
TopBoundary().mark(markers, 1) # marking to 1 while all other facets are 0
form = 1*ds(1, domain=mesh)
ans = assemble(form)   # circumnference
assert abs(ans - Width) < 1E-13
print ans

---

Aha! I see now, since I redefined marks, I needed to re-assign 1 just to the top boundary. Then, I assemble the form on the domain, using only the marked region ds(1) on the whole mesh ds(1, domain=mesh).

Thanks!