Formulation of finite element matrices for a 
2d elastic square element [-1,1]×[-1,1]

Clear["Global`*"]; Off[General::spell, General::spell1];

Shape functions and derivatives

(Subscript[ϕ, 1])  =  (1/4) (1 + x) (1 + 
              y);
(Subscript[ϕ, 2])  =  (1/4) (1 - x) (1 + 
              y);
(Subscript[ϕ, 3])  =  (1/4) (1 - x) (1 - 
              y);
(Subscript[ϕ, 4])  =  (1/4) (1 + x) (1 - 
              y);
dϕ1x  =  D[(Subscript[ϕ, 1]), x];
dϕ1y  =  D[(Subscript[ϕ, 1]), y];
dϕ2x  =  D[(Subscript[ϕ, 2]), x];
dϕ2y  =  D[(Subscript[ϕ, 2]), y];
dϕ3x  =  D[(Subscript[ϕ, 3]), x];
dϕ3y  =  D[(Subscript[ϕ, 3]), y];
dϕ4x  =  D[(Subscript[ϕ, 4]), x];
dϕ4y  =  D[(Subscript[ϕ, 4]), y];

Element displacement interpolation matrix

H = (({{(Subscript[ϕ, 1]), 0, (Subscript[ϕ, 2]), 0, (Subscript[ϕ, 3]), 0, (Subscript[ϕ, 4]), 0}, {0, (Subscript[ϕ, 1]), 0, (Subscript[ϕ, 2]), 0, (Subscript[ϕ, 3]), 0, (Subscript[ϕ, 4])}}));
MatrixForm[
      Transpose[H]]

Element nodal displacement vector

Clear[U]
UU  =  
        Table[(Subscript[U, i, a]), {a, 4}, {i, 2}] // 
          Flatten;
MatrixForm[UU]

Element strain-operator matrix

B = (({{dϕ1x, 0, dϕ2x, 0, dϕ3x, 0, dϕ4x, 0}, {0, dϕ1y, 0, dϕ2y, 0, dϕ3y, 0, dϕ4y}, {dϕ1y, dϕ1x, dϕ2y, dϕ2x, dϕ3y, dϕ3x, dϕ4y, dϕ4x}}));
MatrixForm[
      B]

Compute a few strain fields from given nodal displacements

Nodal coordinates

(Subscript[XY, 1]) =  {1, 1};
(Subscript[XY, 2]) =  {-1, 1};
(Subscript[XY, 3]) =  {-1, -1};
(Subscript[XY, 4]) =  {1, -1};
XYn  =  Table[(Subscript[XY, i]), {i, 4}]

XYp  =  
      Append[XYn, XYn[[1]]]

Nodal displacements

a = 0.25;
(Subscript[UV, 1]) =  a {1, 0};
(Subscript[UV, 2]) =  a {1, 0};
(Subscript[UV, 3]) =  a {2, 0};
(Subscript[UV, 4]) =  a {0, -3};
UVn  =  Table[(Subscript[UV, i]), {i, 4}]

UVp  =  
      Append[UVn, UVn[[1]]]

U  =  
      Flatten[UVn]
  

    XYd  =  XYn  +  UVn
  

    u  =  H . U // Simplify // Chop

ϵ  =  B . U // 
          Simplify // Chop

Show[
        Graphics[{Thickness[0.01], RGBColor[1, 0, 0], Line[XYp], 
            RGBColor[0, 0, 1], Line[XYp + UVp]}, 
          AspectRatio -> Automatic, Axes -> True]];

Elastic moduli (Plane strain)

CC = (EE/((1 - (ν^2))))(({{1, ν, 0}, {ν, 1, 0}, {0, 0, (((1 - ν))/2)}}));
MatrixForm[
      CC]

K  =  
        Integrate[Transpose[B] . CC . B, {x, -1, 1}, {y, -1, 1}] // 
          Simplify;
MatrixForm[K]