How do I implement a thin domain approximation?

variational formulation
1

I have a problem with a Poisson equation which is in essence 2 large domains separated by a very thin domain (green) and a source term on a boundary (thick line).

ex

The problem is the thin domain is very small so it’s quite expensive to mesh.

What I would like to do is to replace the thin domain with a boundary condition. Impose continuity of the current across it, but not the voltage and say there’s a voltage drop d*Jr where d is the width and r is the resistance per unit length.

I’m not quite sure how to implement this in fenics and any suggestions would be appreciated.

Thanks,
Aaron

2 Likes
2

Hi,

I would suggest defining your thin domain as a submesh of lower dimension, and imposing your condition on this interface using a Lagrange multiplier. You can have a look to the dedicated mixed-dimensional features developed in FEniCS, in the branch cecile/mixed-dimensional, or to the multiphenics library.

The mixed-dimensional branch is also available as a Docker container ceciledc/fenics_mixed_dimensional. In particular, the Docker container includes a demo with a Lagrange multiplier defined at the interface between subdomains. This post might also be relevant.

1 Like
3

Hi, let me mention an alternative to previous answers which avoids the multiscale formulation. As such you can readily implement it in standard FEniCS. In addition, compared to the suggested mortaring approach you will get a better representation of flux/current. Moreover, standard preconditioners will work (for the multiscale preconditioners see here for mortaring and here mixed mortaring approach). Unless you have a use for the Lagrange multiplier further in your solver, as is done eg. here, the only downside is that a larger problem is to be solved.

from dolfin import *

# Let Gamma be the interface between Omega_1, Omega_2 and consider
# the problem
# 
#   -div(kappa_i*grad(u_i)) = f_i in Omega_i
#
# such that on Gamma
#
#   kappa_0*dot(grad(u_0), n_0) + kappa_1*dot(grad(u_1), n_1) = 0, (1)
#   u_0 - u_1 = J.                                                  (2)
#
# Here n_i is the outer normal of Omega_i. Condition (1) is the flux/current
# continuity while (2) forces jumps in the potentials/voltages.
#
# Let us introduce fluxes sigma_i = kappa_i*grad(u_i) as explicit unknowns.
# By (1) the fluxes are continuous accross the interfaces. With the new
# unknowns the system reads
#
#   kappa^{-1}_i sigma_i - grad(u_i) = 0   in Omega_i (3)
#   -div(sigma_i) = f_i                    in Omega_i
#
# while the interface conditions are
#
#   dot(sigma_0,  n_0) + dot(sigma_1, n_1) = 0  (4)
#   u_0 - u_1 = J.
#
# We can enforce (3) by construction if sigma, such that sigma|_Omega_i = sigma_i
# is sought after in H-div conforming (vector valued) finite element space. 
# Let then tau be a test function from such a space. Then integrating (3)
# on each subdomain leads to
#
#   inner(kappa^{-1}_i*sigma_i, tau_i)*dx(i) + inner(div(tau_i), u_i)*dx(i) +
#
# an interface term -inner(u_i, dot(tau_i, n_i))*dS(1). However, (4) for
# tau by construction the equations can be added yielding the interface
# term
# 
#    -inner((u_0 - u1), dot(tau, n_0))*dS(1).
#
# Plugging u_0 - u_1 = J, the interface term contributes to the right
# hand side. Let's implement the idea:

mesh = UnitSquareMesh(32, 32)
left = CompiledSubDomain('x[0] < 0.5+DOLFIN_EPS')
thin = CompiledSubDomain('near(x[0], 0.5)')

# Conductivities for subdomains
kappa_1, kappa_2 = Constant(1), Constant(2)
# Forcing on subdomains
f_1, f_2 = Constant(1), Constant(2)
# Interface forcing
J = Constant(3)


tdim = mesh.topology().dim()
cell_f = MeshFunction('size_t', mesh, tdim, 2)
left.mark(cell_f, 1)
# Interface
facet_f = MeshFunction('size_t', mesh, tdim-1, 0)
thin.mark(facet_f, 1)

cell = mesh.ufl_cell()
# Flux and potential
elm = MixedElement([FiniteElement('Raviart-Thomas', cell, 1),
                    FiniteElement('Discontinuous Lagrange', cell, 0)])
W = FunctionSpace(mesh, elm)

sigma, u = TrialFunctions(W)
tau, v = TestFunctions(W)

dx = Measure('dx', domain=mesh, subdomain_data=cell_f)
dS = Measure('dS', domain=mesh, subdomain_data=facet_f)
n = FacetNormal(mesh)

a = inner((1./kappa_1)*sigma, tau)*dx(1) + inner((1./kappa_2)*sigma, tau)*dx(2)
a += inner(div(tau), u)*dx() + inner(div(sigma), v)*dx()

L = inner(-f_1, v)*dx(1) + inner(-f_2, v)*dx(2)
L += inner(avg(J), dot(avg(tau), n('+')))*dS(1)

wh = Function(W)
solve(a == L, wh)

sigma, u = wh.split(deepcopy=True)

File('u.pvd') << u
File('sigma.pvd') << sigma
1 Like