Animating Boundary Elements

Question

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["NDSolve`FEM`"]
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, 
       0, -Y/(2*(1 + \[Nu]))}}.Inactive[Grad][
      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,
        0, -Y/(2*(1 + \[Nu]))}}.Inactive[Grad][
      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 - 
              2*\[Nu])*(1 + \[Nu]))), 0}}.Inactive[Grad][
      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

Answer 1

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[
      Developer`ToPackedArray[#["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.