I am looking for a mathematical explanation for the singularity caused by a Dirichlet boundary condition partially imposed at a boundary.

For instance $$ \nabla^2u=0 ~ \text{in}~\Omega $$ where $\Omega$ is a rectangle of 10 by 4. Dirichlet boundary condition $u=0$ is imposed at $ y=0 ~\text{and}~ y = 4 ~\text{and}~ x < 2 $ and a flux is applied at the boundary $x=10$

Plotting the flux, you can see the singularity where the Dirichlet boundary condition ends

enter image description here

I plotted this with the following FEniCS code

from fenics import *

mesh = RectangleMesh(Point(0, 0), Point(10, 4), 400, 160, diagonal='crossed')

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

a = inner(grad(u), grad(v))*dx

def MyDirichlet(x, on_boundary):
    return on_boundary and (x[0] < 2.0) and (x[0] > 0.0)

class MyNeumann(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 10.0)

subdomain = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
myneumann = MyNeumann()
myneumann.mark(subdomain, 1)
ds = Measure('ds', domain=mesh, subdomain_data=subdomain)

bc = DirichletBC(V, Constant(0.0), MyDirichlet)

L = Constant(1.0)*v*ds(1)

u_sol = Function(V)

solve(a==L, u_sol, bcs=bc)
File("sing_example.pvd") << u_sol
Vflux = VectorFunctionSpace(mesh, 'CG', 1)
File("sing_flux.pvd") << project(grad(u_sol), Vflux)

  • $\begingroup$ Just out of curiosity, is there any particular reason to use this weird configuration of boundary conditions for Laplace equation? $\endgroup$ – Alone Programmer Jan 3 '20 at 0:56
  • $\begingroup$ No specific reason. $\endgroup$ – balborian Jan 3 '20 at 3:10
  • $\begingroup$ So, is this for research or right now just for learning and then moving on to solving a real problem in the future? $\endgroup$ – Alone Programmer Jan 3 '20 at 3:37
  • $\begingroup$ Both, the real problem is a convection diffusion problem where the advection field enters the right side of the rectangle and exists through the left side. The idea of placing the Dirichlet boundary conditions this way is to keep them away from the right side (the inlet). I have not thought of an alternative way to pose this real problem, but in any case I would like to know why this phenomenon happens. $\endgroup$ – balborian Jan 3 '20 at 15:59
  • 2
    $\begingroup$ As far as I know, you need to have at least some sort of boundary condition specified on every part of the boundary. Since you do not specify any conditions on $x=0$ or $(x>2)(x<10)$ I have no idea what Fenics will do. It probably depends on the details of solve. $\endgroup$ – knl Feb 1 '20 at 21:15

It looks to me that you don't have a well posed problem, so who knows how the fenics is handling that. I know that you can get some odd results for ill-posed poisson problems when the solver doesn't diverge, but instead returns garbage. This seems to be on of those cases. Relatedly you can have certain fixed point iterations that solve linear systems of singular matrices without diverging, but any inverse of a singular matrix is clearly meaningless.

  • $\begingroup$ It’s not an answer really... It could be just a comment probably. $\endgroup$ – Alone Programmer Jan 3 '20 at 0:57
  • 2
    $\begingroup$ Fair. I figured it does sort of answer his question as it points out that his problem is Ill posed and the details of how it handles that are up to the software and we can't answer it. $\endgroup$ – EMP Jan 3 '20 at 6:47

As others have mentioned, the boundary-conditions are not clear. But such singularities in the gradient do occur in the neighborhood of a discontinuity in boundary conditions. Here's a mathematical explanation. Take the following prototypical boundary condition: $u = 0$ for $y = 0, \ x > 0$ and $u = \pi$ for $y=0, \ x < 0$. Then a solution of the Laplace equation for $y\geq 0$, which obeys the boundary-conditions and is bounded at infinity is simply $u = \phi$, where $\phi$ is the angle from the positive $x$ axis. The gradient is then $\nabla u = \frac{1}{\rho}\mathbf{1}_\phi$, where $\rho$ is the distance from the origin and $\mathbf{1}_\phi$ is a unit-vector in the $\phi$ direction, and it is singular near the origin. In this example, the flux is normal to the $x$ axis boundary, so either one of the Dirichlet conditions (but not both) may be readily replaced by a Neumann condition without altering the result.

To put it in more intuitive terms, the potential at $y = 0^+$ must jump discontinuously between two different values, therefore the gradient must be infinite near the jump.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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