# How to DensityPlot3D the solution of a PDE?

I solve a Dirichlet problem on a parallelepiped in 3D

p = Parallelepiped[{0, 0, 0}, {{1, 0, 0}, {1, 1, 0}, {0, 1, 1}}];
sol = NDSolve[{Laplacian[u[x, y, z], {x, y, z}] == 0,
DirichletCondition[u[x, y, z] == x^4 + y^2 + z, True]}, u[x, y, z], {x, y, z} \[Element] p]


{{u[x,y,z]->InterpolatingFunction[Domain: {{0.,2.},{0.,2.},{0.,1.}} Output: scalar][x,y,z]}}

I'd like to vizualize that solution. Following the documentation (see the Applications section), I try

DensityPlot3D[ u[x, y, z] /. sol, {x, 0, 2}, {y, 0, 2}, {z, 0, 1}]


The above command works, but produces an empty plot without any error messages. My next attempt is

Clear[x,y,z,u];DensityPlot3D[Evaluate[u[x, y, z] /. sol], {x, 0, 2}, {y, 0, 2}, {z, 0, 1}]


with an unsatisfactory result as well as

Clear[x,y,z,u];DensityPlot3D[Evaluate[u[x, y, z] /. sol], {x, 0, 2}, {y, 0, 2}, {z, 0, 1},
PlotPoints -> 20, WorkingPrecision -> 12]


The slices of sol are plotted well, for example,

Plot3D[(u[x, y, z] /. sol) /. z -> 0.2, {x, 0, 2}, {y, 0, 2}]


• Clear[x, y, z, u]; sol1 = NDSolveValue[{Laplacian[u[x, y, z], {x, y, z}] == 0, DirichletCondition[u[x, y, z] == x^4 + y^2 + z, True]}, u[x, y, z], {x, y, z} \[Element] p];DensityPlot3D[u[x, y,z] /. sol1, {x, 0, 2}, {y, 0, 2}, {z, 0, 1}]  also produces an empty plot. Mar 3, 2021 at 15:51
• this is because you should use NDSolve here if you want to use it a as rule. Mar 4, 2021 at 7:10

This seems to work?

Clear[x, y, z, u,p,U];
p = Parallelepiped[{0, 0, 0}, {{1, 0, 0}, {1, 1, 0}, {0, 1, 1}}];
U = NDSolveValue[{Laplacian[u[x, y, z], {x, y, z}] == 0,
DirichletCondition[u[x, y, z] == x^4 + y^2 + z, True]},
u, {x, y, z} \[Element] p];


Then

U[x, y, z] // Evaluate, {x, y,
z} \[Element] p, Contours -> 5, MeshFunctions -> (Norm[{#1, #2, \
#3}] &)]


Update

As requested by OP, you can find the range of values via:

   NMaximize[U[x, y, z], {x, y, z} \[Element] p] // Quiet


{12.713,{x->1.80957,y->1.28142,z->0.30654}}

   NMinimize[U[x, y, z], {x, y, z} \[Element] p] // Quiet


{1.07235,{x->0.996752,y->0.101257,z->0.0664789}}

So that

  ContourPlot3D[U[x, y, z] // Evaluate, {x, y, z} \[Element] p,
Contours -> (val = {1,2,5}),
MeshFunctions -> (Norm[{#1, #2, #3}] &),
ContourStyle -> ColorData["Heat"] /@ (Range[5]/5.),
PlotLegends -> val]


• Thank you. You first plot is a duplication of my first plot. Which surfaces are drawn in your second plot? Can you demonstrate more such surfaces? The behavior of the solution near the origin is not clear to me. Unfortunately, I cannot accept your answer in its present form. Mar 4, 2021 at 7:13
• +1. However, there are still deficiences: some surfaces are partly out of p and the levels of U are not known. Mar 4, 2021 at 7:32
• ContourPlot3D[U[x, y, z] // Evaluate, {x, y, z} \[Element] p, Contours -> {1, 2, 5}] suits better to me. I'll be wating some time for other answers (maybe, through ListDensityPlot3D). Mar 4, 2021 at 7:44