# Re-entrant Corner Problem¶

Here we will consider the classic re-entrant corner problem, \begin{align*} -\Delta u &= f, && \text{in } \Omega, \\ u &= g, && \text{on } \partial\Omega, \end{align*} where the domain is given using polar coordinates, \begin{gather*} \Omega = \{ (r,\varphi)\colon r\in(0,1), \varphi\in(0,\Phi) \}~. \end{gather*} For the boundary condition $$g$$, we set it to the trace of the function $$u$$, given by \begin{gather*} u(r,\varphi) = r^{\frac{\pi}{\Phi}} \sin\big(\frac{\pi}{\Phi} \varphi \big) \end{gather*}

[1]:

import matplotlib
matplotlib.rc( 'image', cmap='jet' )
import numpy
import matplotlib.pyplot as pyplot
from dune.fem.plotting import plotPointData as plot
import dune.grid as grid
import dune.fem as fem
from dune.fem.space import lagrange as solutionSpace
from dune.alugrid import aluConformGrid as leafGridView
from dune.fem.function import integrate, uflFunction
from ufl import *
from dune.ufl import DirichletBC

# set the angle for the corner (0<angle<=360)
cornerAngle = 320.

# use a second order space
order = 2


We first define the domain and set up the grid and space. We need this twice - once for a computation on a globally refined grid and once for an adaptive one so we put the setup into a function:

