Animating Boundary Elements

I am trying to visualise the eigenvectors of a thick disk. I manage to calculate them and look at one but fail to look at more -probably due to using too much memory. I start by making a mesh and then use NDEigensystem.

Needs["NDSolveFEM"]
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[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[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[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[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[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,
v[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[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[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 - 2*\[Nu])*(1 + \[Nu])))}}.Inactive[
Grad][w[x, y, z], {x, y, z}], {x, y, z}]}
r = 0.3; (* radius of cylinder *)
L = 0.2; (* length of cylinder *)
Y = 10^6 ;(* Ratio of modulus of elasticity to density *)
\[Nu] = 0.2; (* Poission ratio *)

mesh = ToElementMesh[Cylinder[{{0, 0, 0}, {0, 0, L}}, r],
MaxCellMeasure -> .5];
mesh["Wireframe"]


Here I calculate the eigenvalues and eigenvectors. The first 6 are zero (approximatly) because they are rigid body modes.

{vals, vecs} =
NDEigensystem[
stressOperator[Y, \[Nu]], {u, v, w}, {x, y, z} \[Element] mesh, 14];
freqs = 1/(2 \[Pi]) Sqrt[Abs[Chop[vals]]];
TableForm[freqs]


Now I have the problem of animating the eigenvectors. The following module does this but is wrong because I am animating every node. Only the boundary nodes can be seen and thus only these should be animated. If I make too many animations using this DynamicModule then Mathematica freezes.

    ClearAll[animateEigenVectors];
animateEigenVectors[{freqs_, vecs_}, n_] := DynamicModule[{dmesh},
Animate[
dmesh =
ElementMeshDeformation[mesh, vecs[[n]],
"ScalingFactor" -> 0.01 Cos[t]];
Show[{dmesh[
"Wireframe"[
"ElementMeshDirective" ->
Directive[EdgeForm[Gray], FaceForm[Lighter[Blue]]]]]},
PlotLabel ->
"Natural Frequency = " <> ToString[freqs[[n]]] <>
" Hz  Mode = " <> ToString[n], Boxed -> False,
PlotRange -> 1.5 {{-r, r}, {-r, r}, {0, L}}
],
{t, 0, 2 \[Pi]}, SaveDefinitions -> True]
]

animateEigenVectors[{freqs, vecs}, 7]


(I can't do the fancy animation bit in SE yet.) Is it possible to extract and animate just the boundary elements? Thanks

• Could you provide the code you used to make the mesh? Dec 9, 2015 at 17:00
• @user21 Rats I forgot two lines of code. Now edited in. Thanks for looking at this.
– Hugh
Dec 9, 2015 at 21:31

Try this:

ClearAll[animateEigenVectors];
animateEigenVectors[{freqs_, vecs_}, n_] :=
DynamicModule[{dmesh}, Animate[
(*dmesh=ElementMeshDeformation[mesh,vecs[[n]],
"ScalingFactor"\[Rule]0.01 Cos[t]];*)

scalefactor = 0.01 Cos[t];
move =
scalefactor*
Transpose[
DeveloperToPackedArray[#["ValuesOnGrid"] & /@ vecs[[n]]]];
dmesh =
ToBoundaryMesh["Coordinates" -> mesh["Coordinates"] + move,
"BoundaryElements" -> mesh["BoundaryElements"],
"CheckIntersections" -> False];

Show[{dmesh[
"Wireframe"[
"ElementMeshDirective" ->
Directive[EdgeForm[Gray], FaceForm[Lighter[Blue]]]]]},
PlotLabel ->
"Natural Frequency = " <> ToString[freqs[[n]]] <>
" Hz  Mode = " <> ToString[n], Boxed -> False,
PlotRange -> 1.5 {{-r, r}, {-r, r}, {0, L}}], {t, 0, 2 \[Pi]},
SaveDefinitions -> True]]
`

This generates a boundary mesh only and switches off the expensive intersection checking. Since I could not reproduce the original problem, let me know if this helps to solve the issue you see.