ufl.log.UFLException: Cannot determine geometric dimension from expression

Dear Community

I am running into the following error

ufl.log.UFLException: Cannot determine geometric dimension from expression.

while trying to evaluate the last expression in the code snippet below, in which “M” is defined to be the integral over a subdomain. I have provided only minimal portions of the code below. I could share more details if required.

x = SpatialCoordinate(mesh)
strm = FiniteElement("CG",triangle,2)
vort = FiniteElement("CG",triangle,2)
W = FunctionSpace(mesh,strm*vort)

### This is followed by code segment in which the variational problem is set up

## Compute solution
w = Function(W)
solve(a == L, w, bcs)
## split using deep copy
(psi, omeg) = w.split(True)

w_eval =-1.0/(r*sin(theta))*(psi.dx(0).dx(0)+psi.dx(1).dx(1))+((x[0]*x[0])/(r*x[1]*x[1]))*psi.dx(1)-(x[0]/(x[1]*r*r))*psi.dx(1)

w_chk= ((x[0]*w_eval.dx(0)) + (x[1]*w_eval.dx(1)) - w_eval)
M = assemble(w_chk*ds(1))

Could you please help me out with this issue?

Thank You
Warm Regards

Please make sure that your code reproduces the error message when copy-pasting it into a new file.
You are missing several key definitions, such as a mesh (can be a built in one such as a UnitSquareMesh) and the variables theta and r. Note that you can remove the line solve(a==L, w, bcs), as the code is only dependent on the split of w, not the actual values.

For more information about how to easily get your questions answered, see: Read before posting: How do I get my question answered?

Thank you for the prompt response. Given below is the complete code with the mesh definitions and the definitions of “r” and “theta”, and executing it would reproduce the error mentioned above.

from dolfin import *
from mshr import *
import math
import numpy as np

Re = 5.
Ri = 1.
rect = Rectangle(Point(0., 0.), Point(Re, Re))
num_segments=200
domain = Circle(Point(0., 0.), Re, num_segments) - Circle(Point(0., 0.), Ri, num_segments)

#creating semi-circular domain
domain = (domain - Rectangle(Point(0., 0.), Point(Re, -Re))
                - Rectangle(Point(0., 0.), Point(-Re, -Re)))

mesh_resolution = 128
mesh = generate_mesh(domain, mesh_resolution)
x = SpatialCoordinate(mesh)

strm = FiniteElement("CG",triangle,2)
vort = FiniteElement("CG",triangle,2)
W = FunctionSpace(mesh,strm*vort)

r = Expression("sqrt(x[0]*x[0]+x[1]*x[1])", degree=2)
theta = Expression("atan2(x[1],x[0])", degree=2)
U = Constant(1.0)

## Constructing known solution for \psi and vorticity and applying on boundaries
psi_mfg = Expression("U*sin(theta)*sin(theta)*(0.5*r*r-0.75*r+(0.25/r))",r=r,
                      theta=theta,U=U,degree=2)

omeg_mfg = Expression("-1.5*U*sin(theta)/(r*r)",r=r,
                      theta=theta,U=U,degree=2)

### Defining boundary conditions
### and marking subdomains
class Inside(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1.5*Ri * (1 - np.cos(np.pi / num_segments))
        return near(sqrt(x[0]**2 + x[1]**2), Ri, tol) and on_boundary
class Outside(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1.5*Re * (1 - np.cos(np.pi / num_segments))
        return near(sqrt(x[0]**2 + x[1]**2), Re, tol) and on_boundary

sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
sub_domains.set_all(0)
Inside().mark(sub_domains, 1)
Outside().mark(sub_domains, 2)

## Assembling the BCs
inner_circ_psi  = DirichletBC(W.sub(0), psi_mfg, sub_domains, 1)
outer_circ_psi  = DirichletBC(W.sub(0), psi_mfg, sub_domains, 2)
inner_circ_omeg = DirichletBC(W.sub(1), omeg_mfg, sub_domains, 1)
outer_circ_omeg = DirichletBC(W.sub(1), omeg_mfg, sub_domains, 2)
bcs = [inner_circ_psi, outer_circ_psi, inner_circ_omeg, outer_circ_omeg]
### Defining test and trial functions
(psi, omeg) = TrialFunctions(W)
(u, v) = TestFunctions(W)
## Assembling bilinear equation
a1 = (inner(grad(omeg),grad(v))-(Constant(1.0)/(r*sin(theta)))*v*(omeg.dx(1))
      +(omeg*v/(r*r*sin(theta)*sin(theta))))
a2 = (inner(grad(psi),grad(u))+(Constant(1.0)/(r*sin(theta)))*u*(psi.dx(1))
     -(omeg*r*sin(theta)*u))
a = (a1+a2)
f = Constant(0.0)
g = Constant(0.0)
a = (a1+a2)*dx
L = f*v*dx+g*u*dx
## Compute solution
w = Function(W)
solve(a == L, w, bcs)
# split using deep copy
(psi, omeg) = w.split(True)
w_eval =-1.0/(r*sin(theta))*(psi.dx(0).dx(0)+psi.dx(1).dx(1))+((x[0]*x[0])/(r*x[1]*x[1]))*psi.dx(1)-(x[0]/(x[1]*r*r))*psi.dx(1)

w_chk= ((x[0]*w_eval.dx(0)) + (x[1]*w_eval.dx(1)) - w_eval)
M = assemble(w_chk*ds(1))

@dokken : I think the issue seems to arise when I attempt to calculate derivatives of “w_eval”, as done in the expression for “w_chk”. The following expression, which is based on “w_eval” and does not need its derivatives, evaluates just fine.

M = assemble(w_eval*sin(theta)*sin(theta)*ds(1))

Increasing the order of the finite element from 2 to 4 does not seem to solve the problem. I expected it would because “w_chk” implicitly evaluates the third order derivative of “psi”.

Just wanted to add this bit of information, in case it helps in tracing the root of the problem.

Thank You

You should use spatial coordinates instead of expressions, as you cannot really differentiate through the Expression class. Please note that the code below shows a minimal example, where one use a built in mesh, avoid the definition of code that is not required to produce the issue.

from dolfin import *
from ufl.classes import Atan2

mesh = UnitSquareMesh(10, 10)
x = SpatialCoordinate(mesh)

strm = FiniteElement("CG",triangle,2)
vort = FiniteElement("CG",triangle,2)
W = FunctionSpace(mesh,strm*vort)

r = sqrt(x[0]*x[0] + x[1]*x[1])
theta = Atan2(x[1], x[1])
U = Constant(1.0)

## Constructing known solution for \psi and vorticity and applying on boundaries
psi_mfg = U*sin(theta)*sin(theta)*(0.5*r*r-0.75*r+(0.25/r))

omeg_mfg = -1.5*U*sin(theta)/(r*r)

w = Function(W)
# split using deep copy
(psi, omeg) = w.split(True)
(psi, omeg) = w.split(True)
w_eval =-1.0/(r*sin(theta))*(psi.dx(0).dx(0)+psi.dx(1).dx(1))+((x[0]*x[0])/(r*x[1]*x[1]))*psi.dx(1)-(x[0]/(x[1]*r*r))*psi.dx(1)

w_chk= ((x[0]*w_eval.dx(0)) + (x[1]*w_eval.dx(1)) - w_eval)
M = assemble(w_chk*ds(1))
1 Like

Thank you for describing, with an example, the construction of a minimal example. It is quite helpful!

Warm Regards