I do not have a complete answer as of now (time constraints), but here is what I would do:
import dolfinx
import numpy as np
from mpi4py import MPI
import basix
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 2, 2)
V = dolfinx.fem.FunctionSpace(mesh, ("N1curl", 2))
tdim = mesh.topology.dim
b_el = V.element.basix_element
mesh.topology.create_connectivity(tdim-1, tdim)
mesh.topology.create_connectivity(tdim, tdim-1)
f_to_c = mesh.topology.connectivity(tdim-1, tdim)
c_to_f = mesh.topology.connectivity(tdim, tdim-1)
exterior_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
cell_facet_pairs = np.zeros((len(exterior_facets), 2), dtype=np.int32)
for i, facet in enumerate(exterior_facets):
cell = f_to_c.links(facet)[0]
l_facets = c_to_f.links(cell)
l_index = np.argwhere(l_facets == facet)[0, 0]
cell_facet_pairs[i] = [cell, l_index]
interpolation_points = b_el.x
# First derivatives of basis functions for the coordinate element (entity_dim, entity_index, (num_points, num_basis_functions))
basis_functions = []
b_cell = basix.cell.string_to_type(mesh.topology.cell_name())
cmap_element = basix.create_element(basix.finite_element.string_to_family(
"Lagrange", mesh.topology.cell_name()), b_cell,
mesh.geometry.cmap.degree, basix.LagrangeVariant.equispaced)
i_map = {}
j = 0
for dim, entity_points in enumerate(interpolation_points):
dim_values = []
for index, points in enumerate(entity_points):
dphi = cmap_element.tabulate(1, points)[1:, :, :, 0]
dim_values.append(dphi)
for p in range(points.shape[0]):
i_map[(dim, index, p)] = j
j += 1
basis_functions.append(dim_values)
# Verify map
ip = V.element.interpolation_points()
for key in i_map.keys():
assert np.allclose(
ip[i_map[key]], interpolation_points[key[0]][key[1]][key[2]])
gdim = mesh.geometry.dim
x = mesh.geometry.x[:, :gdim]
normals = basix.cell.facet_normals(b_cell)
def compute_interpolation_jacobian(mesh, cell, entity_dim, entity_index, point_index):
x_dofs = mesh.geometry.dofmap.links(cell)
cell_coords = x[x_dofs]
dphi_cell = basis_functions[entity_dim][entity_index][:, point_index, :]
J = np.dot(cell_coords.T, dphi_cell.T)
K = np.linalg.inv(J)
return K
def compute_normal(mesh, cell, facet_index):
"""
Push forward facet normal using covariant Piola transform
"""
num_points = interpolation_points[tdim-1][facet_index].shape[0]
n = np.zeros((num_points, mesh.geometry.dim))
for i in range(num_points):
K = compute_interpolation_jacobian(mesh, cell, tdim-1, facet_index, i)
n[i] = np.dot(K.T, normals[facet_index, :]) / \
np.linalg.norm(np.dot(K.T, normals[facet_index, :]))
return n
# Create unique array, (num_unique_cells * num_interpolations_points_per_cell, value_size*block_size)
for (cell, facet) in cell_facet_pairs:
cell_normals = compute_normal(mesh, cell, facet)
integration_points_on_facet = interpolation_points[tdim-1][facet]
num_points = integration_points_on_facet.shape(0)
# Now map each entry (tdim-1, facet, point_i) to the interpolation points for a given cell
# Build array for input to C++ interpolation into appropriate function from V to use for BC
# https://github.com/FEniCS/dolfinx/blob/main/python/dolfinx/fem/function.py#L374-L377
If you have any questions or comments, let me know.