Solve PDE with discontinuous Dirichlet condition


I want to solve the following PDE with Dirichlet condition:
\begin{cases} -\nabla \cdot (\sigma \nabla u)=0 & \text{in }\Omega\\ u = f & \text{on }\Gamma\\ u = 0 & \text{on } \partial \Omega \backslash \Gamma\end{cases}

Where \Omega is the unit disk. The function f is assumed strictly increasing or decreasing which makes the function along the boundary composed of 0 and f discontinuous. For our theory it is important that the function along the boundary is solely increasing or decreasing, which I find very hard to achieve in my implementation (as there always appears a smoothing at the discontinuity).


Is there a way in FEniCS to achieve a solution to the above PDE with an actual jump in the boundary condition, while maintaining interior regularity (being able to take derivatives in the interior)? Or are there any recommendations on what would be the best solution?

Progress so far

I tried using three different types of basis functions: “DG” (with 0th or 1st order polynomials) and “CG” with 1st order polynomials. I feel “DG” with 0th order polynomials is closest to what I want at the boundary, as this seems the only way to enforce a strict discontinuity at the boundary. However, then I am not able to take derivatives in the interior, as I seem to get no contributions from the terms in the variational equation, where derivatives of u or the test function enter. The other two possibilities work way better in the interior, however, I feel I violate my assumption of a strict discontinuity at the boundary.

Illustration of my progress

Zoom in on discontinuity for solution when using “DG” with 0th order polynomials:

Code used for creating figure

from dolfin import *
import matplotlib.pyplot as plt
from dolfin import plot as dolfinplot
from mshr import generate_mesh, Circle
import numpy as np
import ufl

def plotV(vec,**kwargs):
        """ dolfin.plot-wrapper """
        fn = Function(V,vec)
        return dolfinplot(fn,**kwargs)

#Define Boundary
class DirichletBoundary(SubDomain):
  def inside(self, x, on_boundary):
    return (abs((x[0]**2+x[1]**2)-1)  < 1e-3)

#Create mesh
C = Circle(Point(0,0),1)
msh = generate_mesh(C,100)

#Refine mesh at the boundary
class Border(SubDomain):
    def inside(self, x, on_boundary):
        return (abs((x[0]**2+x[1]**2)-1)  < 3e-2)

Border = Border()

nRef = 4

for i in range(nRef):
    cell_markers = MeshFunction('bool',msh,msh.topology().dim())

    msh = refine(msh,cell_markers)

#Define Function space
V = FunctionSpace(msh, 'DG', 0)
#V = FunctionSpace(msh, 'DG', 1)
#V = FunctionSpace(msh, 'CG', 1)

xy = V.tabulate_dof_coordinates().reshape((-1,2))

#Define conductivity
sig = Function(V)
sig.vector()[:] = 1 + 5*np.exp(-20*(np.power(xy[:,0]+1/2,2)+np.power(xy[:,1],2))) + 5*np.exp(-20*(np.power(xy[:,0],2)+np.power(xy[:,1]+1/2,2))) + 5*np.exp(-50*(np.power(xy[:,0]-1/2,2)+np.power(xy[:,1]-1/2,2)))

# Define test and trial functions
u = TrialFunction(V)
v = TestFunction(V)

# Define normal vector and mesh size
n = FacetNormal(msh)
h = CellDiameter(msh)
h_avg = (h('+') + h('-'))/2

# Define the source term f and Dirichlet term u0
class MyExpression1(UserExpression):
        def eval(self, value, x):
            if -(1.0/8.0)*ufl.atan(1)*4.0 < ufl.atan_2(x[1],x[0]) and ufl.atan_2(x[1],x[0]) < (6.0/8.0)*ufl.atan(1)*4.0-0.1:
                value[0] = (3.0/5.0)*ufl.sqrt(ufl.atan_2(x[1],x[0])+(1.0/8.0)*ufl.atan(1)*4.0)
            else :
                value[0] = 0

u0 = MyExpression1(degree=0)
f = Constant(0)

# Mark facets of the mesh
boundaries = MeshFunction('size_t',msh,msh.topology().dim()-1)
DirichletBoundary().mark(boundaries, 1)

# Define outer surface measure aware of Dirichlet and Neumann boundaries
ds = Measure('ds', domain=msh, subdomain_data=boundaries)

# Define parameters
alpha = 4.0
gamma = 8.0

a = dot(grad(v), sig*grad(u))*dx \
   - dot(avg(sig*grad(v)), jump(u, n))*dS \
   - dot(jump(v, n), avg(sig*grad(u)))*dS \
   + alpha/h_avg*dot(jump(v, n), jump(u, n))*dS \
   - dot(sig*grad(v), u*n)*ds(1) \
   - dot(v*n, sig*grad(u))*ds(1) \
   + (gamma/h)*v*u*ds(1)
L = v*f*dx - u0*dot(sig*grad(v), n)*ds(1) + (gamma/h)*u0*v*ds(1)

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

plot_settings = {
        'cmap': plt.get_cmap('RdBu'),

# Plot solution
h = plotV(u.vector(),**plot_settings)

# Zoom in on discontinuity at the boundary
h = plotV(u.vector(),**plot_settings)
x1, x2, y1, y2 = -0.659995, -0.60221, 0.74591, 0.78522

If you want to use weak BC enforcement via Nitsche’s method, as in your MWE, the smoothing effect of the interior PDE will always “compete” with the discontinuous boundary data. However, you can prioritize the enforcement of the BC by increasing the penalty parameter gamma.

The case of discontinuous Dirichlet data is actually a nontrivial mathematical setting, because you can’t even assume that the solution is in H^1(\Omega). (If it were, the trace on \partial\Omega would need to be in H^{1/2}(\partial\Omega), which does not permit jump discontinuities.) This thesis derives an error estimate by considering the boundary discontinuity to be smoothed over an O(h) length scale and applying the standard H^1-conforming formulation (Theorem 2.2.1). I didn’t read all the fine print, but I believe it would amount in practice to directly interpolating the discontinuous data and applying the BC strongly, since the O(h) smoothing could fit between adjacent boundary nodes.

If you want to use a degree-0 DG space for u, you might also consider a mixed Poisson formulation with flux in a discrete H(\operatorname{div}) space, as in the demo here.

1 Like

Thank you for your quick reply and valuable input! I was not aware of the fact, that I could use “gamma” to prioritize enforcement of the BC.

I was aware that I could not obtain H^1-regularity for the solution over the whole of \Omega. However, away from the boundary I should be able to control the smoothness of u by the conductivity \sigma. But this seems rather difficult to achieve numerically. The reference to the thesis seems a really good way to quantify the error of a smoothed representation of the boundary condition. I will look further into that!

Actually I was hoping to avoid a degree-0 DG space for u, so when possible I will focus on the other options.