Implement Morley element for plate bending problems

[ahojukka5]

Let's implement Morley element for plate bending problems. I'm not sure if this has any practical usage in daily calculations but it's basis functions seems to be quite interesting and giving us some ideas what kind of cumbersome basis functions we should be able to automatically generate using FEMBasis.jl.

[ahojukka5]

Basis functions (in principle):

var('x,y')
p = a[1] + a[2]*x + a[3]*y + a[4]*x**2 + a[5]*x*y + a[6]*y**2
x1, y1 = 0, 0
x2, y2 = 1, 0
x3, y3 = 0, 1
eq1 = p.subs({x: x1, y: y1})
eq2 = p.subs({x: x2, y: y2})
eq3 = p.subs({x: x3, y: y3})
X1 = Matrix([x1, y1])
X2 = Matrix([x2, y2])
X3 = Matrix([x3, y3])
t1 = (X2-X1)/(X2-X1).norm()
t2 = (X3-X2)/(X3-X2).norm()
t3 = (X1-X3)/(X1-X3).norm()
Q = Matrix([[0, 1], [-1, 0]])
n1 = Q*t1
n2 = Q*t2
n3 = Q*t3
dp = Matrix([p.diff(x), p.diff(y)])
eq4 = dp.dot(n1).subs({x: Rational(1,2)*(x1+x2), y: Rational(1,2)*(y1+y2)})
eq5 = dp.dot(n2).subs({x: Rational(1,2)*(x2+x3), y: Rational(1,2)*(y2+y3)})
eq6 = dp.dot(n3).subs({x: Rational(1,2)*(x3+x1), y: Rational(1,2)*(y3+y1)})

coeffs = [a[1], a[2], a[3], a[4], a[5], a[6]]
eqs1 = [eq1-1, eq2, eq3, eq4, eq5, eq6]
eqs2 = [eq1, eq2-1, eq3, eq4, eq5, eq6]
eqs3 = [eq1, eq2, eq3-1, eq4, eq5, eq6]
eqs4 = [eq1, eq2, eq3, eq4-1, eq5, eq6]
eqs5 = [eq1, eq2, eq3, eq4, eq5-1, eq6]
eqs6 = [eq1, eq2, eq3, eq4, eq5, eq6-1]
sols = [solve(eq, coeffs) for eq in [eqs1, eqs2, eqs3, eqs4, eqs5, eqs6]]
N = [p.subs(sol) for sol in sols]

How to implement?