7
$\begingroup$

Here is a wooden board, with dimensions shown on the picture below. How we can use Mathematica's newly build-in finite element analysis features to show the different modes of its vibrations. Assuming the board is made of spruce.

enter image description here

$\endgroup$
  • 4
    $\begingroup$ Once again, interesting question but you are missing the fact that you should at least try something. $\endgroup$ – Öskå Aug 14 '14 at 18:32
  • 3
    $\begingroup$ @Mr.Wizard Would you also think that "please show me how to find the ground state of the 2d Kitaev model with Mathematica's linear algebra tools" is of reasonable scope? If not, why not? (in fact, mine is a much, much less technically involved question--unless we do include the spruce bit). $\endgroup$ – acl Aug 14 '14 at 18:53
  • 5
    $\begingroup$ @Mr.Wizard, before you close this, give me some time - I might get a Reversal if I wait a little longer.... $\endgroup$ – user21 Aug 14 '14 at 19:02
  • 1
    $\begingroup$ At the least the OP could provide (links to) the equations describing the vibration, or do we have to do the whole (home)work ourselves? $\endgroup$ – Sjoerd C. de Vries Aug 14 '14 at 22:49
  • 1
    $\begingroup$ Have a look at the example "Solve a Wave Equation in 2D" here. It shows how you can solve the 2D wave equation over a non-trivial region. For instance, if I understand your diagram correctly, your region would be RegionDifference[Rectangle[{0, 0}, {40, 76}], RegionUnion[Rectangle[{40/2 - 5/2, 76/2 - 68/2}, {40/2 + 5/2, 76/2 + 68/2}], Rectangle[{40/2 - 36/2, 76/2 - 5/2}, {40/2 + 36/2, 76/2 + 5/2}]]]. You would need to define an initial condition to start off the 2D wave solution. $\endgroup$ – Stephen Luttrell Aug 14 '14 at 22:51
66
$\begingroup$

The only reason I am attempting to answer this is to perhaps get a Reversal badge. There you go...

We will go slowly and this answer is the basis for what comes next. Let's start with two dimensions. You'll see why.

We create a rectangular region:

Needs["NDSolve`FEM`"]
mesh = ToElementMesh[FullRegion[2], {{0, 5}, {0, 1}}, "MeshOrder" -> 1, "MaxCellMeasure" -> 0.005];

We use the plane stress equations of a material sheet. For every material you need to set the Young's modulus and Poisson's ratio.

planeStress = {Inactive[
       Div][{{0, -((Y*ν)/(1 - ν^2))}, {-(Y*(1 - ν))/(2*(1 \
- ν^2)), 0}}.Inactive[Grad][v[t, x, y], {x, y}], {x, y}] + 
     Inactive[
       Div][{{-(Y/(1 - ν^2)), 
         0}, {0, -(Y*(1 - ν))/(2*(1 - ν^2))}}.Inactive[Grad][
        u[t, x, y], {x, y}], {x, y}], 
    Inactive[
       Div][{{0, -(Y*(1 - ν))/(2*(1 - ν^2))}, {-((Y*ν)/(1 \
- ν^2)), 0}}.Inactive[Grad][u[t, x, y], {x, y}], {x, y}] + 
     Inactive[
       Div][{{-(Y*(1 - ν))/(2*(1 - ν^2)), 
         0}, {0, -(Y/(1 - ν^2))}}.Inactive[Grad][
        v[t, x, y], {x, y}], {x, y}]} /. {Y -> 10^3, ν -> 33/100};

Make a coupled, time dependent PDE and constrain the region to not move at the left boundary.

pde2D = D[{u[t, x, y], v[t, x, y]}, t] + planeStress == {0, 0};
(* held fixed at left *)
bcs = DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, x==0];

What follows next, is pretty much the same from the other post.

