# To obtain eigenvalues and eigenvectors from FEM solution

I have a model which solve a 2D coupled elasticity problem. Now I'm trying to get the eigenvalues and eigenvector from that solution. But, following the very first step to get stiffness and load matrices, I faced up an error. Could someone help me?

My code:

Needs["NDSolveFEM"];

G = 6.894745 10^9;
E1 = 26.25 G; E2 = 1.49 G; G12 = 1.04 G; nu12 = 0.28; nu21 = (E2*nu12)/E1; t = 0.0050 .0254;
Son = {{1/E1, -nu12/E1, 0}, {-nu21/E2, 1/E2, 0}, {0, 0, 1/G12}};
Qon = Inverse[Son];
Q11 = Qon[[1, 1]]; Q12 = Qon[[1, 2]]; Q22 = Qon[[2, 2]]; Q66 = Qon[[3, 3]];
U1 = (3 Q11 + 3 Q22 + 2 Q12 + 4 Q66)/8; U2 = (Q11 - Q22)/2;
U3 = (Q11 + Q22 - 2 Q12 - 4 Q66)/8; U4 = (Q11 + Q22 + 6 Q12 - 4 Q66)/8;
U5 = (Q11 + Q22 - 2 Q12 + 4 Q66)/8;

alpha = 0.0; a = 1; b = 1; d = a Cos[alpha] + b Sin[alpha];
u0 = Rationalize[0.001];
angle1 = 30; angle0 = 0;
angles = {{angle0, angle1}, {-angle0, -angle1}, {angle0, angle1}, {-angle0, -angle1},
{angle0, angle1}, {-angle0, -angle1}, {-angle0, -angle1}, {angle0, angle1},
{-angle0, -angle1}, {angle0, angle1}, {-angle0, -angle1}, {angle0, angle1}};

num = Dimensions[angles][[1]]; h = num*t; pos = Table[0, num + 1];
pos[[1]] = -h/2; For[i = 2, i <= num + 1, i++, pos[[i]] = pos[[i - 1]] + t];
ξA = {0, 0, 0, 0, 0}; ξB = {0, 0, 0, 0, 0}; ξD = {0, 0, 0, 0, 0};

For[i = 1, i <= num, i++,

T0 = angles[[i, 1]]; T1 = angles[[i, 2]];

func[s_] := ((2/d) (T1 - T0) s + T0) (π/180);
theta[x_, y_] := alpha + func[x Cos[alpha] + y Sin[alpha]];

zA = pos[[i + 1]] - pos[[i]]; zB = pos[[i + 1]]^2 - pos[[i]]^2;
zD = pos[[i + 1]]^3 - pos[[i]]^3;
V1 = Cos[2 theta[x, y]]; V2 = Sin[2 theta[x, y]];
V3 = Cos[4 theta[x, y]]; V4 = Sin[4 theta[x, y]];

ξa = {1, V1, V2, V3, V4} zA; ξb = {1, V1, V2, V3, V4} zB;
ξd = {1, V1, V2, V3, V4} zD;

ξA = ξA + ξa; ξB = ξB + ξb; ξD = ξD + ξd;];

mU = {{U1, U2, 0, U3, 0}, {U4, 0, 0, -U3, 0}, {U1, -U2, 0, U3, 0}, {0,0, U2/2, 0, U3},
{0, 0, U2/2, 0, -U3}, {U5, 0, 0, -U3, 0}};

mA = mU.ξA; mB = (mU.ξB)/2; mD = (mU.ξD)/3;

A11[x_, y_] = mA[[1]]; A12[x_, y_] = mA[[2]]; A16[x_, y_] = mA[[4]];
A22[x_, y_] = mA[[3]]; A26[x_, y_] = mA[[5]]; A66[x_, y_] = mA[[6]];

D11[x_, y_] = mD[[1]]; D12[x_, y_] = mD[[2]]; D16[x_, y_] = mD[[4]];
D22[x_, y_] = mD[[3]]; D26[x_, y_] = mD[[5]]; D66[x_, y_] = mD[[6]];

Nx[x_, y_] = A11[x, y] D[u[x, y], {x, 1}] + A12[x, y] D[v[x, y], {y, 1}] +
A16[x, y] (D[u[x, y], {y, 1}] + D[v[x, y], {x, 1}]);
Ny[x_, y_] = A12[x, y] D[u[x, y], {x, 1}] + A22[x, y] D[v[x, y], {y, 1}] +
A26[x, y] (D[u[x, y], {y, 1}] + D[v[x, y], {x, 1}]);
Nxy[x_, y_] = A16[x, y] D[u[x, y], {x, 1}] + A26[x, y] D[v[x, y], {y, 1}] +
A66[x, y] (D[u[x, y], {y, 1}] + D[v[x, y], {x, 1}]);

PDEs = {D[Nx[x, y], {x, 1}] + D[Nxy[x, y], {y, 1}], D[Ny[x, y], {y, 1}] + D[Nxy[x, y], {x,
1}]};
gammaD = {DirichletCondition[{u[x, y] == -u0}, x == a/2],
DirichletCondition[{u[x, y] == 0}, x == 0],
DirichletCondition[{v[x, y] == 0}, y == 0],
DirichletCondition[{v[x, y] == 0}, y == b/2]};
omega = Rectangle[{0, 0}, {a/2, b/2}];
mesh = ToElementMesh[omega, MaxCellMeasure -> 0.001];
{U, V} = NDSolveValue[{PDEs == {0, 0}, gammaD}, {u, v}, {x, y} ∈ mesh];

{state} = NDSolveProcessEquations[{PDEs == {0, 0}, gammaD}, {u, v}, {x, y} ∈ mesh,
Method -> {"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement"}}}];


I'm getting the following error:

NDSolveProcessEquations::ivone: Boundary values may only be specified for one independent variable. Initial values may only be specified at one value of the other independent variable.

Set::shape: Lists {state} and NDSolveProcessEquations[<<1>>] are not the same shape.

• Have a look at the documentation: ref/message/NDSolve/ivone – user21 Apr 21 '20 at 6:58
• @user21, thank you for answering me! I already looked at the documentation. What I don't understand is that using NDSolveValue works like a charm, but when try to use NDSolveProcessEquations, I get that error. I really don't know how to change my BCs to solve this problem, for what I understood from documentation that I should do. Can you, please, help me to understand the reason of that error? – Professor P. Cosmo Klunk Apr 21 '20 at 13:16
• @user21, do you know another way to get the FEM matrices to evaluate eigenvalues and eigenvectors? – Professor P. Cosmo Klunk Apr 21 '20 at 13:21
• Try NDEigensystem, or this and this. – user21 Apr 22 '20 at 11:10