Outward-facing normal on interior boundary

Question

I'm still trying to get a grasp on FEniCS and it seems like similar questions have been asked previously:

https://fenicsproject.org/qa/10737/compute-normal-derivative-on-boundary
https://fenicsproject.org/qa/11723/boundary-condition-that-depends-on-facet-normal
https://fenicsproject.org/qa/12074/calculation-functions-containing-normal-vector-boundary

but I'm still struggling with the concept of how to define the outward-facing normal on an interface. I've been experimenting with the following code that implements a DG method for a circular domain immersed in a square domain. I'd like to enforce a jump condition on the interior boundary (of the circle) that depends on the outward-facing normal on that interface. The idea would be to solve the Poisson problem at each time point, update the forced jump condition, then use the updated jump condition to solve the Poisson problem again. Is there an efficient way to do this?

from fenics import *
import os
import numpy as np
import matplotlib.pyplot as plt

# Import the mesh and determine relevant boundaries/subdomains
ofilename = '2D-DG_cell_fine.xml'
mesh = Mesh(ofilename)
subdomains0 = MeshFunction("size_t", mesh, "%s_physical_region.xml"%ofilename.split('.')[0])
boundaries0 = MeshFunction("size_t", mesh, "%s_facet_region.xml"%ofilename.split('.')[0])

# Rename subdomain numbers from GMSH
subdomains = CellFunction("size_t", mesh, 0)
subdomains.set_all(0)
subdomains.array()[subdomains0.array()==17] = 1
subdomains.array()[subdomains0.array()==18] = 2

# Rename boundary numbers from GMSH
boundaries = FacetFunction("size_t", mesh, 0)
boundaries.set_all(0)
boundaries.array()[boundaries0.array()==14] = 1
boundaries.array()[boundaries0.array()==15] = 2
#boundaries.array()[boundaries0.array()==16] = 3

interfaces = FacetFunction("size_t", mesh, 0)
interfaces.set_all(0)
interfaces.array()[boundaries0.array()==16] = 1

ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
dS = Measure('dS', domain=mesh, subdomain_data=interfaces)
dx = Measure('dx', domain=mesh, subdomain_data=subdomains)

V = FunctionSpace(mesh, 'DG', 1)
u = TrialFunction(V)
v = TestFunction(V)

n = FacetNormal(mesh)
h = CellSize(mesh)
h_avg = (h('+') + h('-'))/2

f = Constant(0.0)
s0 = 1*1e-6
s1 = 1*1e-6
dt = 0.1e-9
sm = 1e-9
e0 = 12e-12
dm = 10e-6
Cm = e0/dm*1e-12
Rm = sm/dm*1e-12

alpha = 1e4
gamma = alpha

D1= 0.
D2 = 0.

forced_jump = 0.00

for i in range(2):
    D1 = 80

    ## Deal with volume integrals (stiffness matrix)
    a = dot(grad(v), s0*grad(u))*dx(1) \
       + dot(grad(v), s1*grad(u))*dx(2)

    ## Deal with internal element interface conditions on non-boundary elements
    a += -dot(jump(u, n), avg(s0*grad(v)))*dS(0) \
       - dot(jump(v, n), avg(s0*grad(u)))*dS(0) \
       + (alpha/h_avg)*dot(jump(v, n), jump(u, n))*dS(0) \
       - dot(jump(v, n), avg(sm*grad(u)))*dS(1) \
       + (alpha/h_avg)*dot(jump(v,n), jump(u, n))*dS(1)

    ## Deal with boundary conditions
    a += - dot(s0*grad(v), u*n)*ds(1) \
       - dot(v*n, s0*grad(u))*ds(1) \
       + (gamma/h)*v*u*ds(1) \
       - dot(s0*grad(v), u*n)*ds(2) \
       - dot(v*n, s0*grad(u))*ds(2) \
       + (gamma/h)*v*u*ds(2) \

    ## Deal with source term  
    L = v*f*dx(1) + v*f*dx(2)

    ## Deal with Dirichlet boundary conditions explicitly
    L += - D1*dot(s0*grad(v), n)*ds(1) \
        + (gamma/h)*D1*v*ds(1) \
        - D2*dot(s0*grad(v), n)*ds(2) \
        + (gamma/h)*D2*v*ds(2) \
        - forced_jump*avg(dot(sm*grad(v), n))*dS(1) \
        + (gamma/h_avg)*forced_jump*jump(v)*dS(1)


    u0 = Function(V)
    solve(a == L, u0)

    ## THE ISSUE IS HERE--HOW DO I DEFINE THE OUTWARD FACING NORMAL ON A SURFACE???
    forced_jump = dt/Cm*(-s1*dot(grad(u0('+')),n('+')) - 1/Rm*forced_jump) + forced_jump

Comments

I wanted to update this with the jump condition I got to work on the internal boundary.

forced_jump = dt*(dot(grad(u0('-')), n('-')) - jump(u0)) + jump(u0)

and for the relevant terms in the linear form

    L += - inner(forced_jump*n('+'), avg(grad(v)))*dS(3) \
    + (gamma/h_avg)*inner(forced_jump, jump(v))*dS(3) 

Thanks again for all the help!

Answer 1

Read Chapter 30 of the FEniCS book, Modeling evolving discontinuities.

https://fenicsproject.org/book/

Comments

The documentation makes it appear that the cell_domains variable must be passed to the assembler--is that correct? How does that happen? Does this get passed to the assembler as a keyword argument to a function or is this something that gets set a different way? Thanks!

---

Orientation will work if the form has a volume integral, where you specify the orienting cell function as subdomain_data. Note, that you can always include such a term as L += Constant(0)*inner(v)*dx(...) where L would be your dx-free form.

---

So something like this should account for orientation if I add it to the LHS?
L += Constant(0)*v*dx(subdomain_data=subdomains)
Or is this something I have to add to each volume integral?

---

It should be added only to forms where you want to control normal orientation.

---

I've been working on this problem for a week and I'm wondering if there is a quick way to plot the direction the '+' and '-' refer to for a given edge/boundary. Is this something that can be done easily, even if this is normally hidden? Thanks!