# Efficient matrix assembly over subdomain

**URL:** <https://fenicsproject.discourse.group/t/efficient-matrix-assembly-over-subdomain/7000>\
**Category:** Uncategorized\
**Created:** [November 5, 2021, 5:28pm UTC](https://fenicsproject.discourse.group/t/efficient-matrix-assembly-over-subdomain/7000 "2021-11-05T17:28:16Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![ChristinaNasika](https://avatars.discourse-cdn.com/v4/letter/c/0ea827/32.png) [@ChristinaNasika](https://fenicsproject.discourse.group/u/ChristinaNasika)\
**Post date:** [November 5, 2021, 5:28pm UTC](https://fenicsproject.discourse.group/t/efficient-matrix-assembly-over-subdomain/7000/1 "2021-11-05T17:28:16Z")

</div>

Hello,  
I have a question about how dolfin assembles a stiffness matrix over a subdomain.

I am timing the assembly of a stiffness matrix on a mesh that has 2800 elements.  
Then I define a subdomain on that mesh, with only 10 elements.  
I time the assembly of the same matrix form over the 10 element subdomain.  
It is only about 7 times faster than the full assembly, even though it should only have to loop through 10 cells.  
Does somebody know why that is?

Thank you

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [November 5, 2021, 6:42pm UTC](https://fenicsproject.discourse.group/t/efficient-matrix-assembly-over-subdomain/7000/2 "2021-11-05T18:42:31Z")

</div>

As you have not supplied any code, it is hard to pinpoint why this is the case for you.  
I think one thing your are not taking into account is the setup cost of the sparse matrices.  
See the following example which illustrates run-times where you only iterate over a small sub-set of your mesh.

```python
from dolfin import *

N = 400
frac = 0.02
mesh = UnitSquareMesh(N, N)
V = FunctionSpace(mesh, "CG", 1)

class domain(SubDomain):
    def inside(self, x, on_boundary):
        return x[0]< frac

mf = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
domain().mark(mf, 1)

dx_c = Measure("dx", domain=mesh, subdomain_data=mf)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(grad(u), grad(v))*dx_c
a_form = Form(a)
A = assemble(a_form)

a_sub = inner(grad(u), grad(v))*dx_c(1)
a_sub_form = Form(a_sub)
A_sub = assemble(a_sub_form)

import timeit

runs = 100
full = timeit.timeit(lambda: assemble(a_form, tensor=A), number=runs)

sub = timeit.timeit(lambda: assemble(a_sub_form, tensor=A_sub), number=runs)

print(sub/full, frac, sub, frac)

```

For this example, this returns:  
`0.028362406769913293 0.02 0.20799705399986124 0.02`  
If you reduce the fraction even further, for instance to (0.001) you will observe that we do not get optimal performance any longer, as we get:  
`0.01412187111184388 0.001 0.1051290079994942 0.001`  
This is because there are certain costs associated with assembling a matrix. Assembling over only 10 cells is very little, and the setup cost is dominating the run-time.