8
$\begingroup$

I am slightly puzzled over the example "Solve PDEs over 3D Regions" provided by Wolfram. The example shows how to solve Laplace's equation over a space shuttle mesh. However, if I run the code in 11.0.1.0 the solution produced appears to be completely homogeneous:

mr = BoundaryDiscretizeGraphics[
ExampleData[{"Geometry3D", "SpaceShuttle"}]];
uif = NDSolveValue[{Inactive[Laplacian][u[x, y, z], {x, y, z}] == 0, 
DirichletCondition[u[x, y, z] == 1, z <= -1.3]}, u, {x, y, z} ∈ mr];

Then

vals = uif["ValuesOnGrid"];
{Max[vals], Min[vals]}

produces the output {1.,1.}, suggesting that the temperature field is homogeneous. Somehow, ElementMeshSurfacePlot3D[uif, Boxed -> False, ViewPoint -> {0, -4, 2}] still produces the beautiful coloured plot of the documentation. This confuses me.

Here is a partially successful modified example with slightly altered Dirichlet boundary conditions:

uif = NDSolveValue[{Inactive[Laplacian][u[x, y, z], {x, y, z}] == 0, 
DirichletCondition[u[x, y, z] == 1, -1 <= x <= 1], 
DirichletCondition[u[x, y, z] == -1, -5 <= x <= 4]}, 
u, {x, y, z} ∈ mr];

The ElementMeshSurfacePlot3D suggests that a plausible solution was computed:

With[{vals = uif["ValuesOnGrid"]}, 
 Legended[ElementMeshSurfacePlot3D[uif, ViewPoint -> {0, -4, 2}, 
   Axes -> True, AxesLabel -> {X, Y, Z}], 
  BarLegend[{ColorFunction /. Options[ElementMeshSurfacePlot3D] // 
     First, MinMax@vals}]]]

enter image description here

DensityPlot[uif[X, Y, 0.4], {X, -7, 7}, {Y, -4, 4}, PlotPoints -> 70, 
  PlotLegends -> Automatic, FrameLabel -> {X, Y, Z}, 
  PlotLabel -> "Slice at Z=0.4"] // Quiet

enter image description here

However, a more careful look suggests that my Dirichlet condition is not exactly met at X=0:

Plot[uif[X, 0, 0.4], {X, -7, 7}, AxesLabel -> {X, u}, 
  PlotLabel -> "Slice at Y=0, Z=0.4"] // Quiet

enter image description here

I would really appreciate your help fully understanding this example, as I also intend to solve FEM problems on BoundaryMeshRegions imported from STLs.

$\endgroup$
  • 2
    $\begingroup$ I'm also confused by the example from Wolfram. First problem : how is it possible that the Dirichlet condition DirichletCondition[u[x, y, z] == 1, z <= -1.3] is outside mr ??? $\endgroup$ – andre314 Aug 20 '17 at 16:49
  • 1
    $\begingroup$ The first part is most likely a typo. The equation should have been a Poisson equation; since the visualization kind of looks right, it was not noticed. The fact that it kind of looks right is because of a numerical curiosity. A perfect example that garbage in does not always mean garbage out.... The second issue you raise, I can not reproduce in 11.1 and I no longer have a 11.0 around to check. I get a flat line at 1. $\endgroup$ – user21 Aug 21 '17 at 17:32
  • $\begingroup$ @andre, what do you mean with outside? $\endgroup$ – user21 Aug 21 '17 at 17:40
  • $\begingroup$ @user21 I mean : -1.3 is under the lowest point of the shuttle $\endgroup$ – andre314 Aug 21 '17 at 17:43
  • 1
    $\begingroup$ The Example page is updated now. Thanks for the feedback. $\endgroup$ – user21 Sep 11 '17 at 19:58
7
$\begingroup$

You question has two parts. The example on the page you reference was not quite correct (EDIT- It's fixed now.) There are many ways to fix this. One idea is to have two distinct DirichletConditions like this:

mr = BoundaryDiscretizeGraphics[
   ExampleData[{"Geometry3D", "SpaceShuttle"}]];
uif = NDSolveValue[{Inactive[Laplacian][u[x, y, z], {x, y, z}] == 0, 
    DirichletCondition[u[x, y, z] == 1, z <= -1.3], 
    DirichletCondition[u[x, y, z] == 0, z >= 1]}, 
   u, {x, y, z} \[Element] mr];
Needs["NDSolve`FEM`"]
ElementMeshSurfacePlot3D[uif, Boxed -> False, ViewPoint -> {0, -4, 2}]

Better

Another way to fix this is to use a different PDE, e.g. use Poisson's equations in stead.

The fact that ElementMeshSurfacePlot3D did show a visualization that seemed valid is because the values computed do have slight variations. You will see that if you copy and paste the values in the notebook. What is surprising to me too is that this noise does seem to have a 'direction' that looks like a correct flow... Perfect example where garbage in does not result in garbage out. I'll try to get the web page updated.

For the second point: Maybe I am missing the point but the coordinate {0,0,0.4} is inside the region and not on the boundary.

emsf = ElementMeshSurfacePlot3D[uif, Boxed -> False, 
   ViewPoint -> {0, -4, 2}];
Show[emsf, Graphics3D[{Red, Ball[{0, 0, 0.4}, 0.5]}], 
 PlotRange -> {Automatic, {0, 5}, {-2, 5}}]

That, I think, explains why the point does not match the value on the boundary.

enter image description here

$\endgroup$
  • $\begingroup$ This is perfect, thank you, explains everything. $\endgroup$ – Alexander Erlich Aug 22 '17 at 15:20

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.