Hello everyone,
I am experimenting with the MixedFunctionSpace, and I am trying to couple different functions spaces. For now, I consider a simple diffusion-convection problem in a square, as preparation for more complex problems. Currently, I have a vector field u_1 that is defined on the whole square, while some scalar function u_2 only exists for x > 0.5. Here is my current code:
from dolfin import *
# Dirichlet boundary conditions
class DirichletBoundary1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary
class DirichletBoundary2(SubDomain):
def inside(self, x, on_boundary):
return abs(x[0] - 1.0) < DOLFIN_EPS and on_boundary
class DirichletBoundary3(SubDomain):
def inside(self, x, on_boundary):
return abs(x[0] - 0.5) < DOLFIN_EPS and on_boundary
# Create meshes
n = 16
mesh = UnitSquareMesh(n, n)
marker = MeshFunction("size_t", mesh, 2)
class RightDomain(SubDomain):
def inside(self, x, on_boundary):
return x[0] >= 0.5
r_domain = RightDomain()
r_domain.mark(marker, 1)
submesh1 = MeshView.create(marker, 0)
submesh2 = MeshView.create(marker, 1)
# Create function spaces
W1 = VectorFunctionSpace(mesh, "Lagrange", 1)
W2 = FunctionSpace(submesh2, "Lagrange", 1)
# Define the mixed function space W = W1 x W2
W = MixedFunctionSpace(W1, W2)
# Define boundary conditions
vector_fn = Expression(("x[0]", "0.0"), degree=0)
bc1 = DirichletBC(W1, vector_fn, "on_boundary")
bc2 = DirichletBC(W2, 1.0, "x[0] > 1 - DOLFIN_EPS && on_boundary")
bc3 = DirichletBC(W2, 0.5, "x[0] < 0.5 + DOLFIN_EPS && on_boundary")
bcs = [bc1, bc2, bc3]
# Define mixed-domains variational problem
(v1,v2) = TestFunctions(W)
u = Function(W)
u1 = u.sub(0)
u2 = u.sub(1)
dx2 = Measure("dx", domain=W.sub_space(1).mesh())
dx1 = Measure("dx", domain=W.sub_space(0).mesh(), subdomain_data=marker)
F1 = inner(u1, v1)*dx1(0) - inner(vector_fn, v1)*dx1(0) \
+ inner(u1, v1)*dx1(1) - u2*inner(Constant((1.0, 0.0)), v1) * dx2
F2 = inner(grad(u2), grad(v2))*dx2 #+ inner(grad(u2), u1)*v2*dx2
F = F1 + F2
# Compute solution
solve(F == 0, u, bcs)
s_1, s_2 = u.split()
File("test/vector.pvd") << s_1
File("test/diffu.pvd") << s_2
My problem now is, that I can not couple both functions. I can use the function u_1 together with the measure dx_2 without any problems (The part in F2 currently commented out). So the coupling from the whole domain down to the subdomain works fine.
But the other way around gives either errors or the newton-solver does not converge.
If I use:
- u2*inner(Constant((1.0, 0.0)), v1) * dx2
The newton-solver does only one “correct” iteration step and afterwards the residual stays constant (around 10^{-4}). If I use:
- u2*inner(Constant((1.0, 0.0)), v1) * dx1(1)
I get the error:
*** Error: Unable to successfully call PETSc function 'MatSetValuesLocal'.
*** Reason: PETSc error code is: 63 (Argument out of range).
*** Where: This error was encountered inside /build/dolfin-lgM8f6/dolfin-2019.2.0~git20220407.d29e24d/dolfin/la/PETScMatrix.cpp.
Is it possible to somehow couple the subdomain function u_2 together with the test functions of u_1? The problem itself should be well posed (solution without the convection in F2 is just u_1 = (x[0], 0) and u_2=x[0]). Most likely, I have to define some different measure? Any help would be appreciated, thanks!