Hi, I am trying to solve a non-linear three-dimensional problem coupling elasticity and solvent diffusion. The state variables, defined over a three-dimensional body, are the displacement vector u, the pressure p and the solvent concentration c. To enforce the boundary conditions I have to define additional degrees of freedom (corresponding to a Lagrange multiplier) on the boundary. So, in the end, I have a mixed function space consisting of three spaces (for u, c and p) over the three-dimensional body and one function space over the boundary. When I run the script (see below) I get this error:
Sub elements must live in the same top level domain. Traceback (most
recent call last): File "swelling_gel_equilibrium_clean.py", line
14, in
W = MixedFunctionSpace([V,Q,R,M]) File "/Applications/FEniCS.app/Contents/Resources/lib/python2.7/site-packages/dolfin/functions/functionspace.py",
line 502, in init
element = ufl.MixedElement([V.ufl_element() for V in spaces]) File
"/Applications/FEniCS.app/Contents/Resources/lib/python2.7/site-packages/ufl/finiteelement/mixedelement.py",
line 64, in init
error("Sub elements must live in the same top level domain.") File
"/Applications/FEniCS.app/Contents/Resources/lib/python2.7/site-packages/ufl/log.py",
line 154, in error
raise self._exception_type(self._format_raw(message)) ufl.log.UFLException: Sub elements must live in the same top level
domain.
My guess is that I cannot define a mixed function space combining function spaces over the 3D domain and function spaces over the boundary: is it right? Is there a workaround for this problem, in case this feature is not supported?
Thanks a lot! AL
from dolfin import *
import math
import ufl
mesh = UnitCubeMesh(1,1,1)
V = VectorFunctionSpace(mesh, 'Lagrange', 2) # displacement
Q = FunctionSpace(mesh,'Lagrange',1) # pressure
R = FunctionSpace(mesh,'Lagrange',2) # bulk concentration
boundary_mesh = BoundaryMesh(mesh, "exterior")
M = FunctionSpace(boundary_mesh,'Lagrange',2) # Lagrange multiplier
W = MixedFunctionSpace([V,Q,R,M])
# Test and trial functions
vqrm = TestFunction(W) # Test function in the mixed space
v,q,r,tcs = split(vqrm)
dupccs = TrialFunction(W) # Trial function in the mixed space (solution increment)
upccs = Function(W) # Function in the mixed space to store the solution
u,p,c,cs = split(upccs) # Function in each subspace to write the functional
# nonlinear effective diffusion coefficient
def g(c):
return D*((2*chi-1.)*Omega*Jo*c-1.)/pow((1.+Omega*Jo*c),3)
def dmudc(c):
return -(Rg*T)*((2.*chi-1.)*Omega*Jo-1./c)/pow((1.+Omega*Jo*c),3)
# Parameters
Omega = 6.023e-5
D = 1E-1
Gdry = 4E5
Rg = 8.314
T = 293.0
lamo = 1.001
Jo = lamo**3
cod = (lamo**3-1.0)/Omega
co = cod/Jo
chi = 0.2
po = Gdry/lamo
muo = Rg*T*(math.log(Omega*cod/(1+Omega*cod))+1/(1+Omega*cod)+ chi/pow((1+Omega*cod),2)) \
+Omega*po
mus = muo
mutop = -10000
# Kinematics
I = Identity(3)
F = I + grad(u)
J = det(F)
# Constitutive equations
S = (Gdry/lamo)*F-p*cofac(F)
mu = Rg*T*(ufl.ln(Omega*Jo*c/(1+Omega*Jo*c)) + 1./(1.+Omega*Jo*c)+ \
chi/pow((1.+Omega*Jo*c),2))+Omega*p
h = g(c)*grad(c)+(-c*D/(Rg*T))*Omega*grad(p)
dmudp = Omega
dmu = mu-mus
# Boundary conditions
def left_boundary(x, on_boundary): # x = 0
return on_boundary and abs(x[0]) < DOLFIN_EPS
def back_boundary(x, on_boundary): # y = 0
return on_boundary and abs(x[1]) < DOLFIN_EPS
def bottom_boundary(x, on_boundary): # z = 0
return on_boundary and abs(x[2]) < DOLFIN_EPS
bcl = DirichletBC((W.sub(0)).sub(0), Constant(0.), left_boundary) # u = 0 on x = 0
bcb = DirichletBC((W.sub(0)).sub(1), Constant(0.), back_boundary) # v = 0 on y = 0
bcbo = DirichletBC((W.sub(0)).sub(2), Constant(0.), bottom_boundary) # w = 0 on z = 0
bc = [bcl,bcb,bcbo]
# WEAK FORM OF BALANCE EQUATIONS
a = inner(S,grad(v))*dx + q*(J-(1./Jo+Omega*c))*dx + inner(h,grad(r))*dx + inner(dmu,tcs)*ds \
+ inner(dmudc(c)*cs,r)*ds + inner(dmudp*cs,q)*ds
Jac = derivative(a,upccs,dupccs)
problem = NonlinearVariationalProblem(a, upccs, bc, Jac)
solver = NonlinearVariationalSolver(problem)
# Initial guess for unknowns
c = interpolate(Constant(co), R)
p = interpolate(Constant(po), Q)
cs = interpolate(Constant(co),M)
# solve non-linear problem
solver.solve()
# plot solvent concentration
w = Function(V)
w = upccs.split()[2]
plot(w)
interactive()
