UFLException: Cannot determine geometric dimension from expression

variational formulation
1

Hi community,

I am struggling to figure the error I get when I try to do ufl.lhs(F) on my weak form. Apparently, I am not writing it in the right UFL way. In case it helps, I am trying to do 2D electromagnetic wave propagation on complex dolfinx 0.6.0. Following is an MWE that reproduces the error. Thanks a lot for your help!

import dolfinx, petsc4py, gmsh, ufl, mpi4py, dolfinx_mpc
from dolfinx.io import gmshio
import numpy as np

#%% Input parameters

fov_x = 0.5
fov_y = 1
wlen = 0.25
U0=1
theta = np.pi/4*1

#%% geometry
def create_geom(fov_x, fov_y):
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0) # disable output messages
    gmsh.clear()
    gdim = 2
    model_rank = 0

    occ = gmsh.model.occ

    fov = occ.add_rectangle(-fov_x/2, -fov_y/2, 0, fov_x, fov_y)
    sub = occ.add_rectangle(-fov_x/2, -fov_y/2, 0, fov_x, fov_y/2)
    # fov = occ.addBox()
    # sub = occ.addBox(-fov_x/2, -fov_y/2, 0, fov_x, fov_y/2, 0)
    occ.synchronize()
    dom = occ.fragment([(gdim, fov)], [(gdim, sub)])
    occ.synchronize()

    # number all domains
    all_doms = gmsh.model.getEntities(gdim)
    for j, dom in enumerate(all_doms):
        gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1)  # create the main group/node

    # number all boundaries
    all_edges = gmsh.model.getEntities(gdim - 1)
    for j, edge in enumerate(all_edges):
        gmsh.model.addPhysicalGroup(edge[0], [edge[1]], j + 1)  # create the main group/node

    gmsh.model.mesh.generate(gdim)
    model_rank = 0
    mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, model_rank, gdim)
    gmsh.finalize()
    return mesh, ct, ft

mesh, ct, ft = create_geom(fov_x, fov_y)

degree = 2
curl_el = ufl.FiniteElement("N1curl", mesh.ufl_cell(), degree)
Omega = dolfinx.fem.FunctionSpace(mesh, curl_el)
U = ufl.TrialFunction(Omega)
Ut = ufl.TestFunction(Omega)

# wavevector in incidence and transmission media
k0 = 2*np.pi/wlen
kxi = k0*np.sin(theta)

kb_vec = ufl.as_vector((kxi, 0, 0)) # vector representing FB phase
U_3d = ufl.as_vector((U[0], U[1], 0))
Ut_3d = ufl.as_vector((Ut[0], Ut[1], 0))

F = -k0**2*ufl.inner(U, Ut)*ufl.dx # this is ok
F += ufl.inner(ufl.cross(1j*kb_vec, (ufl.curl(U_3d) + ufl.cross(1j*kb_vec, U_3d))), Ut_3d)*ufl.dx # this is problematic
F += ufl.inner((ufl.curl(U_3d)+ufl.cross(1j*kb_vec, U_3d)), ufl.curl(Ut_3d))*ufl.dx # this is also problematic

a = ufl.lhs(F) # Error!
L = ufl.rhs(F)
2

I figured the problem. The ufl.curl does not like to take ufl.as_vector() fields.

The following code fixes it

# Model an inclined plane wave using periodic boundary conditions. First we consider p polarization i.e. vector case
# MWE prepared for seeking help from forum
# FIXED

import dolfinx, petsc4py, gmsh, ufl, mpi4py, dolfinx_mpc
from dolfinx.io import gmshio
import numpy as np

#%% Input parameters

fov_x = 0.5
fov_y = 1
wlen = 0.25
U0=1
theta = np.pi/4*1

#%% geometry
def create_geom(fov_x, fov_y):
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0) # disable output messages
    gmsh.clear()
    gdim = 2
    model_rank = 0

    occ = gmsh.model.occ

    fov = occ.add_rectangle(-fov_x/2, -fov_y/2, 0, fov_x, fov_y)
    sub = occ.add_rectangle(-fov_x/2, -fov_y/2, 0, fov_x, fov_y/2)
    # fov = occ.addBox()
    # sub = occ.addBox(-fov_x/2, -fov_y/2, 0, fov_x, fov_y/2, 0)
    occ.synchronize()
    dom = occ.fragment([(gdim, fov)], [(gdim, sub)])
    occ.synchronize()

    # number all domains
    all_doms = gmsh.model.getEntities(gdim)
    for j, dom in enumerate(all_doms):
        gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1)  # create the main group/node

    # number all boundaries
    all_edges = gmsh.model.getEntities(gdim - 1)
    for j, edge in enumerate(all_edges):
        gmsh.model.addPhysicalGroup(edge[0], [edge[1]], j + 1)  # create the main group/node

    gmsh.model.mesh.generate(gdim)
    model_rank = 0
    mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, model_rank, gdim)
    gmsh.finalize()
    return mesh, ct, ft

mesh, ct, ft = create_geom(fov_x, fov_y)

degree = 2
curl_el = ufl.FiniteElement("N1curl", mesh.ufl_cell(), degree)
Omega = dolfinx.fem.FunctionSpace(mesh, curl_el)
U = ufl.TrialFunction(Omega)
Ut = ufl.TestFunction(Omega)

# wavevector in incidence and transmission media
k0 = 2*np.pi/wlen
kxi = k0*np.sin(theta)

kb_vec = ufl.as_vector((kxi, 0, 0)) # vector representing FB phase
U_3d = ufl.as_vector((U[0], U[1], 0))
Ut_3d = ufl.as_vector((Ut[0], Ut[1], 0))

def curl_2d(u):
    '''calculate 3 component curl of two component vector field u'''
    return ufl.as_vector((0,0,u[1].dx(0)-u[0].dx(1)))

F = -k0**2*ufl.inner(U, Ut)*ufl.dx # this is ok
F += ufl.inner(curl_2d(U), curl_2d(Ut))*ufl.dx # this is also problematic

a = ufl.lhs(F)
L = ufl.rhs(F)