Ffcx error using quadrature elements


I’m trying to update some code from version 0.6 to 0.7.2. I’ve implemented the obvious changes between the versions, but ffcx still refuses to compile the form. The error is when I sum an integration over the domain and the facets. The 2 parts of the form compile correctly when taken in isolation.

here is the code stripped of all the unnecessary :

import basix
from ufl import (Mesh,
                 SpatialCoordinate,as_matrix,as_vector,as_tensor,inner,inv,sqrt,cross, action,sym)
from basix.ufl import element, quadrature_element

cell = "hexahedron"
# Coordinate Element for Mesh
elementV = element("Lagrange", cell, 1, shape=(3,))  
mesh = Mesh(elementV) 
normal = FacetNormal(mesh)

# Lagrangian spaces
V = FunctionSpace(mesh, elementV)

# Quadrature space for stress
VQV_ele = quadrature_element(cell, degree=2, scheme="default", value_shape=(6,))
VQV = FunctionSpace(mesh, VQV_ele)

# Get points and weights for custom integration measure
points, weights = basix.make_quadrature( basix.cell.string_to_type(cell), 2, basix.quadrature.string_to_type("default"))
metadata = {"quadrature_rule": "custom", "quadrature_points": points, "quadrature_weights": weights}

dx = Measure("dx", domain=mesh, metadata=metadata)

# Define Functions for linear problem

v = TestFunction(V)

# Define quadrature functions for the weak form
q_sigma = Coefficient(VQV)

# strain measure function (used for strain increment in weak form)
def eps(v):
        e = sym(grad(v))
        return as_vector([e[0, 0], e[1, 1], e[2, 2], 2 * e[1, 2], 2 * e[0, 2], 2 * e[0, 1]] )

# Load application on single facet
Pmax = Constant(mesh) #Max load

# Variational Form
L_form = -inner(eps(v), q_sigma) * dx + Pmax * dot(normal, v) * ds
forms = [L_form]

here is the error:

root@bbb62a0311a2:~# ffcx test.py
WARNING:py.warnings:/usr/local/lib/python3.10/dist-packages/ufl/utils/sorting.py:84: UserWarning: Applying str() to a metadata value of type ndarray, don't know if this is safe.
  warnings.warn(f"Applying str() to a metadata value of type {type(value).__name__}, "

Traceback (most recent call last):
  File "/usr/local/bin/ffcx", line 8, in <module>
  File "/usr/local/lib/python3.10/dist-packages/ffcx/main.py", line 67, in main
    code_h, code_c = compiler.compile_ufl_objects(
  File "/usr/local/lib/python3.10/dist-packages/ffcx/compiler.py", line 102, in compile_ufl_objects
    ir = compute_ir(analysis, object_names, prefix, options, visualise)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/ir/representation.py", line 197, in compute_ir
    irs = [_compute_integral_ir(fd, i, analysis.element_numbers, integral_names, finite_element_names,
  File "/usr/local/lib/python3.10/dist-packages/ffcx/ir/representation.py", line 197, in <listcomp>
    irs = [_compute_integral_ir(fd, i, analysis.element_numbers, integral_names, finite_element_names,
  File "/usr/local/lib/python3.10/dist-packages/ffcx/ir/representation.py", line 492, in _compute_integral_ir
    integral_ir = compute_integral_ir(itg_data.domain.ufl_cell(), itg_data.integral_type,
  File "/usr/local/lib/python3.10/dist-packages/ffcx/ir/integral.py", line 84, in compute_integral_ir
    mt_table_reference = build_optimized_tables(
  File "/usr/local/lib/python3.10/dist-packages/ffcx/ir/elementtables.py", line 361, in build_optimized_tables
    t = get_ffcx_table_values(quadrature_rule.points, cell,
  File "/usr/local/lib/python3.10/dist-packages/ffcx/ir/elementtables.py", line 124, in get_ffcx_table_values
    entity_points = map_integral_points(points, integral_type, cell, entity)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/ir/representationutils.py", line 92, in map_integral_points
    assert points.shape[1] == tdim - 1

can anybody help?
I’m using the latest stable version [dolfinx 0.7.2; ffcx 0.7.0]

See: Create a DOLFINx form for a quadrature element - #3 by mscroggs

This might have been fixed in Check for custom quadrature for each integral, not each form by mscroggs · Pull Request #622 · FEniCS/ffcx · GitHub. Can you try with a docker image dolfinx/dolfinx:nightly?

yes, the nightly version seems to have fixed the issue.

Thanks both for the help.