I'm studying the dealii tutorial number 4,5 and I understand the workflow. I've also been able to find the EOC by using manufactured solution where $f$ is a smooth r.h.s. and $\alpha(x)$ smooth too.

Now I'm interested in checking the EOC when $\alpha(x)$ is discontinuous (link to Doxygen docs of dealii)

Let's say I have $\Omega = [-1,1]\times[-1,1]$ and I assume $u(x,y)=\sin(\pi x ) \sin(\pi y)$ to be the solution of

$$- \nabla \cdot (\alpha(x) \nabla{u}) = f$$ with homogeneous Dirichlet b.cs and $\alpha(x)$ defined as in the link. The forcing term I obtain is:

$$f(x,y)= \begin{cases} 20 \cdot 2\pi^2 \sin(\pi x)\cos(\pi y) \\ 2 \pi^2 \sin(\pi x)\cos(\pi y) \end{cases}$$

Looking at the EOC, it seems that my numerical solution is not converging to the true solution $u(x,y)$ (I obtain an EOC which decreases from 0.23 to 0.02...) So I'm wondering if this approach is not okay.

  • How can I construct a manufactured solution when the $\alpha(x)$ is not continuous?


Using that $\alpha(x)$, I could set $\nabla{u}=[1,20]^T$ on $B_{1/2}(0,0)$ and $\nabla{u}=[20,40]^T$ outside $B_{1/2}(0,0)$.

This implies that $\alpha(x,y) \nabla{u}$ is continuous, actually it's constant, hence I obtain $f=0$. Also, $$u(x,y)= x+20y$$ in $B_{1/2}(0,0)$, and $$u(x,y)=20x+40y$$ outside $B_{1/2}(0,0)$. With this argument, I always have a jump in $\nabla{u}$ due to $\alpha(x,y)$.

However, I'm still facing issues in showing the convergence, and I think this may be due to the fact that I have $f=0$ as forcing term.


Following the comment from @MaximUmansky, let $\alpha \nabla{u} = \vec{f}$ such that $$-\operatorname{div}(\vec{f}) = f_x + f_y= \sin(x)+\cos(y)$$

so $\vec{f}=[\cos(x), - \sin(y)]^T$

With $\alpha$ as above, we need $\nabla{u} = [\frac{\cos(x)}{20}, -\frac{\sin(y)}{20}]$ in $B_{1/2}(0,0)$

and $\nabla{u} = [\cos(x),-\sin(y)]$ outside $B_{1/2}(0,0)$.

So I have, in $B_{1/2}(0,0)$, $$u(x,y)=\frac{1}{20}(\sin(x)+\cos(y))$$

and $$u(x,y)=\sin(x) +\cos(y)$$ outside $B_{1/2}(0,0)$.

  • $\begingroup$ grad(u) has to be discontinuous to make the product $\alpha$ grad(u) smooth. For example, in 1D make $\alpha$=1 in x$\in [0,0.5]$ and $\alpha=2$ in x$\in$ [0.5,1], and $\partial_x$u=2 in the left half-domain, $\partial_x$u=1 in the right half-domain. $\endgroup$ Apr 10, 2021 at 14:23
  • $\begingroup$ @MaximUmansky Thanks, I see now why $\alpha \nabla u$ needs to be continuous (I need to take its divergence), but I can't see how I can do this in practice. I mean, $u$ is what I have to find with my solver! What should I change in my problem above? $\endgroup$
    – FEGirl
    Apr 10, 2021 at 15:14
  • $\begingroup$ I think I got your point: $u_x = 2$ implies $u(x)=2x$ on the left-half domain, while on the right-half we have $u(x)=x$. I'm trying to work this out in the 2D case @MaximUmansky $\endgroup$
    – FEGirl
    Apr 10, 2021 at 15:42
  • $\begingroup$ @MaximUmansky I worked out a possible way in my edit, but the r.h.s $f$ turns out to be $0$ now, and I'm having some problems in showing the convergence. Do you have any hint, maybe another simple solution with a non-trivial $f$? $\endgroup$
    – FEGirl
    Apr 10, 2021 at 16:30
  • $\begingroup$ I mean, I'd like to find a discontinuous $f$ as r.h.s. @MaximUmansky $\endgroup$
    – FEGirl
    Apr 10, 2021 at 17:57

2 Answers 2


@cfdlab's answer gives a general way to construct solutions for discontinuous coefficients, but here is a slightly more theoretical perspective to it as well. All of this can be understood in 1d, so imagine all of the functions below to be functions of just one argument $x$ for a moment.

