Mesh connectivity arrays, vectorized computaitons

mesh
1

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)
2

You might first look into whether your estimator can be expressed in UFL, as a Form, then assembled, 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/

1 Like
3

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())
4

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])
2 Likes