(*We use NDSolve as a pre-processor:*)
{state} = 
  NDSolve`ProcessEquations[{pde2D, bcs, u[0, x, y] == 0, 
    v[0, x, y] == 0}, {u, v}, {t, 0, 1}, {x, y} \[Element] mesh, 
   Method -> {"PDEDiscretization" -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"FiniteElement"}}}];

(*Extract the finite element data:*)
femdata = state["FiniteElementData"];

initBCs = femdata["BoundaryConditionData"];
methodData = femdata["FEMMethodData"];
initCoeffs = femdata["PDECoefficientData"];

(*Set up the solution and variable data:*)
vd = methodData["VariableData"];
nr = ToNumericalRegion[mesh];
sd = NDSolve`SolutionData[{"Space" -> nr, "Time" -> 0.}];

(*Discretize the PDE and the boundary conditions:*)
discretePDE = DiscretizePDE[initCoeffs, methodData, sd];
discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd];
(*Extract the system matrices:*)
load = discretePDE["LoadVector"];
stiffness = discretePDE["StiffnessMatrix"];
damping = discretePDE["DampingMatrix"];

(*Deploy the boundary conditions:*)
DeployBoundaryConditions[{load, 
  stiffness, damping}, discreteBCs]

(*Set the number of X smallest eigen values we would like to compute \
but ignore the Dirichlet positions.*)

