# Why is the solution of NDSolveValue asymmetric when the regions it's solving over are symmetric?

I am creating a mesh over a 3D region and then using NDSolveValue to try and solve for a function over that mesh. The structure is fairly simple, just two polygon squares inside of a big cube. I define the squares and big cube with a set of polygons, create a MeshRegion from them, then use ToBoundaryMesh on that. I then pass that boundary mesh to ToElementMesh and use a mesh refinement function around the two squares. This all seems to work pretty well.

Here is a link to a .txt file of the mesh, because it is 3MB (~24k mesh elements) and I probably can't paste it here.

You can turn it into a usable mesh by doing

mesh=ToExpression@Import@"Path\\to\\file.txt";


myMesh["Wireframe"] lets us see the boundary mesh:

And this slice lets you see how the mesh is a lot denser near the squares:

mesh["Wireframe"["MeshElement" -> "MeshElements",
PlotRange -> {{0, 100}, {0, 100}, {49.5, 50.5}}]]


Now I just want to solve for the function u (it's essentially Laplace's equation, where one of the squares is at V = 10 and the other is at V = 0). features is just two polygons that the squares were created from:

features={Polygon[{{85/2, 46.7, 46.7}, {85/2, 46.7, 53.3}, {85/2, 53.3,
53.3}, {85/2, 53.3, 46.7}}],
Polygon[{{115/2, 46.7, 46.7}, {115/2, 46.7, 53.3}, {115/2, 53.3,
53.3}, {115/2, 53.3, 46.7}}]};
sol = NDSolveValue[{D[u[x, y, z], x, x] + D[u[x, y, z], y, y] +
D[u[x, y, z], z, z] == 0,
DirichletCondition[u[x, y, z] == 10,
RegionMember[First@features, {x, y, z}]],
DirichletCondition[u[x, y, z] == 0.0,
RegionMember[Last@features, {x, y, z}]]},
u, {x, y, z} \[Element] mesh];


So I get the solution, sol, and then plot it to look at it. Similar to the mesh slice above, I plot the solution in slices, so I can see its behavior:

Show[SliceDensityPlot3D[
sol[x, y, z], {{z == 50}}, {x, 0, 100}, {y, 0, 100}, {z, 0, 100},
ColorFunction -> "Rainbow", PlotRange -> {All, All, All},
PlotLegends -> Automatic, ImageSize -> Large,
AxesLabel -> {"x", "y", "z"}, ViewPoint -> {0, 0, 25},
PerformanceGoal -> "Quality", ClippingStyle -> Black],
Graphics3D[{Opacity[0], First@features, Opacity[0], Last@features},
ImageSize -> Large]]


You can see that there's this "blockiness" that seems rougher than the mesh. However, what worries me more is that there's an asymmetry to this, where there really shouldn't be any (because the whole thing is symmetric).

If I set z = 50 and use a regular DensityPlot instead, it's actually worse:

DensityPlot[sol[x, y, 50], {x, 0, 100}, {y, 0, 100},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
ClippingStyle -> Black, ImageSize -> Large,
PlotRange -> {All, All, All}, AxesLabel -> {"x", "y"}]


The blockiness I might attribute to needing a finer mesh, but why am I seeing this fairly prominent asymmetry?

• It would really, really, really help if your posts could be a bit more copy and paste friendly. In essence I would like to be able to copy your complete code from this site and paste it into a notebook and run the stuff and see exactly what you show here. For example: "Path\\to\\file.txt" - If you absolutely have to use an external file then the it should be Import[https://gist.githubusercontent.com/anonymous/afcf7d999dafd644500ccfa8e3ce1f73/raw/bd30ceb2411700059018dac7ad78fe1eb78ab3c2/gistfile1.txt]. For these problems it makes it easier if the mesh generation is part of the post. Commented Jun 9, 2016 at 23:39

If you look at what RegionMember[First@features, {x, y, z}] actually finds you see something like this (after a small re-write)

ply = {Polygon[{{85/2, 46.7, 46.7}, {85/2, 46.7, 53.3}, {85/2, 53.3,
53.3}, {85/2, 53.3, 46.7}}],
Polygon[{{115/2, 46.7, 46.7}, {115/2, 46.7, 53.3}, {115/2, 53.3,
53.3}, {115/2, 53.3, 46.7}}]};
features = RegionMember[#, {x, y, z}] & /@ ply;
p1 = Pick[mesh["Coordinates"],
Function[{x, y, z}, Evaluate[features[[1]]]] @@@
mesh["Coordinates"], True]
{{42.5, 46.7, 46.7}, {42.5, 53.3, 53.3}, {42.5, 53.3,
48.35}, {42.5, 46.7, 51.65}}

Show[
mesh["Wireframe"["MeshElement" -> "MeshElements",
PlotRange -> {{42, 43}, {46, 54}, {46, 54}}]]
, Graphics3D[
Polygon[{{85/2, 46.7, 46.7}, {85/2, 46.7, 53.3}, {85/2, 53.3,
53.3}, {85/2, 53.3, 46.7}}]]
, Graphics3D[{Red, PointSize[0.1], Point[p1]}]
]


It'd advise against using RegeionMember in boundary conditions. In the way you have set it it's also very expensive to compute because for each time a DirichletCondition is used it needs to run the whole RegionMember setup.

If you use something like:

p1 = Pick[mesh["Coordinates"], Function[{x, y, z}, Evaluate[x == 42.5 && 46 <= y <= 54 && 46 <= z <= 54]] @@@ mesh["Coordinates"], True]

Show[
mesh["Wireframe"["MeshElement" -> "MeshElements",
PlotRange -> {{42, 43}, {46, 54}, {46, 54}}]]
, Graphics3D[
Polygon[{{85/2, 46.7, 46.7}, {85/2, 46.7, 53.3}, {85/2, 53.3,
53.3}, {85/2, 53.3, 46.7}}]]
, Graphics3D[{Red, PointSize[0.05], Point[p1]}]
]


sol = NDSolveValue[{Laplacian[u[x, y, z], {x, y, z}] == 0,
DirichletCondition[u[x, y, z] == 10,
x == 42.5 && 46 <= y <= 54 && 46 <= z <= 54],
DirichletCondition[u[x, y, z] == 0.0,
x == 57.5 && 46 <= y <= 54 && 46 <= z <= 54]},
u, {x, y, z} \[Element] mesh];

ContourPlot[sol[x, y, 50], {x, 0, 100}, {y, 0, 100}, PlotRange -> All,
ColorFunction -> "Rainbow"]


For the problem at hand, if you look at

Plot[sol[50, y, 50], {y, 0, 100}]


it looks to me as if the ContourPlot is overly sensitive to the data. If you want to avoid that you could add a face in the boundary mesh at x==50 that should help it.