EDIT: I reported this to Mathematica support. Will update later.

I just posted on this issue, but I found a more elementary example which I think will make the issue more transparent.

Apparently, some boundary points are ignored/missed when solving a PDE using FEM, at least for the following case.

Here is almost the simplest mesh you can imagine. We start by defining some functions to set up a first order quad mesh of the unit square.

MakeCoords[Nx_, Ny_] := Module[{i, j}, Flatten[Table[N@{i/Nx, j/Ny}, {j, 0, Ny}, {i, 0, Nx}], 1]]
MakeTuples[Nx_, Ny_] := Module[{i, j, i1, i2, i3, i4, if},
  if[i_, j_] := i + (j - 1) (Nx + 1);
   Table[i1 = if[i, j]; i2 = if[i + 1, j]; i3 = if[i + 1, j + 1];i4 = if[i, j + 1];
    {i1, i2, i3, i4}, {j, 1, Ny}, {i, 1, Nx}], 1]

Here then is a 4x2 mesh:

ONx = 4; ONy = 2;
meshO = ToElementMesh["Coordinates" -> MakeCoords[ONx, ONy], "MeshElements" -> {QuadElement[MakeTuples[ONx, ONy]]}];

enter image description here

We run into trouble as before by attempting to solve Laplace's equation with Dirichlet boundary conditions, where part of this is enforced using PeriodicBoundaryCondition:

