Coupled equations on intersection of domains

[fewagner]

Hi,

I am trying to implement two equations, that are defined on two different domains, but the domains are overlapping on a certain subdomain. On this overlapping subdomain, the two equations are coupled. I tried to write this as a minimal example with a Poisson equation:

In code, I would try to split the coupled terms and the regular terms, maybe like that:

@fem.BilinearForm
def a(u, v, w):
    poisson = dot(grad(u), grad(v))
    return poisson

@fem.BilinearForm
def a_coupling(u, v, w):
    coupling = - w['kappa'] * u * v
    return coupling

@fem.LinearForm
def l(v, w):
    x, y = w.x  # global coordinates
    f = np.sin(np.pi * x)
    return f * v

@fem.LinearForm
def l_coupling(v, w):
    x, y = w.x  # global coordinates
    coupling = w['kappa'] * w['u_other']
    return coupling * v

and define the two meshes like that:

points = np.array([[0, 1, 2, 3, 0, 1, 2, 3],  # x 
                   [0, 0, 0, 0, 1, 1, 1, 1],  # y
                  ])  # mesh.p
cells = np.array([[0, 5, 1, 6, 2, 7], 
                  [1, 4, 2, 5, 3, 6], 
                  [4, 1, 5, 2, 6, 3]
                 ])  # mesh.t, topology connection betw points

mesh0 = fem.MeshTri(points, cells)

points = np.array([[0, 1, 0, 1, 0, 1, 0, 1],  # x 
                   [0, 0, 1, 1, 2, 2, 3, 3],  # y
                  ])  # mesh.p
cells = np.array([[0, 1, 2, 3, 4, 5], 
                  [1, 2, 3, 4, 5, 6], 
                  [2, 3, 4, 5, 6, 7]])  # mesh.t, topology connection betw points

mesh1 = fem.MeshTri(points, cells)

However, I am not really sure how to proceed from there. I would want a basis that is defined on both of the meshes, but then build the assemblies only on the subdomains - and add up all the assemblies when calling solve and condense? I would appreciate any advice or pointers on helpful pieces in the documentation, examples, and previous discussions.

[kinnala]

Do the meshes have to be separate or is it OK to generate one mesh for the union of the domains? The latter will greatly simplify the solver as you do not need to find the overlapping elements and their intersections one by one.

[fewagner]

It is OK to generate one mesh for the unions of the domains, like that

mesh = mesh0 + mesh1

However, I would then need to constrain the non-coupled assemblies a and l on the individual subdomains of mesh0 and mesh1? And the coupled terms to the intersection of the grids?

[kinnala]

That summation is probably not the best idea as it might lead to nontrivial bugs if you are not completely certain that all overlapping elements are equal in shape and the node locations are equal up to the floating point accuracy.

I would simply generate a mesh for the union, whatever shape it is, and tag the elements belonging to subdomain 1 and subdomain 2. Then create separate basis objects for assembly: basis1 = Basis(unionmesh, elem, elements='subdomain1') and basis2 = Basis(unionmesh, elem, elements='subdomain2').

The resulting matrices have some extra zero rows/cols corresponding to the ”unused” degrees of freedom but it is easy enough to slice off those rows if you know how the rows/cols are ordered.

[fewagner]

Thank you already for your response, this sounds like a much better approach! I made much progress but am facing now a behavior that makes me doubt my understanding of the basis on subdomains. I am posting my code below, its executable and commented:

# define forms for equations: 
# (grad u0, grad v) + kappa*(u0, v) = (f,v) + kappa*(u1, v)
# (grad u1, grad v) + kappa*(u1, v) = kappa*(u0, v)

@fem.BilinearForm
def a(u, v, _):
    poisson = dot(grad(u), grad(v))
    return poisson

@fem.BilinearForm
def coupling(u, v, w):
    coupling = w['kappa'] * u * v
    return coupling

@fem.LinearForm
def l(v, _):
    f = 1. 
    return f * v

# main mesh L shaped [0,3]x[0,3] without [1,3]x[1,3]
points = np.array([[0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 0, 1],  # x 
                   [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 3, 3],  # y
                  ])  # mesh.p
cells = np.array([[0, 5, 1, 6, 2, 7, 4, 5, 8, 9], 
                  [1, 4, 2, 5, 3, 6, 5, 8, 9, 10], 
                  [4, 1, 5, 2, 6, 3, 8, 9, 10, 11]
                 ])  # mesh.t, topology connection betw points
mesh = fem.MeshTri(points, cells)
mesh = mesh.with_subdomains({
    'mesh0': lambda x: x[0] <= 1.,
    'mesh1': lambda x: x[1] <= 1.,
    'intersect': lambda x: np.logical_and(x[0] <= 1., x[1] <= 1.),
}).refined(3)

# define basis for individual domains and intersection

Vh0 = fem.Basis(mesh, fem.ElementTriP1(), elements='mesh0') 
Vh1 = fem.Basis(mesh, fem.ElementTriP1(), elements='mesh1')
Vhi = fem.Basis(mesh, fem.ElementTriP1(), elements='intersect')
Vh = fem.Basis(mesh, fem.ElementTriP1())

# build assemblies

A0 = a.assemble(Vh0)
A1 = a.assemble(Vh1)
b0 = l.assemble(Vh0)
G = coupling.assemble(Vhi, kappa=1.)

# get dof for boundary conditions
 
D = Vh.get_dofs()
D = np.concatenate((D.all(), D.all() + D.N))

# build coupled system

M = block_array([[A0+G, -G], [-G, A1+G]], format='csr')
c = np.concatenate((b0, Vh.zeros()))

# solve

x = fem.solve(*fem.condense(M, c, D=D)) # this printes dgstrf info 394 

x0 = x[:Vh.N]
x1 = x[Vh.N:]

