3
$\begingroup$

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["NDSolve`FEM`"];

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

\[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[[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[{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},
    gammaD,
    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.

Improve this question
$\endgroup$
3
  • 2
    $\begingroup$ 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? $\endgroup$
    – Henrik Schumacher
    Commented Apr 13, 2020 at 20:10
  • 1
    $\begingroup$ 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. $\endgroup$
    – Nasser
    Commented Apr 13, 2020 at 20:13
  • $\begingroup$ 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. $\endgroup$
    – Ponyboy Curtis
    Commented Apr 13, 2020 at 20:29

1 Answer 1

Reset to default
5
$\begingroup$

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["NDSolve`FEM`"];

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

\[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[[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[{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];
Improve this answer
$\endgroup$
5
  • $\begingroup$ 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? $\endgroup$
    – Ponyboy Curtis
    Commented Apr 14, 2020 at 12:20
  • $\begingroup$ @DiegoMagela Make sure the PDE you put into FullSimplify are all rational numbers. In other words no reals! $\endgroup$
    – user21
    Commented Apr 14, 2020 at 12:33
  • 1
    $\begingroup$ @DiegoMagela I dumped the code in the post. $\endgroup$
    – user21
    Commented Apr 14, 2020 at 12:50
  • $\begingroup$ 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? $\endgroup$
    – Ponyboy Curtis
    Commented Apr 14, 2020 at 13:07
  • 1
    $\begingroup$ @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. $\endgroup$
    – user21
    Commented Apr 14, 2020 at 13:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.