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["NDSolve`FEM`"];
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,
load, stiffness, damping},
(* 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,
0, -Y/(2*(1 + ν))}}.Inactive[Grad][
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 -
2*ν)*(1 + ν))), 0}}.Inactive[Grad][
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 -
2*ν)*(1 + ν)))}}.Inactive[Grad][
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} =
NDSolve`ProcessEquations[{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 = NDSolve`SolutionData[{"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*)
load = discretePDE["LoadVector"];
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*)
DeployBoundaryConditions[{load, stiffness, damping}, discreteBCs];
{load, stiffness, damping}
]
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}];
TableForm[Reverse@res[[1]], TableHeadings -> {Automatic, None}]
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.
