# Finite Element Method fails to give symmetric stiffness matrix for small element size

In a previous post here a method for getting the matrices from a finite element analysis were developed. I am now using this method and it works for large elements but not for small elements. My overall objective is to get the matrices so that I can do an eigenvalue calculation with them.

The start is a large module that processes the mesh to make the matrices.

Needs["NDSolveFEM"];

ClearAll[MakeMatrices];
MakeMatrices[mesh_, modE_, pr_] :=
Module[{stressOperator, t, u, v, w, x, y, z, Y, ν, pde3D, bcs,
state, femdata, initBCs, methodData, initCoeffs, vd, nr, sd,
discretePDE, discreteBCs,

(* Stress operator for 3D time dependant elasticity *)
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,
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 -
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 -
w[t, x, y, z], {x, y, z}], {x, y, z}]};

(* Set the modulus of elasticity and Poission ratio *)
Y = modE;
ν = pr;

(* 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.*)

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} =
NDSolveProcessEquations[{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 = NDSolveSolutionData[{"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*)
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*)
]


I now make a mesh. Note that I have a "MaxCellMeasure" of 0.005

L = 1;
a = 0.5;
tk = 0.3;
coords = {{0, 0, 0}, {L, 0, 0}, {L, a, 0}, {0, a, 0}, {0, 0, tk}, {L,
0, tk}, {L, a, tk}, {0, a, tk}};
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, 7, 8}}]}];
Show[bmesh["Wireframe"],
bmesh["Wireframe"["MeshElement" -> "PointElements",
"MeshElementIDStyle" -> Red]]]
mesh = ToElementMesh[bmesh, "MeshOrder" -> 1,
"MaxCellMeasure" -> 0.005];
mesh["Wireframe"]


I now use my module to get the matrices

{load, stiffness, damping} = MakeMatrices[mesh, 100, 1/3];


The damping matrix is nice and symmetric but the stiffness matrix is not.

damping == Transpose[damping]
stiffness == Transpose[stiffness]

(*
True

False
*)


However if I repeat with a larger "MaxCellMeasure" 0f 0.01

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


I get symmetric matrices and I can go on to calculate my eigenvalues

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


So what is happening? How do I make a fine mesh? Thanks for any help.

Edit I have also noted that if you run the problem in Help concerning the Swingng Beam - Transient coupled PDEs which may be found here then the stiffness matrix is also not symmetrix. Investigations suggested that the loss of symmetry is more than just due to numerical noise. This is a two dimensional problem so the issues seems widespread.

Edit 2 Thanks to user21 we have a partial answer below and it seems that some lack of symmetry is not too much of a problem and may just be a rounding issue. However it can be a significant problem sometimes. Taking the geometry from your previous problem

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})}]}];
mesh = ToElementMesh[bmesh, "MeshOrder" -> 1, "MaxCellMeasure" -> 10];

{load, stiffness, damping} = MakeMatrices[mesh, 100, 1/3];

SymmetricMatrixQ[stiffness]

(* False  *)


The eigen solver now fails

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

(*
Eigensystem::nosymh: The input matrix SparseArray[Automatic,{3207,3207},0.,<<1>>] should be real symmetric or complex Hermitian.
*)


This is actually where I met the problem first. I have been using this example as a test case. I then tried to make a minimum working example which you have analysed in your answer. From the example immediately above the loss of symmetry can be so significant that it stops the eigen calculation. So is the question now along the lines of how much loss of symmetry can be tolerated? Also this example did work in previous versions so something has happened. Is there a workaround to this problem? Thanks once again for your help.

This is a numerical precision thing with Equal.

mesh = ToElementMesh[bmesh, "MeshOrder" -> 1,
"MaxCellMeasure" -> 0.005];
{load, stiffness, damping} = MakeMatrices[mesh, 100, 1/3];
damping == Transpose[damping]
stiffness == Transpose[stiffness]

True
False


When we look at the actual values they do not seem too large:

MinMax[SparseArray[stiffness - Transpose[stiffness]]["NonzeroValues"]]

{-7.10543*10^-15, 7.10543*10^-15}


Also the eigenvalues seem fine:

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

res2 = Eigensystem[{stiffness, damping}, -10, Method -> {"Arnoldi"}];

res1[[1, ;; -7]] - res2[[1, ;; -7]]
{-5.330321073415689*^-9, -4.699813871411607*^-10, \
-4.128878572373651*^-9, -8.99490260053426*^-10}


The remaining 6 eigenvalues are rigid body modes and thus close to zero.

Note also that

SymmetricMatrixQ[stiffness]
True


Update

For the second case you can either use the Arnoldi method (FEAST takes for ever for these type of problems):

AbsoluteTiming[
res = Eigensystem[{stiffness, damping}, -10, Method -> {"Arnoldi"}];]


which will give a warning about the damping matrix not being positive definite but should be OK. Or, if you know that you stiffness matrix is symmetric, you can test the symmetry property to a certain Tolerance like so

SymmetricMatrixQ[stiffness, Tolerance -> 10^(-10)]
True


In case that test gives true I'd then just for the stiffness matrix to be symmetric:

stiffness = (stiffness + Transpose[stiffness])/2;

• Thank you for looking at this. For this example the loss of symmetry does not stop the eigen calculation. I have added an example to the original question (see Edit 2) where the loss of symmetry does stop the eigen calculation. Please could you look at that? How much symmetry can be lost and is there a work around? Thanks
– Hugh
Jun 6, 2017 at 9:37