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

Errors
1

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

1 Like
2

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?

3

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))
4

@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

5

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
6

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

Warm Regards

7

Hello. I get the same ERROR:

Traceback (most recent call last):
File “cube_non_New.py”, line 148, in
a = K*dot(sigma(u_n) * grad(u), grad(v)) * dx
File “cube_non_New.py”, line 95, in sigma
E = -grad(u_value)
File “/usr/lib/python2.7/dist-packages/ufl/operators.py”, line 322, in grad
return Grad(f)
File “/usr/lib/python2.7/dist-packages/ufl/differentiation.py”, line 124, in new
dim = f.geometric_dimension()
File “/usr/lib/python2.7/dist-packages/ufl/core/expr.py”, line 304, in geometric_dimension
error(“Cannot get geometric dimension from an expression with no domains!”)
File “/usr/lib/python2.7/dist-packages/ufl/log.py”, line 151, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Cannot get geometric dimension from an expression with no domains!
this is my code:

from dolfin import *
import numpy as np
#I just make a simple cube without any domain by UnitCubeMesh(6,6,6)
mesh = Mesh("cube_mesh.xml") 
# Define the conductivity function
sigma0 = 5.5E+7
alpha = 0.5

def sigma(u):
    # Define the function space for sigma
    W = FunctionSpace(mesh, "CG", 1)
    
    # Define the function sigma
    #sigma = Function(W)
    sigma_values = np.zeros(W.dim())
    
    mx=np.loadtxt("data6/Mx_0002.txt").T
    my=np.loadtxt("data6/My_0002.txt").T
    mz=np.loadtxt("data6/Mz_0002.txt").T

    d2v = dof_to_vertex_map(V)
    Mx = Function(V); Mx.vector().set_local(mx[d2v])
    My = Function(V); My.vector().set_local(my[d2v])
    Mz = Function(V); Mz.vector().set_local(mz[d2v])
    
    u_function = interpolate(u,W)
           
    # Compute sigma at each vertex of the mesh
    for vertex in vertices(mesh):
        # Get the coordinates of the vertex
        x = vertex.point().x()
        y = vertex.point().y()
        z = vertex.point().z()
            
        # Evaluate the function u at the vertex
        u_value = u_function(x, y, z)
            
        # Calculate E (gradient of u)
        E = -grad(u_value)
            
        # Define M as a vector at the vertex
        M = as_vector((Mx(x,y,z), My(x,y,z), Mz(x,y,z)))
            
        # Calculate dot products
        dot_ME = dot(M, E)
        dot_EE = dot(E, E)
        dot_MM = dot(M, M)
            
        # Calculate cos(theta)
        costheta = dot_ME / sqrt(dot_EE * dot_MM)
            
        # Evaluate sigma_vertex to obtain a scalar value
        sigma_vertex_value = sigma_vertex.evaluate()
        
        # Calculate sigma at the vertex
        sigma_vertex_scalar = sigma0 / (1 + alpha * costheta**2) * sigma_vertex_value
    
        # Assign sigma value to the function sigma
        sigma_values = np.zeros(W.dim())
        
        for vertex in vertices(mesh):
            # Set the scalar value of sigma_vertex for each vertex
            sigma_values[vertex.index()] = sigma_vertex_scalar
    
        # Create a Function object to represent the sigma values
    sigma_function = Function(W)
    sigma_function.vector().set_local(sigma_values)
    
    sigma.interpolate(sigma_function)
       
    return sigma

as you can see I take the information about vector M from external files, In Mx, My, Mz you just can see information about coordinate in each vertex. I want to calculate the sigma at each vertex of the mesh depending on the vector M, since logically the angle will affect the sigma at each vertex. Could you please give any suggestion?
In my perception I have to create FunctionSpace for it, Interpolate it and Evaluate the Function.

8

