White spots in the plot

I/O
1

Hello everybody,
since I’m pretty new to Fenics I’m trying to figure out subdomains and Neumann Conditions.
Doing so I came across a strange effect in my Plot of a pretty basic Laplace equation.

Screenshot from 2022-04-07 21-05-09

The above plot results from the following example code

from dolfin import *
from fenics import *
from mshr import *
import matplotlib.pyplot as plt

g = Constant(5000320.0)

space = Rectangle(Point(0,0), Point(1,1))
cylinder = Circle(Point(0.5, 0.5), 0.01)
domain = space - cylinder

mesh = generate_mesh(domain, 100 )
 
Q = FunctionSpace(mesh, 'Lagrange', 1)

u  = TrialFunction(Q)
v  = TestFunction(Q)

aL  = inner(grad(u), grad(v))*dx
LL = g * v *ds

u = Function(Q)

solve(aL == LL, u) 

plot(mesh, title='the mesh')
plt.show()

Laplace=plot(u, title='show laplace')
plt.colorbar(Laplace)
plt.show()

As you can see there are these white spots in the plot and I have no idea what they are from as the mesh covers these areas and there are no error messages.
So what are they and from what do they result?

Addtionally the Neumann condition doesnt seem to be applied right since the results seem to be uneven even with them applied to all edges.

As a note g is intentionally so large, both to show the effect and for the program this all stems from.

Thank you for your help!

2

I don’t think your problem is well posed. You have insufficient data to yield a unique solution. You’re essentially plotting arbitrary numbers.

For a singular system see, for example, this demo.

3

Thank you for the fast answer and sorry that the reply took so long.
Anyways it got me looking in how to solve the problem but sadly I don’t normally work with Differential Equations so my “solutions” might be faulty.

So experimentally I changed g from a constant to g=Expression("sin(x[0])", degree=2) similar to the demo and added an f*v*dx with f=0 to the LL so that I still get a Lapalce equation. The result still got blank spots, even if it should give the code as much data to solve as the equation in the demo and if I ain’t wrong it should even reduce the amount of unknowns that need to be solved compared to the demo.

Thinking about it since the Neumann Condition is applied to all edges it should provide Border Conditions all around similar to the ft01_poisson.py tutorial granted that one got Dirichlet Conditions on the edges but since it doesn’t solve the problem I need to look further.

Theoretically I got a DirichletBC for that problem as well but there are some problems with implementing that I need to look further into or open up another topic about.

Sadly I still don’t know what the white spots are from since if the program giving out arbitrary numbers that get plotted why are a few filled in and not all or why aren’t all NaN?

4

To expand on @nate’s answer, there are two basic well-posedness issues that come up with pure Neumann BCs. The first is compatibility of the Neumann data with the volume source term. If you have the interior PDE

-\nabla\cdot(\nabla u) = f\text{ ,}

integrating both sides over \Omega and applying the divergence theorem on the left implies that

\int_{\partial\Omega}-\nabla u\cdot\mathbf{n} \,d\partial\Omega = \int_\Omega f\,d\Omega\text{ .}

If you also have the Neumann BC of \nabla u\cdot\mathbf{n} = g, then you need

\int_{\partial\Omega}-g\,d\partial\Omega = \int_\Omega f\,d\Omega\text{ ,}

i.e., you cannot choose g and f independently from each other. Physically, this is just saying that the flux through the boundary (g) must be balanced by the rate of production on the interior (f) to have a steady equilibrium solution. The second issue is that, even with compatible problem data, you can arbitrarily add a constant to u and get another solution, since only derivatives of u appear in both the PDE and the BC. This arbitrary constant corresponds physically to the initial distribution of u in \Omega, before the problem reached a steady state. These two issues correspond to the two classes of singular linear problems: in the first case, there are no solutions, and in the second, there are infinitely-many solutions. In either case, you’re not guaranteed to get meaningful output from the linear solver.

The demo that @nate linked eliminates the singularity at the linear algebraic level. You may find it easier to physically interpret the demo here, which introduces an auxiliary constraint to make u unique, enforced by a Lagrange multiplier that corrects f by a constant, to be compatible with g.

3 Likes