I am trying to solve the following linear elastic problem:

The number pairs next to each node are the coordinates in meters. The charges are in Newtons.
The part of the codes that probably don't have errors:
ClearAll("Global`*")
(* Load the FEM package *)
<< NDSolve`FEM`
(* 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,
"MeshElements" -> {QuadElement(element)}
);
(* 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 *);
(* Loads *)
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_ {aa}) mathbf {I} big) $
$ boldsymbol { sigma} = displaystyle frac {E} {1- nu ^ 2} big ((1- nu) boldsymbol { varepsilon} + nu ( varepsilon_ {xx} + varepsilon_ {aa}) mathbf {I} big) $
The mesh, for future reference:

The part where I'm not so sure:
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, which converts nodal charges into distributed charges:
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 into pressure multiplied by 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 not valid for the tensor
Compile arguments :: argcompten: the comparison, LessEqual, is not valid
for tensor arguments. Compile :: argcompten: The comparison, LessEqual,
It is not valid for tensor arguments. General :: stop: additional departure from
Compile :: argcompten will be deleted during this calculation.
