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!
#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(Enu / ((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(xx + yy + 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(xx + yy + 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 mutemp + lmbdatr(temp)*Identity(g)
#Variational Form
F = inner(sigma(u), grad(tau))dx + alphagrad§taudx - inner(grad§, grad(v))*dx - (f/c)vdx
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()