# Solving PDE with prescribed rectangular mesh

Question: I want to solve a PDE using NDEigensystem. The region is a rectangle Rectangle[{0,-70}, {20, 70}]. How can I make a mesh that satisfies the following two conditions?

(1) The mesh is finer near $$y \approx 0$$.

(2) The line $$y=0$$ should be the boundary of the mesh.

Motivation: The PDE that I want to solve can be written in the following command: First, I set the mesh as below:

mesh = ToElementMesh[Rectangle[{0, -70}, {20, 70}],
MeshRefinementFunction ->
Function[{vertices, area},
Block[ {x, y}, {x, y} = Mean[vertices];
If[-10 < y < 10, area > 0.1, area > 10]]] ]


and then write down the PDE:

NDEigensystem[{{3*psi3[x, y] - 10*psi2[x, y]*Sign[y] - Derivative[0, 1][psi2][x, y] -
I*Derivative[1, 0][psi2][x, y], 3*psi4[x, y] - 10*psi1[x, y]*Sign[y] + Derivative[0, 1][psi1][x, y] -
I*Derivative[1, 0][psi1][x, y], 3*psi1[x, y] - 10*psi4[x, y]*Sign[y] - Derivative[0, 1][psi4][x, y] -
I*Derivative[1, 0][psi4][x, y], 3*psi2[x, y] - 10*psi3[x, y]*Sign[y] + Derivative[0, 1][psi3][x, y] -
I*Derivative[1, 0][psi3][x, y]}, PeriodicBoundaryCondition[psi1[x, y], x == 0,
TransformationFunction[{{1, 0, 20}, {0, 1, 0}, {0, 0, 1}}]], PeriodicBoundaryCondition[psi2[x, y], x == 0,
TransformationFunction[{{1, 0, 20}, {0, 1, 0}, {0, 0, 1}}]], PeriodicBoundaryCondition[psi3[x, y], x == 0,
TransformationFunction[{{1, 0, 20}, {0, 1, 0}, {0, 0, 1}}]], PeriodicBoundaryCondition[psi4[x, y], x == 0,
TransformationFunction[{{1, 0, 20}, {0, 1, 0}, {0, 0, 1}}]], DirichletCondition[psi1[x, y] == 0,
0 < x < 20 && (y == 70 || y == -70)], DirichletCondition[psi2[x, y] == 0,
0 < x < 20 && (y == 70 || y == -70)], DirichletCondition[psi3[x, y] == 0,
0 < x < 20 && (y == 70 || y == -70)], DirichletCondition[psi4[x, y] == 0,
0 < x < 20 && (y == 70 || y == -70)]}, {psi1, psi2, psi3, psi4}, Element[{x, y}, mesh], 10]


The single important fact about the PDE is that the differential operator contains $$\mathrm{sgn}(y)$$, hence it is not continuous at $$y=0$$. This is why I want the mesh boundary to include $$y=0$$. The mesh I created is visualized as follows:

As wanted, the mesh is finer. However, $$y=0$$ is not the boundary of the mesh.

The result of one solution obtained by solving the PDE is plotted as below:

In the above plot, the horizontal direction is the $$x$$-direction. Note that the solution is localized near $$y\approx 0$$. This is why I want to make the mesh finer near $$y\approx 0$$. I don't believe the result is the good representation of the solution, since the graph looks jagged, which I don't expect since the differential operator has a translation symmetry in $$x$$-direction.

• Line y=0 is not shown in the mesh. What boundary condition supposed to be on the line y=0? Commented Mar 20, 2021 at 21:46

You can achieve your objectives with an anisotropic quad mesh that you construct manually. There are no diffusive terms in your operators, so NDEigensystem will complain about the system being convection dominated. Therefore, you may want to take a closer look at your system of equations and see if there is an opportunity to add diffusive terms.

# Helper functions

The following helper functions will help construct an anisotropic Quad mesh (note that not all functions will be used).

