# Plot FEM-result NDSolveValue in the FEM-mesh

In a simple pde-problem, which I solve in a predefined mesh using NDSolveValue I would like to plot the result in the mesh NDSolve uses:

<< NDSolveFEM
netz = ToElementMesh[Rectangle[], "MeshElementType" -> TriangleElement, MaxCellMeasure -> {"Length" -> .1}]

u = NDSolveValue[{Laplacian[\[Psi][x, y], {x, y}] ==1 + \[Psi][x, y] + NeumannValue[0, x == 1] ,\[Psi][0, y] == 0, \[Psi][x, 0] == 0, \[Psi][x, 1] == 0}, \[Psi] ,Element[{x, y}, netz], Method -> "FiniteElement"]


If I plot the result u[x,y]

Plot3D[u[x, y], {x, 0, 1}, {y, 0, 1}, Mesh -> All]


I get a foursided mesh even though the solution mesh

u["ElementMesh"]["Wireframe"]


is triangular. How can I force Plot3D to show the right mesh?

• A good point. This is easy to make with the ContourPlot such as Show[{ContourPlot[u[x, y], {x, 0, 1}, {y, 0, 1}], netz["Wireframe"] }] but this is not what you are asking about. Dec 15, 2017 at 14:53
• @Alexei Boulbitch: Thanks, something like Plot3D[u[x, y] , Element[{x, y}, netz]] would help Dec 15, 2017 at 16:20

You can use:

Plot3D[u[x, y], {x, y} \[Element] u["ElementMesh"], Mesh -> All,
MaxRecursion -> 0]


In version 11.0 there was a bug that you can work around with something like:

Plot3D[u[x, y], {x, y} \[Element]
MeshRegion[u["ElementMesh"]]["MakeLinear"], Mesh -> All,
MaxRecursion -> 0]

• Thank you, this is the answer I'm looking for. In Mathematica v11.0.1.0 I get the error Plot3D::idomdim: {x,y}\[Element]u[ElementMesh] does not have a valid dimension as a plotting domain. Which is the version you use? Dec 28, 2017 at 12:43
• @UlrichNeumann, I am using 11.2 and I vaguely remember that there was such a bug. See if the edit helps you. Dec 28, 2017 at 13:00
• That's it, thank you! Dec 28, 2017 at 13:18

Plot3D will resample points. If you want to maintain the mesh structure of netz, you could specify z-coords using u function. I specify "MeshOrder" -> 1 to get triangles:

<< NDSolveFEM
netz = ToElementMesh[Rectangle[],
"MeshElementType" -> TriangleElement,
MaxCellMeasure -> {"Length" -> .1}, "MeshOrder" -> 1]

u = NDSolveValue[{Laplacian[\[Psi][x, y], {x, y}] ==
1 + \[Psi][x, y] + NeumannValue[0, x == 1], \[Psi][0, y] ==
0, \[Psi][x, 0] == 0, \[Psi][x, 1] == 0}, \[Psi],
Element[{x, y}, netz], Method -> "FiniteElement"]

triangles = u["ElementMesh"]["MeshElements"][[1, 1]];
coords = u["ElementMesh"]["Coordinates"];

Graphics3D[
GraphicsComplex[{##, u[##]} & @@@
coords, {EdgeForm[{GrayLevel[0.2]}],
Directive[Specularity[GrayLevel[1], 3], RGBColor[
0.880722, 0.611041, 0.142051],
Lighting -> {{"Ambient", RGBColor[
0.30100577, 0.22414668499999998,
0.090484535]}, {"Directional", RGBColor[
0.2642166, 0.18331229999999998, 0.04261530000000001],
ImageScaled[{0, 2, 2}]}, {"Directional", RGBColor[
0.2642166, 0.18331229999999998, 0.04261530000000001],
ImageScaled[{2, 2, 2}]}, {"Directional", RGBColor[
0.2642166, 0.18331229999999998, 0.04261530000000001],
ImageScaled[{2, 0, 2}]}}], Polygon[triangles]}],
BoxRatios -> {1, 1, .4}]

• Thank you for your answer. Trying to understand the look at the 3D-Plot shows a regular mesh structure with holes inside??? In triangles I would have expected some triples of knots , but I found sublists of length 6. Dec 15, 2017 at 15:46
• Now I got it! The holes disappear with Polygon[triangles[[All,{1,2,3}]] ] . But the solution obtained in this way is unfortunately only a workaround. What I'm looking for is the projection of the triangle mesh onto the surface u[x,y] Dec 15, 2017 at 16:11
• did you specify "MeshOrder" -> 1 in netz? Dec 15, 2017 at 18:43
• Not in my first attempt... Dec 16, 2017 at 12:45