Tabulate dof coordinates of N1curl space

Hi,

in dolfin, it was possible to tabulate the dof coordinates of a first- and second-order N1curl space. In dolfinx, running

with XDMFFile(comm, "testmesh.xdmf", "r") as xdmf:
    m = xdmf.read_mesh(name="mesh")
ele = element("N1curl", m.basix_cell(), 1)
V = dfx.fem.functionspace(m, ele)
V.tabulate_dof_coordinates()

results in the following error message:

File ~/miniconda3/envs/fenicsxpg/lib/python3.11/site-packages/dolfinx/fem/function.py:859 in tabulate_dof_coordinates
return self._cpp_object.tabulate_dof_coordinates() # type: ignore

RuntimeError: Cannot evaluate dof coordinates - this element does not have pointwise evaluation.

Is there any alternative/workaround? I thought about calculating the dof coordinates manually, starting with the geometric relations to a “Lagrange” space.

A N1curl space doesn’t really have dof coordinates, as each dof is associated with an entity, i.e. a facet, edge or cell.
What you can get is the coordinates of the interpolation points (which gives you multiple points for the integration across a facet, with


from mpi4py import MPI
import dolfinx
import ufl
import numpy as np
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 1, 1, cell_type=dolfinx.mesh.CellType.triangle)

V = dolfinx.fem.functionspace(mesh, ("N1curl", 1))

x = ufl.SpatialCoordinate(mesh)

int_coord_evaluator = dolfinx.fem.Expression(x, V.element.interpolation_points())


num_cells = mesh.topology.index_map(mesh.topology.dim).size_local
coords = int_coord_evaluator.eval(mesh, np.arange(num_cells, dtype=np.int32))

for i in range(num_cells):
    print(i, coords[i].reshape(-1, mesh.geometry.dim))

returning

0 [[1.  0.5]
 [0.5 0.5]
 [0.5 0. ]]
1 [[0.5 1. ]
 [0.5 0.5]
 [0.  0.5]]

For a higher order space (degree 2), you get

0 [[1.         0.21132487]
 [1.         0.78867513]
 [0.21132487 0.21132487]
 [0.78867513 0.78867513]
 [0.21132487 0.        ]
 [0.78867513 0.        ]
 [0.33333333 0.16666667]
 [0.83333333 0.66666667]
 [0.83333333 0.16666667]]
1 [[0.21132487 1.        ]
 [0.78867513 1.        ]
 [0.21132487 0.21132487]
 [0.78867513 0.78867513]
 [0.         0.21132487]
 [0.         0.78867513]
 [0.16666667 0.33333333]
 [0.66666667 0.83333333]
 [0.16666667 0.83333333]]

where more than one point is associated with each degree of freedom.
They are related through

V.element.basix_element.interpolation_matrix

which tells you how values at each point in the reference cell should be combined to compute the dual basis evaluation