Crash when trying to obtain coordinates of quadrature points

Errors
2

You cannot get the dof coordinates of a quadrature space in that way.
This is because they are not “well defined” finite elements, they do not have basis functions defined in the whole domain.

Consider the following code:

import numpy as np
import basix
from ffcx.element_interface import create_quadrature
from IPython import embed
import ufl
from dolfinx import fem, mesh

from mpi4py import MPI

msh = mesh.create_rectangle(comm=MPI.COMM_WORLD,
                            points=((0.0, 0.0), (1.0, 1.0)), n=(32, 16),
                            cell_type=mesh.CellType.triangle,)

QF_ele = ufl.FiniteElement("Quadrature", ufl.triangle,
                           degree=2, quad_scheme="default")
QF = fem.FunctionSpace(msh, QF_ele)
x = ufl.SpatialCoordinate(msh)
points, _ = create_quadrature(
    QF_ele.cell().cellname(), QF_ele.degree(), QF_ele.quadrature_scheme())
x_expr = fem.Expression(x, points)
num_cells = msh.topology.index_map(msh.topology.dim).size_local
vals = x_expr.eval(np.arange(num_cells, dtype=np.int32))
vals = vals.reshape(num_cells, points.shape[0], points.shape[1])
dofmap = QF.dofmap.list

dof_coordinates = np.zeros(
    (num_cells*points.shape[0], points.shape[1]), dtype=np.float64)
for i in range(num_cells):
    for j, dof in enumerate(dofmap.links(i)):
        print(f"Cell {i}, dof {dof}, coordinate {vals[i,j]}")
        dof_coordinates[dof] = vals[i, j]

print(dof_coordinates)
1 Like