# Finite element boundary breaking

I am trying to model a cantilever beam which can vibrate. On the left the beam is clamped. I am largely following user21 here and also the example in help. I start by doing a static beam which works then I work out the vibration modes (eigenvectors) of a beam with no boundary constraint, this works, and then it all goes wrong when I try to do the eigenvectors of the clamped beam. I start with some helper modules.

Needs["NDSolveFEM"];
(*plane stress equations two versions with and without time *)
ClearAll[planeStress];
planeStress::usage =
"planeStress[u,v,x,y,Y,\[Nu]] (no time) or planeStress[u,v,t,x,y,Y,\
\[Nu]] (with time) Y is modulus of elasticity and \[Nu] is Poission \
ratio. ";
planeStress[u_, v_, x_, y_,
Y_, \[Nu]_] := {Inactive[
Div][{{0, -((Y*\[Nu])/(1 - \[Nu]^2))}, {-(Y*(1 - \[Nu]))/(2*(1 - \
\[Nu]^2)), 0}}.Inactive[Grad][v[x, y], {x, y}], {x, y}] +
Inactive[
Div][{{-(Y/(1 - \[Nu]^2)),
0}, {0, -(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))}}.Inactive[Grad][
u[x, y], {x, y}], {x, y}],
Inactive[Div][{{0, -(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))}, \
{-((Y*\[Nu])/(1 - \[Nu]^2)), 0}}.Inactive[Grad][u[x, y], {x, y}], {x,
y}] + Inactive[
Div][{{-(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2)),
v[x, y], {x, y}], {x, y}]}
planeStress[u_, v_, t_, x_, y_,
Y_, \[Nu]_] := {Inactive[
Div][{{0, -((Y*\[Nu])/(1 - \[Nu]^2))}, {-(Y*(1 - \[Nu]))/(2*(1 - \
\[Nu]^2)), 0}}.Inactive[Grad][v[t, x, y], {x, y}], {x, y}] +
Inactive[
Div][{{-(Y/(1 - \[Nu]^2)),
0}, {0, -(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))}}.Inactive[Grad][
u[t, x, y], {x, y}], {x, y}],
Inactive[Div][{{0, -(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))}, \
u[t, x, y], {x, y}], {x, y}] +
Inactive[
Div][{{-(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2)),
v[t, x, y], {x, y}], {x, y}]}

