Assembling scalar over a single cell

On Wed, 01 May 2013 13:26:11 -0000
Jesse Chan <email address hidden> wrote:
> New question #227944 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/227944
>
> I'm trying to assemble a quadratic functional cell-wise using a
> computed solution, something akin to evaluating
>
> J = inner(grad(u),grad(u))*dx
> energy = assemble(J, "a specific cell")
>
> Is there a way to integrate over just a specific cell?
>
> Thanks
>

cell_domains = CellFunction('size_t', mesh)

for c in cells(mesh):
    cell_domains[c] = ...

# Compute energy functional on cells with cell_domains[c]=42
energy = assemble(inner(grad(u), grad(u))*dx(42),
                  cell_domains=cell_domains)

Thanks Jan - I'm a bit newer to Fenics, however. Can you elaborate a bit on the meaning of some of the code? For example

- What is 'size_t' referring to?
- What should be put in cell_domains[c] = ...?
- Is dx(42) implying that the integration is done over the cell at the 42nd index?

On Wed, 01 May 2013 20:46:05 -0000
Jesse Chan <email address hidden> wrote:
> Question #227944 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/227944
>
> Jesse Chan posted a new comment:
> Thanks Jan - I'm a bit newer to Fenics, however. Can you elaborate a
> bit on the meaning of some of the code? For example
>
> - What is 'size_t' referring to?

From technical point of view it is type which can hold arbitrary memory
adress on considered platform. From practical point of view it is
unsigned int 32-bit/64-bit (depending on platform).
http://en.wikipedia.org/wiki/Size_t#stddef.h

It enables DOLFIN to use maximum indexing range on arbitrary platform.

> - What should be put in cell_domains[c] = ...?

There you need to feed up cell_domains with data you want. For example
if you want to pick cell by its indices you can do:

for c in cells(mesh):
   cell_domains[c] = c.index()

Or if you want to pick say small cells:

for c in cells(mesh):
   if c.inradius() <= 0.042:
       cell_domains[c] = 42
   else
       cell_domains[c] = 666

> - Is dx(42) implying that the integration is done over the cell at
> the 42nd index?

That depends what is in cell_domains. It only says that integration is
performed over cells where cell_domains takes value 42.

Thanks Jan. This was very helpful, and it allowed me to implement what I needed successfully.

For context, I am computing an error indicator over elements. I just found the adaptive poisson demo, and am using it to look for ways to speed this up. Thanks!

You might also want to try the following "trick" that takes advantage of the fact that each DG0 basis function is 1 on a certain cell and zero elsewhere:

DG0 = FunctionSpace(mesh, "DG", 0)
w = TestFunction(DG0)
M = w*inner(grad(u), grad(u))*dx
b = assemble(M)

Now b will be a vector containing the integral over each cell for each cell.

Marie, that appears to be much more efficient in evaluating error indicators than what I was doing. Thank you very much.

Out of curiosity, with the vector that is returned by assemble(M), is there a simple way to sort the contents in descending order? Something similar to sort(...) in Matlab? I'd like to take the indices of the largest entries.

On Thu, 02 May 2013 15:56:23 -0000
Jesse Chan <email address hidden> wrote:
> Question #227944 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/227944
>
> Jesse Chan posted a new comment:
> Marie, that appears to be much more efficient in evaluating error
> indicators than what I was doing. Thank you very much.

Good trick - bypasses form recompilation for every cell.

>
> Out of curiosity, with the vector that is returned by assemble(M), is
> there a simple way to sort the contents in descending order?
> Something similar to sort(...) in Matlab? I'd like to take the
> indices of the largest entries.
>

Extract numpy array:
v = assemble(...)
v_numpy = v.array()

There will certainly be this functionality in numpy...