# Shape Mismatch in Mixed Function Space

**URL:** <https://fenicsproject.discourse.group/t/shape-mismatch-in-mixed-function-space/1931>\
**Category:** Uncategorized\
**Created:** [November 19, 2019, 5:14pm UTC](https://fenicsproject.discourse.group/t/shape-mismatch-in-mixed-function-space/1931 "2019-11-19T17:14:47Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![andrew\_59](https://avatars.discourse-cdn.com/v4/letter/a/5f9b8f/32.png) [@andrew\_59](https://fenicsproject.discourse.group/u/andrew_59)\
**Post date:** [November 19, 2019, 5:14pm UTC](https://fenicsproject.discourse.group/t/shape-mismatch-in-mixed-function-space/1931/1 "2019-11-19T17:14:48Z")

</div>

Hi everyone,

I have a problem with a mixed function space problem, I am attempting to solve the coupled linear elastic/ poroelastic fluid equations across a shell (photograph attached).

In the weak form on line 105 I get a “Can’t add expressions with different shapes” error, and I just can’t figure out what the problem is.

Any help at all would be much appreciated!

 ![pic1](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/1X/b9e5ac4e5af2d5d2a7166eb9cc8779c043a29fd3.jpeg)

#Import some things  
import fenics  
from mshr import \*  
from dolfin import \*  
import matplotlib.pyplot as plt  
from mshr.cpp import Sphere  
import numpy as np

#Parameters  
E = 584 #pa  
nu = 0.35 #poissons ratio  
mu = Constant(E / (1.0 + nu)) # = 2G  
lmbda = Constant(E_nu / ((1.0 + nu)_(1.0 - 2.0\*nu))) #= 2Gv/1-2v

skull\_radius = float(10)  
vent\_radius = float(5)  
resolution = int(32)  
p\_skull = 100  
p\_vent = 1  
u\_skull = Constant((0.0,0.0,0.0))  
u\_vent = Constant((0.1, 0.0 ,0.0))  
traction\_vent = Constant((0.0,0.0,0.0))  
f = Constant((1, 1, 1)) #Source term for water production  
c = 1 #kw  
alpha = Constant((1 ,1 ,1))

#Create a mesh  
center = dolfin.Point()  
sphere = Sphere(center, skull\_radius) - Sphere(center, vent\_radius)  
mesh = generate\_mesh(sphere, resolution)

#Create function space for pressure and displacement

cgp = FiniteElement(‘Lagrange’, mesh.ufl\_cell(), 1) #Set a 1d element  
cgu = VectorElement(‘Lagrange’, mesh.ufl\_cell(), 1) #Set a 3d element  
element = MixedElement([cgp, cgu])  
W = FunctionSpace(mesh, element)

u = TrialFunction(W)  
p = TrialFunction(W)

tau = TestFunction(W)  
v = TestFunction(W)

#Dirichlet BCs  
pr, dis = W.split()

class BoundarySkull(SubDomain):  
def inside(self, x, on\_boundary):  
x = x[0]  
y = x[1]  
z = x[2]  
r = sqrt(x_x + y_y + z\*z)  
return on\_boundary and r \< skull\_radius

class BoundaryVent(SubDomain):  
def inside(self, x, on\_boundary):  
x = x[0]  
y = x[1]  
z = x[2]  
r = sqrt(x_x + y_y + z\*z)  
return on\_boundary and r \> vent\_radius

skullBoundary = BoundarySkull()  
ventBoundary = BoundaryVent()

#Fails on boundary conditions  
bcu1 = DirichletBC(dis, u\_vent, ventBoundary)  
bcu2 = DirichletBC(dis, u\_skull, skullBoundary)

bcu = [bcu1, bcu2]

bcp1 = DirichletBC(pr, p\_vent, ventBoundary)  
bcp2 = DirichletBC(pr, p\_skull, skullBoundary)

bcp = [bcp1, bcp2]

bcs = [bcu, bcp]

#Functions for stress and strain tensors

g = u.geometric\_dimension()  
#new form  
def epsilon(u):  
return 0.5\*(grad(u) + np.transpose(grad(u)))  
def sigma(u):  
temp = epsilon(u)  
return mu_temp + lmbda_tr(temp)\*Identity(g)

#Variational Form

F = inner(sigma(u), grad(tau))_dx + alpha_grad§_tau_dx - inner(grad§, grad(v))\*dx - (f/c)_v_dx

a = lhs(F)  
l = rhs(F)

#Assemble matrices  
A1 = assemble(A1)  
[bc.apply(A1) for bc in bcs]

b1 = assemble(L1)  
[bc.apply(b1) for bc in bcs]

u = Function(W)

#Attempt to solve  
solve(A1, u\_.vector(), b1, “bicgstab”, “hypre\_amg”)  
solve(a==l, u, bcs)

(U,P) = u.split()