ClearAll[finiteElementData];
finiteElementData::usate =
"finiteElementData[start] Extract the finite element data";
finiteElementData[state_] :=
Module[{femdata, initBCs, methodData, initCoeffs, vd, sd,
discretePDE, discreteBCs, load, stiffness, damping, spans, dimDM},
femdata = state["FiniteElementData"];
initBCs = femdata["BoundaryConditionData"];
methodData = femdata["FEMMethodData"];
initCoeffs = femdata["PDECoefficientData"];
vd = methodData["VariableData"];
sd = NDSolveSolutionData[{"Space" -> ToNumericalRegion[mesh]}];
(*Discretize the PDE and the boundary conditions:*)
discretePDE = DiscretizePDE[initCoeffs, methodData, sd];
discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd];
(*Extract the system matrices:*)
stiffness = discretePDE["StiffnessMatrix"];
damping = discretePDE["DampingMatrix"];
spans =
Span @@@
Transpose[{Most[# + 1], Rest[#]} &[methodData["IncidentOffsets"]]];
dimDM = Dimensions[discreteBCs["DirichletMatrix"]];
{load, stiffness, damping, methodData, spans, dimDM}
]

ClearAll[solveEigensystem];
solveEigensystem::usage =
"solveEigensystem[n,mesh,data] where n is the number of \
eigenvalues, and data is the output from finiteElementData";
solveEigensystem[n_,
mesh_, {load_, stiffness_, damping_, methodData_, spans_, dimDM_}] :=
Module[{nDiri, numEigen, mm, res, eigenValues, eigenVectors,
eigenVectorsIF},
nDiri = If[Length[#] > 0, First[#], 0] &[dimDM];
numEigen = n + nDiri;
mm = LinearSolve[damping, stiffness];
res = Eigensystem[mm, -numEigen];
res = Reverse /@ res;
eigenValues = res[[1, nDiri + 1 ;; Abs[numEigen]]];
eigenVectors = res[[2, nDiri + 1 ;; Abs[numEigen]]];
eigenVectors =
Transpose[
DeveloperToPackedArray[eigenVectors[[All, #]] & /@ spans], {2, 1,
3}];
eigenVectorsIF = Table[{}, {n}, {Length[spans]}];
Do[eigenVectorsIF[[i, j]] =
NDSolveFEMElementMeshInterpolation[{mesh},
eigenVectors[[i, j]]], {i, n}, {j, Length[spans]}];
{eigenValues, eigenVectorsIF}
]


Here is the beam

(* make mesh *)
\[CapitalOmega] = ImplicitRegion[True, {x, y}];
mesh = ToElementMesh[\[CapitalOmega], {{0, 5}, {0, 0.25}},
"MaxCellMeasure" -> 0.005];
mesh["Wireframe"] Starting with the static case -all is good

(* Set up boundary conditions and PDE*)
bcs = DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, x == 0];
pde2D = planeStress[u, v, x, y, 10^3, 33/100] == {0,
NeumannValue[0.1, x == 5]};
{state} =
NDSolveProcessEquations[{pde2D, bcs}, {u, v}, {x, y} \[Element]
mesh,
Method -> {"PDEDiscretization" -> {"FiniteElement"}}];


Extract the data, solve and plot.

{load, stiffness, damping, methodData, spans, dimDM} =
finiteElementData[state];
uif = ElementMeshInterpolation[{mesh}, solution[[spans[]]]];
vif = ElementMeshInterpolation[{mesh}, solution[[spans[]]]];
dmesh = ElementMeshDeformation[mesh, {uif, vif}];
Show[{
mesh["Wireframe"],
dmesh["Wireframe"[
"ElementMeshDirective" -> Directive[EdgeForm[Red], FaceForm[]]]]}] Next for vibration of the beam with no supports

(* Set up boundary conditions and PDE *)
bcs = Sequence[];
pde2D = {D[u[t, x, y], t], D[v[t, x, y], t]} +
planeStress[u, v, t, x, y, 10^3, 33/100] == {0, 0};
{state} =
NDSolveProcessEquations[{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"}}}];

fed = finiteElementData[state];
{eigenValues, eigenVectorsIF} = solveEigensystem[8, mesh, fed];


The theoretical values can be calculated assuming a slender beam. The first three are zero corresponding to rigid body modes i.e. moving horizontally, vertically and rocking with no flexure. The next four are transverse flexural vibration. The last value corresponds to longitudinal vibration.

theory = {0, 0, 0, 22/L^2 Sqrt[(Y d^2)/(12 1)],
61.7/L^2 Sqrt[(Y d^2)/(12 1)], 121/L^2 Sqrt[(Y d^2)/(12 1)],
200/L^2 Sqrt[(Y d^2)/(12 1)], \[Pi]/L Sqrt[Y/1.]} /. {Y -> 10^3,
d -> 0.25, L -> 5};
TableForm[Transpose[{Sqrt[Abs[eigenValues]], theory}],
TableHeadings -> {Automatic, { "Calculated", "Theory"}}] The eigenvectors plot as

Column[Table[uif = eigenVectorsIF[[n, 1]];
vif = eigenVectorsIF[[n, 2]];
dmesh =
ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 10];
Show[{
mesh["Wireframe"],
dmesh[
"Wireframe"[
"ElementMeshDirective" ->
Directive[EdgeForm[Red], FaceForm[]]]]}],
{n, 8}]] The first three mode shapes need normalizing to give the three rigid body modes corresponding to each coordinate. Now things go wrong - I put in the boundary conditions

bcs = DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, x == 0];
pde2D = {D[u[t, x, y], t], D[v[t, x, y], t]} +
planeStress[u, v, t, x, y, 10^3, 33/100] == {0, 0};
{state} =
NDSolveProcessEquations[{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"}}}];
fed = finiteElementData[state];
{eigenValues, eigenVectorsIF} = solveEigensystem[5, mesh, fed];


The theoretical solutions are as follows (note that the fourth solution is for longitudinal vibration)

 theory = {3.52 Sqrt[(Y d^2)/(12 L^4)], 22 Sqrt[(Y d^2)/(12 L^4)],
61.7 Sqrt[(Y d^2)/(12 L^4)], \[Pi] /2 Sqrt[Y/L^2],
121 Sqrt[(Y d^2)/(12 L^4)]} /. {Y -> 10^3, d -> 0.25, L -> 5.};
TableForm[Transpose[{Sqrt[Abs[eigenValues]], theory}],
TableHeadings -> {Automatic, { "Calculated", "Theory"}}] Column[Table[uif = eigenVectorsIF[[n, 1]];
vif = eigenVectorsIF[[n, 2]];
dmesh =
ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 10];
Show[{
mesh["Wireframe"],
dmesh[
"Wireframe"[
"ElementMeshDirective" ->
Directive[EdgeForm[Red], FaceForm[]]]]}],
{n, 5}]] The first mode frequency and mode shape are wrong. The solution is badly distorted at the left hand end. If we try with a mesh that is fine at the left hand end we get

Lx = 5; Ly = 0.25;
inc = 0.01;
pts = Flatten[
Join[{{{0, 0}, {Lx, 0}, {Lx, Ly}, {0, Ly}},
Table[{0, y}, {y, Ly - inc, inc, -inc}]}], 1];
boundarymesh = ToBoundaryMesh[Polygon[pts]];
mesh = ToElementMesh[boundarymesh, "MaxCellMeasure" -> 0.005];
mesh["Wireframe"] (* Set up boundary conditions and PDE *)
bcs = DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, x == 0];
pde2D = {D[u[t, x, y], t], D[v[t, x, y], t]} +
planeStress[u, v, t, x, y, 10^3, 33/100] == {0, 0};
{state} =
NDSolveProcessEquations[{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"}}}];
fed = finiteElementData[state];
{eigenValues, eigenVectorsIF} = solveEigensystem[5, mesh, fed];
TableForm[Transpose[{Sqrt[Abs[eigenValues]], theory}],
TableHeadings -> {Automatic, { "Calculated", "Theory"}}] Column[Table[uif = eigenVectorsIF[[n, 1]];
vif = eigenVectorsIF[[n, 2]];
dmesh =
ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 10];
Show[{
mesh["Wireframe"],
dmesh[
"Wireframe"[
"ElementMeshDirective" ->
Directive[EdgeForm[Red], FaceForm[]]]]}],
{n, 5}]] Again the first mode is wrong. It is clear that the left hand end elements are badly distorted. What is going wrong? Thanks

I am going to show you where the code that you refer to goes south and how to fix that. However, I'd like to stress that using NDEigensystem does this correctly; and it does more that is not considered here and the other post. A second point I'd like to make is that one can compare Euler beams with what is done here, but one should be aware that the theories are different and a discrepancy in the result is expected. So, it's OK if the results do not match 100%. Let's look at NDEigensystem first.

Define a plane stress operator:

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


{vals, funs} =
NDEigensystem[{ps}, {u, v}, {x, y} \[Element]
Rectangle[{0, 0}, {5, 0.25}], 8];


This gives:

theory = {0, 0, 0, 22/L^2 Sqrt[(Y d^2)/(12 1)],
61.7/L^2 Sqrt[(Y d^2)/(12 1)], 121/L^2 Sqrt[(Y d^2)/(12 1)],
200/L^2 Sqrt[(Y d^2)/(12 1)], \[Pi]/L Sqrt[Y/1.]} /. {Y -> 10^3,
d -> 0.25, L -> 5};
TableForm[Transpose[{Sqrt[Abs[vals]], theory}], The second case is:

bcs = DirichletCondition[{u[x, y] == 0, v[x, y] == 0}, x == 0];
{vals, funs} =
NDEigensystem[{ps, bcs}, {u, v}, {x, y} \[Element]
Rectangle[{0, 0}, {5, 0.25}], 5];


Checking:

theory = {3.52 Sqrt[(Y d^2)/(12 L^4)], 22 Sqrt[(Y d^2)/(12 L^4)],
61.7 Sqrt[(Y d^2)/(12 L^4)], \[Pi]/2 Sqrt[Y/L^2],
121 Sqrt[(Y d^2)/(12 L^4)]} /. {Y -> 10^3, d -> 0.25, L -> 5.};
TableForm[Transpose[{Sqrt[Abs[vals]], theory}], Visualizing:

Needs["NDSolveFEM"]
mesh = funs[[1, 1]]["ElementMesh"];
Column[Table[uif = funs[[n, 1]];
vif = funs[[n, 2]];
dmesh =
ElementMeshDeformation[mesh, {uif, vif}, "ScalingFactor" -> 0.1];
Show[{mesh["Wireframe"],
dmesh["Wireframe"[
"ElementMeshDirective" ->
Directive[EdgeForm[Red], FaceForm[]]]]}], {n, 5}]] What went wrong in the code above. In essence you need to remove the DirichletCondition from the system of equations. You'd notice this if you use Eigensystem with the method Arnoldi option. This gives you a hint with an error message saying that the matrices are not symmetric. Here is the code modified to do that:

    ClearAll[finiteElementData];
finiteElementData::usate =
"finiteElementData[start] Extract the finite element data";
finiteElementData[state_] :=
Module[{femdata, initBCs, methodData, initCoeffs, vd, sd,
discretePDE, discreteBCs, load, stiffness, damping, spans, dimDM},
femdata = state["FiniteElementData"];
initBCs = femdata["BoundaryConditionData"];
methodData = femdata["FEMMethodData"];
initCoeffs = femdata["PDECoefficientData"];
vd = methodData["VariableData"];
sd = NDSolveSolutionData[{"Space" -> ToNumericalRegion[mesh]}];
(*Discretize the PDE and the boundary conditions:*)

discretePDE = DiscretizePDE[initCoeffs, methodData, sd];
discreteBCs =
DiscretizeBoundaryConditions[initBCs, methodData, sd];
(*Extract the system matrices:*)

stiffness = discretePDE["StiffnessMatrix"];
damping = discretePDE["DampingMatrix"];
"ConstraintMethod" -> "Remove"];
spans =
Span @@@
Transpose[{Most[# + 1], Rest[#]} &[methodData["IncidentOffsets"]]];
dimDM = Dimensions[discreteBCs["DirichletMatrix"]];
{load, stiffness, damping, methodData, spans, dimDM, discreteBCs}
]

ClearAll[solveEigensystem];
solveEigensystem::usage =
"solveEigensystem[n,mesh,data] where n is the number of \
eigenvalues, and data is the output from finiteElementData";
solveEigensystem[n_,
mesh_, {load_, stiffness_, damping_, methodData_, spans_, dimDM_,
discreteBCs_}] :=
Module[{nDiri, numEigen, mm, res, eigenValues, eigenVectors,
eigenVectorsIF},
nDiri = If[Length[#] > 0, First[#], 0] &[dimDM];
(*numEigen=n+nDiri;*)
numEigen = n;
(*mm=LinearSolve[damping,stiffness];
res=Eigensystem[mm,-numEigen,Method->"Arnoldi"];*)

res = Eigensystem[{stiffness, damping}, -numEigen,
Method -> "Arnoldi"];
res = Reverse /@ res;
eigenValues = res[[1, 1 ;; Abs[numEigen]]];
eigenVectors = res[[2, 1 ;; Abs[numEigen]]];

eigenVectors =
NDSolveFEMDirichletValueReinsertion[#, discreteBCs] & /@
eigenVectors;

eigenVectors =
Transpose[
DeveloperToPackedArray[eigenVectors[[All, #]] & /@ spans], {2, 1,
3}];
eigenVectorsIF = Table[{}, {n}, {Length[spans]}];
Do[eigenVectorsIF[[i, j]] =
NDSolveFEMElementMeshInterpolation[{mesh},
eigenVectors[[i, j]]], {i, n}, {j, Length[spans]}];
{eigenValues, eigenVectorsIF}
]


What it does it uses the method option "ConstraintMethod" -> "Remove" the eliminate the Dirichlet boundary positions and the returns the data structure that contains the Dirichlet information. It's better the use Eigensystem with the Arnoldi method to get the speed in place of the LinearSolve variant. A bit of shifting the amount of the eigenvalues and vecors that need to be computed is done and then the Dirichlet conditions are then reinserted again.

• (1) Many thanks, once again, to user21 for the full and helpful response. My code now works and the alternative of NDEigensystem is possibly where I should be looking. I completely agree that my theory results (included for comparison above) is only approximate and I should have made this clear in my question. The Bernoulli beam theory does not include Poission ratio effects, shear deformation and rotary inertia. The finite element method should be better in this respect. The new option value "Remove" is what was needed. Finding these options is a challenge. Thanks again – Hugh Nov 22 '15 at 11:30
• (2) I will investigate NDEigensystem further. For the examples you give am I correct in noting that the assumption made is that the mass matrix is not only diagonal but the identity matrix? I will investigate circumstances where this is not the case and see if I can work things out. A look at the mass/damping matrix suggested that it was not diagonal. An alternative assumption would be that a uniform density of 1.0 is associated with the inertia term. – Hugh Nov 22 '15 at 12:01
• @Hugh, no the system matrices are not lumped. The system matrices from the code above are the same as what NDEigensystem uses. You can verify that by inspecting the system matrices of NDEigensystem via the documented "VectorNormalization" method option. Now, the fact that "Remove" (and a few other other constraint deployment methods) are not documented is that they do not yet work across the board. Once they do, they will be documented. Promise. – user21 Nov 23 '15 at 8:34