# Error when trying to solve elastic PDEs using FEM: “Compile::argcompten: The comparison, LessEqual, is invalid for tensor arguments.”

I'm trying to solve the following linear elastic problem:

The pairs of numbers beside each node is the coordinates in meters. The loads are in Newtons.

The part of codes that probably have no mistakes:

ClearAll["Global*"]

(* Load the FEM package *)
<< NDSolveFEM

(* Define the mesh *)
node = {
{0., 0.}, {10., 2.}, {20., 0.},
{2., 6.}, {5., 8.}, {17., 6.}
};
element = {
{1, 2, 5, 4},
{2, 3, 6, 5}
};
mesh = ToElementMesh[
"Coordinates" -> node,
];

(* Constants *)
e = 2.*^11 (* Pa, Young's modulus *);
ν = 0.3    (* dimensionless, Poisson's ratio *);
ρ = 7.8*^3 (* kg m^-3, density *);
h = 0.1    (* m, thickness *);
g = 9.8    (* m s^-2, gravitational acceleration *);

f = {
{0.,     0.   },
{0.,     0.   },
{1.*^4,  0.   },
{0.,    -1.*^4},
{0.,    -1.*^4},
{0.,     0.   }
} (* N, nodal load *);
b = -ρ g (* N m^-3, body force *);

(* Strain (2x2) from displacement (1x2) *)
ε[u_] := Grad[u[x, y], {x, y}];

(* Strain from stress (2x2) *)
σ2ε[σ_] := 1/e ((1 + ν) σ - ν Total@Diagonal[σ] IdentityMatrix[2]);

(* Stress from displacement *)
σ[u_] := e/(1 - ν^2) ((1 - ν) ε[u] + ν Total@Diagonal[ε[u]] IdentityMatrix[2]);


The last three definitions are:

$$\boldsymbol{\varepsilon} = \nabla \boldsymbol{u}$$

$$\boldsymbol{\varepsilon} = \displaystyle\frac{1}{E}\big[(1+\nu)\boldsymbol{\sigma} - \nu(\sigma_{xx} + \sigma_{yy})\mathbf{I}\big]$$

$$\boldsymbol{\sigma} = \displaystyle\frac{E}{1-\nu^2}\big[(1-\nu)\boldsymbol{\varepsilon} + \nu(\varepsilon_{xx} + \varepsilon_{yy})\mathbf{I}\big]$$

The mesh, for future reference:

The part where I'm not so sure about:

The PDE operator, which is basically $$\nabla\cdot\boldsymbol{\sigma}$$ :

(* PDE operator *)
op = Div[σ[{ux[x, y], uy[x, y]}], {x, y}];
opx = Flatten[op][[1]] (* x component *);
opy = Flatten[op][[2]] (* y component *);


Neumann boundary condition, converting the nodal loads to distributed loads:

nbc := Module[
{
line = Line[node[[{#1, #2}]]],
l, sum1, sum2, n, σ
},

l = ArcLength[line];
sum1 = Total[node[[#1]]];
sum2 = Total[node[[#2]]];
n = RotationMatrix[90 Degree].((node[[#2]] - node[[#1]])/l);
σ = Total[
f[[{#1, #2}]] {x + y - sum1, -x - y + sum2} / (l h (sum2 - sum1))
]
IdentityMatrix[2];

{
NeumannValue[(n.σ2ε[σ])[[1]], {x, y} ∈ line] (* x component *),
NeumannValue[(n.σ2ε[σ])[[2]], {x, y} ∈ line] (* y component *)
}
] &;
ΓN = Total@{
nbc @@ {4, 5} (* Edge 4-5 *),
nbc @@ {6, 3} (* Edge 6-3 *)
};


This converts the nodal load to pressure times the shape function of that node, and then finds the normal component.

Dirichlet boundary condition:

ΓD =
{
DirichletCondition[{ux[x, y] == 0, uy[x, y] == 0}, {x, y} == node[[1]]],
DirichletCondition[uy[x, y] == 0, {x, y} == node[[2]]],
DirichletCondition[uy[x, y] == 0, {x, y} == node[[3]]]
};


Solving the PDEs

uHat = NDSolveValue[
{opx == ΓN[[1]], opy + b == ΓN[[2]], ΓD},
{ux[x, y], uy[x, y]},
{x, y} ∈ mesh
];
`

which is where the error occurs:

Compile::argcompten: The comparison, LessEqual, is invalid for tensor arguments. Compile::argcompten: The comparison, LessEqual, is invalid for tensor arguments. Compile::argcompten: The comparison, LessEqual, is invalid for tensor arguments. General::stop: Further output of Compile::argcompten will be suppressed during this calculation.

• Why are you using a two-element mesh? – Alex Trounev Nov 27 '19 at 0:47
• @AlexTrounev It's a course assignment that specifies this mesh. I wanted to try using the FEM package which lead to this problem here. I ended up hard coding the algorithms instead. – Yukai Qian Nov 28 '19 at 0:32