nDiri = If[Length[#] > 0, First[#], 0] &[
   Dimensions[discreteBCs["DirichletMatrix"]]];

numEigenToCompute = 5;
numEigen = numEigenToCompute + nDiri;

Now, things become hard. We solve the eigensystem.

(*Solve the eigen system: this is how you should do it*)
res = Eigensystem[{stiffness, damping}, -numEigen];

You will need patience. Play with the "Arnoldi" method and the shift. (left as an exercise). As a bad alternative you can play with

(* this may be a bit faster but is the dark side... *)
(*
mm=LinearSolve[damping,stiffness];
res=Eigensystem[mm,-numEigen];
*)

Further down I use "FEAST" as a solve.

For each eigenvalue we now have two eigenvectors. Once in the x-direction and one in the y-direction. So we post-process:

res = Reverse /@ res;
eigenValues = res[[1, nDiri + 1 ;; Abs[numEigen]]];
eigenVectors = res[[2, nDiri + 1 ;; Abs[numEigen]]];
(*res=Null;*)
inciOffs = methodData["IncidentOffsets"];
 spans = MapThread[Span, {Most[inciOffs] + 1, Rest[inciOffs]}];
 eigenVectors = 
  Transpose[
   Developer`ToPackedArray[eigenVectors[[All, #]] & /@ spans], {2, 1, 
    3}];
 eigenVectorsIF = Table[{}, {numEigenToCompute}, {Length[spans]}];
 Do[
  eigenVectorsIF[[i, j]] = 
   NDSolve`FEM`ElementMeshInterpolation[{mesh}, 
    eigenVectors[[i, j]]] 
  , {i, numEigenToCompute}, {j, Length[spans]}];
res = {eigenValues, eigenVectorsIF};

And visualize the first 5 eigenvectos in the x-direction:

Show[NDSolve`FEM`ElementMeshPlot3D[res[[2, #, 1]]["Coordinates"][[1]],
     NDSolve`FEM`ElementMeshDirective -> 
     Directive[EdgeForm[Gray], FaceForm[]]], 
   NDSolve`FEM`ElementMeshPlot3D[res[[2, #, 1]]], Boxed -> False, 
   Axes -> False, ImageSize -> 600] & /@ Range[numEigenToCompute];

Visualize the y-direction:

Show[NDSolve`FEM`ElementMeshPlot3D[res[[2, #, 2]]["Coordinates"][[1]],
     NDSolve`FEM`ElementMeshDirective -> 
     Directive[EdgeForm[Gray], FaceForm[]]], 
   NDSolve`FEM`ElementMeshPlot3D[res[[2, #, 2]]], Boxed -> False, 
   Axes -> False, ImageSize -> 600] & /@ Range[numEigenToCompute];

And the "breathing modes" enlarged by a factor. (I just invented that name - it might mean something else)

fact = 5;
Show[{
    NDSolve`FEM`ElementMeshPlot3D[res[[2, 2, 1]]["Coordinates"][[1]], 
     NDSolve`FEM`ElementMeshDirective -> 
      Directive[EdgeForm[Gray], FaceForm[]]],
    NDSolve`FEM`ElementMeshPlot3D[
     NDSolve`FEM`ElementMeshInterpolation[
      res[[2, #, 1]]["Coordinates"], 
      fact*Sqrt[Total[#["ValuesOnGrid"]^2 & /@ res[[2, #]]]]]]
    }, Boxed -> False, ImageSize -> 600] & /@ Range[numEigenToCompute]

enter image description here enter image description here enter image description here enter image description here enter image description here

The 3D case. Unfortunately, my local super-computer center is closed, but here is how to do it. Create a mesh and enlarge the features for now to create not too many elements.

base = {0, 0, 0};
h1 = 5;
h2 = 5;
w1 = 40;
l1 = 76;
cw1 = 5;
cl1 = 68;
cw2 = 36;
cl2 = 5;
offset1 = base + {(w1 - cw1)/2, (l1 - cl1)/2, 0};
offset2 = base + {(w1 - cw2)/2, (l1 - cl2)/2, 0};
offset3 = base + {(w1 - cw1)/2, (l1 - cl2)/2, 0};
ClearAll[rect]
rect[base_, w_, l_, h_] := {base + {0, 0, h}, base + {w, 0, h}, 
  base + {w, l, h}, base + {0, l, h}}
coords = ConstantArray[{0., 0., 0.}, 4 + 4 + 12 + 12];
coords[[{1, 2, 3, 4}]] = rect[base, w1, l1, 0];
coords[[{5, 6, 7, 8}]] = rect[base, w1, l1, h1];
coords[[{9, 10, 15, 16}]] = rect[offset1, cw1, cl1, h1];
coords[[{19, 12, 13, 18}]] = rect[offset2, cw2, cl2, h1];
coords[[{20, 11, 14, 17}]] = rect[offset3, cw1, cl2, h1];
coords[[20 + Range[12]]] = ({0, 0, h2} + #) & /@ 
   coords[[8 + Range[12]]];
bmesh = ToBoundaryMesh["Coordinates" -> coords, 
  "BoundaryElements" -> {QuadElement[{{1, 2, 3, 4}, {1, 2, 6, 5}, {2, 
       3, 7, 6}, {3, 4, 8, 7}, {4, 1, 5, 8},
      {5, 6, 10, 9}, {6, 12, 11, 10}, {6, 7, 13, 12}, {7, 15, 14, 
       13}, {7, 8, 16, 15}, {8, 18, 17, 16}, {8, 5, 19, 18}, {5, 9, 
       20, 19},
      Sequence @@ ({{9, 10, 11, 20}, {11, 12, 13, 14}, {14, 15, 16, 
           17}, {17, 18, 19, 20}, {20, 11, 14, 17}} + 12),
      Sequence @@ (Partition[Join[Range[9, 20]], 2, 1, 
          1] /. {i1_, i2_} :> {i1, i2, i2 + 12, i1 + 12})
      }]}]

If you want to visualize the boundary structure:

Show[
 bmesh["Wireframe"],
 bmesh["Wireframe"["MeshElement" -> "PointElements", 
   "MeshElementIDStyle" -> Red]]
 ];

Create the mesh:

mesh = ToElementMesh[bmesh, "MeshOrder" -> 1, 
   "MaxCellMeasure" -> 10];
mesh["Wireframe"]

enter image description here

Here is the PDE stress operator in 3D:

stressOperator[
  Y_, ν_] := {Inactive[
     Div][{{0, 0, -((Y*ν)/((1 - 2*ν)*(1 + ν)))}, {0, 0, 
       0}, {-Y/(2*(1 + ν)), 0, 0}}.Inactive[Grad][
      w[t, x, y, z], {x, y, z}], {x, y, z}] + 
   Inactive[
     Div][{{0, -((Y*ν)/((1 - 2*ν)*(1 + ν))), 
       0}, {-Y/(2*(1 + ν)), 0, 0}, {0, 0, 0}}.Inactive[Grad][
      v[t, x, y, z], {x, y, z}], {x, y, z}] + 
   Inactive[
     Div][{{-((Y*(1 - ν))/((1 - 2*ν)*(1 + ν))), 0, 
       0}, {0, -Y/(2*(1 + ν)), 0}, {0, 
       0, -Y/(2*(1 + ν))}}.Inactive[Grad][
      u[t, x, y, z], {x, y, z}], {x, y, z}], 
  Inactive[Div][{{0, 0, 0}, {0, 
       0, -((Y*ν)/((1 - 
              2*ν)*(1 + ν)))}, {0, -Y/(2*(1 + ν)), 
       0}}.Inactive[Grad][w[t, x, y, z], {x, y, z}], {x, y, z}] + 
   Inactive[
     Div][{{0, -Y/(2*(1 + ν)), 
       0}, {-((Y*ν)/((1 - 2*ν)*(1 + ν))), 0, 0}, {0, 0, 
       0}}.Inactive[Grad][u[t, x, y, z], {x, y, z}], {x, y, z}] + 
   Inactive[
     Div][{{-Y/(2*(1 + ν)), 0, 
       0}, {0, -((Y*(1 - ν))/((1 - 2*ν)*(1 + ν))), 0}, {0,
        0, -Y/(2*(1 + ν))}}.Inactive[Grad][
      v[t, x, y, z], {x, y, z}], {x, y, z}], 
  Inactive[Div][{{0, 0, 0}, {0, 
       0, -Y/(2*(1 + ν))}, {0, -((Y*ν)/((1 - 
              2*ν)*(1 + ν))), 0}}.Inactive[Grad][
      v[t, x, y, z], {x, y, z}], {x, y, z}] + 
   Inactive[
     Div][{{0, 0, -Y/(2*(1 + ν))}, {0, 0, 
       0}, {-((Y*ν)/((1 - 2*ν)*(1 + ν))), 0, 0}}.Inactive[
       Grad][u[t, x, y, z], {x, y, z}], {x, y, z}] + 
   Inactive[
     Div][{{-Y/(2*(1 + ν)), 0, 0}, {0, -Y/(2*(1 + ν)), 0}, {0,
        0, -((Y*(1 - ν))/((1 - 2*ν)*(1 + ν)))}}.Inactive[
       Grad][w[t, x, y, z], {x, y, z}], {x, y, z}]}

And the 3D PDE. no boundary conditions, it seems an unconstraint analysis is wanted.

(* choose your Y and ν -- no idea what the values for spruce are - 
   is it wet, dry, old, ...? *)
pde3D = D[{u[t, x, y, z], v[t, x, y, z], w[t, x, y, z]}, t] + 
    stressOperator[100, 1/3] == {0, 0, 0};
(* unconstraint? Yes! *)
bcs = Sequence[];

If you want a constraint analysis then you'd have to set DirichletCondition on the boundary (bcs = DirichletCondition[{u[x,y,z]==0,v[x,y,z]==0,w[x,y,z]==0},True])

Now, we use

{state} = 
  NDSolve`ProcessEquations[{pde3D, bcs, u[0, x, y, z] == 0, 
    v[0, x, y, z] == 0, w[0, x, y, z] == 0}, {u, v, w}, {t, 0, 
    1}, {x, y, z} \[Element] mesh, 
   Method -> {"PDEDiscretization" -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"FiniteElement"}}}];

