# Understanding PeriodicBoundaryConditions

Every thing works fine in a simple example with periodic boundary condition u[ 2,y]==u[0,y] from documentation of PeriodicBoundaryConditions

Ω = Rectangle[{0, 0}, {2, 1}];
pde = -Laplacian[u[x, y], {x, y}] ==If[1.25 <= x <= 1.75 && 0.25 <= y <= 0.5,1., 0.];
ΓD =DirichletCondition[u[x, y] == 0, (y == 0 || y == 1) && 0 < x < 2];

pbc = PeriodicBoundaryCondition[u[x, y], x == 0,TranslationTransform[{  2, 0}]];
ufun = NDSolveValue[{pde, pbc, ΓD},u, {x, y} ∈ Ω];
ContourPlot[ufun[x, y], {x, y} ∈ Ω,ColorFunction -> "TemperatureMap", AspectRatio -> Automatic]


But if I modify the periodic boundary conditions slightly from x==0, translation +2 to x==2,translation -2, expecting the same result(!)

pbc = PeriodicBoundaryCondition[u[x, y], x == 2,TranslationTransform[{  -2, 0}]];
ufun = NDSolveValue[{pde, pbc, ΓD},u, {x, y} ∈ Ω];
ContourPlot[ufun[x, y], {x, y} ∈ Ω,ColorFunction -> "TemperatureMap", AspectRatio -> Automatic]


the solution changes significantly!

What's wrong here(Mathematica v11.0.1)?

Thanks!

