Finite Element Code for Stokes Problem using P2P1 Quad Elements

Question

I am trying to write some MATLAB code to solve fluid flow problems using the finite element method. I am starting with the relatively-simple steady Stokes problem. I am using the quadrilateral P2P1 element, and following Reddy and Gartling 2010. As far as I understand, the viscous matrix (K(1:18,1:18) in my code below) should be symmetric positive definite for the linear Stokes problem, but mine is symmetric but not positive definite! I checked every single line and cannot still figure out why it's not working!

Could anyone please take a look at the code and please tell me if anything is wrong!

I would definitely appreciate any additional advice/resources, thank you.

clc; clear;

syms r s mu rho x1 x2 x3 x4 x5 x6 x7 x8 x9 y1 y2 y3 y4 y5 y6 y7 y8 y9

% Velocity Shape functions (Reddy and Gartling P65):
% h > ψ, r > ξ, s > η, H > Ψ:
h1 = ( (1-r)*(1-s)*(-r-s-1)+(1-r^2)*(1-s^2) )/4;
h2 = ( (1+r)*(1-s)*(+r-s-1)+(1-r^2)*(1-s^2) )/4;
h3 = ( (1+r)*(1+s)*(+r+s-1)+(1-r^2)*(1-s^2) )/4;
h4 = ( (1-r)*(1+s)*(-r+s-1)+(1-r^2)*(1-s^2) )/4;
h5 = ( 2*(1-r^2)*(1-s^1)-(1-r^2)*(1-s^2) )/4;
h6 = ( 2*(1+r^1)*(1-s^2)-(1-r^2)*(1-s^2) )/4;
h7 = ( 2*(1-r^2)*(1+s^1)-(1-r^2)*(1-s^2) )/4;
h8 = ( 2*(1-r^1)*(1-s^2)-(1-r^2)*(1-s^2) )/4;
h9 = ( 4*(1-r^2)*(1-s^2) )/4;
H =  [h1; h2; h3; h4; h5; h6; h7; h8; h9];

Hr = diff(H,r); Hs = diff(H,s);

% Coordinate Interpolation (Local to Natural Coor.) (Reddy and Gartling P165)
x_n = [x1; x2; x3; x4; x5; x6; x7; x8; x9];
y_n = [y1; y2; y3; y4; y5; y6; y7; y8; y9];

% The Jacobian operator and its determinant (Reddy and Gartling P69)
J = [Hr.'; Hs.']*[x_n y_n];
iJ = J^-1; dJ = det(J);
% Derivatives of H Shape Functions:
Hxy = iJ*[Hr.'; Hs.']; Hx = transpose(Hxy(1,:)); Hy = transpose(Hxy(2,:));

% Pressure Basis Functions (Reddy and Gartling P65):
HP = [(1-r)*(1-s)/(4); (1+r)*(1-s)/(4); (1+r)*(1+s)/(4); (1-r)*(1+s)/(4)];

% Building the FE Equations (Reddy and Gartling P166):
K11 = mu*(Hx*Hx.')*dJ;
K22 = mu*(Hy*Hy.')*dJ;
K12 = mu*Hx*Hy.'*dJ;
K21 = mu*Hy*Hx.'*dJ;
Q1  = Hx*HP.'*dJ;
Q2  = Hy*HP.'*dJ;
Zp  = zeros(max(size(HP)), max(size(HP)));

K = [2*K11+K22,  K12,       -Q1; ...
     K21,        K11+2*K22, -Q2; ...
    -Q1.',      -Q2.',       Zp];

Answer 1

I realized the shape functions h5, h6, h7, and h8 in the textbook were wrong. They are missing "2" in the second term. They should be as follows:

h5 = (2*(1-r^2)*(1-s^1)-2*(1-r^2)*(1-s^2) )/4;
h6 = (2*(1+r^1)*(1-s^2)-2*(1-r^2)*(1-s^2) )/4;
h7 = (2*(1-r^2)*(1+s^1)-2*(1-r^2)*(1-s^2) )/4;
h8 = (2*(1-r^1)*(1-s^2)-2*(1-r^2)*(1-s^2) )/4;

It makes me sad that such a highly praised book has such a silly error.