# Assembling scalar over a single cell

I'm trying to assemble a quadratic functional cell-wise using a computed solution, something akin to evaluating

J = inner(grad(

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:/

>

> I'm trying to assemble a quadratic functional cell-wise using a

> computed solution, something akin to evaluating

>

> J = inner(grad(

> energy = assemble(J, "a specific cell")

>

> Is there a way to integrate over just a specific cell?

>

> Thanks

>

cell_domains = CellFunction(

for c in cells(mesh):

cell_domains[c] = ...

# Compute energy functional on cells with cell_domains[c]=42

energy = assemble(

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:/

>

> 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://

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:

else

> - 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:/

>

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