Compute element-wise integral on subdomain

Question

I want to know if there is a way to get the element-wise value of the integral of a functional. The following code gives me the value of the integral over the whole boundary subdomain.

f = inner(nabla_grad(u), nabla_grad(u))*ds(1)
E = assemble(f)

But I want to know if there is a way to get it by element. That is, if the boundary subdomain (marked as "1") has "n" edges, then I want the value of each of the "n" integrals, that when added together match the value of the integral over the whole boundary subdomain (shown above).

The only solution that I can think of is to split the subdomain into "n" different subdomains and perform the integration. However I haven't figured out how to do it, so any help would be appreciated.

Answer 1

The master branch now has an assemble_local function that does what you want. It takes a form and a cell. You will have to build FEniCS yourself to use it, or wait for the next release. Building FEniCS is quite easy if you use the scripts that are described on the web page (fenics-install.sh), at least I was able to do it without too many problems.

You can most likely also split the domain into subdomains and use the current version of FEniCS. I guess you only need one subdomain and then change which cell is marked in a loop. I have not used subdomains much, so I would start by looking to see if there is a tutorial that uses subdomains.

Comments

Thanks for the answer Tormod

I am getting closer now, however the code is still not working. I am doing something wrong with the cells, maybe you can tell:

All_edges = SubsetIterator(boundary_parts, 1)
E_t = 0.0
for cell in All_edges:
     f = inner(nabla_grad(u), nabla_grad(u))*ds(1)
     E = assemble_local(f, cell)
     E_t += E
print E_t

Which returns the TypeError:

** TypeError: in method 'assemble_local', argument 2 of type 'dolfin::Cell const &' **

So clearly "cell" is not a cell, I have also tried:

All_edges = SubsetIterator(boundary_parts, 1)
E_t = 0.0
for edge in All_edges:
     cell = Cell(mesh, edge.index())
     f = inner(nabla_grad(u), nabla_grad(u))*ds(1)
     E = assemble_local(f, cell)
     E_t += E
print E_t

Which returns zeros at the first computed values and errors in the next ones presumably because the index is out of range.

** Error: Unable to create mesh entity.
** Reason: Mesh entity index 4987 out of range [0, 4938] for entity of dimension 2.

Which is confusing because the iterator said there were more entities:

** Iterating over subset, found 16 entities out of 7556. **

So am I confusing cells with facets?, Instead of "Cell" I have also tried using "Facet" but then TypeError comes back.
Any help would be great!

---

You need to pass it a Cell object. Just get the index of the cell connected to the facet and pass this cell, Cell(mesh, cell_index). For cells with multiple ds facets (multiple external facets) you will have to use a FacetFunction to mark the facet you are interested in. It will probably be quite a lot faster than integrating globally, but you will need the marking function.

As far as I know all assembly routines in dolfin takes a cell and iterates over interior, interior and exterior facets and points, you cannot assemble over a facet on its own

---

Thanks a lot Tormod!
It works as expected!

Posting here my final solution:

D = mesh.topology().dim()
# Build connectivity between facets and cells
mesh.init(D-1,D) 
# Initialize total 
E_t = 0.0
for edge in SubsetIterator(boundary_parts, 1):
    cell = Cell(mesh, edge.entities(D)[0])
    f = inner(nabla_grad(u), nabla_grad(u))*ds(1)
    E = assemble_local(f, cell)
    E_t += E
print E_t

If there are facets with 2 or more cells assigned the following line should be changed:

cell = Cell(mesh, edge.entities(D)[0])