PDE
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).
Question
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
plt.ion()
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())
cell_markers.set_all(False)
Border.mark(cell_markers,True,False)
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
plt.figure(1).clear()
h = plotV(u.vector(),**plot_settings)
plt.gca().axis('off')
plt.colorbar(h)
#plt.savefig('DG0sol')
# Zoom in on discontinuity at the boundary
plt.figure(7).clear()
h = plotV(u.vector(),**plot_settings)
plt.gca().axis('off')
x1, x2, y1, y2 = -0.659995, -0.60221, 0.74591, 0.78522
plt.xlim(x1,x2)
plt.ylim(y1,y2)
plt.colorbar(h)
#plt.savefig('DG0zoom')
}