I am experimenting with a FETI-based approach (inspired by this example) in order to solve the mixed formulation of the Poisson problem (e.g. here). To this end, I have created a mesh consisting of two parts:
rect = MoveTo(0, .0).Rectangle(1, .5).Face()
rect.faces.name='Omega_1'
shape = rect
rect = MoveTo(0, .5).Rectangle(1, .5).Face()
rect.faces.name='Omega_2'
shape = Glue([shape, rect])
The finite elements for the function \sigma should be Raviart-Thomas:
fes = HDiv(mesh, order=0, RT=True, dirichlet=[])
In order to ensure continuity, it must be imposed that
\int_{\Gamma_{12}} \mu (\sigma_1 - \sigma_2) ds = 0,
where \Gamma_{12} is the interface between the domains, \sigma_{1/2} are the trial functions, and \mu is some test function. In the original FETI paper, this is implemented via an equation B_1 \sigma_1 + B_2 \sigma_2= 0 , where the B_i extract the required DOFs (turning the linear system into a saddle-point one) on the linear algebra level. The DOFs that have to be matched for the RT elements are the normals of the edges in the interface \Gamma_{12}.
The approach in the example above is to use the test functions of a boundary FES as values of \mu. Specifically, they set
feslam = None
for inter in mesh.Boundaries('.*').Split():
doms = inter.Neighbours(VOL).Split()
if len(doms) == 2:
feslami = Compress(H1(mesh, definedon=inter))
feslam = feslam * feslami if feslam else feslami
and construct the bilinear form corresponding to B based on this FES.
The problem is however, that this approach does not work for HDiv. When constructing HDiv on the intersection \Gamma_{12} there are zero free DOFs remaining, and hence no additional constraints are imposed.
My question then is, how do I match the normals of the RT elements along the interface between the meshes? Which subdomain and test functions are required?