Ne(:,1)= (1-xsi)*(1-eta)/4;
        Ne(:,2)= (1+xsi)*(1-eta)/4;
        Ne(:,3)= (1+xsi)*(1+eta)/4;
        Ne(:,4)= (1-xsi)*(1+eta)/4;
        Ne2D =[Ne(:,1) 0    Ne(:,2) 0   Ne(:,3) 0   Ne(:,4) 0;
                0   Ne(:,1) 0    Ne(:,2) 0   Ne(:,3) 0   Ne(:,4)]



gradN1=[dN1dxsi;dN1deta];
        gradN2=[dN2dxsi;dN2deta];
        gradN3=[dN3dxsi;dN3deta];
        gradN4=[dN4dxsi;dN4deta];
        %
        %   Calculate gradients of all basis functions
        %   ==> Be-matrix
Be = [invJ*gradN1  invJ*gradN2     invJ*gradN3     invJ*gradN4];
Be2D = [Be(1,1)     0       Be(1,2)     0       Be(1,3)     0       Be(1,4)     0;
           0    Be(2,1)        0       Be(2,2)     0        Be(2,3)    0     Be(2,4);
        Be(1,1) Be(2,1) Be(1,2) Be(2,2) Be(1,3) Be(2,3) Be(1,4) Be(2,4)];