• Mathematicaly, The periodic boundary condition u[2,y]==u[0,y] is not enough to have unicity of the solution. For having the unicity , it would be sufficient to add one more condition such as flux conservation between the left and the right side (you can see this as a thermal problem on a ribbon which is wrapped on itself). This is a kind of periodic Neumann condition. But instead of speaking of Neumann, the documentation introduces a new notion (a choice between "source" and "target"). I don't know how to express "flux conservation" ... – andre314 Aug 30 '19 at 13:58
• @user21, drastic problem is that solution with periodic boundary conditions is ambiguous, mathenatically incorrect. – Rodion Stepanov Mar 25 at 12:13
• @RodionStepanov, The above problem behaves according to the quote from the documentation and I tried to explain the reasoning for it in my answer. This is the expected and correct behavior. Are you suggesting that the documentation is not clear enough? What do you mean by 'ill-defined PBC'? Most importantly: What do you expect? Unless you clarify that I can not fix it or improve the documentation. Why don't you ask a question illustrating the issue you have. – user21 Mar 26 at 8:48
• @RodionStepanov, if you can point me to literature that shows how your request can be implemented for FEM (not for FDM, there it's clear and works like what you suggest) than I can implement that. My current understanding is that for FEM it is technically not possible to implement what you suggest/request. I'd be happy to learn otherwise. And let me be honest, several people made similar requests but none could show so far how to implement that for FEM. If you find something let me know. Even better if you can demonstrate it with the low level FEM functions. To be perfectly clear: – user21 Mar 27 at 5:50
• I am not questioning that your suggestion is useful. I am simply stating that the code behaves as documented (= is correct) if that is not the desirable functionality then that is a different matter. Thanks for your feedback. And note that for the 'TensorProductGrid' method PBC work as you request. – user21 Mar 27 at 5:51

Nothing wrong here. This is expected. A periodic boundary condition takes whatever boundary conditions is present (explicitly or implicitly) at the source boundary and projects it to the target boundary. Since this seems to be a source of confusion I have tried to further clarify this in the documentation.

Here is what is documented now.

And here is what will appear as a new possible issues example in a future version (post V12.0)

Periodic boundary conditions relate the solution of a PDE from the source to the target boundary. Boundary conditions present, also implicit ones, at the source will affect the solution at the target.

To exemplify the behavior, consider a time-dependent equation discretized with the finite element method. An initial condition u, implicit Neumann zero boundary conditions on both sides and no PeriodicBoundaryCondition are specified:

ufun = NDSolveValue[{D[u[t, x], t] - D[u[t, x], {x, 2}] == 0,
u[0, x] == Sin[x]}, u, {t, 0, 1}, {x, -\[Pi], \[Pi]},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement"}}]


Visualize the solution at various times:

frames = Table[
Plot[ufun[t, x], {x, -\[Pi], \[Pi]}, PlotRange -> {-1, 1}], {t, 0,
1, 0.1}];
ListAnimate[frames, SaveDefinitions -> True]


Note that at both spatial boundaries the implicit Neumann 0 boundary conditions are satisfied.

When a PeriodicBoundaryCondition is used on a source boundary that has an implicit Neumann 0 boundary condition, then that condition will be mapped to the target boundary.

Following is the solution of the same equation and initial condition as previously and an additional periodic boundary condition that has its source on the left and its target on the right:

ufun = NDSolveValue[{D[u[t, x], t] - D[u[t, x], {x, 2}] == 0,
u[0, x] == Sin[x],
PeriodicBoundaryCondition[u[t, x], x == \[Pi],
Function[X, X - 2 \[Pi]]]}, u, {t, 0, 1}, {x, -\[Pi], \[Pi]},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement"}}]


Visualize the solution at various times:

Note how the solution value at the implicit Neumann 0 boundary condition on the left is mapped to the right.

This is the expected behavior for the finite element method. The tensor product grid method behaves differently, as that method does not have implicit boundary conditions:

ufunTPG =
NDSolveValue[{D[u[t, x], t] - D[u[t, x], {x, 2}] == 0,
u[0, x] == Sin[x], u[t, -\[Pi]] == u[t, \[Pi]]},
u, {t, 0, 1}, {x, -\[Pi], \[Pi]},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid"}}]


Visualize the tensor product grid solution at various times:

frames = Table[
Plot[ufunTPG[t, x], {x, -\[Pi], \[Pi]}, PlotRange -> {-1, 1}], {t,
0, 1, 0.1}];
ListAnimate[frames, SaveDefinitions -> True]


A similar behavior can be achieved with the finite element method by specifying a DirichletCondition on the left and a PeriodicBoundaryCondition:

ufunFEM =
NDSolveValue[{D[u[t, x], t] - D[u[t, x], {x, 2}] == 0,
u[0, x] == Sin[x],
PeriodicBoundaryCondition[u[t, x], x == \[Pi],
Function[X, X - 2 \[Pi]]],
DirichletCondition[u[t, x] == Sin[-\[Pi]], x == -\[Pi]]},
u, {t, 0, 1}, {x, -\[Pi], \[Pi]},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement"}}]


Visualize the difference between the finite element and tensor product grid solutions at various times:

frames = Table[
Plot[ufunFEM[t, x] - ufunTPG[t, x], {x, -\[Pi], \[Pi]},
PlotRange -> {-5 10^-4, 5 10^-4}], {t, 0, 1, 0.1}];
ListAnimate[frames, SaveDefinitions -> True]


Alternatively, a DirichletCondition could be specified at each side.

• Thank you very much for your detailed answer. It looks like the TranslationTransform I used in my examples for two variables make the problems? – Ulrich Neumann Aug 30 '19 at 8:01
• @UlrichNeumann, correct, when you change the direction of the mapping you change which boundary condition is mapped to the target. – user21 Aug 30 '19 at 9:28
• I have to clarify my comment. If I take your first example (ufun with PeriodiceBoundaryConditions) and change it to NDSolveValue[{D[u[t, x], t] - D[u[t, x], {x, 2}] == 0, u[0, x] == Sin[x],PeriodicBoundaryCondition[u[t, x], x == \[Pi], TranslationTransform[{0,-2 Pi }]]}, u, {t, 0,1}, {x, -\[Pi], \[Pi]} , Method -> {"MethodOfLines", "SpatialDiscretization" -> {"FiniteElement"}}] it doesn't evaluate anymore. – Ulrich Neumann Aug 30 '19 at 9:58
• That means in the Method of lines only "spatial" periodic boundary conditions are allowed? – Ulrich Neumann Aug 30 '19 at 10:01
• @UlrichNeumann, I can not quite imagine what the purpose of temporal periodic boundary condition is. Maybe there is a purpose, I just do not know about that. – user21 Aug 30 '19 at 11:10

There is a trick to get true periodic solution, i.e. u(t,x)=u(t,2pi+x) and u'(t,x)=u'(t,2pi+x). For that you have to double x-range and to choose x=0 as "source" for both boundaries.

ufunFEM =
NDSolveValue[{D[u[t, x], t] - D[u[t, x], {x, 2}] == 0,
u[0, x] == Sin[x],
PeriodicBoundaryCondition[u[t, x], x == 2 π,
Function[X, X - 2 π]],
PeriodicBoundaryCondition[u[t, x], x == -2 π,
Function[X, X + 2 π]]}, u, {t, 0, 1}, {x, -2 π, 2 π},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement"}}]

Plot[ufunFEM[1, x], {x, -2 π, 2 π}, PlotRange -> All,
PlotLegends -> Automatic]


This is the same result as obtained by the tensor product grid method

ufunTPG =
NDSolveValue[{D[u[t, x], t] - D[u[t, x], {x, 2}] == 0,
u[0, x] == Sin[x], u[t, -\[Pi]] == u[t, \[Pi]]},
u, {t, 0, 1}, {x, -\[Pi], \[Pi]},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {"TensorProductGrid"}}];

Plot[ufunTPG[1, x] - ufunFEM[1, x], {x, -\[Pi], \[Pi]},
PlotRange -> All, PlotLegends -> Automatic]


For 2D case it works too

Ω = Rectangle[{-2, 0}, {2, 1}];
pde = -Derivative[0, 2][u][x, y] - Derivative[2, 0][u][x, y] ==
If[(1.25 <= x + 2 <= 1.75 || 1.25 <= x <= 1.75) &&
0.25 <= y <= 0.5, 1., 0.];

ufun = NDSolveValue[{
pde,
PeriodicBoundaryCondition[u[x, y], x == -2 && 0 <= y <= 1,
TranslationTransform[{2, 0}]],
PeriodicBoundaryCondition[u[x, y], x == 2 && 0 <= y <= 1,
TranslationTransform[{-2, 0}]],
DirichletCondition[
u[x, y] == 0, (y == 0 || y == 1) && -2 < x < 2]},
u, {x, y} ∈ Ω];
ContourPlot[ufun[x, y], {x, y} ∈ Ω,
ColorFunction -> "TemperatureMap", AspectRatio -> Automatic]


This solution is different from two ones if you choose only on target boundary

Ω1 = Rectangle[{0, 0}, {2, 1}];
ufunR = NDSolveValue[{pde,
PeriodicBoundaryCondition[u[x, y], x == 2 && 0 <= y <= 1,
TranslationTransform[{-2, 0}]],
DirichletCondition[
u[x, y] == 0, (y == 0 || y == 1) && 0 < x < 2]},
u, {x, y} ∈ Ω1];
ufunL = NDSolveValue[{pde,
PeriodicBoundaryCondition[u[x, y], x == 0 && 0 <= y <= 1,
TranslationTransform[{2, 0}]],
DirichletCondition[
u[x, y] == 0, (y == 0 || y == 1) && 0 < x < 2]},
u, {x, y} ∈ Ω1];
Row[ContourPlot[#[x, y], {x, y} ∈ Ω1,
ColorFunction -> "TemperatureMap", AspectRatio -> Automatic,
ImageSize -> 300] & /@ {ufun, ufunR, ufunL}]


In fact there is no need to double numerical domain. Just add some ghost vicinity

Ω2 = Rectangle[{-0.01, 0}, {2 + 0.01, 1}];
ufun = NDSolveValue[{
pde,
PeriodicBoundaryCondition[u[x, y], x == -0.01 && 0 <= y <= 1,
TranslationTransform[{2, 0}]],
PeriodicBoundaryCondition[u[x, y], x == 2 + 0.01 && 0 <= y <= 1,
TranslationTransform[{-2, 0}]],
DirichletCondition[
u[x, y] == 0, (y == 0 || y == 1) && -0.01 < x < 2 + 0.01]},
u, {x, y} ∈ Ω2];
ContourPlot[ufun[x, y], {x, y} ∈ Ω2,
ColorFunction -> "TemperatureMap", AspectRatio -> Automatic]


Let's look at the limit of the ghost points to the original region size. Up until down to 10^-14. things work fine, it's only below that that the solution seems to change.

epsilon = 10^-14.;
pde = -Derivative[0, 2][u][x, y] - Derivative[2, 0][u][x, y] ==
If[(1.25 <= x + 2 <= 1.75 || 1.25 <= x <= 1.75) &&
0.25 <= y <= 0.5, 1., 0.];
\[CapitalOmega]2 = Rectangle[{-epsilon, 0}, {2 + epsilon, 1}];
ufun = NDSolveValue[{pde,
PeriodicBoundaryCondition[u[x, y], x == -epsilon && 0 <= y <= 1,
TranslationTransform[{2, 0}]],
PeriodicBoundaryCondition[u[x, y],
x == 2 + epsilon && 0 <= y <= 1, TranslationTransform[{-2, 0}]],
DirichletCondition[
u[x, y] == 0, (y == 0 || y == 1) && -epsilon < x < 2 + epsilon]},
u, {x, y} \[Element] \[CapitalOmega]2];
ContourPlot[ufun[x, y], {x, y} \[Element] \[CapitalOmega]2,
ColorFunction -> "TemperatureMap", AspectRatio -> Automatic]


Also note that if you use triangle elements you can use epsilon=0:

epsilon = 0;
pde = -Derivative[0, 2][u][x, y] - Derivative[2, 0][u][x, y] ==
If[(1.25 <= x + 2 <= 1.75 || 1.25 <= x <= 1.75) &&
0.25 <= y <= 0.5, 1., 0.];
\[CapitalOmega]2 = Rectangle[{-epsilon, 0}, {2 + epsilon, 1}];
ufun = NDSolveValue[{pde,
PeriodicBoundaryCondition[u[x, y], x == -epsilon && 0 <= y <= 1,
TranslationTransform[{2, 0}]],
PeriodicBoundaryCondition[u[x, y],
x == 2 + epsilon && 0 <= y <= 1, TranslationTransform[{-2, 0}]],
DirichletCondition[
u[x, y] == 0, (y == 0 || y == 1) && -epsilon < x < 2 + epsilon]},
u, {x, y} \[Element] \[CapitalOmega]2,
Method -> {"FiniteElement",
"MeshOptions" -> {"MeshElementType" -> "TriangleElement"}}];
ContourPlot[ufun[x, y], {x, y} \[Element] \[CapitalOmega]2,
ColorFunction -> "TemperatureMap", AspectRatio -> Automatic]


• This seems interesting bug I can not reproduce the results you show. 1) you are missing definitions 2) When I add T=1 I still get different results. 3) It's not clear to me what you are trying to say with the last example. Could you clarify this a bit. Also in a comment to this answer you mentioned that a doubling of the domain is not necessary, could you perhaps sow an example? – user21 Apr 30 at 7:37
• Here I have extended the domain from (2 Pi) to (2 Pi +Pi/4) to obtain Neuman + Dirichlet periodicity. It's the same approach, I think. – andre314 Apr 30 at 8:48
• By the way, I'm working on a more conventional approach to this problem. It consists in modifying the stiffness matrices (and maybe too the other matices of the system). I have already done it in 1D, 2D, but only with InperpolationOrder = MeshOrder =1 (here we have MeshOrder = 2) – andre314 Apr 30 at 8:55
• More preciselly, I work on 2 other approachs. – andre314 Apr 30 at 10:15
• @RodionStepanov Why don't you use this approach to this problem mathematica.stackexchange.com/questions/220631/… ? – Alex Trounev May 20 at 18:54

Beginning of explanations are coming later (2 days ?).

The code below is complete, so one can already evaluate it and enjoy.

Short and quick explanations are already possible in in this chatroom, but the subject is really hudge.

If you see a problem or some possible simplification anywhere, don't hesitate to comment.

It could save me some iterations in the construction of this answer.

Needs["NDSolveFEM"]

domain = Rectangle[{0, 0}, {2, 1}];
pde = -Laplacian[u[x, y], {x, y}] ==
If[1.25 <= x <= 1.75 && 0.25 <= y <= 0.5, 1., 0.];
bcFullDirichlet = DirichletCondition[u[x, y] == 0, True];

pointMarkerFunction =
Compile[{{coords, _Real, 2}, {pMarker, _Integer, 1}},
Block[{x = #1[[1]], y = #1[[2]], autoMarker = #2},
Which[
y == 1 , 3,
True, autoMarker]
] &, {coords, pMarker}]];

mesh50 = ToElementMesh[domain, "MeshElementType" -> "QuadElement"
, "MeshOrder" -> 2, "PointMarkerFunction" -> pointMarkerFunction ];

Show[mesh50["Wireframe"["MeshElement" -> "PointElements"
, "MeshElementMarkerStyle" ->
Directive[Black, FontWeight -> Bold, FontSize -> 6]
, "MeshElementStyle" -> (Directive[AbsolutePointSize[4],
Opacity[.8], #] & /@
{Black, Red, Green, Blue})]]
, Frame -> True]

newMesh00 = ToElementMesh[
"Coordinates" -> mesh50 ["Coordinates"]
, "MeshElements" -> mesh50["MeshElements"]
, "BoundaryElements" -> (mesh50["BoundaryElements"] //
, "PointElements" -> (mesh50["PointElements"] //

vd = NDSolveVariableData[{"DependentVariables",
"Space"} -> {{u}, {x, y}}];
nr = ToNumericalRegion[newMesh00];
sd = NDSolveSolutionData[{"Space"} -> {nr}];
bcdata = InitializeBoundaryConditions[vd, sd, {{bcFullDirichlet}}];
mdata = InitializePDEMethodData[vd, sd];

cdata = NDSolveProcessEquations[{pde, bcFullDirichlet}, u,
Element[{x, y}, domain]
, Method -> {"PDEDiscretization" -> {"FiniteElement",
"MeshOptions" ->
{"MeshElementType" -> QuadElement, "MeshOrder" -> 2}}}] //
RightComposition[
First
, #["FiniteElementData"] &
, #[PDECoefficientData] &
];

discretePDE = DiscretizePDE[cdata, mdata, sd
, "SaveFiniteElements" -> True, "AssembleSystemMatrices" -> True];
{load, stiffness, damping, mass} = discretePDE["SystemMatrices"];

dbc1 = DiscretizeBoundaryConditions[bcdata, mdata, sd
, "Stationary", "PartialBoundaryAssembly" -> {1 }];
dbc3 = DiscretizeBoundaryConditions[bcdata, mdata, sd
, "Stationary", "PartialBoundaryAssembly" -> {3 }];

dbc2 = DiscretizeBoundaryConditions[bcdata, mdata, sd
, "Stationary", "PartialBoundaryAssembly" -> {2}] ;
dbc4 = DiscretizeBoundaryConditions[bcdata, mdata, sd
, "Stationary", "PartialBoundaryAssembly" -> {4}];

stiffness[[dbc2["DirichletRows"]]] =
stiffness[[dbc2["DirichletRows"]]] +
stiffness[[dbc4["DirichletRows"]]];
stiffness[[All, dbc2["DirichletRows"]]] =
stiffness[[All, dbc2["DirichletRows"]]] +
stiffness[[All, dbc4["DirichletRows"]]] ;

stiffnessReduced = stiffness //
Delete[#, List /@ dbc4["DirichletRows"]] & //
(Delete[#, List /@ dbc4["DirichletRows"]] & /@ # &);

Fold[Insert[#1, {0.}, {#2}] &, solution20, dbc4["DirichletRows"]];

NDSolveSetSolutionDataComponent[sd, "DependentVariables",
{sol} = ProcessPDESolutions[mdata, sd];

(* beyond this point : visualization of the solution sol *)
myOptions01 = {ColorFunction -> "TemperatureMap",
AspectRatio -> Automatic
, Frame -> {True, True}, PlotRangePadding -> None
, ImagePadding -> {{0, 0}, {30, 10}}};
myDuplicateImage[image_] :=
Rasterize[image] // ImageAssemble[{{#, #}}] &
myViewOptions = {ViewAngle -> 0.42, ViewCenter -> {0.5, 0.5, 0.5}
, ViewMatrix -> Automatic, ViewPoint -> {0.34, -3.36, -0.12}
, ViewProjection -> Automatic, ViewRange -> All
, ViewVector -> Automatic
, ViewVertical -> {0.00378, -0.037, 1.}};
myStreamContourPlot00[ufun_] :=
Column[{
Plot3D[ufun[x, y], {x, y} \[Element] domain,
ColorFunction -> "TemperatureMap"] //
{Show[#, ViewAngle -> 0.42],
Show[#, Evaluate @ myViewOptions]} & // Row
, ContourPlot[Evaluate @ ufun[x, y]
, Element[{x, y}, domain], Evaluate @ myOptions01] //
myDuplicateImage
, StreamDensityPlot[
Evaluate @ {-Grad[ufun[x, y], {x, y}], ufun[x, y]}
, Element[{x, y}, domain], Evaluate @ myOptions01] //
myDuplicateImage
, DensityPlot[Evaluate[Norm @ Grad[ufun[x, y], {x, y}]]
, Element[{x, y}, domain]
, PlotPoints -> 100, Frame -> False, Evaluate @ myOptions01] //

myDuplicateImage} //
Thread[Labeled[#, {"Overviews", "graphic 1 : Dirichlet periodic"
, "graphic 2 : Neuman periodic (flux direction verification)"
,
"graphic 3 : Neuman periodic (flux intensity verification)"},
Top]] &
, Dividers -> None, Spacings -> {1, 4}] //
Style[#, ImageSizeMultipliers -> {1, 1}] &;

Labeled[myStreamContourPlot00[sol]
, Style["\n\n(Dirichlet & Neuman) periodicity visualization\n\n",
FontSize -> 18, FontWeight -> Bold], Top]


• Thanks Andre for your answer. I need think a bit more about that and I am keen to hear your ideas on this. Just for the partial boundary assembly, I should award a bounty ;-) – user21 Jun 8 at 8:13
• I'm still working on this subject – andre314 Jun 16 at 20:41
• Note : The code above is still right, but I have probably found a simplification (that don't delete/restore some row and columns of the matrices). This simplify the workflow, (above all if the problem is furthermore a transcient problem). I will share all this as soon as possible, but I don't want to iterate my answer too much. – andre314 Jun 20 at 23:59
• Code above still right, but the future simplifications make its understanding probably useless (except for those who are interested in the history of my approach). – andre314 Jun 24 at 19:24

Although I anxiously await Andres' complete write up, I thought that I would post some observations that may help in the investigation of the PeriodicBoundaryCondition. In this case, my initial findings are that a combination @Rodion Stepanov's symmetrized PBC and a triangle mesh lead to more robust results without needing a "Ghost Vicinity".

# Default Element Mesh for Rectangle Domains are Quads.

If we copy Rodion's ghost vicinity example and view the mesh, we see that it is a quad mesh.

pde = -Derivative[0, 2][u][x, y] - Derivative[2, 0][u][x, y] ==
If[(1.25 <= x + 2 <= 1.75 || 1.25 <= x <= 1.75) &&
0.25 <= y <= 0.5, 1., 0.];
Ω2 = Rectangle[{-0.01, 0}, {2 + 0.01, 1}];
ufun = NDSolveValue[{pde,
PeriodicBoundaryCondition[u[x, y], x == -0.01 && 0 <= y <= 1,
TranslationTransform[{2, 0}]],
PeriodicBoundaryCondition[u[x, y], x == 2 + 0.01 && 0 <= y <= 1,
TranslationTransform[{-2, 0}]],
DirichletCondition[
u[x, y] == 0, (y == 0 || y == 1) && -0.01 < x < 2 + 0.01]},
u, {x, y} ∈ Ω2];
ContourPlot[ufun[x, y], {x, y} ∈ Ω2,
ColorFunction -> "TemperatureMap", AspectRatio -> Automatic]
ufun["ElementMesh"]["Wireframe"]


# Using Symmetrized PBC's on a Triangle Mesh Requires No Ghost Vicinity

Before I show the workflow, I will set up a colormap so we can compare to another solver later.

(* Banded ColorMap *)
img = Uncompress[
"1:eJzt2+tP02cUB/\
CjYjQMnYuTYHQzLJItGI2OuWA0EpjG6eI07Vi8IFrgZ630Ai3VNjqeGQgCYyAKdlSBAuVS\
ZSgV5A5ekMWBEFEjYkBxBiUoTofxFvjamu2N/8GS8+KcnHOekzxvPm+\
Pb4ROtnMyERncaa1GoZR2TnS3Xq70vVEj6VWRwXq9whwxyTXwccUlV7hrPHyI3l50dKC5G\
ZWVKCpCdjYOHoTJhN27ERaGDRsQHIyAAPj5wccHnp4vp9Dwx9T3GXUtpvMrqeo7KtlMvyk\
peS/tSyTNYdpuI9nvtKqBvr5MX9ykOffJ8znRGw8a+YjuzqPuhdS6nGq+JcePdCyKfomj+\
AMUk0ERuRR6gtbU0rI2WnCdPh2gac8mTBifPv3p3Ll/+fvfCAz8Y/Xqerm8XKHIi41NF+\
LntDSD1SqVlm6qrl538eKKq1cX9ff7PnkyY2xsIkY/\
wOBs9HyOP5eiKQSnNiJPgUwtEvZjTwp2WbDVjvVOBJ3Dkk749mPmI0x+/\
WIqhrxxez6ufIlzQXCuR0E4sqKRZIY5CdFZCC/AxlMIacJX7Zh/G95DmPoCk8bg9RKz/\
sEnI/AbwqL7WNaH4B6suwZZJ7ZeRmQr1C0w1iO+\
CskVOORAjh0223hB3mjB8eFC673CnFtFRzuLslvtRxrtmc7iDEdJen5JmqU09dfS5MSyJH\
NZYowjQek4sO2ECK0Qm8+I7bVCahTRF4S+\
TZjaxU9dIuG6SOkRGX0ia0BYB4VtWJT8LcqfC+crUTsuml7HN4/ua35sbnqwt/\
GOsfGWoaE7tr5DV3dJU9cSXVunqnEqa8qls/\
aI6twdVZbwqkNhZ1K3OFPDKjMVFRblyXxNWbGhuNxU6Iy31SXktqRY29ItHVnZ3TmHe20Z\
A8VpD06mjJxOYk7MiTkxJ+\
bEnJgTc2JOzIk5MSfmxJyYE3NiTsyJOTEn5sScmBNzYk7MiTkxJ+\
bEnJgTc2JOzIk5MSfmxJyYE3NiTsyJOTEn5sScmBNzYk7MiTkxp/8dJ/\
kMIgrVGlRKrRS1VhsnKSV9oNzDNQwxx/17rOfuZEa1ZPB0Fd/\
o1Dq9PEYRKcndd3qyNSHvLX3436WfTDLo1MY4lU6rMrlm7625LwDd/+nVkmKPSqt89/\
KD3ii9BWHVFNA="];
dims = ImageDimensions[img];
colors = RGBColor[#] & /@
ImageData[img][[IntegerPart@(dims[[2]]/2), 1 ;; -1]];


Now, we will force a triangle mesh using ToElementMesh on the domain and we will not use a ghost vicinity as shown in the following workflow.

Needs["NDSolveFEM"]
{length, height, xc, yc} = {2, 1, 0, 0};
{sx, sy, fx, fy} = {0, 0, length, height};
{ssx, ssy, fsx, fsy} = {1.25, 0.25, 1.75, 0.5};
centersource = Mean[{{ssx, ssy}, {fsx, fsy}}];
srcReg = Rectangle[{ssx, ssy}, {fsx, fsy}];
source = If[ssx <= x <= fsx && ssy <= y <= fsy, 1., 0.];
pde = -\!$$\*SubsuperscriptBox[\(∇$$, $${x, y}$$, $$2$$]$$u[x, y]$$\) -
source == 0;
Ω = Rectangle[{sx, sy}, {fx, fy}];
mesh = ToElementMesh[Ω,
"MeshElementType" -> TriangleElement];
mesh["Wireframe"]
ufun = NDSolveValue[{pde,
PeriodicBoundaryCondition[u[x, y], x == sx && 0 <= y <= 1,
TranslationTransform[{length, 0}]],
PeriodicBoundaryCondition[u[x, y], x == fx && 0 <= y <= 1,
TranslationTransform[{-length, 0}]],
DirichletCondition[
u[x, y] == 0, (y == 0 || y == 1) && sx < x < fx]},
u, {x, y} ∈ mesh];
Plot3D[ufun[x, y], {x, y} ∈ mesh,
ColorFunction -> "TemperatureMap", AspectRatio -> Automatic]
ContourPlot[ufun[x, y], {x, y} ∈ mesh,
ColorFunction -> "TemperatureMap", AspectRatio -> Automatic]
Plot3D[Evaluate@Norm[Grad[ufun[x, y], {x, y}]], {x, y} ∈
mesh, PlotPoints -> 250, ColorFunction -> (Blend[colors, #3] &),
BoxRatios -> {2, 1, 1/2}, PerformanceGoal -> "Quality", Mesh -> None,
Background -> Black]
DensityPlot[
Evaluate@Norm[Grad[ufun[x, y], {x, y}]], {x, y} ∈ mesh,
ColorFunction -> "TemperatureMap", PlotPoints -> All,
AspectRatio -> Automatic]


As you can see, it solves without requiring the any extra padding of the domain. We can see that the flux magnitude is quite jagged. We can fix the solution by provide the appropriate refinement zones at the wall and around the source.

# Mesh Refined Solution

The following workflow will refine the mesh and re-solve the PDE.

(* Shrink source 10% *)
smallSrc =
TransformedRegion[srcReg,
ScalingTransform[0.9 {1, 1}, centersource]];
(* Expand source 10% *)
bigSrc = TransformedRegion[srcReg,
ScalingTransform[1.1 {1, 1}, centersource]];
(* Create a Difference Around the Source Edge *)
diff = RegionDifference[bigSrc, smallSrc];
(* Create mesh refinement function *)
mrf = With[{rmf = RegionMember[diff],
rmfinner = RegionMember[smallSrc]},
Function[{vertices, area},
Block[{x, y}, {x, y} = Mean[vertices];
Which[rmf[{x, y}], area > 0.00005,
rmfinner[{x, y}], area > 0.000125,
True, area > 0.00125]]]];
(* Create and display refined mesh *)
mesh = ToElementMesh[Ω,
"MaxBoundaryCellMeasure" -> 0.01,
"MeshElementType" -> TriangleElement,
MeshRefinementFunction -> mrf];
mesh["Wireframe"]
(* Solve and display solution *)
ufun = NDSolveValue[{pde,
PeriodicBoundaryCondition[u[x, y], x == sx && 0 <= y <= 1,
TranslationTransform[{length, 0}]],
PeriodicBoundaryCondition[u[x, y], x == fx && 0 <= y <= 1,
TranslationTransform[{-length, 0}]],
DirichletCondition[
u[x, y] == 0, (y == 0 || y == 1) && sx < x < fx]},
u, {x, y} ∈ mesh];
Plot3D[ufun[x, y], {x, y} ∈ mesh,
ColorFunction -> "TemperatureMap", AspectRatio -> Automatic]
ContourPlot[ufun[x, y], {x, y} ∈ mesh,
ColorFunction -> "TemperatureMap", AspectRatio -> Automatic]
Plot3D[Evaluate@Norm[Grad[ufun[x, y], {x, y}]], {x, y} ∈
mesh, PlotPoints -> 250, ColorFunction -> (Blend[colors, #3] &),
BoxRatios -> {2, 1, 1/2}, PerformanceGoal -> "Quality", Mesh -> None,
Background -> Black]
DensityPlot[
Evaluate@Norm[Grad[ufun[x, y], {x, y}]], {x, y} ∈ mesh,
ColorFunction -> "TemperatureMap", PlotPoints -> All,
AspectRatio -> Automatic]
`

The flux magnitude results look much less jagged.

# Comparison to Another Solver

I always find it useful to compare the Mathematica results to another solver for a sanity check. In this case, I compare the Mathematica results to Altair's AcuSolve and we see that the results are quite similar. I don't know how general the solution is, but I would recommend using Rodion's symmetrized PBC approach and use Triangle or Tet Elements versus Quads or Hexa as there seems to be negative interaction with setting a PBC.

## COMSOL, AcuSolve, and Mathematica Comparison with Same ColorMap.

For completeness, I am positing a comparison of the simulation results of COMSOL, Altair's AcuSolve, and Mathematica on the same ColorMap to show that these FEM codes all are in agreement.

• See my update to Rodion's answer. It's possible to use a very small distance from the original region for the ghost points and for triangle elements no ghost points are needed at all, like you have shown. I am not sure what causes this, though. But it's very useful to have a comparison to other solver results. – user21 Jun 8 at 8:17