This is a form function generation code and result of quarter singular element mesh. Is it correct?
parm = LinearSolve({f(1), f(2), f(3), f(4), f(5), f(6), f(7),
f(8)}, {u1, u2, u3, u4, u5, u6, u7, u8});
f(n_) = {1, Subscript(x, n ), Subscript(y, n ), Subscript(x, n )^2,
Subscript(x, n )*Subscript(y, n ), Subscript(y, n )^2,
Subscript(y, n )*Subscript(x, n )^2,
Subscript(x, n )*Subscript(y, n )^2};
Subscript(x, 1 ) = 1; Subscript(y, 1 ) = 1;
Subscript(x, 2) = -1; Subscript(y, 2 ) = 1;
Subscript(x, 3 ) = -1; Subscript(y, 3 ) = -1;
Subscript(x, 4 ) = 1; Subscript(y, 4 ) = -1;
Subscript(x, 5 ) = 0; Subscript(y, 5 ) = 1;
Subscript(x, 6 ) = -1; Subscript(y, 6 ) = 0;
Subscript(x, 7 ) = 0.5; Subscript(y, 7 ) = -1;
Subscript(x, 8) = 1; Subscript(y, 8 ) = -0.5;
g = Coefficient(parm, #) &;
formfuction = {g(u1), g(u2), g(u3), g(u4), g(u5), g(u6), g(u7),
g(u8)} . {1, x, y, x^2, x*y, y^2, y*x^2, x*y^2};
shapformfuction = formfuction /. Association({x -> r, y -> s});
Print("N(1,2,3,4,5,6,7,8)=", shapformfuction // MatrixForm)