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,

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})
```