# Error using FEM package

I'm getting the following error:

NDSolveValue: The FEMStiffnessElements operator failed.

I looked for this error in FEM documentation and did not find anything. I'm using Mathematica 12.

It follows my code:

<< NumericalDifferentialEquationAnalysis;
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 (\[Pi]/180);
a = 1; b = 1; d = a Cos[alpha] + b Sin[alpha];
omega = Rectangle[{0, 0}, {a, b}];
mesh = ToElementMesh[omega, MaxCellMeasure -> 0.001];
u0 = 0.01;

angle1 = 10; 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][]; h = num*t; pos = Table[0, num + 1];
pos[] = -h/2;
For[i = 2, i <= num + 1, i++, pos[[i]] = pos[[i - 1]] + t];

\[Xi]A = {0, 0, 0, 0, 0};
\[Xi]B = {0, 0, 0, 0, 0};
\[Xi]D = {0, 0, 0, 0, 0};

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

T0 = angles[[i, 1]] ;
T1 = angles[[i, 2]] ;
func[s_] :=
Simplify@((2.0/d) (T1 - T0) Sqrt[(s - d/2)^2] + T0) (\[Pi]/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]];

\[Xi]a = {1, V1, V2, V3, V4} zA;
\[Xi]b = {1, V1, V2, V3, V4} zB;
\[Xi]d = {1, V1, V2, V3, V4} zD;

\[Xi]A = \[Xi]A + \[Xi]a;
\[Xi]B = \[Xi]B + \[Xi]b;
\[Xi]D = \[Xi]D + \[Xi]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.\[Xi]A; mB = (mU.\[Xi]B)/2; mD = (mU.\[Xi]D)/3;

A11[x_, y_] = mA[]; A12[x_, y_] = mA[]; A16[x_, y_] = mA[];
A22[x_, y_] = mA[]; A26[x_, y_] = mA[]; A66[x_, y_] = mA[];
D11[x_, y_] = mD[]; D12[x_, y_] = mD[]; D16[x_, y_] = mD[];
D22[x_, y_] = mD[]; D26[x_, y_] = mD[]; D66[x_, y_] = mD[];

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}]
};

{
DirichletCondition[{v[x, y] == u0, u[x, y] == 0}, y == 0],
DirichletCondition[{v[x, y] == -u0, u[x, y] == 0}, y == b]
};

omega = Rectangle[{0, 0}, {a, b}];

mesh = ToElementMesh[omega, MaxCellMeasure -> 0.001];

{U, V} =
NDSolveValue[{
PDEs == {0, 0},
DirichletCondition[u[x, y] == 0, x == a/2]
},
{u, v}, {x, y} \[Element] mesh
];


Basically I'm trying to solve a 2D elasticity (a plate under a prescribed displacement) problem above. I already solved it considering a 1D variation of the theta function e all went well. Now I need to solve considering a 2D variation of the theta function theta[x,y]. What I changed it was the function theta, but I'm getting this error.

Does anyone knows the reason of this error and how can I solve it?

UPDATE

It using alpha = 0 (\[Pi]/180) my code runs like a charm. But when I set 45 (\[Pi]/180) I get division-by-zero.

• Okay, when I evaluate your code in a fresh Mathematica kernel, I get already a couple of division-by-zero errors before NDSolveValue does anything. Maybe it would be a good idea to resolve these first? Apr 13 '20 at 20:10
• I also tried it and get same error message as Henrik. Division by zero. It is good to look at ALL the error messages in the console and not just at the last error at the end, because the earlier error messages most likely is what caused the last one you saw. Apr 13 '20 at 20:13
• Really, I'm getting division-by-zero as well, I thought could be due to FEM error. I can't understand that division-by-zero error in my code. Could you help me in that issue? I updated my code due to some typo right now and the error still happening. I noticed that using alpha = 0 (\[Pi]/180) my code runs properly. Apr 13 '20 at 20:29

This is a precision issue. Rationalize all your numbers and use

PDEs = FullSimplify[PDEs];


on the resulting PDE. Then it will work.

Code Dump

This works for 12.1:

<< NumericalDifferentialEquationAnalysis;
Needs["NDSolveFEM"];

G = Rationalize[6.894745 10^9]
E1 = Rationalize[26.25 G]; E2 = Rationalize[1.49 G]; G12 =
Rationalize[1.04 G]; nu12 = Rationalize[0.28]; nu21 = (E2*nu12)/E1;

t = Rationalize[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 = 45 (\[Pi]/180);
a = 1; b = 1; d = a Cos[alpha] + b Sin[alpha];
omega = Rectangle[{0, 0}, {a, b}];
mesh = ToElementMesh[omega, MaxCellMeasure -> 0.001];
u0 = Rationalize[0.01];
angle1 = 10; 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][]; h = num*t; pos = Table[0, num + 1];
pos[] = -h/2;
For[i = 2, i <= num + 1, i++, pos[[i]] = pos[[i - 1]] + t];

\[Xi]A = {0, 0, 0, 0, 0};
\[Xi]B = {0, 0, 0, 0, 0};
\[Xi]D = {0, 0, 0, 0, 0};
For[i = 1, i <= num, i++, T0 = angles[[i, 1]];
T1 = angles[[i, 2]];
func[s_] :=
Simplify@((2/d) (T1 - T0) Sqrt[(s - d/2)^2] + T0) (\[Pi]/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]];
\[Xi]a = {1, V1, V2, V3, V4} zA;
\[Xi]b = {1, V1, V2, V3, V4} zB;
\[Xi]d = {1, V1, V2, V3, V4} zD;
\[Xi]A = \[Xi]A + \[Xi]a;
\[Xi]B = \[Xi]B + \[Xi]b;
\[Xi]D = \[Xi]D + \[Xi]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.\[Xi]A; mB = (mU.\[Xi]B)/2; mD = (mU.\[Xi]D)/3;
A11[x_, y_] = mA[]; A12[x_, y_] = mA[]; A16[x_, y_] = mA[];
A22[x_, y_] = mA[]; A26[x_, y_] = mA[]; A66[x_, y_] = mA[];
D11[x_, y_] = mD[]; D12[x_, y_] = mD[]; D16[x_, y_] = mD[];
D22[x_, y_] = mD[]; D26[x_, y_] = mD[]; D66[x_, y_] = mD[];

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[{v[x, y] == u0, u[x, y] == 0}, y == 0],
DirichletCondition[{v[x, y] == -u0, u[x, y] == 0}, y == b]};
omega = Rectangle[{0, 0}, {a, b}];

mesh = ToElementMesh[omega, MaxCellMeasure -> 0.001];

PDEs = FullSimplify[PDEs];

{U, V} = NDSolveValue[{PDEs == {0, 0}, gammaD,
DirichletCondition[u[x, y] == 0, x == a/2]}, {u,
v}, {x, y} \[Element] mesh];

• I did it and it is not working yet. I'm getting division-by-zero. Did you change my code and run it? Only the FullSimplify[PDEs] took 30+ min in my PC. Is it that time of computing normal? Apr 14 '20 at 12:20
• @DiegoMagela Make sure the PDE you put into FullSimplify are all rational numbers. In other words no reals! Apr 14 '20 at 12:33
• @DiegoMagela I dumped the code in the post. Apr 14 '20 at 12:50
• It works like a charm. Thank you. For others angles, for instance 35 degree, it is taking too long to evaluate, while using your code for 45 degree is very fast. Do you know the reason? Apr 14 '20 at 13:07
• @DiegoMagela, in that case the equations can not be simplified as easy as in the 45 case. Try Simplify in stead of FullSimplify and see if that also helps. Apr 14 '20 at 13:17