(*Import required FEM package*)
Needs["NDSolveFEM"];
(*Define Some Helper Functions For Structured Meshes*)
pointsToMesh[data_] :=
MeshRegion[Transpose[{data}],
Line@Table[{i, i + 1}, {i, Length[data] - 1}]];
unitMeshGrowth[n_, r_] :=
Table[(r^(j/(-1 + n)) - 1.)/(r - 1.), {j, 0, n - 1}]
meshGrowth[x0_, xf_, n_, r_] := (xf - x0) unitMeshGrowth[n, r] + x0
firstElmHeight[x0_, xf_, n_, r_] :=
Abs@First@Differences@meshGrowth[x0, xf, n, r]
lastElmHeight[x0_, xf_, n_, r_] :=
Abs@Last@Differences@meshGrowth[x0, xf, n, r]
findGrowthRate[x0_, xf_, n_, fElm_] :=(*Quiet@*)
Abs@FindRoot[
firstElmHeight[x0, xf, n, r] - fElm, {r, 0.00000001,
100000000/fElm}, Method -> "Brent"][[1, 2]]
meshGrowthByElm[x0_, xf_, n_, fElm_] :=
N@Sort@Chop@meshGrowth[x0, xf, n, findGrowthRate[x0, xf, n, fElm]]
meshGrowthByElm0[len_, n_, fElm_] := meshGrowthByElm[0, len, n, fElm]
flipSegment[l_] := (#1 - #2) & @@ {First[#], #} &@Reverse[l];
leftSegmentGrowth[len_, n_, fElm_] := meshGrowthByElm0[len, n, fElm]
rightSegmentGrowth[len_, n_, fElm_] :=
Module[{seg}, seg = leftSegmentGrowth[len, n, fElm];
flipSegment[seg]]
reflectRight[pts_] :=
With[{rt = ReflectionTransform[{1}, {Last@pts}]},
Union[pts, Flatten[rt /@ Partition[pts, 1]]]]
reflectLeft[pts_] :=
With[{rt = ReflectionTransform[{-1}, {First@pts}]},
Union[pts, Flatten[rt /@ Partition[pts, 1]]]]
extendMesh[mesh_, newmesh_] := Union[mesh, Max@mesh + newmesh]


The following workflow will create a symmetric graded mesh about $$y=0$$ and a uniform mesh in the x-direction.

ht = 140;
len = 20;
ny = 100;
nx = 100;
Print["Mesh along the X-direction"]
rh = pointsToMesh@Subdivide[0, len, nx]
Print["Mesh along the Y-direction"]
rv = pointsToMesh@reflectLeft@leftSegmentGrowth[ht/2, ny/2, ht/2/ny/10]
(*Create 2D mesh with RegionProduct*)
rp = RegionProduct[rh, rv];
(* Build mesh based on region product *)
crd = MeshCoordinates[rp];
inc = Delete[0] /@ MeshCells[rp, 2];
mesh = ToElementMesh["Coordinates" -> crd,
Print["Mesh zoomed in around y=0"]
Show[mesh["Wireframe"], PlotRange -> {{0, 20}, {-5, 5}}]


This mesh is refined and has a line about $$y=0$$.

# Solve and plot solution

(*Solve and plot*)
{vals, funs} =
NDEigensystem[{{3*psi3[x, y] - 10*psi2[x, y]*Sign[y] -
Derivative[0, 1][psi2][x, y] - I*Derivative[1, 0][psi2][x, y],
3*psi4[x, y] - 10*psi1[x, y]*Sign[y] +
Derivative[0, 1][psi1][x, y] - I*Derivative[1, 0][psi1][x, y],
3*psi1[x, y] - 10*psi4[x, y]*Sign[y] -
Derivative[0, 1][psi4][x, y] - I*Derivative[1, 0][psi4][x, y],
3*psi2[x, y] - 10*psi3[x, y]*Sign[y] +
Derivative[0, 1][psi3][x, y] - I*Derivative[1, 0][psi3][x, y]},
PeriodicBoundaryCondition[psi1[x, y], x == 0,
TransformationFunction[{{1, 0, 20}, {0, 1, 0}, {0, 0, 1}}]],
PeriodicBoundaryCondition[psi2[x, y], x == 0,
TransformationFunction[{{1, 0, 20}, {0, 1, 0}, {0, 0, 1}}]],
PeriodicBoundaryCondition[psi3[x, y], x == 0,
TransformationFunction[{{1, 0, 20}, {0, 1, 0}, {0, 0, 1}}]],
PeriodicBoundaryCondition[psi4[x, y], x == 0,
TransformationFunction[{{1, 0, 20}, {0, 1, 0}, {0, 0, 1}}]],
DirichletCondition[psi1[x, y] == 0,
0 < x < 20 && (y == 70 || y == -70)],
DirichletCondition[psi2[x, y] == 0,
0 < x < 20 && (y == 70 || y == -70)],
DirichletCondition[psi3[x, y] == 0,
0 < x < 20 && (y == 70 || y == -70)],
DirichletCondition[psi4[x, y] == 0,
0 < x < 20 && (y == 70 || y == -70)]}, {psi1, psi2, psi3, psi4},
Element[{x, y}, mesh], 10];
Table[Plot3D[
funs[[i, 1]][x, y] Conjugate[funs[[i, 1]][x, y]], {x, y} \[Element]
mesh, PlotRange -> All, PlotLabel -> Chop@vals[[i]], Axes -> True,
AxesLabel -> Automatic], {i, Length[vals]}]


The solution no longer is jagged in the X-direction. However, there are many oscillations in the Y-direction. You may need to consider diffusive terms to dampen these oscillations.

• Thank you so much for your detailed answer. Your answer does answer my question, but it is strange that the solution is not localized near $y\approx 0$. There is a good physical reason that the solution should be localized near $y\approx 0$. I may suspect that this is due to the "convection dominated PDE". How is the result reliable? Commented Mar 21, 2021 at 8:32
• @eigenvalue The result may not be reliable. There is a good chance, but not necessarily, that a convection-dominated problem is numerically unstable. I found a question regarding the semi-classical Pauli operator here 148252 that seems to have some similarity to yours, but it has a diffusive term. If there is an opportunity to add this diffusive term to your operator, it could stabilize the results. Alternatively, you could investigate adding artificial diffusion. Commented Mar 21, 2021 at 17:59