2
$\begingroup$

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

Mathematica graphics

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?

$\endgroup$
  • 2
    $\begingroup$ All your code seems to execute fine on my system (MMA 10.2 Win7-64). I obtain at the end an interpolating function with the full 45 s domain (it takes roughly 9 seconds to do so). I don't really know what to do with that function to check whether it is correct, but at least it seems to run error-free. $\endgroup$ – MarcoB Aug 11 '15 at 20:24
  • $\begingroup$ @MarcoB Thanks for looking at this. I have MMA 10.0.2 and Win7-64. It definitely does not give me answers in over 20 minutes and I can't abort it either. So perhaps it is a version problem. If someone could check it who has 10.0.2 that would be helpful. Thanks. $\endgroup$ – Hugh Aug 12 '15 at 2:54
  • $\begingroup$ Works fine for me in 10.0.2 on Linux and it takes about 11 seconds. Do you perhaps have a variable defined that should not be. If you copy your own code from this post into a fresh notebook and a fresh kernel, does it still not work? $\endgroup$ – user21 Aug 12 '15 at 7:56
  • $\begingroup$ @user21 and MarcoB I have checked yet again and done exactly as you suggested copying from here to a fresh kernel and it still gets stuck. Is this a diagnostic point: the displayed time starts small and produces some negative values before going positive and sticking at 0.00158...? Also if I set the time to less than this value i.e. 0.0015 the code still does not stop. As MMA freezes and I have to exit and restart I can't think of a way of capturing these time values. Thanks for your continued interest. $\endgroup$ – Hugh Aug 12 '15 at 13:09
4
$\begingroup$

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: Mathematica graphics

After restart: Mathematica graphics

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

Thanks for the help

$\endgroup$
  • $\begingroup$ +1 for making the effort to write a detailed self-answer with a nice link $\endgroup$ – LLlAMnYP Aug 21 '15 at 9:28

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.