Projecting subdomain solutions onto the full domain

General
1

Dear all,

I am wondering how to project solutions of 2 subdomains (which divide the full domain into two equal parts) onto the full domain, while preserving the the correct mapping of DoFs from these subdomains onto the full domain. Here is a MWE of what I am trying to say and do:

from dolfin import *
import math
import matplotlib.pyplot as plt

# Creating mesh
mesh_trial = RectangleMesh(Point(-1.0, -1.0), Point(1.0, 1.0), 8, 8,"left")

# Marking different regions in the domain
class Top_1(SubDomain):
    def inside(self,x,on_boundary):
        return x[1] > 1.0 - DOLFIN_EPS and x[0] < 0.0 + DOLFIN_EPS 
    
class Top_2(SubDomain):
    def inside(self,x,on_boundary):
        return x[1] > 1.0 - DOLFIN_EPS and x[0] > 0.0 - DOLFIN_EPS 
    
class Bottom_1(SubDomain):
    def inside(self,x,on_boundary):
        return x[1] < -1.0 + DOLFIN_EPS and x[0] < 0.0 + DOLFIN_EPS 
    
class Bottom_2(SubDomain):
    def inside(self,x,on_boundary):
        return x[1] < -1.0 + DOLFIN_EPS and x[0] > 0.0 - DOLFIN_EPS 
    
class Left(SubDomain):
    def inside(self,x,on_boundary):
        return x[0] < -1.0 + DOLFIN_EPS

class Right(SubDomain):
    def inside(self,x,on_boundary):
        return x[0] > 1.0 - DOLFIN_EPS

class Interface(SubDomain):
    def inside(self,x,on_boundary):
        return near(x[0], 0.0)

# Marking subdomains and interface between them
domain = MeshFunction("size_t", mesh_trial, mesh_trial.geometry().dim(), 0)
for c in cells(mesh_trial):
    domain[c] = c.midpoint().x() > 0.0
    
    
facets = MeshFunction("size_t", mesh_trial, mesh_trial.geometry().dim()-1, 0)
Interface().mark(facets,1)

# Creat submeshes of subdomains and interface from the full domain

mesh_L = MeshView.create(domain, 0)
mesh_R = MeshView.create(domain, 1)
mesh_I = MeshView.create(facets, 1)

# Marking boundaries of subdomains
mf_1 = MeshFunction("size_t", mesh_L, 1)
mf_1.set_all(0)
Top_1().mark(mf_1, 1)
Bottom_1().mark(mf_1, 2)
Left().mark(mf_1, 3)
Interface().mark(mf_1, 4)


mf_2 = MeshFunction("size_t", mesh_R, 1)
mf_2.set_all(0)
Top_2().mark(mf_2, 5)
Bottom_2().mark(mf_2, 6)
Right().mark(mf_2, 7)

# Defining function spaces and intergration measure for different subdomains

# Function space
V1 = FunctionSpace(mesh_L, "CG", 1)
V2 = FunctionSpace(mesh_R, "CG", 1)
Lmp = FunctionSpace(mesh_I,"R",0)

V = MixedFunctionSpace(V1,V2,Lmp)

u1,u2,u_lmp = TrialFunctions(V)
v1,v2,v_lmp = TestFunctions(V)

# Integration measure
dx1 = Measure("dx",domain=V.sub_space(0).mesh())
ds1 = Measure("ds",domain=V.sub_space(0).mesh(),subdomain_data=mf_1)

dx2 = Measure("dx",domain=V.sub_space(1).mesh())
ds2 = Measure("ds",domain=V.sub_space(1).mesh(),subdomain_data=mf_2)

ds_lmp   = Measure("dx", domain=V.sub_space(2).mesh())

# Define variational problem 

kappa1 = Constant(1.0)
kappa2 = Constant(10.0)
miu = Constant(10.0)
g = miu
a = kappa1*inner(grad(u1), grad(v1))*dx1 + kappa2*inner(grad(u2), grad(v2))*dx2
a+= inner(u_lmp, (v1-v2))*ds_lmp +inner(v_lmp,(u1-u2))*ds_lmp 
L = g*v1*ds1(3)


# BCs 
bc_dirich = DirichletBC(V.sub_space(1), Constant((0)), mf_2, 7)


# Compute solution
u = Function(V)
solve(a == L, u, bc_dirich)

# Splitting solution
u1,u2,u_lmp = u.split(deepcopy=True)

# Plotting solutions of subdomains
plt.colorbar(plot(u1)) 
plt.colorbar(plot(u2))

Running this example gives me solutions of two subdomains, which look good and has the the continuity at the interface (which is of course achieved using Lagrange multiplier).

S1
S12000×2000 78.5 KB
S2
S22000×2000 75.3 KB

I cannot however figure out how to project these two subdomain-solutions onto the full domain such that I get one full solution in the full domain.

Could someone please help me with that?