# Why do the shapes below not match?

+1 vote

Hello. I am trying to implement a variational problem involving pressure and one-dimensional displacement. I get the following error from the code:

Shapes do not match: <Grad id=139643460028464> and <Indexed id=139643460375568>.


The error occurs at the lines where the variational problem is defined ("F = p_v_1dx + ...").

Here is the relevant code:

from fenics import *
import numpy as np
import matplotlib.pyplot as plt

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

# 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

# Create mesh and define function space.
mesh = IntervalMesh(100,0,l)
P1 = FiniteElement('P', interval, 1)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)

# Define boundary conditions for pressure.
tol = 1E-14

def upper_boundary(x,on_boundary):
return on_boundary and abs(x[0])<tol

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

bc_p = DirichletBC(V.sub(0),Constant(P_surf),upper_boundary)

# Define boundary conditions for displacements.
def bottom_boundary(x,on_boundary):
return on_boundary and abs(x[0]-l)<tol

der_w_D = Constant(der_w_top)

bc_w = DirichletBC(V.sub(1),Constant(0.0),bottom_boundary)

# Aggregate the Dirichlet boundary conditions in a list.
bc = [bc_p, bc_w]

# Define the test function space.
v_1, v_2 = TestFunctions(V)

# Define functions for solutions at previous and current time steps.
# Previous time step.
# Define the initial conditions.
init_cond = Expression(('P0','sigma_0/Kv*(l-x[0])'), degree = 1,
sigma_0=sigma_0, Kv=Kv, l=l,P0=P0)
u_n = interpolate(init_cond,V)
p_n,w_n = u_n.split()
# Current time step.
u = Function(V)
p_, w_ = split(u)

# Define the variational problem.
- p_n*v_1*dx \

# Change the variational problem into a bilinear form.
a,L = lhs(F),rhs(F)

# Now, declare the current time step's variable as a function.
p_ = Function(V)

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

# Open the file to store results in VTK format.
vtkfile_w = File('Terzaghi_results_FI/pressure.pvd')

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

# Update current time.
t += dt

# Solve the variational problem.
solve(a == L, u, bc)

# Extract the solution.
p_,w_ = u.split()
p_comp = p_.compute_vertex_values(mesh)
w_comp = w_.compute_vertex_values(mesh)

# Write the results to a VTK file for visualization in Paraview.
vtkfile_w << (p_,t)
vtkfile_w << (w_,t)

# Update the previous solution.
u_n.assign(u)

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


Any thoughts on why this may be occurring?

It is easy test out the individual terms: It is the second last term, which is causing the issue.

I suggest replacing

 cm*inner(grad(p_),v_2)*dx


with

 cm*p_.dx(0)*v_2*dx

answered Jan 3, 2017 by FEniCS Expert (12,900 points)
selected Jan 6, 2017 by SM

KristianE, thank you for your response. I appreciate it. The error no longer occurs at the above line. Two follow up questions:

1. How can one check the shapes of the individual terms? I am using
Spyder as my editor, and the variable window does not show functions
such as p_ and w_.
What is the .dx property?

.dx(0) is the first (x) spatial derivative. An alternative solution, is to use a VectorFiniteElement for the velocity. This is probably the better solution, since you would need to that anyway for dimensions higher than 1.