for pre-processing. There may be a waring about no DirichletCondition or no NeumannValue; this is save to ignore in this case.

And now the same as above - it will take a long time.... I did not want to wait. (Also I did not want to think about how to visualize this in 3D... that's for you...)

When you do this do not forget that the result need to be sorted and post processed as in the 2D example above. i.e.

res = Reverse /@ res;
eigenValues = res[[1, nDiri + 1 ;; Abs[numEigen]]];
eigenVectors = res[[2, nDiri + 1 ;; Abs[numEigen]]];

Update:

The Eigensystem solution above takes about 450 seconds on my machine. You can use

AbsoluteTiming[
 res = Eigensystem[{stiffness, damping}, -numEigen, 
    Method -> {"FEAST", "Tolerance" -> 10^-6}];]

to get it down to 45 seconds. Which is a bit better.

Here are the deformations for the eigenmodes 7 to 10 - the first 6 are zero. (I must admit that that I am not sure if that makes sense)

res[[1]]
{0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.011403583383327644`, \
0.01526089137692353`, 0.05661022352859022`, 0.07266104128273859`}

And the visualizations:

MeshRegion[
   ElementMeshDeformation[mesh, res[[2, #]], 
    "ScalingFactor" -> 100]] & /@ Range[7, numEigenToCompute]

enter image description here enter image description here enter image description here enter image description here

When you run the 3D exmaple, please adjust the numEigenToCompute to be appropriate.

Anything else?

$\endgroup$
  • $\begingroup$ @m_goldberg, thanks for the edit. $\endgroup$ – user21 Aug 15 '14 at 6:18
  • $\begingroup$ @Öskå, the problem with Eigensystem as it is used here is that this calls a direct solver for the generalized eigenvalue problem. I need to find a better solution for that. You do not by any change have access to Abaqus/Ansys or the like and could try? $\endgroup$ – user21 Aug 15 '14 at 9:23
  • $\begingroup$ Let us continue this discussion in chat. $\endgroup$ – user21 Aug 15 '14 at 9:27
  • 1
    $\begingroup$ @Mr.Wizard, ;-) yea I did :-P - I am such a scoundrel. Thanks. $\endgroup$ – user21 Aug 19 '14 at 7:51
  • 2
    $\begingroup$ I either don't understand or disagree with your assertion, at least least for the examples given. This is a good answer yet the question has a negative score. The other Q&A also has good answers and it has a very positive score. While I have seen a question suddenly rise in popularity after a good answer is posted I think that is only natural given the additional interest (and views) that the answer can create; this is especially true for a question that is hard to understand but which when illuminated is interesting. $\endgroup$ – Mr.Wizard Aug 19 '14 at 8:25
8
$\begingroup$

This is not another answer but an attempt to repeat the calculations in the other answer in version 11. There are problems at the end for which I need help (edit: now sorted and corrected below). Note that there is a more direct method using NDEigensystem which is discussed here. There are issues with the eigen solver which may take a long time; this is discussed above and in the referenced question.

I am also writing this answer as it provides a workflow for solving a vibration eigenvalue problem where details need to be extracted. I make some comments that may help those used to doing vibration calculations.

Set up the stress operator for a general structural, time dependent systems. (Note that you can't just type in the textbook equations because there needs to be appropriate Inactive commands).

Needs["NDSolve`FEM`"]


ClearAll[stressOperator, u, v, w, t, x, y, z]; 
stressOperator[
  Y_, ν_] := {Inactive[
     Div][{{0, 0, -((Y*ν)/((1 - 2*ν)*(1 + ν)))}, {0, 0, 
       0}, {-Y/(2*(1 + ν)), 0, 0}}.Inactive[Grad][
      w[t, x, y, z], {x, y, z}], {x, y, z}] + 
   Inactive[
     Div][{{0, -((Y*ν)/((1 - 2*ν)*(1 + ν))), 
       0}, {-Y/(2*(1 + ν)), 0, 0}, {0, 0, 0}}.Inactive[Grad][
      v[t, x, y, z], {x, y, z}], {x, y, z}] + 
   Inactive[
     Div][{{-((Y*(1 - ν))/((1 - 2*ν)*(1 + ν))), 0, 
       0}, {0, -Y/(2*(1 + ν)), 0}, {0, 
       0, -Y/(2*(1 + ν))}}.Inactive[Grad][
      u[t, x, y, z], {x, y, z}], {x, y, z}], 
  Inactive[Div][{{0, 0, 0}, {0, 
       0, -((Y*ν)/((1 - 
              2*ν)*(1 + ν)))}, {0, -Y/(2*(1 + ν)), 
       0}}.Inactive[Grad][w[t, x, y, z], {x, y, z}], {x, y, z}] + 
   Inactive[
     Div][{{0, -Y/(2*(1 + ν)), 
       0}, {-((Y*ν)/((1 - 2*ν)*(1 + ν))), 0, 0}, {0, 0, 
       0}}.Inactive[Grad][u[t, x, y, z], {x, y, z}], {x, y, z}] + 
   Inactive[
     Div][{{-Y/(2*(1 + ν)), 0, 
       0}, {0, -((Y*(1 - ν))/((1 - 2*ν)*(1 + ν))), 0}, {0,
        0, -Y/(2*(1 + ν))}}.Inactive[Grad][
      v[t, x, y, z], {x, y, z}], {x, y, z}], 
  Inactive[Div][{{0, 0, 0}, {0, 
       0, -Y/(2*(1 + ν))}, {0, -((Y*ν)/((1 - 
              2*ν)*(1 + ν))), 0}}.Inactive[Grad][
      v[t, x, y, z], {x, y, z}], {x, y, z}] + 
   Inactive[
     Div][{{0, 0, -Y/(2*(1 + ν))}, {0, 0, 
       0}, {-((Y*ν)/((1 - 2*ν)*(1 + ν))), 0, 0}}.Inactive[
       Grad][u[t, x, y, z], {x, y, z}], {x, y, z}] + 
   Inactive[
     Div][{{-Y/(2*(1 + ν)), 0, 0}, {0, -Y/(2*(1 + ν)), 0}, {0,
        0, -((Y*(1 - ν))/((1 - 2*ν)*(1 + ν)))}}.Inactive[
       Grad][w[t, x, y, z], {x, y, z}], {x, y, z}]}

Now make the geometry

 base = {0, 0, 0};
    h1 = 5;
    h2 = 5;
    w1 = 40;
    l1 = 76;
    cw1 = 5;
    cl1 = 68;
    cw2 = 36;
    cl2 = 5;
    offset1 = base + {(w1 - cw1)/2, (l1 - cl1)/2, 0};
    offset2 = base + {(w1 - cw2)/2, (l1 - cl2)/2, 0};
    offset3 = base + {(w1 - cw1)/2, (l1 - cl2)/2, 0};
    ClearAll[rect]
    rect[base_, w_, l_, h_] := {base + {0, 0, h}, base + {w, 0, h}, 
      base + {w, l, h}, base + {0, l, h}}
    coords = ConstantArray[{0., 0., 0.}, 4 + 4 + 12 + 12];
    coords[[{1, 2, 3, 4}]] = rect[base, w1, l1, 0];
    coords[[{5, 6, 7, 8}]] = rect[base, w1, l1, h1];
    coords[[{9, 10, 15, 16}]] = rect[offset1, cw1, cl1, h1];
    coords[[{19, 12, 13, 18}]] = rect[offset2, cw2, cl2, h1];
    coords[[{20, 11, 14, 17}]] = rect[offset3, cw1, cl2, h1];
    coords[[20 + Range[12]]] = ({0, 0, h2} + #) & /@ 
       coords[[8 + Range[12]]];
    bmesh = ToBoundaryMesh["Coordinates" -> coords, 
       "BoundaryElements" -> {QuadElement[{{1, 2, 3, 4}, {1, 2, 6, 5}, {2,
             3, 7, 6}, {3, 4, 8, 7}, {4, 1, 5, 8}, {5, 6, 10, 9}, {6, 12, 
            11, 10}, {6, 7, 13, 12}, {7, 15, 14, 13}, {7, 8, 16, 15}, {8, 
            18, 17, 16}, {8, 5, 19, 18}, {5, 9, 20, 19}, 
           Sequence @@ ({{9, 10, 11, 20}, {11, 12, 13, 14}, {14, 15, 16, 
                17}, {17, 18, 19, 20}, {20, 11, 14, 17}} + 12), 
           Sequence @@ (Partition[Join[Range[9, 20]], 2, 1, 
               1] /. {i1_, i2_} :> {i1, i2, i2 + 12, i1 + 12})}]}];

Show the geometry

Show[bmesh["Wireframe"], 
 bmesh["Wireframe"["MeshElement" -> "PointElements", 
   "MeshElementIDStyle" -> Red]]]

Mathematica graphics

mesh = ToElementMesh[bmesh, "MeshOrder" -> 1, "MaxCellMeasure" -> 10];
mesh["Wireframe"]

Mathematica graphics

Set up the pde and give values for modulus of elasticity and Poisson ratio. Those users familiar to vibration may be surprised that the derivative with respect to time is only first order and not second (Newton's law is second). This is because the equations are converted to first order within the code and we are just trying to get the matrices. For this reason what is normally called the mass matrix is called the damping matrix here.

Y = 100;
ν = 1/3;
pde3D = D[{u[t, x, y, z], v[t, x, y, z], w[t, x, y, z]}, t] + 
    stressOperator[Y, ν] == {0, 0, 0};

Define the boundary conditions

bcs = Sequence[];

The equations are now processed.

{state} = 
  NDSolve`ProcessEquations[{pde3D, bcs, u[0, x, y, z] == 0, 
    v[0, x, y, z] == 0, w[0, x, y, z] == 0}, {u, v, w}, {t, 0, 
    1}, {x, y, z} ∈ mesh, 
   Method -> {"PDEDiscretization" -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"FiniteElement"}}}];

The finite element data are now extracted.

femdata = state["FiniteElementData"];
initBCs = femdata["BoundaryConditionData"];
methodData = femdata["FEMMethodData"];
initCoeffs = femdata["PDECoefficientData"];

The solution and variable data are extracted

vd = methodData["VariableData"];
nr = ToNumericalRegion[mesh];
sd = NDSolve`SolutionData[{"Space" -> nr, "Time" -> 0.}];

The PDE and boundary conditions are discretized

discretePDE = DiscretizePDE[initCoeffs, methodData, sd];
discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd];

The system matrices can now be extracted

load = discretePDE["LoadVector"];
stiffness = discretePDE["StiffnessMatrix"];
damping = discretePDE["DampingMatrix"];

The system matrices must now be modified to account for the boundary conditions. This is done "in place" so there are no new matrices

DeployBoundaryConditions[{load, stiffness, damping}, discreteBCs];

Now we can calculate the eigen system. I look for the smallest 10 frequencies (this is the -10). They are given in order of largest to smallest so they need to be reversed to get them from small to large. These are actually the frequencies squared because of the formulation as a standard eigenvalue problem. Take the square root to get the frequencies and divide by 2 Pi to get them in Hz.

AbsoluteTiming[res = Eigensystem[{stiffness, damping}, -10,
    Method -> {"FEAST", "Tolerance" -> 10^-6}];]

This takes about 25 seconds on my machine. Here are the results

TableForm[Reverse@res[[1]], TableHeadings -> {Automatic, None}]

Mathematica graphics

These agree with the original post so all is good. (This is edited from my original post, see comments, because I forgot to revers.) The first six are zero, as they should be, because there are no constraints so there are six modes of rigid body vibration at zero frequency.

$\endgroup$
  • $\begingroup$ Have you maybe overlooked the reversing of the result (as in the 2D case?) like so: res = Reverse /@ res; eigenValues = res[[1, nDiri + 1 ;; Abs[numEigen]]]; eigenVectors = res[[2, nDiri + 1 ;; Abs[numEigen]]]; I have added a note to that respect in the post. $\endgroup$ – user21 Sep 7 '16 at 14:19
  • $\begingroup$ You are, of course, correct and I had missed the Reverse. I have updated my answer. Thanks once again for your help. $\endgroup$ – Hugh Sep 8 '16 at 17:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.