The problem is that the solution x contains mostly nans, probably everywhere where it is not constrained by the Dirichlet conditions. As a utility check I removed the G from the matrix M, with the same result. Somehow the solve function does not seem to like bases that are constrained to a subdomain? Or am I on a completely wrong trail?

EDIT: I changed to code such that I am not summing meshes, which does not change the outcome.
EDIT2: I added a line after the definition of the Dirichlet conditions, such that they are defined for all blocks in the larger, coupled matrix (didnt change the problem though).
EDIT3: I think what Scipy wants to tell me is, that the matrices are singular. This also makes sense due to the extra zero rows/cols for the subbasis. So probably I need to understand now how to exclude the zero-elements from the matrices?

[kinnala]

Yes, so basically you would slice those unnecessary zero rows/cols. You can find the corresponding dofs using I1 = Vh.get_dofs(elements='mesh1').all() and then you can slice the matrix A0 = A0[I1][:, I1].

[fewagner]

Hey, it took me a little but I wanted to come back to that. Thanks a lot for your help, I managed to get it running. I leave my code for completeness and future reference as an executable script below. Overall it was not trivial, but pretty cool that it worked! :)

#!/usr/bin/env python
# coding: utf-8

# In[1]:


import numpy as np
import matplotlib.pyplot as plt
import skfem as fem
from skfem.helpers import dot, grad
from scipy.sparse import block_array


# In[2]:


# define forms for equations: 
# (grad u0, grad v) + kappa*(u0, v) = f + kappa*(u1, v)
# (grad u1, grad v) + kappa*(u1, v) = kappa*(u0, v)

@fem.BilinearForm
def a(u, v, _):
    poisson = dot(grad(u), grad(v))
    return poisson

@fem.BilinearForm
def coupling(u, v, w):
    coupling = w['kappa'] * u * v
    return coupling

@fem.LinearForm
def l(v, _):
    f = 1. 
    return f * v


# In[3]:


# main mesh L shaped [0,3]x[0,3] without [1,3]x[1,3]
points = np.array([[0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 0, 1],  # x 
                   [0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 3, 3],  # y
                  ])  # mesh.p
cells = np.array([[0, 5, 1, 6, 2, 7, 4, 5, 8, 9], 
                  [1, 4, 2, 5, 3, 6, 5, 8, 9, 10], 
                  [4, 1, 5, 2, 6, 3, 8, 9, 10, 11]
                 ])  # mesh.t, topology connection betw points
mesh = fem.MeshTri(points, cells)
mesh = mesh.with_subdomains({
    'mesh0': lambda x: x[0] <= 1.,
    'mesh1': lambda x: x[1] <= 1.,
    'intersect': lambda x: np.logical_and(x[0] <= 1., x[1] <= 1.),
}).refined(3)


# In[4]:


# define basis

Vh0 = fem.Basis(mesh, fem.ElementTriP1(), elements='mesh0') 
Vh1 = fem.Basis(mesh, fem.ElementTriP1(), elements='mesh1')
Vhi = fem.Basis(mesh, fem.ElementTriP1(), elements='intersect')
Vh = fem.Basis(mesh, fem.ElementTriP1())


# In[5]:


# build assemblies

kappa = 1.

A0 = a.assemble(Vh0)
A1 = a.assemble(Vh1)
b0 = l.assemble(Vh0)
G = coupling.assemble(Vhi, kappa=kappa)


# In[6]:


# get dof for slicing

I0 = Vh.get_dofs(elements='mesh0').all()
I1 = Vh.get_dofs(elements='mesh1').all()
Ii = Vh.get_dofs(elements='intersect').all()


# In[7]:


# get boundary dofs
# mesh 0
x, y = Vh.get_dofs().doflocs[:, I0]
D0 = []
for d,x_,y_ in zip(I0,x,y):
    if np.isclose(x_,0.) or np.isclose(x_,1.) or np.isclose(y_,0.) or np.isclose(y_,3.):
        D0.append(d)
D0 = np.array(D0)

# mesh 1
x, y = Vh.get_dofs().doflocs[:, I1]
D1 = []
for d,x_,y_ in zip(I1,x,y):
    if np.isclose(x_,0.) or np.isclose(x_,3.) or np.isclose(y_,0.) or np.isclose(y_,1.):
        D1.append(d)
D1 = np.array(D1)


# In[8]:


# slice

A0 = A0[I0][:, I0]
A1 = A1[I1][:, I1]
b0 = b0[I0]
G00 = G[I0][:, I0]
G01 = G[I0][:, I1]
G10 = G[I1][:, I0]
G11 = G[I1][:, I1]


# In[9]:


def reduce_dofs(D, I, N):  
    flag = np.zeros(N, dtype=bool)
    flag[D] = True
    flag = flag[I]
    return flag.nonzero()[0]


# In[10]:


D0 = reduce_dofs(D0, I0, Vh.N)
D1 = reduce_dofs(D1, I1, Vh.N)
D = np.concatenate((D0, D1 + A0.shape[0]))


# In[11]:


# build coupled system

M = block_array([[A0+G00, -G01], [-G10, A1+G11]], format='csr')
c = np.concatenate((b0, np.zeros(A1.shape[0])))


# In[12]:


# solve

x = fem.solve(*fem.condense(M, c, D=D))

#x0 = x[:Vh.N]
#x1 = x[Vh.N:]


# In[13]:


x0 = Vh.zeros()
x0[I0] = x[:I0.shape[0]]

x1 = Vh.zeros()
x1[I1] = x[I0.shape[0]:]


# In[14]:


# plot x0

Vh.plot(x0, colorbar=True)
plt.show()

# plot x1

Vh.plot(x1, colorbar=True)
plt.show()


# In[ ]: