System of PDEs defined on different subdomains

Question

Hi

I'm trying to solve the following nonlinear coupled equations

$$ (\nabla u_1,\nabla v_1)_{[0,1]} + (\nabla u_2, \nabla v_2)_{[0.5,1]} = (f(u_1,u_2),v_1)_{[0,1]} + (f(u_1,u_2),v_2)_{[0.5,1]} ,
$$

using the code sniper (testing with f = constant)

from __future__ import print_function
from fenics import *
import sympy as sym

# Create mesh and define function space
mesh = UnitIntervalMesh(50)
P1 = FiniteElement('P', mesh.ufl_cell(), 2)
element = MixedElement([P1,P1])
W = FunctionSpace(mesh, element)

# Defin subdomain
tol = 1e-14
class Omega0(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] <= 0.5 + tol
class Omega1(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] >= 0.5 - tol
subdomains = CellFunction("size_t", mesh)
subdomain0 = Omega0()
subdomain1 = Omega1()
subdomain0.mark(subdomains, 0)
subdomain1.mark(subdomains, 1)

# Define boundary condition
def boundary(x, on_boundary):
    return on_boundary
def boundary_right(x, on_boundary):
    return near(x[0], 1, tol)

bc0 = DirichletBC(W.sub(0), 0, boundary)
bc1 = DirichletBC(W.sub(1), 0, boundary_right)
bcs = [bc0, bc1]

# Define variational problem
v1, v2 = TestFunction(W)
u = Function(W)
u1, u2 = split(u)

dx = Measure('dx', domain=mesh, subdomain_data=subdomains) 
F = dot(grad(u1), grad(v1))*dx - v1*dx\
   +dot(grad(u2), grad(v2))*dx(1) - v2*dx(1)

# Compute solution
solve(F == 0, uu, bcs)

but it doesn't work. The variable $u_1$ is defined over the whole mesh $[0,1]$, whereas $u_2$ is only defined over $[0.5,1]$. I found this using submesh; my problem would then be how to mix mesh and submesh?

Thanks

Comments

What are the boundary conditions you want to enforce? Can you post a minimal working example?

Answer 1

Hi,
I think your equation is linear, so I change your code like this.

from __future__ import print_function
from fenics import *

# Create mesh and define function space
mesh = UnitIntervalMesh(10)
P1 = FiniteElement('P', mesh.ufl_cell(), 2)
element = MixedElement([P1, P1])
W = FunctionSpace(mesh, element)

# Define subdomains
tol = 1e-14

class Omega1(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] >= 0.5 - tol

# Mark subdomains
domain_marker = CellFunction("size_t", mesh)
domain_marker.set_all(0)

subdomain1 = Omega1()
subdomain1.mark(domain_marker, 1)

# Redefine the integration measure on subdomain
dx = Measure('dx', domain = mesh, subdomain_data = domain_marker) 

# Define boundary condition
def boundary(x, on_boundary):
    return on_boundary

def boundary_right(x, on_boundary):
    return (near(x[0], 1, tol) and on_boundary)

bc0 = DirichletBC(W.sub(0), 0, boundary)
bc1 = DirichletBC(W.sub(1), 0, boundary_right)
bcs = [bc0, bc1]


# Define trial and test functions
v1, v2 = split(TestFunction(W))
u1, u2 = split(TrialFunction(W))

# Define variational formulation
F = dot(grad(u1), grad(v1))*dx - v1*dx + dot(grad(u2), grad(v2))*dx(1) - v2*dx(1)
a, L = lhs(F), rhs(F)

# Compute solution
u_ = Function(W)
solve(a == L, u_, bcs)
u_1, u_2 = u_.split()

# Plot information
plot(domain_marker)
plot(u_1, title = 'u_1')
plot(u_2, title = 'u_2')

# Hold plot
interactive()