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:
-
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) -
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()