Example of custom kernel for facet assembly

General
1

I am looking for an example of a custom kernel that assembles a facet contribution (interior or exterior), either in c++ or Python. I found such examples for cells in the demos / tests but not for facets. In particular I am curious about what is done with the local_facet and permutation arguments.

2

For a custom kernel, you might not need them ( that depends on the kernel).
What local_facet is is the local index of a facet (relative to the input cell). This is used to accessed the correct coordinates of the mesh nodes for this facet, and appropriate basis values.

The permutation integer is used to ensure that for dS integrals for two adjacent cells, that the facet quadrature rule is aligned.

Here is an example of manual usage of the ufcx kernel for exterior facets:

from mpi4py import MPI
import numpy as np
import dolfinx
import cffi

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 1, 2, dolfinx.cpp.mesh.CellType.triangle)

V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1, (2, )))
u = dolfinx.fem.Function(V)
u.interpolate(lambda x: (np.sin(2.3*x[0]), x[1]))

def right_facets(x):
    return np.isclose(x[0], 1)


facets = dolfinx.mesh.locate_entities_boundary(mesh, mesh.topology.dim-1, right_facets)

integration_entities = dolfinx.fem.compute_integration_domains(dolfinx.fem.IntegralType.exterior_facet,mesh.topology, facets, mesh.topology.dim-1).reshape(-1, 2)

import ufl

ds = ufl.Measure("ds", domain=mesh, subdomain_data=([(2, integration_entities.flatten())]), subdomain_id=2)
v = ufl.TestFunction(V)
L = dolfinx.fem.form(ufl.inner(u, v) * ds)



x = mesh.geometry.x
v = mesh.geometry.dofmap

dtype= np.float64
entity_local_index = np.array([0], dtype=np.intc)
perm = np.array([0], dtype=np.uint8)
geometry = np.zeros((3, 3), dtype=x.dtype)
constants = np.zeros(1, dtype=dtype)

b_local = np.zeros(V.dofmap.dof_layout.num_dofs*V.dofmap.bs, dtype=dtype)

mesh.topology.create_entity_permutations()
perms = mesh.topology.get_facet_permutations()

# Assume all cells has the same number of facets

num_facets_per_cell = dolfinx.cpp.mesh.cell_num_entities(mesh.topology.cell_type, mesh.topology.dim-1)

kernel = getattr(L.ufcx_form.form_integrals[0], f"tabulate_tensor_{np.dtype(dtype).name}")
ffi = cffi.FFI()
b_func = dolfinx.fem.Function(V)
all_coeffs = dolfinx.fem.assemble.pack_coefficients(L._cpp_object)[(dolfinx.fem.IntegralType.exterior_facet, 2)].reshape(integration_entities.shape[0], -1)
coeffs = np.zeros(all_coeffs.shape[1], dtype=dtype)

for i, (cell, local_facet_index) in enumerate(integration_entities):
    for j in range(v.shape[1]):
        geometry[j] = x[v[cell, j], :]
    entity_local_index[0] = local_facet_index
    perm[0] = perms[cell*num_facets_per_cell + local_facet_index]
    coeffs[:] = all_coeffs[i]
    b_local.fill(0.0)
    kernel(
        ffi.from_buffer(b_local),
        ffi.from_buffer(coeffs),
        ffi.from_buffer(constants),
        ffi.from_buffer(geometry),
        ffi.from_buffer(entity_local_index),
        ffi.from_buffer(perm),
    )
    local_dofs = V.dofmap.cell_dofs(cell)
    for j in range(3):
        for k in range(V.dofmap.bs):
            b_func.x.array[local_dofs[j]*V.dofmap.bs+k] += b_local[j*V.dofmap.bs+k]


b_ref = dolfinx.fem.assemble_vector(L)
np.testing.assert_allclose(b_func.x.array, b_ref.array)
2 Likes