First, you can't just choose any $u$ and $\alpha$ to obtain a function $f(x)$. That's because in general, if $u$ is continuous and $\alpha$ is discontinuous, then $\alpha \nabla u = \alpha\frac{du}{dx}$ is a discontinuous function and if you take the negative divergence of it (we're in 1d, so the divergence is again just $\frac{d}{dx}$), you end up with a function $f$ that is not, in fact, a function but a distribution: It has delta functions in all of those places where $\alpha \nabla u$ is discontinuous. One can, in theory, solve the PDE you are interested in with such right hand sides $f$, but not easily with finite elements and for sure not if you use quadrature, because the way you integrate the right hand side terms doesn't see these delta functions.

Secondly, the Poisson equation models some physical behavior. In general, this would be a diffusive process of some sort where $j=-\alpha\nabla u$ is a flux and the divergence of the flux equals the source terms $f$. Fluxes are generally thought to be continuous functions, though not necessarily smooth: they can have kinks, and so their derivatives can be discontinuous. That means that if $\alpha$ is discontinuous, then $\nabla u$ must be discontinuous as well, in such a way that $-\alpha\nabla u$ is continuous again. What this implies is that you can't just choose $\alpha$ and $u$ independently. At its heart, that is why @cfdlab in their answer above starts by prescribing $\alpha$ and the flux (rather than $\alpha$ and $u$) and then deriving what $u$ needs to be.


If you want to use the coefficient of the form $$ \alpha = \begin{cases} \alpha_1 & r < 0.5 \\ \alpha_2 & r > 0.5 \end{cases} $$ it is useful to work in polar coordinates. You build two pieces of the solution $$ u_1(r,\theta), \qquad r < 0.5 $$ and $$ u_2(r,\theta), \qquad r > 0.5 $$ Then at the interface you ensure solution and flux continuity $$ u_1(0.5, \theta) = u_2(0.5, \theta), \qquad \alpha_1 \partial_r u_1(0.5,\theta) = \alpha_2 \partial_r u_2(0.5,\theta) $$ Maybe you can try to build a solution of the form $$ u_i(r,\theta) = f_i(r) g(\theta), \qquad i=1,2 $$ You can first try to build a purely radial solution, which should be easy to do, it will be just like a 1d problem.

Note that the discontinuity is on a circle and to observe good convergence you must make a mesh that aligns with this circle.

If you know fenics, you can see below example where the mesh is made in gmsh https://github.com/cpraveen/fembook/tree/master/fenics/step-10

  • $\begingroup$ Thanks @cfdlab for your answer. I tried to follow your argument, and tried to build $u_1$, $u_2$, by setting $$\partial_r ( \alpha \partial_r u)=r$$ with $\alpha$ discontinuous as in my link, and I obtained $u^1_r(r)=\frac{1}{40} r^2$ and $u^2_r(r)=\frac{r^2}{2}$. Imposing the continuity at the interface I found $u(r)=\frac{r^3}{120} + \frac{13}{320}$ and $u_2(r)=\frac{r^3}{3}$. You can see that it's continuous at $r=0.5$, but I'm a bit afraid that maybe I should have written the original equation in polar coordinates. What do you think? $\endgroup$
    – FEGirl
    Apr 11, 2021 at 23:07
  • $\begingroup$ Yes, write PDE in polar coordinates, find a solution and then write solution in Cartesian coordinates. $\endgroup$
    – cfdlab
    Apr 12, 2021 at 4:31
  • $\begingroup$ Okay, so if I look at a purely radial solution I may consider $$- \frac{1}{r} \partial_r ( \alpha r \partial_r u) = 1$$ By setting $\alpha r \partial_r u= f$ I got $ \partial_r f=-r$, and finally $u_1(r)=- \frac{r^2}{80} - \frac{19}{320}$ and $u_2(r)=- \frac{r^2}{4}$, where the constant has been chosen to ensure continuity at $r=0.5$. Of course now I'll use the fact $r = \sqrt{x^2+y^2}$. Is this what you meant? $\endgroup$
    – FEGirl
    Apr 12, 2021 at 9:15
  • $\begingroup$ Looks correct. This is a quadratic solution which an FEM can compute exactly. Try a more nonlinear solution. $\endgroup$
    – cfdlab
    Apr 12, 2021 at 11:16
  • $\begingroup$ You mean to start with something like $\alpha r \partial_r u = \sin(r)$ and then get $u$ ? @cfdlab $\endgroup$
    – FEGirl
    Apr 12, 2021 at 13:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.