Swinging Beam example from Help. Not Understanding/Not Working

Question

To help me understand finite element analysis with Mathematica I have been reading the Finite Element Method User Guide in Help and am stuck on the example in Coupled PDEs. Here there is a static beam problem and a swinging beam problem.

I have taken the code from the static beam problem and this works nicely

 Needs["NDSolve`FEM`"]

\[CapitalOmega] = ImplicitRegion[True, {x, y}];
mesh = ToElementMesh[\[CapitalOmega], {{0, 5}, {0, 1}}, 
   "MaxCellMeasure" -> 0.1];
vd = NDSolve`VariableData[{"DependentVariable", 
     "Space"} -> {{u, v}, {x, y}}];
sd = NDSolve`SolutionData["Space" -> ToNumericalRegion[mesh]];

bcDu0 = DirichletCondition[u[x, y] == 0, x == 0];
bcDv0 = DirichletCondition[v[x, y] == 0, x == 0];
bcNL = NeumannValue[-1, x == 5];
initBCs = 
  InitializeBoundaryConditions[vd, sd, {{bcDu0}, {bcDv0, bcNL}}];

diffusionCoefficients = 
  "DiffusionCoefficients" -> {{{{-(Y/(1 - \[Nu]^2)), 
        0}, {0, -((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2)))}}, {{0, -((
         Y \[Nu])/(1 - \[Nu]^2))}, {-((Y (1 - \[Nu]))/(
         2 (1 - \[Nu]^2))), 
        0}}}, {{{0, -((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2)))}, {-((
         Y \[Nu])/(1 - \[Nu]^2)), 
        0}}, {{-((Y (1 - \[Nu]))/(2 (1 - \[Nu]^2))), 
        0}, {0, -(Y/(1 - \[Nu]^2))}}}} /. {Y -> 10^3, \[Nu] -> 33/100};
initCoeffs = 
  InitializePDECoefficients[vd, sd, {diffusionCoefficients}];


methodData = InitializePDEMethodData[vd, sd];
discretePDE = DiscretizePDE[initCoeffs, methodData, sd];
split = {#[[1]] + 1 ;; #[[2]], #[[2]] + 1 ;; #[[3]]} &[
   methodData["IncidentOffsets"]];
discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd];
{load, stiffness, damping, mass} = discretePDE["SystemMatrices"];

DeployBoundaryConditions[{load, stiffness}, discreteBCs];

solution = LinearSolve[stiffness, load];
uif = ElementMeshInterpolation[{mesh}, solution[[split[[1]]]]];
vif = ElementMeshInterpolation[{mesh}, solution[[split[[2]]]]];
dmesh = ElementMeshDeformation[mesh, {uif, vif}];

I can plot results as follows using

 mesh["Wireframe"]
ContourPlot[uif[x, y], {x, y} \[Element] mesh, 
 ColorFunction -> "Temperature", AspectRatio -> Automatic]
ContourPlot[vif[x, y], {x, y} \[Element] mesh, 
 ColorFunction -> "Temperature", AspectRatio -> Automatic]
Show[{
  mesh["Wireframe"],
  dmesh["Wireframe"[
    "ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}]

to give

All this is good but it goes wrong for me when I move onto the swinging beam which modifies the above code. I can't get this to run even if I copy directly from Help. Here is the relevant set-up code from Help

    massCoefficients = "MassCoefficients" -> {{1, 0}, {0, 1}};
initCoeffs = 
  InitializePDECoefficients[vd, 
   sd, {diffusionCoefficients, massCoefficients}];
discretePDE = DiscretizePDE[initCoeffs, methodData, sd];
{load, stiffness, damping, mass} = discretePDE["SystemMatrices"];
rayleighDamping = 0.1*mass + 0.04*stiffness;
DeployBoundaryConditions[{load, stiffness, rayleighDamping, mass}, 
  discreteBCs];
dof = methodData["DegreesOfFreedom"];
init = dinit = ConstantArray[0, {dof, 1}];
sparsity = 
  ArrayFlatten[{{mass["PatternArray"], 
     mass["PatternArray"]}, {rayleighDamping["PatternArray"], 
     rayleighDamping["PatternArray"]}}];

It is the next bit of the solution code that gets stuck. The time monitor stops around 0.00158 so I have tried running the code to 0.0015 rather than 45 as in Help but it still gets stuck.

 Dynamic["time: " <> ToString[CForm[currentTime]]]
AbsoluteTiming[
 tif = NDSolveValue[{
    mass.u''[ t] + rayleighDamping.u'[ t] + stiffness.u[ t] == load
    , u[ 0] == init, u'[ 0] == dinit}, u, {t, 0, 0.0015}
   , Method -> {"EquationSimplification" -> "Residual"}
   , Jacobian -> {Automatic, Sparse -> sparsity}
   , EvaluationMonitor :> (currentTime = t;)
   ]]

Are there any suggestions for how to make this work?

Answer 1

It seems that my version of Mathematica was corrupted. I did a Shift+Control while starting Mathematica (see here for good instructions) and the code now works perfectly. Thanks to user21 and MarcoB for confirming that the code works. Wolfram Support also confirmed that the code works and suggested a restart.

Out of interest I compared a benchmark from before the Shift+Control to after the restart. The benchmark is found by running

    Needs["Benchmarking`"]
    BenchmarkReport[]

Before restart:

After restart:

As you can see I have gone from anti-penultimate to top-of-the-heap!

Thanks for the help