Assembling scalar over a single cell

Asked by Jesse Chan on 2013-05-01

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

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Jesse Chan
Solved:
2013-05-02
Last query:
2013-05-02
Last reply:
2013-05-01
Jan Blechta (blechta) said : #1

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)

Jesse Chan (jchan985) said : #2

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?

Jan Blechta (blechta) said : #3

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.

Jesse Chan (jchan985) said : #4

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!

Marie Rognes (meg-simula) said : #5

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.

Jesse Chan (jchan985) said : #6

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.

Jan Blechta (blechta) said : #7

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...