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

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.
F = p_*v_1*dx + c*dt*inner(grad(p_),grad(v_1))*dx - c*dt*der_p_D*v_1*ds \
    - p_n*v_1*dx \
    + inner(grad(w_),grad(v_2))*dx + cm*inner(grad(p_),v_2)*dx-der_w_D*v_2*ds

# 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')

# 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.

    # Update the progress bar.

Any thoughts on why this may be occurring?

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

1 Answer

+2 votes
Best answer

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

I suggest replacing



answered Jan 3, 2017 by KristianE 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_.
  2. Regarding your expression above, is p_.dx(0) equivalent to grad(p_)?
    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.

Thank you for your answer. The code now works.