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?
