# Want to get better quality out of a DensityPlot3D[]

The plot shown below was obtained from a similar program mentioned in the following post.

3D Elastic waves in a glass

How can I have a contour bar (color contours) along with the plot?

The code shown below generated the plot, but the colors are not good and the plot is neither continuous and nor uniform. It is too blurred. How I can change it to show high-resolution graphics? Also, how I can put labels on the axes; for example, $$x$$, $$y$$, and $$z$$?

Could anyone help me?

Needs["NDSolveFEM"]
L = .14; L1 = .01; r1 = .085/2; r2 = .055/
2; del = .006;(*cg=3962 m/s,3980,5100,5640*);
reg = RegionUnion[
ImplicitRegion[(r2 + (r1 - r2) (-L1)/(L - L1))^2 <=
x^2 + y^2 <= (r2 + (r1 - r2) (-L1)/(L - L1) + del)^2 &&
0 <= z <= L, {x, y, z}],
ImplicitRegion[(r2 + (r1 - r2) (-L1)/(L - L1))^2 <=
x^2 + y^2 <= (r2 + (r1 - r2) (-L1)/(L - L1) + del)^2 &&
0 <= z <= L, {x, y, z}]];
param = {Y -> 169*10^9, \[Nu] -> 0.3}; rho = 2332;
ClearAll[stressOperator];
stressOperator[
Y_, \[Nu]_] := {Inactive[
Div][{{0, 0, -((Y*\[Nu])/((1 - 2*\[Nu])*(1 + \[Nu])))}, {0, 0,
0}, {-Y/(2*(1 + \[Nu])), 0, 0}}.Inactive[Grad][
w[t, x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{0, -((Y*\[Nu])/((1 - 2*\[Nu])*(1 + \[Nu]))),
0}, {-Y/(2*(1 + \[Nu])), 0, 0}, {0, 0, 0}}.Inactive[Grad][
v[t, x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{-((Y*(1 - \[Nu]))/((1 - 2*\[Nu])*(1 + \[Nu]))), 0,
0}, {0, -Y/(2*(1 + \[Nu])), 0}, {0,
u[t, x, y, z], {x, y, z}], {x, y, z}],
Inactive[
Div][{{0, 0, 0}, {0,
0, -((Y*\[Nu])/((1 -
2*\[Nu])*(1 + \[Nu])))}, {0, -Y/(2*(1 + \[Nu])),
0}}.Inactive[Grad][w[t, x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{0, -Y/(2*(1 + \[Nu])),
0}, {-((Y*\[Nu])/((1 - 2*\[Nu])*(1 + \[Nu]))), 0, 0}, {0, 0,
0}}.Inactive[Grad][u[t, x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{-Y/(2*(1 + \[Nu])), 0,
0}, {0, -((Y*(1 - \[Nu]))/((1 - 2*\[Nu])*(1 + \[Nu]))),
0}, {0, 0, -Y/(2*(1 + \[Nu]))}}.Inactive[Grad][
v[t, x, y, z], {x, y, z}], {x, y, z}],
Inactive[
Div][{{0, 0, 0}, {0,
0, -Y/(2*(1 + \[Nu]))}, {0, -((Y*\[Nu])/((1 -
v[t, x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{0, 0, -Y/(2*(1 + \[Nu]))}, {0, 0,
0}, {-((Y*\[Nu])/((1 - 2*\[Nu])*(1 + \[Nu]))), 0,
0}}.Inactive[Grad][u[t, x, y, z], {x, y, z}], {x, y, z}] +
Inactive[
Div][{{-Y/(2*(1 + \[Nu])), 0, 0}, {0, -Y/(2*(1 + \[Nu])),
0}, {0, 0, -((Y*(1 - \[Nu]))/((1 -
w[t, x, y, z], {x, y, z}], {x, y, z}]};

{vals, funs} =
NDEigensystem[
stressOperator[169*10^9, 0.3] +
rho {D[u[t, x, y, z], {t, 2}], D[v[t, x, y, z], {t, 2}],
D[w[t, x, y, z], {t, 2}]} == {0, 0, 0}, {u, v, w},
t, {x, y, z} \[Element] reg, 15];
Abs[vals]/(2 Pi)
DensityPlot3D[
Re[funs[[15, 3]][x, y, z]], {x, y, z} ∈ reg,
ColorFunction -> "Rainbow",
Boxed -> False,
PlotLabel -> Row[{"f = ", Abs[vals[[15]]]/2/Pi}],
BoxRatios -> Automatic, PlotPoints -> 50]

• Hard to help you when you don't give us the code for funs[[15, 3]] Commented Dec 10, 2020 at 1:36
• Increase the value of PlotPoints and add PlotLegends -> Automatic Commented Dec 10, 2020 at 1:41
• Thanks. I provided the code in my main question.
– iman
Commented Dec 10, 2020 at 5:20

First, although it will not make a difference, you will want to look at your region. You take a RegionUnion of two identical implicit regions. Just one implicit region would suffice.

Next, if you look at the answers in the post that you linked, they suggest that the default mesher may not capture features well for cylindrical geometry. You can see that that it indeed applies to your case as well by executing the following function.

funs[[15, 3]]["ElementMesh"]["Wireframe"]


If you want to make pretty images, it is best to start with a mesh that captures the sharp features. You can follow one of the answers' approaches, or you can use the structured mesh approach that I will provide below.

# Structured hex mesh

## Helper functions

We import the required FEM package and define some helper functions that will allow us to create an annular structured hex mesh.

(*Import required FEM package*)
Needs["NDSolveFEM"];
pointsToMesh[data_] :=
MeshRegion[Transpose[{data}],
Line@Table[{i, i + 1}, {i, Length[data] - 1}]];
rp2DMesh[rh_, rv_, marker_] := Module[{sqr, crd, inc, msh, mrkrs},
sqr = RegionProduct[rh, rv];
crd = MeshCoordinates[sqr];
inc = Delete[0] /@ MeshCells[sqr, 2];
mrkrs = ConstantArray[marker, First@Dimensions@inc];
msh = ToElementMesh["Coordinates" -> crd,
]
rp3DMesh[r2D_, rv_, marker_ : 1] := Module[{rp, crd, inc, msh, mrkrs},
rp = RegionProduct[r2D, rv];
crd = MeshCoordinates[rp];
inc = Delete[0] /@ MeshCells[rp, 3];
mrkrs = ConstantArray[marker, First@Dimensions@inc];
msh = ToElementMesh["Coordinates" -> crd,
"MeshElements" -> {HexahedronElement[inc, mrkrs]}]
]
annularMesh[r_, th_, rh_, rv_, marker_] :=
Module[{r1, r2, th1, th2, m1, crd, melms, newcrd, em, mr},
{r1, r2} = r;
{th1, th2} = th;
m1 = rp2DMesh[rh, rv, marker];
crd = m1["Coordinates"];
melms = m1["MeshElements"];
newcrd = {#1 Cos[#2], #1 Sin[#2]} & @@@ ({r1 + (r2 - r1) #1,
th1 + (th2 - th1) #2} & @@@ crd);
em = ToElementMesh["Coordinates" -> newcrd, "MeshElements" -> melms];
mr = MeshRegion@em;
<|"Mesh" -> em, "Region" -> mr|>
]


## Construct structured annular cylindrical hex mesh

The following code will construct a nice structured annular cylindrical hex mesh. You can course and or refine the mesh if desired—the parameters used in this mesh required about 16 GB of memory to solve.

subs = {L -> .14, L1 -> .01, r1 -> .085/2, r2 -> .055/2, del -> .006};
rinner = (r2 + (r1 - r2) (-L1)/(L - L1)) /. subs;
router = (r2 + (r1 - r2) (-L1)/(L - L1) + del) /. subs;
l = L /. subs;
nAngular = 90;
theta = 360 Degree;
afrac = theta/(360 Degree);
rv = pointsToMesh[Subdivide[0, 1, nAngular afrac]];
m = annularMesh[{rinner, router}, {0 Degree, theta}, rh, rv, 1];
mesh = rp3DMesh[m["Region"],
pointsToMesh[Subdivide[0, L, 40] /. subs]];
bmesh = ToBoundaryMesh[mesh];
groups = bmesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp
bmesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]]


As you can see, the structured mesh captures these sharp features at the end caps.

# Solve the Eigensystem

To solve the Eigensystem, you need to make a slight modification to NDEigensystem, as shown below:

{vals, funs} =
NDEigensystem[
stressOperator[169*10^9, 0.3] +
rho {D[u[t, x, y, z], {t, 2}], D[v[t, x, y, z], {t, 2}],
D[w[t, x, y, z], {t, 2}]} == {0, 0, 0}, {u, v, w},
t, {x, y, z} ∈ mesh, 15];


# Visualizations

## DensityPlot3D

DensityPlot3Dis probably not the most efficient way to view objects with thin walls for three reasons. First, DensityPlot3D requires much memory versus other possible visualization techniques. Second, resolving thin regions becomes doubly difficult. Third, it can take a long time. The following code produces a density plot that looks superior to that in the OP, but in my opinion, it is not a great visualization.

DensityPlot3D[Re[funs[[15, 3]][x, y, z]], {x, y, z} ∈ mesh,
ColorFunction -> "Rainbow", Boxed -> False,
PlotLabel -> Row[{"f = ", Abs[vals[[15]]]/2/Pi}],
BoxRatios -> Automatic,
RegionFunction ->
Function[{x, y, z,
f}, (x^2 + y^2 <= router^2 || x^2 + y^2 >= rinner^2) &&
0 <= z <= l], PlotPoints -> 200, PlotLegends -> Automatic,
AxesLabel -> {"x", "y", "z"}]


## SliceDensityPlot3D

SliceDensityPlot3Dshould be the function to create clip planes and iso-surfaces. As you can see, if you put exact values for the curved surfaces, the results are not very good.

Show[SliceDensityPlot3D[
Re[funs[[15, 3]][x, y, z]], {x^2 + y^2 == rinner^2,
x^2 + y^2 == router^2, z == 0, z == l}, {x, y, z} ∈ mesh,
ColorFunction -> "Rainbow", Axes -> False,
PlotLabel -> Row[{"f = ", Abs[vals[[15]]]/2/Pi}],
BoxRatios -> Automatic, PlotPoints -> {100, 100, 40},
PlotLegends -> Automatic], Boxed -> False]


If you put a slight offset in for each of the radii, the results are much better. However, if you look closely, you can see a slight offset at the end caps.

Show[SliceDensityPlot3D[
Re[funs[[15, 3]][x, y, z]], {x^2 + y^2 == 1.01 rinner^2,
x^2 + y^2 == 0.99 router^2, z == 0, z == l}, {x, y, z} ∈
mesh, ColorFunction -> "Rainbow", Axes -> False,
PlotLabel -> Row[{"f = ", Abs[vals[[15]]]/2/Pi}],
BoxRatios -> Automatic, PlotPoints -> {100, 100, 40},
PlotLegends -> Automatic], Boxed -> False]


The advantage of this method is that it is much faster than a volumetric rendering. The visualization has a much higher resolution and better contrast on the boundary surfaces than the volumetric rendering.

## Custom plotting

One advantage of using a structured mesh is selecting nodal points where the interpolation functions are interpreted. First, it will create some helper functions to create a variety of DensityPlot's.

### Helper functions

(*Helper functions for various ListDensityPlots*)
eps = 0.000001;
epsf = eps/1000;
fPos[x_, v_] := Flatten@Position[x, _?(Abs[# - v] < eps &), 1]
Clear[ldprfn]
ldprfn[ridx_, cf_ : "Rainbow", opac_ : 0.75] :=
Module[{cd, ids, data, app, ldpg, ldpo, ldp, gc, gco, cylo, cyl,
vcpos},
cd = If[StringQ[cf], ColorData[cf], cf];
ids = ridxfn[ridx];
data = {θs[[ids]], zs[[ids]], zvals[[ids]]} // Transpose;
app = MapAt[-# &, data[[fPos[θs[[ids]], Pi]]], {All, 1}];
data = Join[data, app];
ldpg = ListDensityPlot[data, ColorFunction -> GrayLevel,
PlotRange -> zmnmx];
vcpos =
Position[ldpg, HoldPattern[VertexColors -> _]]~Join~{2} // Flatten;
ldpo = ReplacePart[ldpg, vcpos -> (cd /@ ldpg[[## & @@ vcpos]])];
ldp = ReplacePart[ldpg,
vcpos -> (Opacity[opac, #] & /@ (cd /@ ldpg[[## & @@ vcpos]]))];
gco = First@Cases[ldpo, _GraphicsComplex, Infinity];
gc = First@Cases[ldp, _GraphicsComplex, Infinity];
cylo = Graphics3D[
gco /. {θ_, z_} /;
AllTrue[{θ, z},
DeveloperRealQ] :> {uniqRs[[ridx]] Cos[θ],
uniqRs[[ridx]]  Sin[θ], z}];
cyl = Graphics3D[
gc /. {θ_, z_} /;
AllTrue[{θ, z},
DeveloperRealQ] :> {uniqRs[[ridx]] Cos[θ],
uniqRs[[ridx]]  Sin[θ], z}];
<|"Data" -> data, "DensityPlot" -> ldp, "DensityPlotOpaque" -> ldpo,
"CylDensityPlot" -> cyl, "CylDensityPlotOpaque" -> cylo|>
]
(*Z level ListDensityPlot by index*)
Clear[ldpzfn]
ldpzfn[zidx_, cf_ : "Rainbow", opac_ : 0.75] :=
Module[{cd, ids, data, dataxy, app, lp, ldpg, ldpo, ldp, gc, gco,
cylo, cyl, vcpos, z},
cd = If[StringQ[cf], ColorData[cf], cf];
ids = zidxfn[zidx];
data = {xs[[ids]], ys[[ids]], zvals[[ids]]} // Transpose;
ldpg = ListDensityPlot[data, ColorFunction -> GrayLevel,
PlotRange -> zmnmx,
RegionFunction ->
Function[{x, y, z}, 1.0001*rinner <= Norm[{x, y}] <= router]];
vcpos =
Position[ldpg, HoldPattern[VertexColors -> _]]~Join~{2} // Flatten;
ldpo = ReplacePart[ldpg, vcpos -> (cd /@ ldpg[[## & @@ vcpos]])];
ldp = ReplacePart[ldpg,
vcpos -> (Opacity[opac, #] & /@ (cd /@ ldpg[[## & @@ vcpos]]))];
gco = First@Cases[ldpo, _GraphicsComplex, Infinity];
gc = First@Cases[ldp, _GraphicsComplex, Infinity];
cylo = Graphics3D[
gco /. {x_, y_} /; AllTrue[{x, y}, DeveloperRealQ] :> {x, y,
uniqZs[[zidx]]}];
cyl = Graphics3D[
gc /. {x_, y_} /; AllTrue[{x, y}, DeveloperRealQ] :> {x, y,
uniqZs[[zidx]]}];
<|"Data" -> data, "DensityPlot" -> ldp, "DensityPlotOpaque" -> ldpo,
"CylDensityPlot" -> cyl, "CylDensityPlotOpaque" -> cylo|>
]


### Collect solution data and defined functions pick off radial and Z coordinates

(*Eigen mode*)
mode = 15;
(*Get coordinates and solution data*)
crd = mesh["Coordinates"];
{xs, ys, zs} = crd // Transpose;
zvals = Re[funs[[mode, 3]]["ValuesOnGrid"]];
zmnmx = MinMax[zvals];
(*Create substitutions to map between cylindrical and Cartesian \
coordinates*)
cart2cyl = {x_, y_, z_} ->
CoordinateTransform["Cartesian" -> "Cylindrical", {x, y, z}];
cyl2cart = {r_, θ_, z_} ->
CoordinateTransform[
"Cylindrical" -> "Cartesian", {r, θ, z}];
cylcrds = crd /. cart2cyl;
{rs, θs, zs} = cylcrds // Transpose;
(*Use nearest functions to pick off selected radial and Z coordinates*)
rnf = Nearest[rs -> "Index"];
uniqRs = Union[Round[rs, epsf]];
ridxfn = rnf[uniqRs[[#]], {All, eps}] &;
znf = Nearest[zs -> "Index"];
uniqZs = Union[Round[zs, epsf]];
zidxfn = znf[uniqZs[[#]], {All, eps}] &;
(*A custom highly banded color map*)
img = Uncompress[Import["https://pastebin.com/raw/X1RTREmc"]];
dims = ImageDimensions[img];
colors2 =
RGBColor[#] & /@
ImageData[img][[IntegerPart@(dims[[2]]/2), 1 ;; -1]];


### Create cylindrical bounding surface with custom color map

We will now conduct tests on our workflow using a custom color map with a highly banded structure. The following will show the bounding surfaces with both opaque and translucent coloring.

cf = (Blend[colors2, #] &);
(*cf = "Rainbow";*)
ldri = ldprfn[1, cf];
ldro = ldprfn[5, cf];
ldzb = ldpzfn[1, cf];
ldzt = ldpzfn[uniqZs // Length, cf];
Show[{ldri["CylDensityPlotOpaque"], ldro["CylDensityPlotOpaque"],
ldzb["CylDensityPlotOpaque"], ldzt["CylDensityPlotOpaque"]},
Boxed -> False]
Show[{ldri["CylDensityPlot"], ldro["CylDensityPlot"],
ldzb["CylDensityPlot"], ldzt["CylDensityPlot"]}, Boxed -> False]


These results appear to work as intended in the highly banded color map gives a nice banded structure without resorting to mesh lines or a contour plot.

### Two-dimensional ListDensityPlot

Fundamentally, the basic approach is to take a two-dimensional ListDensityPlot, converted into a GraphicsComplex and perform a 2D to 3D coordinate transformation.

You can see going from the inner wall to the outer wall (left to right) the ListDensityPlot in 2D.

ldprs = ldprfn[#, (Blend[colors2, #] &)] & /@ Range[uniqRs // Length];
GraphicsRow[#["DensityPlotOpaque"] & /@ ldprs, ImageSize -> 1000]


### Peeling the onion

With this approach, it is relatively straightforward to create a Manipulate function that will allow one to peel back the outer layers to view the inner layers. An example is shown below:

(* Surface Association for CheckBoxBar *)
rsurfs = Range[uniqRs // Length];
surfaces =
ldprs[[#]]["CylDensityPlotOpaque"] & /@ rsurfs];
surfaces =
Append[surfaces, {"top" -> ldzt["CylDensityPlotOpaque"],
"bot" -> ldzb["CylDensityPlotOpaque"]}];
(* Visualization *)
Manipulate[Show[{
{Append[choices, {"top", "bot"}] /. surfaces}},
Boxed -> False], {{choices, rsurfs}, rsurfs, CheckboxBar},
ControlPlacement -> Top]


• Dear Tim Thanks so much for your help.I encountered this error. Could you please help or send the total program? The version of my software is 12.Thanks Set::shape: Lists {vals,funs} and NDEigensystem[{stressOperator[169000000000,0.3]+rho (u^(2,0,0,0))[t,x,y,z],stressOperator[169000000000,0.3]+rho (v^(2,0,0,0))[t,x,y,z],stressOperator[169000000000,0.3]+rho (w^(2,0,0,0))[t,x,y,z]}=={0,0,0},{u,v,w},t,{x,y,z}[Element]mesh,15] are not the same shape.
– iman
Commented Dec 15, 2020 at 18:42
• @iman I will check it out. I am also planning to update the answer with an additional approach. Commented Dec 15, 2020 at 19:33
• @iman I started a new session of Mathematica and copied the code from this post into a new notebook. It executed fine. I am on version 12.1.1 (Windows). If you start a new session and execute the following NotebookOpen["https://pastebin.com/raw/ikbhK3wY"] , I think it should work because there should not be that big of a difference in the versions. Commented Dec 15, 2020 at 21:41
• Dear Tim, Thanks so much. Could you please contact me through my email?. I have a personal question that prefer to not stated here. Best
– iman
Commented Dec 15, 2020 at 22:56
• @iman You might not want to post your email directly. You could try Compress and return the output string like In[1]:= Compress["[email protected]"] Out[1]= "1:eJxTTMoPChZiYGAoLU4tykvMTXVIz03MzNFLzs8FAG49CME=". Commented Dec 16, 2020 at 2:01