Expressions with different shapes - Mixed Spaces

[Kaa]

Hello all,

I am following a finite element course (they use FEniCS) “High Performance Finite Element Modeling”. Some of the codes are out dated:

# FEM function spaces and functions
V = VectorFunctionSpace(mesh, "CG", 1); 
Q = FunctionSpace(mesh, "CG", 1);
W = V * Q;   # outdated part
h = CellSize(mesh);
(v, q) = TestFunctions(W);
w = Function(W); 
(u, p) = (as_vector((w[0], w[1])), w[2]); 
u0 = Function(V)

I wrote elements and function spaces in a different way, now I have another problem UFLException: Can't add expressions with different shapes. It originates from elements or function spaces. I spent some time on it -i think it has a simple solution but couldn’t find something similar, and asked on course’s discussion board has no replies- so I couldn’t figure it out. I have the following code:

# Generate mesh
mresolution=30
mesh = generate_mesh(Rectangle(Point(G[0], G[2]), Point(G[1], G[3])) - Circle(Point(.5, .5), .1), mresolution)     

# Elements
el = FiniteElement("CG", mesh.ufl_cell(), degree=1);     
vel = VectorElement("CG", mesh.ufl_cell(), degree=1);     
element = MixedElement([vel, el]);     

# Function spaces
V = VectorFunctionSpace(mesh, vel, degree=1);       
W = FunctionSpace(mesh, element);      
h = 2*Circumradius(mesh);
(v, q) = TestFunctions(W);
w = Function(W);      
(u, p) = (as_vector((w[0], w[1])), w[2]);      
u0 = Function(V)      

The rest is the same w/ the original example. At the end I have this error

--> 109     um = theta*u + (1.0-theta)*u0

UFLException: Can’t add expressions with different shapes.

Since theata = 0.5, I think there is a problem with u and u0.

What I did and failed?

I thought since we have this

 (u, p) = (as_vector((w[0], w[1])), w[2]);  

u have (w[0], w[1]) . So the dimensions of u0 (degree of V is 1) should match with the dimension of u (which should be 2 since we have w[0] and w[1]) and I thought degree of V should be 2 instead of 1 but it didn’t work.

So I have 2 questions:

1. How can I solve this problem?
2. What is happening in (u, p) = (as_vector((w[0], w[1])), w[2]);

Thanks for following along !

Note1: I don’t know much about FEniCS, so I wasn’t able to post a minimal working example from the existed code.
Note2: I am new to finite elements and coding so I might not be using the correct terminology.

[volkerk]

I am not sure what you are trying to do, and you do not provide a proper minimal working example, but I guess the problem is that you define the function space V as a VectorFunctionSpace of a VectorElement, which means that your function u0 is a matrix valued function. Try

V = FunctionSpace(mesh, vel)

and see if this gets you what you want.

[Kaa]

@volkerk

V = FunctionSpace(mesh, vel)

Hi,
Thank you for your reply. I tried but it asked for degree TypeError: VectorFunctionSpace() missing 1 required positional argument: 'degree'

This is the whole code:

import matplotlib.pyplot as plt;
from dolfin import *; 
from mshr import *; 
from IPython.display import display, clear_output; 
import time
import logging; 
logging.getLogger('FFC').setLevel(logging.WARNING)

# Compact plot utility function
def plot_compact(u, t, stepcounter): 
    if stepcounter % 5 == 0:
        uEuclidnorm = project(sqrt(inner(u, u)), Q); 
        ax.cla(); 
        fig = plt.gcf();
        fig.set_size_inches(16, 2)
       
       	# Plot norm of velocity
        plt.subplot(1, 2, 1);
        mplot_function(uEuclidnorm); 
        plt.title("Velocity")
        
        if t == 0.: 
           plt.colorbar(); 
           plt.axis(G)
        
        plt.subplot(1, 2, 2);
        
        if t == 0.:
           plt.triplot(mesh2triang(mesh)); 
           plt.title("Mesh") # Plot mesh
    display(pl)

    plt.suptitle("Navier-Stokes t: %f" % (t));
    plt.tight_layout();
    clear_output(wait=True); 

############## Important code starts here ##############

# Generate domain 
XMIN = 0.; 
XMAX = 4.; 
YMIN = 0; 
YMAX = 1.; 
G = [XMIN, XMAX, YMIN, YMAX]; 
eps = 1e-5

# Generate mesh
mresolution=30
mesh = generate_mesh(Rectangle(Point(G[0], G[2]), Point(G[1], G[3])) - Circle(Point(.5, .5), .1), mresolution)      ##############################

# Elements
el = FiniteElement("CG", mesh.ufl_cell(), degree=1);     
vel = VectorElement("CG", mesh.ufl_cell(), degree=1);     
element = MixedElement([vel, el]);     

# Function spaces
V = VectorFunctionSpace(mesh, vel, degree=1);      #related w/ u0 
W = FunctionSpace(mesh, element);                  #related w/ u
h = 2*Circumradius(mesh);
(v, q) = TestFunctions(W);
w = Function(W);      
(u, p) = (as_vector((w[0], w[1])), w[2]);          # u
u0 = Function(V)                                   # u0 

# Inlet velocity
uin = Expression(("4*(x[1]*(YMAX-x[1]))/(YMAX*YMAX)", "0."), YMAX=YMAX, degree=1) 

# Mark regions for boundary conditions
om = Expression("x[0] > XMAX - eps ? 1. : 0.", XMAX=XMAX, eps=eps, degree=1) 
im = Expression("x[0] < XMIN + eps ? 1. : 0.", XMIN=XMIN, eps=eps, degree=1)
nm = Expression("x[0] > XMIN + eps && x[0] < XMAX - eps ? 1. : 0.", XMIN=XMIN, XMAX=XMAX, eps=eps, degree=1)

# Timestep, viscosity and stabilization parameters
k = 0.1; 
nu = 1e-6; 
d = .2*h**(3./2.) 

# Time interval and penalty parameter
t, T = 0., 10.; 
gamma = 10*1./h 

pl, ax = plt.subplots(); 

# Initialize time stepping
stepcounter = 0; 
timer0 = time.clock()

# Time stepping method
theta = 0.5 # 0.5 - Midpoint rule, 1.0 - Implicit Euler, 0.0 - Explicit Euler

# Time-stepping loop
while t < T: 
    # Weak residual of stabilized FEM for Navier-Stokes eq.
    um = theta*u + (1.0-theta)*u0               #####ERROR#####
    
    # Navier-Stokes equations in weak residual form
    r = ((inner((u - u0)/k + grad(p) + grad(um)*um, v) + nu*inner(grad(um), grad(v)) + div(um)*q)*dx +
        gamma*(om*p*q + im*inner(u - uin, v) + nm*inner(u, v))*ds + # Weak boundary conditions
        d*(inner(grad(p) + grad(um)*um, grad(q) + grad(um)*v) + inner(div(um), div(v)))*dx) # Stabilization
    
    # Solve the Navier-Stokes PDE (one timestep)
    solve(r==0, w)  

    # Visualize the solution
    plot_compact(u, t, stepcounter) # Plot all quantities (see implementation above)
    
    # Shift to next timestep
    t += k; 
    stepcounter += 1; 
    u0 = project(u, V); 

print(f"elapsed CPU time: ", (time.clock() - timer0))  

[volkerk]

As I told you in my last answer, you should use FunctionSpace instead of VectorFunctionSpace for your space V. With that your code does not produce this error message.

[Kaa]

Ok my bad thanks for pointing out.