We first define the grid for this domain (vertices are the origin and 4 equally spaced points on the unit sphere starting with (1,0) and ending at (cos(cornerAngle), sin(cornerAngle))

Next we define the model together with the exact solution.

[2]:

def setup():
vertices = numpy.zeros((8, 2))
vertices[0] = [0, 0]
for i in range(0, 7):
vertices[i+1] = [numpy.cos(cornerAngle/6*numpy.pi/180*i),
numpy.sin(cornerAngle/6*numpy.pi/180*i)]
triangles = numpy.array([[2,1,0], [0,3,2], [4,3,0],
[0,5,4], [6,5,0], [0,7,6]])
domain = {"vertices": vertices, "simplices": triangles}
gridView.hierarchicalGrid.globalRefine(2)
space = solutionSpace(gridView, order=order, storage="numpy")

from dune.fem.scheme import galerkin as solutionScheme
u = TrialFunction(space)
v = TestFunction(space)
x = SpatialCoordinate(space.cell())

# exact solution for this angle
Phi = cornerAngle / 180 * pi
phi = atan_2(x[1], x[0]) + conditional(x[1] < 0, 2*pi, 0)
exact = dot(x, x)**(pi/2/Phi) * sin(pi/Phi * phi)

# set up the scheme
laplace = solutionScheme([a==0, DirichletBC(space, exact, 1)], solver="cg",
parameters={"newton.linear.preconditioning.method":"jacobi"})
uh = space.interpolate([0], name="solution")
return uh, exact, laplace


We will start with computing the $$H^1$$ error on a sequence of globally refined grids

Note that by using fem.globalRefine instead of hierarchicalGrid.globalRefine we can prolongate discrete functions to the next level. The second argument can also be a list/tuple of discrete functions to prolong. With this approach we optain a good initial guess for solving the problem on the refined grid.

[3]:

uh, exact, laplace = setup()
errorGlobal = []
dofsGlobal  = []
for count in range(10):
laplace.solve(target=uh)
error = numpy.sqrt(integrate(uh.space.gridView, h1error, 5))
errorGlobal += [error]
dofsGlobal  += [uh.space.size]
print(count, ": size=", uh.space.gridView.size(0), "error=", error)
fem.globalRefine(1,uh)

WARNING (ignored): Could not open file 'alugrid.cfg', using default values 0 < [balance] < 1.2, partitioning method 'ALUGRID_SpaceFillingCurve(9)'.

You are using DUNE-ALUGrid, please don't forget to cite the paper:
Alkaemper, Dedner, Kloefkorn, Nolte. The DUNE-ALUGrid Module, 2016.


Created parallel ALUGrid<2,2,simplex,conforming> from input stream.

0 : size= 32 error= 0.15076788548358533
1 : size= 68 error= 0.1376904360341084
2 : size= 144 error= 0.10101929787844428
3 : size= 296 error= 0.09271236005573762
4 : size= 608 error= 0.06812963617257639
5 : size= 1232 error= 0.06261875789608652
6 : size= 2496 error= 0.04604780053712485
7 : size= 5024 error= 0.042351345493640294
8 : size= 10112 error= 0.03115388858628561
9 : size= 20288 error= 0.028661875392923558


Theory tells us that \begin{align*} \int_\Omega |\nabla(u-u_h)|^2 \leq \sum_K \eta_K, \end{align*} where on each element $$K$$ of the grid the local estimator is given by \begin{align*} \eta_K = h_K^2 \int_K |\triangle u_h|^2 + frac{1}{2}\sum_{S\subset \partial K} h_S \int_S [\nabla u_h]^2. \end{align*} Here $$[\cdot]$$ is the jump in normal direction over the edges of the grid.

We compute the elementwise indicator by defining a bilinear form \begin{align*} \eta(u,v) = \int_\Omega h^2 |\triangle u_h|^2 v + \int_{I_h} h_S [\nabla u_h]^2 \{v\}, \end{align*} where $$\{\cdot\}$$ is the average over the cell edges and $$[\cdot]$$ the jump. With $$h$$ and $$h_S$$ we denote local grid spacings and with $$I_h$$ the set of all facets in the grid. This bilinear form can be easily written in UFL and by using it to define a discrete operator $$L$$ from the second order Lagrange space into a space containing piecewise constant functions we have $$L[u_h]|_{K} = \eta_K$$.

[4]:

uh, exact, laplace = setup()

from dune.fem.space import finiteVolume as estimatorSpace
from dune.fem.operator import galerkin as estimatorOp

fvspace = estimatorSpace(uh.space.gridView)
estimate = fvspace.interpolate([0], name="estimate")

u = TrialFunction(uh.space.as_ufl())
v = TestFunction(fvspace)
hT = MaxCellEdgeLength(fvspace.cell())
he = MaxFacetEdgeLength(fvspace.cell())('+')
n = FacetNormal(fvspace.cell())
estimator_ufl = hT**2 * (div(grad(u)))**2 * v * dx +        he * inner(jump(grad(u)), n('+'))**2 * avg(v) * dS
estimator = estimatorOp(estimator_ufl)
tolerance = 0.05


Created parallel ALUGrid<2,2,simplex,conforming> from input stream.


WARNING (ignored): Could not open file 'alugrid.cfg', using default values 0 < [balance] < 1.2, partitioning method 'ALUGRID_SpaceFillingCurve(9)'.


Let us solve over a loop (solve,estimate,mark) and plot the solutions side by side.

[5]:

h1error = uflFunction(uh.space.gridView, name="h1error", order=uh.space.order, ufl=dot(grad(uh - exact), grad(uh - exact)))
fig = pyplot.figure(figsize=(10,10))
count = 0
errorVector    = []
estimateVector = []
dofs           = []
while True:
laplace.solve(target=uh)
if count%9 == 8:
plot(uh, figure=(fig, 131+count//9), colorbar=False)
error = numpy.sqrt(integrate(uh.space.gridView, h1error, 5))
estimator(uh, estimate)
eta = numpy.sqrt( sum(estimate.dofVector) )
dofs           += [uh.space.size]
errorVector    += [error]
estimateVector += [eta]
if count%3 == 2:
print(count, ": size=", uh.space.gridView.size(0), "estimate=", eta, "error=", error)
if eta < tolerance:
break
marked = fem.doerflerMark(estimate,0.6,layered=0.1)
fem.adapt(uh) # can also be a list or tuple of function to prolong/restrict
count += 1
plot(uh, figure=(fig, 131+2), colorbar=False)

2 : size= 38 estimate= 0.8627288159354651 error= 0.09255491144826526
5 : size= 60 estimate= 0.6093716434396423 error= 0.06487578048081612
8 : size= 74 estimate= 0.47671784411747725 error= 0.03909897004286079
11 : size= 112 estimate= 0.3431156830700155 error= 0.028194496294179996
14 : size= 174 estimate= 0.2031333515022477 error= 0.01625868861215105
17 : size= 316 estimate= 0.11427960643459709 error= 0.010442204473154349


Let’s have a look at the center of the domain:

[6]:

pyplot.close('all')
fig = pyplot.figure(figsize=(35,15))
plot(uh, figure=(fig, 131), xlim=(-0.5, 0.5), ylim=(-0.5, 0.5),
gridLines="white", colorbar={"shrink": 0.75}, linewidth=2)
plot(uh, figure=(fig, 132), xlim=(-0.1, 0.1), ylim=(-0.1, 0.1),
gridLines="white", colorbar={"shrink": 0.75}, linewidth=2)
plot(uh, figure=(fig, 133), xlim=(-0.02, 0.02), ylim=(-0.02, 0.02),
gridLines="white", colorbar={"shrink": 0.75}, linewidth=2)


We can also visualize the grid levels used

[7]:

fig = pyplot.figure(figsize=(35,15))
from dune.fem.function import levelFunction
levels = levelFunction(uh.space.gridView)
plot(levels, figure=(fig, 131), xlim=(-0.5, 0.5), ylim=(-0.5, 0.5),
gridLines="white", colorbar={"shrink": 0.75}, linewidth=0.5)
plot(levels, figure=(fig, 132), xlim=(-0.1, 0.1), ylim=(-0.1, 0.1),
gridLines="white", colorbar={"shrink": 0.75}, linewidth=0.5)
plot(levels, figure=(fig, 133), xlim=(-0.02, 0.02), ylim=(-0.02, 0.02),
gridLines="white", colorbar={"shrink": 0.75}, linewidth=0.5)


Finally, let us compare the globally refined solution and the adaptive one plotting number of degrees of freedom versus the error and the estimator

[8]:

pyplot.close('all')
pyplot.loglog(dofsGlobal,errorGlobal,label="H^1 error (global refine)")

[8]:

<matplotlib.legend.Legend at 0x7f1594e17880>