# Mesh connectivity arrays, vectorized computaitons

**URL:** <https://fenicsproject.discourse.group/t/mesh-connectivity-arrays-vectorized-computaitons/6545>\
**Category:** mesh\
**Created:** [September 10, 2021, 4:13pm UTC](https://fenicsproject.discourse.group/t/mesh-connectivity-arrays-vectorized-computaitons/6545 "2021-09-10T16:13:05Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![beschlcc](https://avatars.discourse-cdn.com/v4/letter/b/91b2a8/32.png) [@beschlcc](https://fenicsproject.discourse.group/u/beschlcc)\
**Post date:** [September 10, 2021, 4:13pm UTC](https://fenicsproject.discourse.group/t/mesh-connectivity-arrays-vectorized-computaitons/6545/1 "2021-09-10T16:13:05Z")

</div>

Hi everybody,

working with a-posteriori error estimates, I am interested in edge-cell connectivity to compute gradient fluxes over edges and edge-vertex connectivity to compute edge normals fast.

Is there a way to avoid a for loop to obtain these connectivities as arrays to work with them in vectorized computations?

Currently I am doing it as in the code example below, separating by boundary/non-boundary edges, but it slows down the vectorized a-posteriori error estimation quite bit.

Any help is appreciated! Thank you very much.

```
import fenics as fn
import numpy as np
from time import time

N = 400
mesh = fn.UnitSquareMesh(N,N)

time_init = 0
time_cells = 0
time_append = 0

start = time()
mesh.init(1,2) # Initialise mesh for edgewise computations
mesh.init(1,0)
top = mesh.topology()
c12 = top(1,2)    
c10 = top(1,0)
num_edges = mesh.num_facets()
time_init = time() - start

inner_edge_cells = []
outer_edge_cells = []

inner_edge_vertices = []
outer_edge_vertices = []

for ind_edge in range(num_edges):  
     start = time()
     cells = c12(ind_edge) # get cell/cells adjacent to edge
     time_cells += time() - start

     if(len(cells)==1):
         start = time()
         outer_edge_cells.append(cells[0])
         outer_edge_vertices.append(c10(ind_edge)) # get vertices adjacent to edge
         time_append += time() - start
     else:
         start = time()
         inner_edge_cells.append([cells[0], cells[1]]) 
         inner_edge_vertices.append(c10(ind_edge)) # get vertices adjacent to edge
         time_append += time() - start

outer_edge_vertices_arr = np.array(outer_edge_vertices)
inner_edge_vertices_arr = np.array(inner_edge_vertices)

inner_edge_cells_arr = np.array(inner_edge_cells)
outer_edge_cells_arr = np.array(outer_edge_cells)        
    
print("Time init: ", time_init)
print("Time cells: ", time_cells)
print("Time append: ", time_append)
```

---

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [September 10, 2021, 4:34pm UTC](https://fenicsproject.discourse.group/t/mesh-connectivity-arrays-vectorized-computaitons/6545/2 "2021-09-10T16:34:20Z")

</div>

You might first look into whether your estimator can be expressed in UFL, as a `Form`, then `assemble`d, following the standard FEniCS workflow. See, e.g., this archived thread from the old Q&A forum:

[https://fenicsproject.org/qa/1949/a-posteriori-error-estimator-face-residuals/](https://fenicsproject.org/qa/1949/a-posteriori-error-estimator-face-residuals/)

---

<div class="post-metadata">

**Author:** ![beschlcc](https://avatars.discourse-cdn.com/v4/letter/b/91b2a8/32.png) [@beschlcc](https://fenicsproject.discourse.group/u/beschlcc)\
**Post date:** [September 10, 2021, 10:07pm UTC](https://fenicsproject.discourse.group/t/mesh-connectivity-arrays-vectorized-computaitons/6545/3 "2021-09-10T22:07:56Z")

</div>

Hi, thank you very much for your fast response.  
I have been working with assemble\_local before. Unfortunately it is too slow for what I am doing, because I am solving stochastic PDEs and need to compute the error estimates for e.g. 10000 samples. That’s why I am working on a vectorization of the error estimates.

When using assemble only, which might be faster, as in the example you sent, and for the MWE below, my Spyder Kernel dies when accessing tensor=cell\_residual.vector().

I am using dolfin 2019.1.0 with Python 3.8.3 and Spyder 4.1.3 in case it is relevant.

```
import fenics as fn

N = 20
mesh = fn.UnitSquareMesh(N,N)
# h = fn.CellSize(mesh) # does not exist anymore, alternative function?
n = fn.FacetNormal(mesh)

Vcg1 = fn.FunctionSpace(mesh, 'CG', 1)
Vdg0 = fn.FunctionSpace(mesh, 'DG', 0)

cell_residual = fn.Function(Vdg0)

u_expr = fn.Expression("x[0]*x[1]", degree = 1)
u = fn.project(u_expr, Vcg1)

residual = fn.jump(fn.grad(u), n)*fn.dS
fn.assemble(residual, tensor=cell_residual.vector())
```

---

<div class="post-metadata">

**Author:** ![beschlcc](https://avatars.discourse-cdn.com/v4/letter/b/91b2a8/32.png) [@beschlcc](https://fenicsproject.discourse.group/u/beschlcc)\
**Post date:** [September 13, 2021, 1:39pm UTC](https://fenicsproject.discourse.group/t/mesh-connectivity-arrays-vectorized-computaitons/6545/4 "2021-09-13T13:39:14Z")

</div>

Hi, here is a little update on how I now obtain the vectorized mesh connectivity data without a for loop:

```
import fenics as fn
import numpy as np

N = 400
mesh = fn.UnitSquareMesh(N,N)

# Initialise mesh for edgewise computations
mesh.init(1,2) 
mesh.init(1,0)
top = mesh.topology()
c12 = top(1,2)    
c10 = top(1,0)
num_edges = mesh.num_facets()

# Obtain boundary markers
boundary = fn.CompiledSubDomain('on_boundary')
boundary_function = fn.MeshFunction('size_t', mesh, dim=1) 
boundary_function.set_all(0)
boundary.mark(boundary_function, 1)
boundary_markers = boundary_function.array()

# Obtain edge vertex connectivity
edge_vertices_array = np.fromiter(c10(), dtype=int).reshape([num_edges, 2])
outer_edge_vertices_arr = edge_vertices_array[boundary_markers==1]
inner_edge_vertices_arr = edge_vertices_array[boundary_markers==0]

# Obtain edge cell connectivity
boundary_markers_cells_2 = np.insert(boundary_markers, np.where(boundary_markers==0)[0], 0)
edge_cells_array = np.fromiter(c12(), dtype=int)
outer_edge_cells_arr = edge_cells_array[boundary_markers_cells_2==1]
inner_edge_cells_arr = edge_cells_array[boundary_markers_cells_2==0].reshape([num_edges-np.count_nonzero(boundary_function.array()), 2])
```