Why is this not a valid UFL expression?

General
1

I am currently trying to port Automatic calibration of damping layers in finite element time domain simulations , written by Steven Vandekerckhove and Garth Wells, from FEniCS to FEniCSx and I am struggling.

To limit my question, I’d love to know why qhat in the MWE below is not a valid UFL expression. It creates a numpy array which I thought was a valid UFL expression.

from dolfinx import fem, mesh
from mpi4py import MPI
from basix.ufl import element, mixed_element
from ufl import (
    FacetNormal,
    TestFunctions,
    TrialFunctions,
    as_vector,
    as_matrix,
    avg,
    dx,
    grad,
    inner,
    lhs,
    rhs,
)
import numpy as np

M = 6
domain = mesh.create_unit_square(MPI.COMM_WORLD, M, M)

DG  = fem.functionspace(domain, ("DG", 0))          # this is for the CML layers
U   = element("DG", domain.basix_cell(), 2, shape=(2, )) # velocity
V   = element("DG", domain.basix_cell(), 2, shape=(3, )) # velocity
X = fem.functionspace(domain, mixed_element([U,V]))

# Set up test and trial functions (and vectors)
m1, n1 = TestFunctions(X)
v1, f1 = TrialFunctions(X) # velocity, deformation

# set material parameters
c11 = 85.e9 					# [Pa]
c12, c22 = 25.e9, 85.e9 		# [Pa]
c16, c26, c66 = 0, 0, 30.e9     # [Pa]
rho = 2500                      # [kg/m^3]

# Wave parameters
C  = np.sqrt(c11/rho)
C2 = np.sqrt(c66/rho)

q = as_vector((v1[0], v1[1], f1[0], f1[1], f1[2]))
l = as_vector((m1[0], m1[1], n1[0], n1[1], n1[2]))

Ax = as_matrix([[0, 0, -c11, -c12, -c16],
                [0, 0, -c16, -c26, -c66],
                [-1./rho, 0, 0, 0, 0],
                [0, 0, 0, 0, 0],
                [0, -1./rho, 0, 0, 0]])

Ay = as_matrix([[0, 0, -c16, -c26, -c66],
                [0, 0, -c12, -c22, -c26],
                [0, 0, 0, 0, 0],
                [0, -1./rho, 0, 0, 0],
                [-1./rho, 0, 0, 0, 0]])

G = as_matrix([[0., 0, 0, 0, 0],
               [0, 0., 0, 0, 0],
               [0, 0, 1., 0, 0],
               [0, 0, 0, 1., 0],
               [0, 0, 0, 0, 1.]])
dt = 4.e-8          # Time step size
n  = FacetNormal(domain) # Facet normals
sx = 0.1
q0 = fem.Function(X)

qold = (q + q0)/2.
qdot = (q - q0)/dt

#Provide difference function
def dif(l):
    return l("-")-l("+")

# Define fluxes on interior and exterior facets
qhat    = n[0]("-")*avg(Ax*qold) + n[1]("-")*avg(Ay*qold) + C*0.5*dif(qold)
qhatbnd = (n[0]*Ax + n[1]*Ay)*(G*qold)

F = inner(qdot, l)*dx \
  - (inner(Ax*qold, l.dx(0)) + inner(Ay*qold, l.dx(1))-inner(sx*qold, l))*dx \
  + inner(qhat, dif(l))*dS \
  + inner(qhatbnd, l)*ds\
  - inner(source, l)*dx

This gives me the following error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[3], line 80
     75 qhat    = n[0]("-")*avg(Ax*qold) + n[1]("-")*avg(Ay*qold) + C*0.5*dif(qold)
     76 qhatbnd = (n[0]*Ax + n[1]*Ay)*(G*qold)
     78 F = inner(qdot, l)*dx \
     79   - (inner(Ax*qold, l.dx(0)) + inner(Ay*qold, l.dx(1))-inner(sx*qold, l))*dx \
---> 80   + inner(qhat, dif(l))*dS \
     81   + inner(qhatbnd, l)*ds\
     82   - inner(source, l)*dx

File /dolfinx-env/lib/python3.12/site-packages/ufl/operators.py:209, in inner(a, b)
    204 def inner(a, b):
    205     """Take the inner product of a and b.
    206 
    207     The complex conjugate of the second argument is taken.
    208     """
--> 209     a = as_ufl(a)
    210     b = as_ufl(b)
    211     if a.ufl_shape == () and b.ufl_shape == ():

File /dolfinx-env/lib/python3.12/site-packages/ufl/constantvalue.py:517, in as_ufl(expression)
    515     return ComplexValue(expression)
    516 else:
--> 517     raise ValueError(
    518         f"Invalid type conversion: {expression} can not be converted to any UFL type."
    519     )

ValueError: Invalid type conversion: [Sum(Indexed(Sum(ComponentTensor(Product(Indexed(ComponentTensor(Product(FloatValue(0.5), ... MESSAGE TRUNCATED BY USER

I would really apprecaite it if someone could help me on this. My guess is that the problem is in my formulation of the mixed element space. The correct way to do this seems to have changed a lot in recent years, and I’m having trouble figuring out what to do here.

Any help is greatly apprecaiated. I’m approching a deadline (PhD student) and I’m exhausting my resources trying to figure this out. Thank you!

VERSION INFORMATION

either tutorial compatable docker images

2

Looks like numpy is trying to do some sneaky type conversion in its overloading of the __rmul__ operator. Try the following where I explicitly cast C to a float.

qhat    = n[0]("-")*avg(Ax*qold) + n[1]("-")*avg(Ay*qold) + float(C)*0.5*dif(qold)
1 Like
3

Thank you @nate. That allows me to keep going.

I’m glad I reached out because that would have been a very difficult bug for me to find. Could you explain a bit about your thought process for debugging that so I learn a bit?

Thanks again.

4

See, e.g., 3. Data model — Python 3.13.1 documentation

These functions are how the operators *, /, +, -, etc., are implemented. Somewhere numpy is implementing one of these operators, converting one of the arguments which is a UFL object to a numpy array, and then giving you that result. This is clearly not what you want. It’s best to stick with UFL functions/classes where possible for defining finite element formulations, and explicitly converting numpy types to what you need accordingly.

1 Like
5

Thank you very much for the help.