Please format your code with 3xand ensure that it is reproducible, ie. use a built in mesh, and artificial data forMx,My, Mz`.

```python
# add code here
```
9

I submit everything here

10

I would strongly encourage you to try to make the example independent of Google drive.
Reasons for this is:

  1. permissions to files might change over the years, meaning that the post loses its value to the community
  2. storage is not persistent, you might remove the files at any time, with the same loss to the community.
  3. by not making the effort to reduce the problem to a minimal reproducible example, you place the efforts of generalizing your problem to the community.
11

Got it. Files TXT just has info about coordinates, I will submit it on GitHub then. The minimal reproducible example is solving nonlinear equation using my sigma function. This is my code:

from dolfin import *
from fenics import *
import numpy as np

# Create/download mesh and define function space
mesh = UnitCubeMesh(3, 3, 3)
plot(mesh, interactive=True)
file_mesh = File("cube_mesh.xml")
file_mesh << mesh
V = FunctionSpace(mesh, "CG", 1)
VV = VectorFunctionSpace(mesh, "CG", 1)

# Define boundary conditions.
tol = 1e-14
def Dirichlet_boundary0(x, on_boundary):
    return 0-tol<=x[0]<=tol
def Dirichlet_boundary1(x, on_boundary):
    return 1-tol<=x[0]<=1+tol
bc0 = DirichletBC(V, Constant(0), Dirichlet_boundary0)
bc1 = DirichletBC(V, Constant(10), Dirichlet_boundary1)
bc = [bc0, bc1]

# Define variational problem 
u = TrialFunction(V)
v = TestFunction(V)
u_n = interpolate(Expression('1.0 - 2*(x[0])', degree=1), V) # previous known u

sigma0 = 5.5E+7
alpha = 0.5

# define sigma
def sigma(u):
    # Define the function space for sigma
    W = FunctionSpace(mesh, "CG", 1)
    sigma = Function(W)
    # Define the function sigma
    sigma_values = np.zeros(W.dim())
    
     # Load Mx, My, and Mz from files, TXT files just contain info about coordinates of vector in vertex
     mx = np.loadtxt("data6/Mx_0002.txt").T
     my = np.loadtxt("data6/My_0002.txt").T
     mz = np.loadtxt("data6/Mz_0002.txt").T
    
    
    # Assign loaded values to Function objects (retern the same level to loop with mx)
     d2v = dof_to_vertex_map(V)
     Mx = Function(V); Mx.vector().set_local(mx[d2v])
     My = Function(V); My.vector().set_local(my[d2v])
     Mz = Function(V); Mz.vector().set_local(mz[d2v])
 
    u_function = project(u,W)
    geometric_dim = u_function.geometric_dimension()
    print("Geometric dimension of the function:", geometric_dim)       
    sigma_vertex = Function(W)
    # Compute sigma at each vertex of the mesh
    for vertex in vertices(mesh):
        # Get the coordinates of the vertex
        x = vertex.point().x()
        y = vertex.point().y()
        z = vertex.point().z()
        u = Function(W)
        # Define M as a vector at the vertex
        M = as_vector((Mx(x,y,z), My(x,y,z), Mz(x,y,z)))
        E = project( -grad(u), VectorFunctionSpace(mesh, 'CG', 1))   
        # Calculate dot products
        dot_ME = dot(M, E)
        dot_EE = dot(E, E)
        dot_MM = dot(M, M)
            
        # Calculate cos(theta)
        costheta = dot_ME / sqrt(dot_EE * dot_MM)
            
        # Evaluate sigma_vertex to obtain a scalar value
        sigma_vertex_value = sigma_vertex(x, y, z)
        
        # Calculate sigma at the vertex
        sigma_vertex_scalar = (sigma0 / (1 + alpha * costheta**2)) * sigma_vertex_value
        
        # Assign sigma value to the function sigma
        sigma_values = np.zeros(W.dim())
        
        for vertex in vertices(mesh):
            # Set the scalar value of sigma_vertex for each vertex
            sigma_values[vertex.index()] = sigma_vertex_scalar
    
        # Create a Function object to represent the sigma values
    sigma_function = Function(W)
    sigma_function.vector().set_local(sigma_values)
    print("Function set local done:", sigma_function)
    
    sigma.interpolate(sigma_function) 
    return sigma

# PICARD ITERATION
a = dot(sigma(u_n) * grad(u), grad(v)) * dx # doesn't work when sigma(u)
f = Constant(0.0)
L = f * v * dx

# Picard iteration
u = Function(V)
eps = 1.0
tol = 1.0e-5
iter = 0
maxiter = 500
while eps > tol and iter < maxiter:
    iter += 1
    solve(a == L, u, bc) # for Picard
    du = u.vector().array() - u_n.vector().array()
    eps = np.linalg.norm(du, ord=np.Inf)
    print('iter=%d: norm=%g' % (iter, eps))
    u_n.assign(u)

# Plot results
plot(u, interactive=True)
plot(grad(u), interactive=True)`
12

actually after my change (so code above doesn’t show ‘Cannot determine geometric dimension from expression’), but i get this error:
*** Error: Unable to solve linear system using PETSc Krylov solver. *** Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_NANORINF, residual norm ||r|| = -nan).
without Sigma this code works perfect. So I’m sure that made a mistake while write Sigma function. Could you give me some idea how I can do it?

13

Your Picard iteration is certainly wrong. The bilinear form a does not get updated at every iteration, it will just use whatever value u_n had a the very first iteration.

A simple workaround is to define a inside the for loop. A more appropriate workaround would be to define sigma_function outside of the sigma function, and update it at every picard iteration.

It’s unlikely anyone will be able to actually run your code, because it uses an ancient version of dolfin. With that being said, if I were you I’d check if the values in the text file result in sigma always positive: if not, the linearized problem may not even be coercive.