1
$\begingroup$

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:

enter image description here

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}}]]

enter image description here

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]]

enter image description here

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"}]

enter image description here

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

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – user21 Jun 9 '16 at 23:39
1
$\begingroup$

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]}]
 ]

enter image description here

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]}]
 ]

enter image description here

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"]

enter image description here

For the problem at hand, if you look at

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

enter image description here

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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