differential equations – PeriodicBoundaryConditions misses points with coarse mesh?

I’m noticing that some boundary points are “missed” by PeriodicBoundaryCondition in the case below. This points qualify in the sense that they satisfy the prescribed condition (predicate), but are nevertheless not identified by NDSolveValue.

The issue becomes a real problem only when I use a coarse mesh, and then no boundary points at all are found that satisfy the given condition.

To motivate this question, note that using a mesh with only a few points in one independent variable is useful for example when solving test cases with symmetries, or problems where the solution is believed to vary slowly in one dimension.

For the example, we solve Laplace’s equation in a cylinder. Using the command ToElementMesh(Cylinder()) is not useful here because the mesh is not coarse enough. So I turn to the package MeshTools. This can be installed in Mathematica 12.0+ with ResourceFunction("GitHubInstall")("c3m-labs", "MeshTools"). An “extruded” mesh is nice for my purposes, because I know exactly where all the mesh points are and can really precise control of my boundary conditions — or so I thought! Here’s a coarse cylindrical mesh:

N1 = 8;
DMesh = DiskMesh({0, 0}, 1, N1);
N2 = 4;
mesh = ExtrudeMesh(DMesh, 2 Pi, N2);

cylindrical mesh using MeshTools

Then I solve Laplace’s equation with Dirichlet boundary conditions.. The problem occurs when I do this using PeriodicBoundaryCondition to target the points with phi > Pi (using the fact that our solution will be symmetric about phi == Pi):

ReflectPhi = Function({#((1)), #((2)), 2 Pi - #((3))});
{uf} = NDSolveValue({Laplacian(u(x, y, phi), {x, y, phi}) == 0, 
   DirichletCondition(u(x, y, phi) == Cos(phi) y^2, Or(0 <= phi <= Pi, phi == 2 Pi)), 
   PeriodicBoundaryCondition(u(x, y, phi), Pi < phi < 2 Pi, ReflectPhi)}, {u}, Element({x, y, phi}, mesh))

The error produced is

NDSolveValue: No places were found on the boundary where (Pi)<phi<2
(Pi) was True, so PeriodicBoundaryCondition… will effectively be

And NDSolveValue fails. This example is a bit artificial because I could just do the whole boundary condition with a single DirichletCondition but it reveals the issue. As mentioned, the error is not encountered when I increase the resolution, e.g. N2 = 8. In that case, places were found where the predicate is satisfied — but now I wonder if Mathematica still “missing” some boundary points..?

To finish, let us just demonstrate there are points that satisfy the condition for N2 = 4. These are points with phi == 3 Pi/2.

bCoords = Map(mesh("Coordinates")((#)) &, DeleteDuplicates@Flatten@mesh("BoundaryElements")((1, 1)));
Cases(bCoords, _?(Pi < #((3)) < 2 Pi &)) // Length


Indeed there are 32 points corresponding to the boundary points of DMesh, which have been “extruded” to the phi == 3 Pi/2 location.