{uf} = NDSolveValue[{Laplacian[u[x, y], {x, y}] == 0, 
DirichletCondition[u[x, y] == (x - 1/2)^2, Or[x == 1, x <= 0.5]], 
PeriodicBoundaryCondition[u[x, y], 0.5 < x < 1, {1 - #[[1]], #[[2]]} &]}, {u}, Element[{x, y}, meshO]]

NDSolveValue fails and complains:

NDSolveValue: No places were found on the boundary where 0.5<x<1 was True, so PeriodicBoundaryCondition[u,0.5<x<1,{1-#1[1],#1[2]}&] will effectively be ignored

Mathematica says that the predicate is not satisfied on any boundary points. But as we know, there are exactly two boundary points that satisfy 0.5 < x < 1, namely (0.75, 0) and (0.75, 1). For some reason, there is trouble with this specification of boundary conditions. If a finer mesh is used, the error goes away, but does the problem itself? Are points on the boundary lost?

Any ideas? If one needs to implement mixed boundary conditions involving some PeriodicBoundaryConditions, is there a way to do it, to avoid this potential problem?

Here is another example that might be related.

meshO = ToElementMesh[ImplicitRegion[True, {{x, 0, 1}, {y, 0, 1}}],"MaxCellMeasure" -> 0.5];

simple unit square mesh with triangle elements

{uf} = NDSolveValue[{Laplacian[u[x, y], {x, y}] == 0, 
   DirichletCondition[u[x, y] == (x - 1/2)^2, True]}, {u}, Element[{x, y}, meshO]]
Plot3D[uf[x, y], Element[{x, y}, meshO]]

plot with Dirichlet

{uf} = NDSolveValue[{Laplacian[u[x, y], {x, y}] == 0, 
   DirichletCondition[u[x, y] == (x - 1/2)^2, x <= 0.5], 
   PeriodicBoundaryCondition[u[x, y], x > 0.5, {1 - #[[1]], #[[2]]} &]}, {u}, Element[{x, y}, meshO]]
Plot3D[uf[x, y], Element[{x, y}, meshO]]

plot using method with periodic bcs

Like the case with the custom mesh, you can clearly see something is going wrong with implementing the periodic boundary conditions for the these "intermediate" boundary points (0.75, 0), (0.75, 1), etc. Maybe it's related...

  • 1
    $\begingroup$ OK, I looked at the code. This is not a bug per se. It behaves as expected. PBC are similar to NeumannValues in that they operate on edges in 2D and surfaces in 3D. So the single mesh point in not enough for the PBC to work. While it may be workable for this specific case I am fairly confident that when I designed this that was a pretty good reason that PBCs behave like NeumannValues. Now, what is a bug is the lousy error message. I will make sure this information gets into the documentation such that your effort was not in vain! So, it's a bug but more of a documentation issue. $\endgroup$
    – user21
    Sep 22, 2020 at 14:34
  • $\begingroup$ And thanks a lot for reporting it. $\endgroup$
    – user21
    Sep 22, 2020 at 14:34
  • $\begingroup$ And thanks to you for looking into it. I had a worry about how PBC might work with so few points, but I figured at least it should behave OK if the boundary points are symmetric. By symmetric I mean an exact correspondence between source and target points, such that DiscretizeBoundaryConditions produces a DirichletMatrix with precisely two entries on each row, one +1 and the other -1. Then deploy might function as intended even with very few points, as my example below. Do you think this makes sense? $\endgroup$
    – Will.Mo
    Sep 23, 2020 at 4:21

1 Answer 1


I'm afraid most people won't be interested in this issue, because it may seem obscure, and I suspect the neglect of a few points in the boundary may not have a very big effect for well-conditioned problems.

However, I believe it is potentially useful to work this out, for problems where precise control over BCs is needed. Hope it will help some people.

I found a workaround for this issue, although I do not know if it will always work, and I think there might need to be some official "fix" from Mathematica.

To recap, we want to enforce mixed boundary conditions including Dirichlet and periodic conditions, but some boundary points are missed when the desired BCs are discretized (during the call to DiscretizeBoundaryConditions).

One clue about what is happening: note that if we simplify the conditions slightly, so the Periodic boundary condition is inclusive of the upper bound, then it works fine and all the correct boundary coordinates are identified:

{uf} = NDSolveValue[{Laplacian[u[x, y], {x, y}] == 0, 
   DirichletCondition[u[x, y] == (x - 1/2)^2, x <= 0.5], 
   PeriodicBoundaryCondition[u[x, y], 
    x > 0.5, {1 - #[[1]], #[[2]]} &]}, {u}, Element[{x, y}, meshO]]

Note how DirichletCondition only targets x <= 0.5, while PeriodicBoundaryCondition includes all x > 0.5, including x == 1. Although it is an equivalent problem, it is not the way we want to solve it -- the point was to be able to choose the predicates freely, which is needed for more difficult problems. But the success of this gives a hint that the problem occurs when PeriodicBoundaryCondition is dealing with exclusive intervals, e.g. 0.5 < x < 1. It could not find the x == 0.75 point in that case.

So to work around this behavior, we can do the boundary conditions in two separate steps and combine them at the end. Here is the mesh we want to work with:

ONx = 4; ONy = 2;
meshO = ToElementMesh["Coordinates" -> MakeCoords[ONx, ONy], 
   "MeshElements" -> {QuadElement[MakeTuples[ONx, ONy]]}];

Here are separated boundary conditions (yes, the periodic BCs include x==1 but we will trim off the extra points later manually):

DirichletFcn[x_, y_] := (x - 1/2)^2
bcD = {DirichletCondition[u[x, y] == DirichletFcn[x, y], 
    Or[x == 1, x <= 0.5]]};
bcP = {PeriodicBoundaryCondition[u[x, y], 
    0.5 < x <= 1, {1 - #[[1]], #[[2]]} &]};

We use FEM programming to continue.

vd = NDSolve`VariableData[{"DependentVariables", 
     "Space"} -> {{u}, {x, y}}];
sd = NDSolve`SolutionData[{"Space" -> ToNumericalRegion[meshO]}];
dofd = 1; dofi = 2;
Cu = Table[
   DiscreteDelta[k - l], {i, dofd}, {j, dofd}, {k, dofi}, {l, dofi}];
coefficients = {"DiffusionCoefficients" -> Cu};
initCoeffs = InitializePDECoefficients[vd, sd, coefficients];
initBCsD = InitializeBoundaryConditions[vd, sd, bcD] ;
initBCsP = InitializeBoundaryConditions[vd, sd, bcP] ;

These steps are all well documented, but we are doing two calls to InitializeBoundaryConditions instead of the usual one. Also note that the final command produces a warning from Mathematica about lack of Dirichlet conditions and non-uniqueness. We are not worried about that; it will be well-posed when we assemble all the BCs together in the end. Continuing:

methodData = 
  InitializePDEMethodData[vd, sd, Method -> {"FiniteElement"}];
discretePDE = DiscretizePDE[initCoeffs, methodData, sd];
{load, stiffness, damping, mass} = discretePDE["SystemMatrices"];
discreteBCsD = 
  DiscretizeBoundaryConditions[initBCsD, methodData, sd];
discreteBCsP = DiscretizeBoundaryConditions[initBCsP, methodData, sd];

Again, there are two calls to DiscretizeBoundaryConditions; normally there is only one. We now have the two BCs in two separate DiscretizedBoundaryConditionData objects, and we can combine them. The problem is that the periodic boundary conditions as we defined them conflict with the Dirichlet condition -- they both include all the x==1 boundary points. Our strategy is to defer to the Dirichlet condition, wherever a conflict occurs. Then we will have succeeded in implementing our specific BCs.

Continuing, we will have to extract the part of the periodic BCs that targets points not present in the Dirichlet condition. These points (indexed 4 and 14 as can be found by inspecting meshO["Coordinates"]) can be visualized as follows:

DirichletCoords = 
  Map[meshO["Coordinates"][[#]] &, discreteBCsP["DirichletRows"]];
KeepCoords = Map[meshO["Coordinates"][[#]] &, {4, 14}];
 Graphics[{PointSize[Large], Red, Point[DirichletCoords]}], 
 Graphics[{PointSize[Large], Blue, Point[KeepCoords]}]]

nodes to keep and nodes to discard

We want to keep the blue ones and discard the red ones. This is done with the following code. First we populate all the discrete BCs data from the automatically generated Dirichlet data:

diriMat = discreteBCsD["DirichletMatrix"];
diriRows = discreteBCsD["DirichletRows"];
diriVals = discreteBCsD["DirichletValues"];
dof = Length[meshO["Coordinates"]];

Then we will add to this data the non-conflicting part of the periodic BC data:

CdiriRows = discreteBCsP["DirichletRows"];(* "candidate DiriRows" *)
CdiriMat = discreteBCsP["DirichletMatrix"];
CdiriVals = discreteBCsP["DirichletValues"];
For[i = 1, i <= Length@CdiriRows, i++,
 If[Not[MemberQ[diriRows, CdiriRows[[i]]]],
   AppendTo[diriRows, CdiriRows[[i]]];
   AppendTo[diriMat, CdiriMat[[i]]];
   AppendTo[diriVals, CdiriVals[[i]]];

Now we define a new DiscretizedBoundaryConditionData object:

lmdof = Length@
discreteBCs = 
 DiscretizedBoundaryConditionData[{SparseArray[{}, {dof, 1}], 
   SparseArray[{}, {dof, dof}], diriMat, diriRows, 
   diriVals, {dof, 0, lmdof}}, 1];

This is the hacked discretized BC data. It is just the Dirichlet data with extra rows in the matrices coming from the periodic boundary condition data, were the target is not present in the list of Dirichlet targets, discreteBCsD["DirichletRows"].

The rest is just the usual steps:

DeployBoundaryConditions[{load, stiffness}, discreteBCs];
solution = LinearSolve[stiffness, load];
NDSolve`SetSolutionDataComponent[sd, "DependentVariables", 
{uf} = ProcessPDESolutions[methodData, sd];
Plot3D[uf[x, y], Element[{x, y}, meshO]]

visualize the solution


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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