Problem with internal facet (MeshFunction) in adaptive solution

Dear all,

I’m facing an issue with adaptive meshing for a nonlinear problem with a Neumann boundary condition, defined on an internal facet. Apparently, there is some trouble with the MeshFunction in the adaptive procedure. I’d be really glad if you could help me with this problem.

Here’s a minimal example code:

from dolfin import *
bounds = [-20., 20.]
class Middle(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.)

mesh = IntervalMesh(200, bounds[0], bounds[1])
V = FunctionSpace(mesh, "CG", 1)
phi = TrialFunction(V)
v = TestFunction(V)
u_ = Function(V)
middle = Middle()
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
middle.mark(boundaries, 1)
dS = Measure("dS", domain=mesh, subdomain_data=boundaries)

Neumann_bc = (exp(-2*phi+1.) - exp(phi))
a = inner(grad(phi), grad(v))*dx
L = (exp(-2*phi)-exp(phi))*v*dx + Neumann_bc("+")*v("+")*dS(1) + Neumann_bc("-")*v("-")*dS(1)

F = a-L
F = action(F, u_)
J = derivative(F, u_, phi)
goal_functional = u_*dx()
problem = NonlinearVariationalProblem(F, u_, [], J)
solver = AdaptiveNonlinearVariationalSolver(problem, goal_functional)
solver.solve(1e-10)

And this is the error message that I get:
*** Error: Unable to adapt mesh function.
*** Reason: Unable to extract information about parent mesh entities.
*** Where: This error was encountered inside adapt.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.1.0
*** Git changeset: unknown

Thank you in advance!

Out of curiosity, what is the “internal” boundary condition mathematically?

As far as I can tell, Neumann_bc is continuous (as phi is continuous), and v is continuous, this term should be zero, or am I missing something?

Physically, the Neumann boundary condition shall represent a surface charge density, which results in a kink in the electric potential.

To be honest, I can’t really answer your question precisely — I will try anyway: phi and v are continuous, that’s right. I think the point is that phi does have a finite value at the location of the internal facet. In the variational formulation, the surface charge density (Neumann_bc) kind of “compensates” for the boundary terms of the form nabla phi * v. The discontinuous character of the boundary condition is represented by the finite value of the surface charge density. The charge density (per volume) is, indeed, something like a Dirac delta function in that point; but this distribution is already “integrated out” in my formulation to give this finite surface density — in any other point but this particular one, the surface charge density would be zero…

Does that help?

I think I get your point, but to be certain, could you write down the strong form of the boundary condition?

I played around with your problem earlier, and it seemed like the exp-function is what is causing the issue. Im am not certain why yet, but if you supply the strong form I’ll have another go.

1 Like

Sure! So the problem is
-\frac{\partial^2 \phi}{\partial x^2} = e^{-2\phi}-e^{\phi}

With the Neumann boundary conditions, defined in the middle (x=0):
\left.\frac{\partial\phi}{\partial x}\right|_{x=+0} = -\frac12\left(e^{1-2\phi}-e^\phi\right)
\left.\frac{\partial\phi}{\partial x}\right|_{x=-0} = \frac12\left(e^{1-2\phi}-e^\phi\right)

And with the Dirichlet boundary conditions
\phi(-20) = 0
\phi(20) = 0

Sorry for the very informal way of writing this down, I am not trained in this field…
by “+0” and “-0” I mean the sides to the right/left of the internal facet.

BTW: Solving this problem with a NonlinearVariationalSolver does not fail, so to me it appears to be specifically an issue of the adaptive procedure.

Thanks a lot for your help!

Seems like your variational formulation is ok (up to a constant, i think it each term should be multiplied by 0.5).
I’ve been able to reproduce your issue without any adaptive problem:

from dolfin import *
bounds = [-20., 20.]
class Middle(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.)

mesh = IntervalMesh(200, bounds[0], bounds[1])
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
Middle().mark(boundaries, 1)

mesh = adapt(mesh)
md = mesh.data()

mf = adapt(boundaries, mesh)

I’ll see if I can get to the root of the issue.

1 Like

Unfortunately, parent facets only seems to be implemented for 2D:


Thus changing your problem to:
from dolfin import *
parameters["refinement_algorithm"] = "plaza_with_parent_facets";

bounds = [-20., 20.]
class Middle(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.)

mesh = RectangleMesh(Point(bounds[0], bounds[0]),Point(bounds[1], bounds[1]), 40, 40)

would work. I do not know how to extend the parent_facet refinement to 1D ;(

1 Like

Thanks a lot for your help! You’re right about the factor 0.5; I thought I was sure that this factor was inherent to the specifiers ("+" and “-”), but now I checked it again, on a simpler problem, and apparently I got that wrong…

For the problem with the MeshFunction, I found a workaround by implementing an adaptive procedure like it is done in the adaptive-Poisson demo on Bitbucket:

So now, I can “rebuild” the MeshFunction in each refinement iteration. I also have to rebuild the variational problem, including function space, function, trial function, test function and Dirichlet boundary conditions. It’s not elegant, but it works. And apparently, there is no straightforward way to do it in 1D for now, if I got that right?

1 Like

Your sketched non-elegant solution is the only solution as far as I can tell. Note that in general the plus and minus restriction is not consistent, see: Integrating over an interior surface for more details. In your case it doesn’t matter as the boundary condition on each side ends up being the same.