This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

How Can One Create Test Functions for a MixedElement?

+1 vote

I am trying to create a fully-implicit solver for a system of two PDEs. I am using a mixed element consisting of a vector element for displacement and a scalar-valued element for pressure. I have been looking at the examples in the FEniCS tutorial (the Navier-Stokes example and the advection-reaction diffusion equations), but when I run my code, I get the following error:

Can only integrate scalar expressions. The integrand is a tensor
expression with value shape (2,) and free indices with labels ().

The relevant code:

from fenics import *
import numpy as np

# Parameters
c = 1.25
l = 5.0 
P_surf = 0.0 
sigma_0 = 3.0 
gamma = 0.65 
P0 = sigma_0*gamma 
Kv = 10.0 
cm = 5.0 
# Derivative of pressure at the bottom surface.
der_p_bottom = 0.0
# Derivative of displacement at the top surface.
der_w_top = cm*gamma*sigma_0

# The initial time.
t = 0.0
# The final time
T = 2.0
# The number of time steps to take.
num_steps = 40
# The step size
dt = (T-t)/num_steps
# Width
width = 3.0
# Thickness
thickness = 2.0

# Create mesh and define function space.
# Mesh for the 2-D case.
mesh = RectangleMesh(Point(0.0,0.0),Point(l,width),100,10)
# Use a vector-valued element for displacement.
P1 = VectorElement('P', triangle, 1)
# Use a scalar-valued element for pressure.
P2 = FiniteElement('P', triangle, 1)
# Create a mixed element to solve the system using the fully-implicit method.
element = MixedElement([P1, P2])
# Create the function space corresponding to this mixed element.
V = FunctionSpace(mesh,element)
# Create the test functions.
v_1, v_2 = TestFunctions(V)
# Create the trial functions.
ans = Function(V)
w_, p_ = split(ans)

# Define boundary conditions for pressure.
upper_boundary = 'on_boundary && near(x[0],0)'

# An expression for the derivative of P at the bottom surface.
der_p_D = Expression(('near(x[0],l) ? der_p_bottom : 0.0','0.0'),degree = 1, \
    der_p_bottom=der_p_bottom, l=l)

# Set the pressure at the porous medium's top to 0 (second subsystem of V).
bc_p = DirichletBC(V.sub(1),Constant(P_surf),upper_boundary)

# Define boundary conditions for displacements.
bottom_boundary = 'on_boundary && near(x[0],%s)'%(l)

der_w_D = Expression(('near(x[0],0.0) ? der_w_top : 0.0','0.0'),degree=0, \
    der_w_top=der_w_top)

# Set the displacement at the porous medium's bottom to 0 (first dimension of
# the first subsystem of the FunctionSpace V).
bc_w = DirichletBC(V.sub(0).sub(0),Constant(0.0),bottom_boundary)

# Define functions for solutions at previous and current time steps.
# Previous time step.
# Define the initial conditions.
init_cond = Constant(P0)
p_n = interpolate(init_cond,V.sub(1).collapse())           

# Define the variational problem.
# Equation 1
Eq1 = p_*v_1*dx - p_n*v_1*dx \
    +c*dt*inner(grad(p_),grad(v_1))*dx \
    -c*dt*der_p_D*v_1*ds 

# Equation 2
Eq2 = inner(grad(w_),grad(v_2))*dx \
    -der_w_D*v_2*ds \
    +cm*p_.dx*v_2*dx

# Add the equations to make a monolithic scheme.
F = Eq1 + Eq2

# Create the progress bar.
progress = Progress('Time-stepping')
set_log_level(PROGRESS)

# Form and solve the linear system.
for n in range(num_steps):

    # Update current time.
    t += dt

    # Solve the variational problem.
    solve(F == 0, ans, bc)

    # Check on the solution.
    w_,p_ = ans.split()
    p_comp = p_.compute_vertex_values(mesh)
    w_comp = w_.compute_vertex_values(mesh)

    # Update the previous solution.
    p_n.assign(p_)

    # Update the progress bar.
    progress.update(t/T)

The error occurs at the line "-cdtder_p_Dv_1ds" in defining Eq1. Have the test functions been made properly? If yes, what could be the problem in this case?

asked Jan 12, 2017 by SM FEniCS Novice (630 points)

2 Answers

+2 votes

v_1 is a vector test function, so you do not get a scalar.

Maybe you need the normal vector on the boundary? You can get it using

 n = FacetNormal(mesh)
answered Jan 12, 2017 by KristianE FEniCS Expert (12,900 points)

Thank you for your answer. No, I did not need the normal vector: I just needed to specify the derivative at that point.

I found out that the problem was that der_p_D and der_w_D yielded vectors when multiplied by the problem's variations (v_1, v_2). I changed the dimensions of der_p_D and der_w_D . When this was done (along with the use of the dot() function instead of * in some places), the issue was resolved.

0 votes

Maybe you need inner products instead of dot products in your form.
e.g. instead of

p_v_1dx

use:

inner(p_,v_1)*dx

answered Jan 16, 2017 by str FEniCS User (1,600 points)
...