Crash when trying to obtain coordinates of quadrature points

Errors
1

I am trying to obtain the coordinates of each quadrature points using tabulate_dofs_coordinates however I get the following error:

[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see https://petsc.org/release/faq/#valgrind and https://petsc.org/release/faq/
[0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run 
[0]PETSC ERROR: to get more information on the crash.
[0]PETSC ERROR: Run with -malloc_debug to check if memory corruption is causing the crash.
Abort(59) on node 0 (rank 0 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0

Here is a minimal example of what I am trying to do.

import ufl
from dolfinx import fem, mesh

from mpi4py import MPI

msh = mesh.create_rectangle(comm=MPI.COMM_WORLD,
                            points=((0.0, 0.0), (2.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)

coords = QF.tabulate_dof_coordinates()
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