NDSolve with interior boundary

Question

The Dirichlet boundary condition of NDSolve, as advertised, should pick up any part of the boundary specified in the "pred" argument and apply the condition specified in the "beqn" argument. But I'm finding that if the boundary has slots, as shown here, it ignores the internal boundaries. Why is it doing this? I saw an answer that involved manually building a mesh, but why should I need to do this?

Specify geometry parameters

x1 = 5; y1 = 1; y2 = 2; y3 = 10; x2 = 10; ϵ = 0.1;

Specify region

region = 
  ImplicitRegion[
    (y - (y1 + (y2 - y1) UnitStep[x - x1]))^2 >= ϵ^2, {{x, 0, x2}, {y, 0, y3}}]

Solve Laplace's equation (cylindrical coordinates) and contour plot. It sees only the boundary conditions on the top and bottom:

slpolar = 
  NDSolve[
    {D[x D[u[x, y], x], x]/x + D[u[x, y], {y, 2}] == 0, 
     u[x, 0] == 30, u[x, y3] == 0, 
     DirichletCondition[u[x, y] == 0, 0 < x < x2 && y > 0]}, 
    u, {x, y} ∈ region]

Answer 1

Thanks. I also found that the problem goes away if I make the slots wide enough. In this example the slots had to be more than exactly 4/3 wide.

Answer 2

ToElementMesh (and DiscretizeRegion) have problems generating the mesh from the region you specified. Using a different specification helps to work around this:

region = RegionDifference[
  RegionDifference[Rectangle[{0, 0}, {10, 10}], 
   Rectangle[{0, 1 - \[Epsilon]}, {5, 1 + \[Epsilon]}]], 
  Rectangle[{5, 2 - \[Epsilon]}, {10, 2 + \[Epsilon]}]]

You can generate the mesh and look at it:

Needs["NDSolve`FEM`"]
m = ToElementMesh[region];
(*m["Wireframe"]*)

And here is the solution:

slpolar = 
  NDSolveValue[{D[x D[u[x, y], x], x]/x + D[u[x, y], {y, 2}] == 0, 
    u[x, 0] == 30, u[x, y3] == 0, 
    DirichletCondition[u[x, y] == 0, 0 < x < x2 && y > 0]}, 
   u, {x, y} \[Element] m];
ContourPlot[slpolar[x, y], {x, y} \[Element] m, PlotRange -> All, 
 ColorFunction -> "TemperatureMap"]