# Assembling scalar over a single cell

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

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:
Last query:
 Revision history for this message Jan Blechta (blechta) said on 2013-05-01: #1

On Wed, 01 May 2013 13:26:11 -0000
Jesse Chan <email address hidden> wrote:
> New question #227944 on DOLFIN:
>
> I'm trying to assemble a quadratic functional cell-wise using a
> computed solution, something akin to evaluating
>
> 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
cell_domains=cell_domains)

 Revision history for this message Jesse Chan (jchan985) said on 2013-05-01: #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?

 Revision history for this message Jan Blechta (blechta) said on 2013-05-01: #3

On Wed, 01 May 2013 20:46:05 -0000
Jesse Chan <email address hidden> wrote:
> Question #227944 on DOLFIN changed:
>
> 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):
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.

 Revision history for this message Jesse Chan (jchan985) said on 2013-05-02: #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!

 Revision history for this message Marie Rognes (meg-simula) said on 2013-05-02: #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)
b = assemble(M)

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

 Revision history for this message Jesse Chan (jchan985) said on 2013-05-02: #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.

 Revision history for this message Jan Blechta (blechta) said on 2013-05-02: #7

On Thu, 02 May 2013 15:56:23 -0000
Jesse Chan <email address hidden> wrote:
> Question #227944 on DOLFIN changed:
>
> 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...