Overlapping coupled domains

variational formulation
1

I’m trying to solve for the voltage on two partly overlapping domains with a leakage resistance in the overlapping parts. I believe it can be described as follows.

\lambda_1^2 \nabla^2 u_1 = \left\{ \begin{array}{} 0, & \text{ on } \Omega_1\setminus\Omega_2 \\ u_1 - u_2 , & \text{ on } \Omega_1 \cap \Omega_2 \\ \end{array} \right.
\lambda_2^2 \nabla^2 u_2 = \left\{ \begin{array}{} 0, & \text{ on } \Omega_2\setminus\Omega_1 \\ u_2 - u_1, & \text{ on } \Omega_1 \cap \Omega_2 \\ \end{array} \right.
\lambda_1^2 \int_{\Omega_1} \nabla u_1 \cdot \nabla v_1 dA + \int_{\Omega_1 \cap \Omega_2} \left( u_1 - u_2 \right) v_1 dA = 0
\lambda_2^2 \int_{\Omega_2} \nabla u_2 \cdot \nabla v_2 dA + \int_{\Omega_1 \cap \Omega_2} \left( u_2 - u_1 \right) v_2 dA = 0

I have solved this using the code below, partly copied from here However, I expect it should be possible to solve this in a more elegant way, or at least to construct the system matrix more directly. Any suggesting on how to improve this code?

import matplotlib.pyplot as plt
from dolfin import *
import numpy as np

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return (x[0] <= 1/3+DOLFIN_EPS)

class Middle(SubDomain):
    def inside(self, x, on_boundary):
        return (x[0] >= 1/3-DOLFIN_EPS) and (x[0] < 2/3+DOLFIN_EPS)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return (x[0] >= 2/3-DOLFIN_EPS)

# Define mesh and subdomains
dim2=False
if dim2:
    mesh = UnitSquareMesh(270, 270)
    P1 = FiniteElement('P', triangle, 1)
else:
    mesh = UnitIntervalMesh(270)
    P1 = FiniteElement('P', interval , 1)
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(0)
Left().mark(domains, 1)
domains2 = MeshFunction("size_t", mesh, mesh.topology().dim())
domains2.set_all(0)
Middle().mark(domains2, 2)
domains3 = MeshFunction("size_t", mesh, mesh.topology().dim())
domains3.set_all(0)
Right().mark(domains3, 3)

# Define new measures associated with the interior domains
dx = Measure("dx", domain=mesh, subdomain_data=domains, subdomain_id=1)
dx2 = Measure("dx", domain=mesh, subdomain_data=domains2, subdomain_id=2)
dx3 = Measure("dx", domain=mesh, subdomain_data=domains3, subdomain_id=3)

elm = MixedElement([P1,P1])
# Define function space and basis functions
V = FunctionSpace(mesh,elm)
u1,u2 = TrialFunctions(V)
v1,v2 = TestFunctions(V)

def left_boundary ( x, on_boundary ):
    return on_boundary and (x[0]==0)
def right_boundary ( x, on_boundary ):
    return on_boundary and (x[0]==1)

bc0 = DirichletBC ( V, [-1,0], left_boundary )
bc1 = DirichletBC ( V, [0,1], right_boundary )

l1 = Constant(5)
l2 = Constant(1)

A1_1, b1_1 = assemble_system(l1*inner(grad(u1),grad(v1))*dx, Constant(0)*v1*dx, bc0)
A1_2, b1_2 = assemble_system(l1*inner(grad(u1),grad(v1))*dx2+(u1-u2)*v1*dx2, Constant(0)*v1*dx2)

A2_2, b2_2 = assemble_system(l2*inner(grad(u2),grad(v2))*dx2+(u2-u1)*v2*dx2, Constant(0)*v2*dx2)
A2_3, b2_3 = assemble_system(l2*inner(grad(u2),grad(v2))*dx3, Constant(0)*v2*dx3, bc1)

A=A1_1+A1_2+A2_2+A2_3
b=b1_1+b1_2+b2_2+b2_3
A.ident_zeros()

w = Function(V)
solve(A,w.vector(),b)

Thanks

2

How about the following?

import matplotlib.pyplot as plt
from dolfin import *
import numpy as np

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return (x[0] <= 1/3+DOLFIN_EPS)

class Middle(SubDomain):
    def inside(self, x, on_boundary):
        return (x[0] >= 1/3-DOLFIN_EPS) and (x[0] < 2/3+DOLFIN_EPS)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return (x[0] >= 2/3-DOLFIN_EPS)

# Define mesh and subdomains
dim2=False
if dim2:
    mesh = UnitSquareMesh(270, 270)
    P1 = FiniteElement('P', triangle, 1)
else:
    mesh = UnitIntervalMesh(270)
    P1 = FiniteElement('P', interval , 1)
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(0)
Left().mark(domains, 1)
Middle().mark(domains, 2)
Right().mark(domains, 3)

# Define new measures associated with the interior domains
dx = Measure("dx", domain=mesh, subdomain_data=domains)

elm = MixedElement([P1,P1])
# Define function space and basis functions
V = FunctionSpace(mesh,elm)
u1,u2 = TrialFunctions(V)
v1,v2 = TestFunctions(V)

def left_boundary ( x, on_boundary ):
    return on_boundary and (x[0]==0)
def right_boundary ( x, on_boundary ):
    return on_boundary and (x[0]==1)

bc0 = DirichletBC ( V, [-1,0], left_boundary )
bc1 = DirichletBC ( V, [0,1], right_boundary )

l1 = Constant(5)
l2 = Constant(1)

a = (
    l1*inner(grad(u1),grad(v1))*dx(1)
    + l1*inner(grad(u1),grad(v1))*dx(2)
    + (u1-u2)*v1*dx(2)
    + l2*inner(grad(u2),grad(v2))*dx(2)
    + (u2-u1)*v2*dx(2)
    + l2*inner(grad(u2),grad(v2))*dx(3)
)
L = Constant(0)*v1*dx + Constant(0)*v2*dx

A, b = assemble_system(a, L, [bc0, bc1])
A.ident_zeros()

w = Function(V)
solve(A,w.vector(),b)
1 Like
3

Yes, that indeed looks cleaner.
Thanks

1 Like