Inconsistent shape with grad

Hi,
I am relatively new to FEniCS. I am solving the poroelasticity problem. To begin with, I have three equations: an equation for displacement and two equations for pressure. A weak formulation is presented below:

I get the following error:

Cannot take symmetric part of rectangular matrix with dimensions (2, 4).
Traceback (most recent call last):
displ_part = (-lmbd * div(u) * div(v) - 2 * mu * inner(epsilon(u), epsilon(v)) + alphaArt * pa * div(v) + alphaVen * pv * div(v)) * dx
line 35, in epsilon
return sym(nabla_grad(val))
ufl.log.UFLException: Cannot take symmetric part of rectangular matrix with dimensions (2, 4).

My code:

lmbd = 504
mu = 216
alphaArt = 0.25
alphaVen = 0.01
arterialK = 1.0e-10
arterialMU = 2.67e-3
venousK = 1.0e-10
venousMU = 2.67e-3
# --------------------------------------------------
# Boundary conditions
# --------------------------------------------------
uExt = Constant((0, 0))
arterialP = 8000
venousP = 650


def epsilon(val):
    return sym(nabla_grad(val))

# --------------------------------------------------
# Construct FE mesh and normal
# --------------------------------------------------
mesh = Mesh("geometry.xml")

subdomains = MeshFunction("size_t", mesh, "geometry_physical_region.xml")
boundaries = MeshFunction("size_t", mesh, "geometry_facet_region.xml")
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
dx = Measure('dx', domain=mesh, subdomain_data=subdomains)
n = FacetNormal(mesh)

# --------------------------------------------------
# Construct FE space
# --------------------------------------------------
P1 = FiniteElement('P', triangle, 2)
P2 = VectorElement('P', triangle, 2)
W = FunctionSpace(mesh, MixedElement([P2, (P1, P1)]))

# --------------------------------------------------
# Define variational problem
# --------------------------------------------------

u = TrialFunction(W)
pa = TrialFunctions(W)
pv = TrialFunction(W)

v = TestFunction(W)
qa = TestFunction(W)
qv = TestFunction(W)



bcs = [
    DirichletBC(W.sub(0), uExt, boundaries, 22),
    DirichletBC(W.sub(1).sub(0), arterialP, boundaries, 22),
    DirichletBC(W.sub(1).sub(1), venousP, boundaries, 22)
]
displ_part = (-lmbd * div(u) * div(v) - 2 * mu * inner(epsilon(u), epsilon(v)) + alphaArt * pa * div(v) + alphaVen * pv * div(v)) * dx
pa_part = ((arterialK / arterialMU) * dot(grad(pa), grad(qa))) * dx
pv_part = ((venousK / venousMU) * dot(grad(pv), grad(qv))) * dx
bound_part = (pa * dot(n, v) + pv * dot(n, v)) * ds(11)
b = displ_part + pa_part + pv_part - bound_part


z = Function(W)
a, L = lhs(b), rhs(b)

solve(a == L, z, bcs)
(u, pa, pv) = z.split()

I also have a number of questions:

  1. Why can’t I unpack the trial and test function as follows: u, pa, pv = TrialFunctions(W); v, qa, qv = TestFunctions(W) ? I get an error: ValueError: not enough values to unpack (expected 3, got 2)

  2. What is the correct way to get a solution in my problem? Am I doing the right thing?:

z = Function(W)
a, L = lhs(b), rhs(b)

solve(a == L, z, bcs)
(u, pa, pv) = z.split()

Hi,
you should not enclose P1,P1 with parenthesis since in this case W is the product of a P2 space and a mixed space which is the product of P1 and P1. Thus TrialFunctions(W) yields two functions, the second one that you can further split in two:

u, p = TrialFunctions(W)
pa, pv = split(p)

So just change to

W = FunctionSpace(mesh, MixedElement([P2, P1, P1]))

u,pa,pv = TrialFunctions(W)
v,qa,qv = TestFunctions(W)

this will also fix the shape error you had earlier

2 Likes

Thanks bleyerj, works great :slightly_smiling_face: