# Code: Lid driven cavity using pressure free velocity form

(Difference between revisions)
 Revision as of 03:06, 29 June 2011 (view source) (Lid-driven cavity code, description to follow.)← Older edit Revision as of 04:19, 29 June 2011 (view source)Newer edit → Line 418: Line 418: num2str(etime(clock,ETstart)/60) ' min']); % Elapsed time from start <<< num2str(etime(clock,ETstart)/60) ' min']); % Elapsed time from start <<< -
-                                                                                                              function [Dm,RowNdx,ColNdx]=DMatW(Xe,Elcon,nn2nft)
-                                                                                                              %DMATW - Returns the affine-mapped element diffusion matrix for the simple cubic Hermite
-                                                                                                              %   basis functions on 4-node straight-sided quadrilateral elements with 3 DOF
-                                                                                                              %   per node using Gauss quadrature on the reference square and row/col indices.

-                                                                                                              %                                                                                                                +                                                             [[PFV diffusion matrix | Diffusion matrix for pressure-free velocity method]]
-                                                                                                              % Cubic-complete, fully-conforming, divergence-free, Hermite basis                                               +
-                                                                                                              %   functions on 4-node rectangular elements with 3 DOF per node using                                           +
-                                                                                                              %   Gauss quadrature on the 2x2 reference square.                                                                +
-                                                                                                              % The assumed nodal numbering starts with 1 at the lower left corner                                             +
-                                                                                                              %   of the element and proceeds counter-clockwise around the element.                                            +
-                                                                                                              % Uses second derivatives of stream function.                                                                    +
-                                                                                                              %                                                                                                                +
-                                                                                                              % Usage:                                                                                                         +
-                                                                                                              %   [Dm,Rndx,Cndx] = DMatW(Xe,Elcon,nn2nft)                                                                      +
-                                                                                                              %   Xe(1,:) -  x-coordinates of corner nodes of element.                                                         +
-                                                                                                              %   Xe(2,:) -  y-coordinates of corner nodes of element.                                                         +
-                                                                                                              %      and shape of the element. It is constant for affine elements.                                             +
-                                                                                                              %   Elcon  - connectivity matrix for this element.                                                               +
-                                                                                                              %   nn2nft - global number and type of DOF at each node                                                          +
-                                                                                                              %                                                                                                                +
-                                                                                                              % Jonas Holdeman, August 2007, revised June 2011                                                                 +

-                                                                                                              % Constants and fixed data                                                                                       +                                                             [[PFV convection matrix | Convection matrix for pressure-free velocity method]]
-                                                                                                              nd = 3;  nd4=4*nd;  ND=1:nd;   % nd = number of dofs per node,                                                   +
-                                                                                                              nn=[-1 -1; 1 -1; 1 1; -1 1];   % defines local nodal order                                                       +

-                                                                                                              % Define 4-point quadrature data once, on first call.                                                            +                                                             [[PFV get pressure | Consistent pressure for pressure-free velocity method]]
-                                                                                                              % Gaussian weights and absissas to integrate 7th degree polynomials exactly.                                     +
-                                                                                                              global GQ4;                                                                                                      +
-                                                                                                              if (isempty(GQ4))       % Define 4-point quadrature data once, on first call.                                    +
-                                                                                                              Aq=[-.861136311594053,-.339981043584856,.339981043584856, .861136311594053]; %Abs                                +
-                                                                                                              Hq=[ .347854845137454, .652145154862546,.652145154862546, .347854845137454]; %Wts                                +
-                                                                                                              GQ4.size=16; GQ4.xa=[Aq;Aq;Aq;Aq]; GQ4.ya=GQ4.xa';                                                               +
-                                                                                                              wt=[Hq;Hq;Hq;Hq]; GQ4.wt=wt.*wt';                                                                                +
-                                                                                                              end                                                                                                              +

-                                                                                                              xa=GQ4.xa; ya=GQ4.ya; wt=GQ4.wt; Nq=GQ4.size;                                                                    +                                                             [[PFV GMRES solver| GMRES solver with ILU preconditioning and Essential BC for pressure-free velocity method]]

-                                                                                                              % -----------------------------------------------                                                                +                                                             [[PFV mesh regrade | grade node spacing for pressure-free velocity method]]
-                                                                                                              global Zs3412d2;                                                                                                 +
-                                                                                                              if (isempty(Zs3412d2)|size(Zs3412d2,2)~=Nq)                                                                      +
-                                                                                                              % Evaluate and save/cache the set of shape functions at quadrature pts.                                          +
-                                                                                                              Zs3412d2=cell(4,Nq);                                                                                             +
-                                                                                                              for k=1:Nq                                                                                                       +
-                                                                                                              for m=1:4                                                                                                        +
-                                                                                                              Zs3412d2{m,k}=D3s(nn(m,:),xa(k),ya(k));                                                                          +
-                                                                                                              end                                                                                                              +
-                                                                                                              end                                                                                                              +
-                                                                                                              end  % if(isempty(*))                                                                                            +
-                                                                                                              % --------------- End fixed data ----------------                                                                +
-                                                                                                                                                                                                                               +
-                                                                                                              Ti=cell(4);                                                                                                      +
-                                                                                                                                                                                                                               +
-                                                                                                              for m=1:4                                                                                                        +
-                                                                                                              Jt=Xe*GBL(nn(:,:),nn(m,1),nn(m,2));   % transpose of Jacobian at node m                                          +
-                                                                                                              JtiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)];   % det(J)*inv(Jt)                                                    +
-                                                                                                              Ti{m}=blkdiag(1,JtiD);                                                                                           +
-                                                                                                              end                                                                                                              +
-                                                                                                                                                                                                                               +
-                                                                                                              % Move Jacobian evaluation inside k-loop for general convex quadrilateral.                                       +
-                                                                                                              % Jt=[x_q, x_r; y_q, y_r];                                                                                       +
-                                                                                                                                                                                                                               +
-                                                                                                              Dm=zeros(nd4,nd4); Sx=zeros(2,nd4); Sy=zeros(2,nd4);   % Pre-allocate arrays                                     +
-                                                                                                                                                                                                                               +
-                                                                                                              for k=1:Nq                                                                                                       +
-                                                                                                                                                                                                                               +
-                                                                                                              Jt=Xe*GBL(nn(:,:),xa(k),ya(k));       % transpose of Jacobian at (xa,ya)                                         +
-                                                                                                              Det=Jt(1,1)*Jt(2,2)-Jt(1,2)*Jt(2,1);  % Determinant of Jt & J                                                    +
-                                                                                                              TL=[Jt(2,2)^2, -2*Jt(2,1)*Jt(2,2), Jt(2,1)^2; ...                                                                +
-                                                                                                              -Jt(1,2)*Jt(2,2), Jt(1,1)*Jt(2,2)+Jt(2,1)*Jt(1,2), -Jt(1,1)*Jt(2,1); ...                                         +
-                                                                                                              Jt(1,2)^2, -2*Jt(1,1)*Jt(1,2),  Jt(1,1)^2]/Det^2;                                                                +
-                                                                                                                                                                                                                               +
-                                                                                                              % Initialize functions and derivatives at the quadrature point (xa,ya).                                          +
-                                                                                                              for m=1:4                                                                                                        +
-                                                                                                              mm=nd*(m-1);                                                                                                     +
-                                                                                                              Ds = TL*Zs3412d2{m,k}*Ti{m};                                                                                     +
-                                                                                                              Sx(:,mm+ND) = [Ds(2,:); -Ds(1,:)];     % [Pyx, -Pxx]                                                             +
-                                                                                                              Sy(:,mm+ND) = [Ds(3,:); -Ds(2,:)];     % [Pyy, -Pxy]                                                             +
-                                                                                                              end  % loop m                                                                                                    +
-                                                                                                                                                                                                                               +
-                                                                                                              Dm = Dm+(Sx'*Sx+Sy'*Sy)*(wt(k)*Det);                                                                             +
-                                                                                                                                                                                                                               +
-                                                                                                              end  % end loop k over quadrature points                                                                         +
-                                                                                                                                                                                                                               +
-                                                                                                              gf=zeros(nd4,1);                                                                                                 +
-                                                                                                              m=0;                                                                                                             +
-                                                                                                              for n=1:4                 % Loop over element nodes                                                              +
-                                                                                                              gf(m+ND)=(nn2nft(1,Elcon(n))-1)+ND;  % Get global freedoms                                                       +
-                                                                                                              m=m+nd;                                                                                                          +
-                                                                                                              end                                                                                                              +
-                                                                                                                                                                                                                               +
-                                                                                                              RowNdx=repmat(gf,1,nd4);      % Row indices                                                                      +
-                                                                                                              ColNdx=RowNdx';               % Col indices                                                                      +
-                                                                                                                                                                                                                               +
-                                                                                                              Dm = reshape(Dm,nd4*nd4,1);                                                                                      +
-                                                                                                              RowNdx=reshape(RowNdx,nd4*nd4,1);                                                                                +
-                                                                                                              ColNdx=reshape(ColNdx,nd4*nd4,1);                                                                                +
-                                                                                                                                                                                                                               +
-                                                                                                              return;                                                                                                          +
-                                                                                                                                                                                                                               +
-                                                                                                              % -------------------------------------------------------------------                                            +
-                                                                                                                                                                                                                               +
-                                                                                                              function P2=D3s(ni,q,r)                                                                                          +
-                                                                                                              % Second derivatives [Pxx; Pxy; Pyy] of simple cubic stream function.                                            +
-                                                                                                              qi=ni(1); q0=q*ni(1); q1=1+q0;                                                                                   +
-                                                                                                              ri=ni(2); r0=r*ni(2); r1=1+r0;                                                                                   +
-                                                                                                              P2=[-.75*qi^2*(r0+1)*q0, 0, -.25*qi*(r0+1)*(3*q0+1); ...                                                         +
-                                                                                                              .125*qi*ri*(4-3*(q^2+r^2)), .125*qi*(r0+1)*(3*r0-1), ...                                                         +
-                                                                                                              -.125*ri*(q0+1)*(3*q0-1); -.75*ri^2*(q0+1)*r0, .25*ri*(q0+1)*(3*r0+1), 0] ;                                      +
-                                                                                                              return;                                                                                                          +
-                                                                                                                                                                                                                               +
-                                                                                                              function G=GBL(ni,q,r)                                                                                           +
-                                                                                                              % Transposed gradient (derivatives) of scalar bilinear mapping function.                                         +
-                                                                                                              % The parameter ni can be a vector of coordinate pairs.                                                          +
-                                                                                                              G=[.25*ni(:,1).*(1+ni(:,2).*r), .25*ni(:,2).*(1+ni(:,1).*q)];                                                    +
-                                                                                                              return;                                                                                                          +
-
+ -
+
-                                                                                                              function [Cm,RowNdx,ColNdx]=CMatW(Xe,Elcon,nn2nft,Vdof)                                                          +
-                                                                                                              %CMATW - Returns the affine-mapped element convection matrix for the simple cubic Hermite                        +
-                                                                                                              %   basis functions on 4-node straight-sided quadrilateral elements with 3 DOF                                   +
-                                                                                                              %   per node using Gauss quadrature on the reference square and row/col indices.                                 +
-                                                                                                              % The columns of the array Vdof must contain the 3 nodal degree-of-freedom                                       +
-                                                                                                              %   vectors in the proper nodal order.                                                                           +
-                                                                                                              % The degrees of freedom in Vdof are the stream function and the two components                                  +
-                                                                                                              %   u and v of the solenoidal velocity vector.                                                                   +
-                                                                                                              % The assumed nodal numbering starts with 1 at the lower left corner of the                                      +
-                                                                                                              %   element and proceeds counter-clockwise around the element.                                                   +
-                                                                                                              %                                                                                                                +
-                                                                                                              % Usage:                                                                                                         +
-                                                                                                              %   [CM,Rndx,Cndx] = CMatW(Xe,Elcon,nn2nft,Vdof)                                                                 +
-                                                                                                              %   Xe(1,:) -  x-coordinates of corner nodes of element.                                                         +
-                                                                                                              %   Xe(2,:) -  y-coordinates of corner nodes of element.                                                         +
-                                                                                                              %   Elcon - this element connectivity matrix                                                                     +
-                                                                                                              %   nn2nft - global number and type of DOF at each node                                                          +
-                                                                                                              %   Vdof  - (3x4) array of stream function, velocity components, and second                                      +
-                                                                                                              %     stream function derivatives to specify the velocity over the element.                                      +
-                                                                                                              %                                                                                                                +
-                                                                                                              % Jonas Holdeman, August 2007, revised  June 2011                                                                +
-                                                                                                                                                                                                                               +
-                                                                                                              % Constants and fixed data                                                                                       +
-                                                                                                              nd = 3;  nd4=4*nd;  ND=1:nd;    % nd = number of dofs per node,                                                  +
-                                                                                                              nn=[-1 -1; 1 -1; 1 1; -1 1];    % defines local nodal order                                                      +
-                                                                                                                                                                                                                               +
-                                                                                                              % Define 5-point quadrature data once, on first call.                                                            +
-                                                                                                              % Gaussian weights and absissas to integrate 9th degree polynomials exactly.                                     +
-                                                                                                              global GQ5;                                                                                                      +
-                                                                                                              if (isempty(GQ5))   % 5-point quadrature data defined? If not, define arguments & weights.                       +
-                                                                                                              Aq=[-.906179845938664,-.538469310105683, .0,               .538469310105683, .906179845938664];                  +
-                                                                                                              Hq=[ .236926885056189, .478628670499366, .568888888888889, .478628670499366, .236926885056189];                  +
-                                                                                                              GQ5.size=25; GQ5.xa=[Aq;Aq;Aq;Aq;Aq]; GQ5.ya=GQ5.xa';                                                            +
-                                                                                                              wt=[Hq;Hq;Hq;Hq;Hq]; GQ5.wt=wt.*wt';                                                                             +
-                                                                                                              end                                                                                                              +
-                                                                                                                                                                                                                               +
-                                                                                                              xa=GQ5.xa; ya=GQ5.ya; wt=GQ5.wt; Nq=GQ5.size;  % Use GQ5 (5x5) for exact integration                             +
-                                                                                                                                                                                                                               +
-                                                                                                              % -----------------------------------------------                                                                +
-                                                                                                              global Zs3412D2c;  global ZS3412c;                                                                               +
-                                                                                                                                                                                                                               +
-                                                                                                              if (isempty(ZS3412c)|isempty(Zs3412D2c)|size(Zs3412D2c,2)~=Nq)                                                   +
-                                                                                                              % Evaluate and save/cache the set of shape functions at quadrature pts.                                          +
-                                                                                                              Zs3412D2c=cell(4,Nq);  ZS3412c=cell(4,Nq);                                                                       +
-                                                                                                              for k=1:Nq                                                                                                       +
-                                                                                                              for m=1:4                                                                                                        +
-                                                                                                              ZS3412c{m,k}= Sr(nn(m,:),xa(k),ya(k));                                                                           +
-                                                                                                              Zs3412D2c{m,k}=D3s(nn(m,:),xa(k),ya(k));                                                                         +
-                                                                                                              end                                                                                                              +
-                                                                                                              end                                                                                                              +
-                                                                                                              end  % if(isempty(*))                                                                                            +
-                                                                                                                                                                                                                               +
-                                                                                                              % --------------- End fixed data ----------------                                                                +
-                                                                                                                                                                                                                               +
-                                                                                                              Ti=cell(4);                                                                                                      +
-                                                                                                              %                                                                                                                +
-                                                                                                              for m=1:4                                                                                                        +
-                                                                                                              Jt=Xe*GBL(nn(:,:),nn(m,1),nn(m,2));                                                                              +
-                                                                                                              JtiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)];   % det(J)*inv(J)                                                     +
-                                                                                                              Ti{m}=blkdiag(1,JtiD);                                                                                           +
-                                                                                                              end                                                                                                              +
-                                                                                                                                                                                                                               +
-                                                                                                              % Move Jacobian evaluation inside k-loop for general convex quadrilateral.                                       +
-                                                                                                              % Jt=[x_q, x_r; y_q, y_r];                                                                                       +
-                                                                                                                                                                                                                               +
-                                                                                                              Cm=zeros(nd4,nd4); Rcm=zeros(nd4,1);                                                                             +
-                                                                                                              S=zeros(2,nd4); Sx=zeros(2,nd4); Sy=zeros(2,nd4);  % Pre-allocate arrays                                         +
-                                                                                                                                                                                                                               +
-                                                                                                              % Begin loop over Gauss-Legendre quadrature points.                                                              +
-                                                                                                              for k=1:Nq                                                                                                       +
-                                                                                                                                                                                                                               +
-                                                                                                              Jt=Xe*GBL(nn(:,:),xa(k),ya(k));         % transpose of Jacobian at (xa,ya)                                       +
-                                                                                                              Det=Jt(1,1)*Jt(2,2)-Jt(1,2)*Jt(2,1);    % Determinant of Jt & J                                                  +
-                                                                                                              Jtd=Jt/Det;                                                                                                      +
-                                                                                                              TL=[Jt(2,2)^2, -2*Jt(2,1)*Jt(2,2),  Jt(2,1)^2; ...                                                               +
-                                                                                                              -Jt(1,2)*Jt(2,2), Jt(1,1)*Jt(2,2)+Jt(2,1)*Jt(1,2), -Jt(1,1)*Jt(2,1); ...                                         +
-                                                                                                              Jt(1,2)^2, -2*Jt(1,1)*Jt(1,2),  Jt(1,1)^2 ]/Det^2;                                                               +
-                                                                                                                                                                                                                               +
-                                                                                                              % Initialize functions and derivatives at the quadrature point (xa,ya).                                          +
-                                                                                                              for m=1:4                                                                                                        +
-                                                                                                              mm=nd*(m-1);                                                                                                     +
-                                                                                                              S(:,mm+ND) = Jtd*ZS3412c{m,k}*Ti{m};                                                                             +
-                                                                                                              Ds = TL*Zs3412D2c{m,k}*Ti{m};                                                                                    +
-                                                                                                              Sx(:,mm+ND) = [Ds(2,:); -Ds(1,:)];     % [Pyx, -Pxx]                                                             +
-                                                                                                              Sy(:,mm+ND) = [Ds(3,:); -Ds(2,:)];     % [Pyy, -Pxy]                                                             +
-                                                                                                              end  % loop m                                                                                                    +
-                                                                                                                                                                                                                               +
-                                                                                                              % Compute the fluid velocity at the quadrature point.                                                            +
-                                                                                                              U = S*Vdof(:);                                                                                                   +
-                                                                                                              % Submatrix ordered by psi,u,v                                                                                   +
-                                                                                                              Cm = Cm + S'*(U(1)*Sx+U(2)*Sy)*(wt(k)*Det);                                                                      +
-                                                                                                              end    % end loop k over quadrature points                                                                       +
-                                                                                                                                                                                                                               +
-                                                                                                              gf=zeros(nd4,1);                                                                                                 +
-                                                                                                              m=0;                                                                                                             +
-                                                                                                              for n=1:4                 % Loop over element nodes                                                              +
-                                                                                                              gf(m+ND)=(nn2nft(1,Elcon(n))-1)+ND;  % Get global freedoms                                                       +
-                                                                                                              m=m+nd;                                                                                                          +
-                                                                                                              end                                                                                                              +
-                                                                                                                                                                                                                               +
-                                                                                                              RowNdx=repmat(gf,1,nd4);      % Row indices                                                                      +
-                                                                                                              ColNdx=RowNdx';               % Col indices                                                                      +
-                                                                                                                                                                                                                               +
-                                                                                                              Cm = reshape(Cm,nd4*nd4,1);                                                                                      +
-                                                                                                              RowNdx=reshape(RowNdx,nd4*nd4,1);                                                                                +
-                                                                                                              ColNdx=reshape(ColNdx,nd4*nd4,1);                                                                                +
-                                                                                                              return;                                                                                                          +
-                                                                                                                                                                                                                               +
-                                                                                                              % ----------------------------------------------------------------------------                                   +
-                                                                                                                                                                                                                               +
-                                                                                                              function P2=D3s(ni,q,r)                                                                                          +
-                                                                                                              % Second derivatives [Pxx; Pxy; Pyy] of simple cubic stream function.                                            +
-                                                                                                              qi=ni(1); q0=q*ni(1); q1=1+q0;                                                                                   +
-                                                                                                              ri=ni(2); r0=r*ni(2); r1=1+r0;                                                                                   +
-                                                                                                              P2=[-.75*qi^2*(r0+1)*q0, 0, -.25*qi*(r0+1)*(3*q0+1); ...                                                         +
-                                                                                                              .125*qi*ri*(4-3*(q^2+r^2)), .125*qi*(r0+1)*(3*r0-1), -.125*ri*(q0+1)*(3*q0-1); ...                               +
-                                                                                                              -.75*ri^2*(q0+1)*r0, .25*ri*(q0+1)*(3*r0+1), 0] ;                                                                +
-                                                                                                              return;                                                                                                          +
-                                                                                                                                                                                                                               +
-                                                                                                              function S=Sr(ni,q,r)                                                                                            +
-                                                                                                              %S  Array of solenoidal basis functions on rectangle.                                                            +
-                                                                                                              qi=ni(1); q0=q*ni(1); q1=1+q0;                                                                                   +
-                                                                                                              ri=ni(2); r0=r*ni(2); r1=1+r0;                                                                                   +
-                                                                                                              % array of solenoidal vectors                                                                                    +
-                                                                                                              S=[ .125*ri*q1*(q0*(1-q0)+3*(1-r^2)), .125*q1*r1*(3*r0-1),    .125*ri/qi*q1^2*(1-q0); ...                        +
-                                                                                                              -.125*qi*r1*(r0*(1-r0)+3*(1-q^2)), .125*qi/ri*r1^2*(1-r0), .125*q1*r1*(3*q0-1)];                                 +
-                                                                                                              return;                                                                                                          +
-                                                                                                                                                                                                                               +
-                                                                                                              function G=GBL(ni,q,r)                                                                                           +
-                                                                                                              % Transposed gradient (derivatives) of scalar bilinear mapping function.                                         +
-                                                                                                              % The parameter ni can be a vector of coordinate pairs.                                                          +
-                                                                                                              G=[.25*ni(:,1).*(1+ni(:,2).*r), .25*ni(:,2).*(1+ni(:,1).*q)];                                                    +
-                                                                                                              return;                                                                                                          +
-
+ -
+
-                                                                                                              function [P,Px,Py] = GetPresW(NumNod,NodNdx,Elcon,nn2nft,Xgrid,Ygrid,Q,EBCPr,nu)                                 +
-                                                                                                              %GETPRESW - Compute continuous simple cubic pressure and derivatives from (simple-cubic)                         +
-                                                                                                              %  velocity field on general quadrilateral grid (bilinear geometric mapping).                                    +
-                                                                                                              %                                                                                                                +
-                                                                                                              % Inputs                                                                                                         +
-                                                                                                              %  NumNod - number of nodes                                                                                      +
-                                                                                                              %  NodNdx - nodal index into Xgrid and Ygrid                                                                     +
-                                                                                                              %  Elcon  - element connectivity, nodes in element                                                               +
-                                                                                                              %  nn2nft - global number and type of (non-pressure) DOF at each node                                            +
-                                                                                                              %  Xgrid  - array of nodal x-coordinates                                                                         +
-                                                                                                              %  Ygrid  - array of nodal y-coordinates                                                                         +
-                                                                                                              %  Q      - array of DOFs for cubic velocity elements                                                            +
-                                                                                                              %  EBCp   - essential pressure boundary condition structure                                                      +
-                                                                                                              %    EBCp.nodn - node number of fixed pressure node                                                              +
-                                                                                                              %    EBCp.val  - value of pressure                                                                               +
-                                                                                                              %  nu - diffusion coefficient                                                                                    +
-                                                                                                              % Outputs                                                                                                        +
-                                                                                                              %  P  - pressure                                                                                                 +
-                                                                                                              %  Px - x-derivative of pressure                                                                                 +
-                                                                                                              %  Py - y-derivative of pressure                                                                                 +
-                                                                                                              % Uses                                                                                                           +
-                                                                                                              %  ilu_gmres_with_EBC - to solve the system with essential/Dirichlet BCs                                         +
-                                                                                                              %  GQ3, GQ4, GQ5  - quadrature rules.                                                                            +
-                                                                                                                                                                                                                               +
-                                                                                                              % Jonas Holdeman,   January 2007, revised June 2011                                                              +
-                                                                                                                                                                                                                               +
-                                                                                                              % Constants and fixed data                                                                                       +
-                                                                                                              nn=[-1 -1; 1 -1; 1 1; -1 1];  % defines local nodal order                                                        +
-                                                                                                              nnd = 4;                      % Number of nodes in elements                                                      +
-                                                                                                              nd = 3;  ND=1:nd;             % Number DOFs in velocity fns (bicubic-derived)                                    +
-                                                                                                              np = 3;                       % Number DOFs in pressure fns (simple cubic)                                       +
-                                                                                                              % Parameters for GMRES solver                                                                                    +
-                                                                                                              GMRES.Tolerance=1.e-9;                                                                                           +
-                                                                                                              GMRES.MaxIterates=8;                                                                                             +
-                                                                                                              GMRES.MaxRestarts=6;                                                                                             +
-                                                                                                              DropTol = 1e-7;                  % Drop tolerence for ilu preconditioner                                         +
-                                                                                                                                                                                                                               +
-                                                                                                              % Define 3-point quadrature data once, on first call (if needed).                                                +
-                                                                                                              % Gaussian weights and absissas to integrate 5th degree polynomials exactly.                                     +
-                                                                                                              global GQ3;                                                                                                      +
-                                                                                                              if (isempty(GQ3))       % Define 3-point quadrature data once, on first call.                                    +
-                                                                                                              Aq=[-.774596669241483, .000000000000000,.774596669241483]; %Abs                                                  +
-                                                                                                              Hq=[ .555555555555556, .888888888888889,.555555555555556]; %Wts                                                  +
-                                                                                                              GQ3.size=9; GQ3.xa=[Aq;Aq;Aq]; GQ3.ya=GQ3.xa';                                                                   +
-                                                                                                              wt=[Hq;Hq;Hq]; GQ3.wt=wt.*wt';                                                                                   +
-                                                                                                              end                                                                                                              +
-                                                                                                              % Define 4-point quadrature data once, on first call (if needed).                                                +
-                                                                                                              % Gaussian weights and absissas to integrate 7th degree polynomials exactly.                                     +
-                                                                                                              global GQ4;                                                                                                      +
-                                                                                                              if (isempty(GQ4))       % Define 4-point quadrature data once, on first call.                                    +
-                                                                                                              Aq=[-.861136311594053,-.339981043584856,.339981043584856, .861136311594053]; %Abs                                +
-                                                                                                              Hq=[ .347854845137454, .652145154862546,.652145154862546, .347854845137454]; %Wts                                +
-                                                                                                              GQ4.size=16; GQ4.xa=[Aq;Aq;Aq;Aq]; GQ4.ya=GQ4.xa';                                                               +
-                                                                                                              wt=[Hq;Hq;Hq;Hq]; GQ4.wt=wt.*wt';     % 4x4                                                                      +
-                                                                                                              end                                                                                                              +
-                                                                                                              % Define 5-point quadrature data once, on first call (if needed).                                                +
-                                                                                                              % Gaussian weights and absissas to integrate 9th degree polynomials exactly.                                     +
-                                                                                                              global GQ5;                                                                                                      +
-                                                                                                              if (isempty(GQ5))   % Has 5-point quadrature data been defined? If not, define arguments & weights.              +
-                                                                                                              Aq=[-.906179845938664,-.538469310105683, .0,               .538469310105683, .906179845938664];                  +
-                                                                                                              Hq=[ .236926885056189, .478628670499366, .568888888888889, .478628670499366, .236926885056189];                  +
-                                                                                                              GQ5.size=25; GQ5.xa=[Aq;Aq;Aq;Aq;Aq]; GQ5.ya=GQ5.xa';                                                            +
-                                                                                                              wt=[Hq;Hq;Hq;Hq;Hq]; GQ5.wt=wt.*wt';     % 5x5                                                                   +
-                                                                                                              end                                                                                                              +
-                                                                                                              % -------------- end fixed data ------------------------                                                         +
-                                                                                                                                                                                                                               +
-                                                                                                              NumEl=size(Elcon,2);            % Number of elements                                                             +
-                                                                                                              [NumNy,NumNx]=size(Xgrid);                                                                                       +
-                                                                                                              NumNod=NumNy*NumNx;             % Number of nodes                                                                +
-                                                                                                              MxVDof=nd*NumNod;               % Max number velocity DOFs                                                       +
-                                                                                                              MxPDof=np*NumNod;               % Max number pressure DOFs                                                       +
-                                                                                                                                                                                                                               +
-                                                                                                              % We can use the same nodal connectivities (Elcon) for pressure elements,                                        +
-                                                                                                              %  but need new pointers from nodes to pressure DOFs                                                             +
-                                                                                                              nn2nftp=zeros(2,NumNod); % node number -> pressure nf & nt                                                       +
-                                                                                                              nf=-np+1;  nt=1;                                                                                                 +
-                                                                                                              for n=1:NumNod                                                                                                   +
-                                                                                                              nf=nf+np;               % all nodes have 3 dofs                                                                  +
-                                                                                                              nn2nftp(:,n)=[nf;nt];   % dof number & type (all nodes type 1)                                                   +
-                                                                                                              end;                                                                                                             +
-                                                                                                                                                                                                                               +
-                                                                                                              % Allocate space for pressure matrix, velocity data                                                              +
-                                                                                                              pMat = spalloc(MxPDof,MxPDof,30*MxPDof);   % allocate P mass matrix                                              +
-                                                                                                              Vdata = zeros(MxPDof,1);       % allocate for velocity data (RHS)                                                +
-                                                                                                              Qp=zeros(MxPDof,1);       % Nodal pressure DOFs                                                                  +
-                                                                                                                                                                                                                               +
-                                                                                                              % Begin essential boundary conditions, allocate space                                                            +
-                                                                                                              MaxPBC = 1;                                                                                                      +
-                                                                                                              EBCp.Mxdof=MxPDof;                                                                                               +
-                                                                                                              % Essential boundary condition for pressure                                                                      +
-                                                                                                              EBCp.dof = nn2nftp(1,EBCPr.nodn(1));  % Degree-of-freedom index                                                  +
-                                                                                                              EBCp.val = EBCPr.val;                         % Pressure Dof value                                               +
-                                                                                                                                                                                                                               +
-                                                                                                              % partion out essential (Dirichlet) dofs                                                                         +
-                                                                                                              p_vec = [1:EBCp.Mxdof]';         % List of all dofs                                                              +
-                                                                                                              EBCp.p_vec_undo = zeros(1,EBCp.Mxdof);                                                                           +
-                                                                                                              % form a list of non-diri dofs                                                                                   +
-                                                                                                              EBCp.ndro = p_vec(~ismember(p_vec, EBCp.dof));	% list of non-diri dofs                                           +
-                                                                                                              % calculate p_vec_undo to restore Q to the original dof ordering                                                 +
-                                                                                                              EBCp.p_vec_undo([EBCp.ndro;EBCp.dof]) = [1:EBCp.Mxdof]; %p_vec';                                                 +
-                                                                                                                                                                                                                               +
-                                                                                                              Qp(EBCp.dof(1))=EBCp.val(1);                                                                                     +
-                                                                                                                                                                                                                               +
-                                                                                                              Vdof = zeros(nd,nnd);             % Nodal velocity DOFs                                                          +
-                                                                                                              Xe = zeros(2,nnd);                                                                                               +
-                                                                                                                                                                                                                               +
-                                                                                                              % BEGIN GLOBAL MATRIX ASSEMBLY                                                                                   +
-                                                                                                              for ne=1:NumEl                                                                                                   +
-                                                                                                              for k=1:4                                                                                                        +
-                                                                                                              ki=NodNdx(:,Elcon(k,ne));                                                                                        +
-                                                                                                              Xe(:,k)=[Xgrid(ki(2),ki(1));Ygrid(ki(2),ki(1))];                                                                 +
-                                                                                                              end                                                                                                              +
-                                                                                                              % Get stream function and velocities                                                                             +
-                                                                                                              for n=1:nnd                                                                                                      +
-                                                                                                              Vdof(ND,n)=Q((nn2nft(1,Elcon(n,ne))-1)+ND); % Loop over local nodes of element                                   +
-                                                                                                              end                                                                                                              +
-                                                                                                              [pMmat,Rndx,Cndx] = pMassMat(Xe,Elcon(:,ne),nn2nftp);     % Pressure "mass" matrix                               +
-                                                                                                              pMat=pMat+sparse(Rndx,Cndx,pMmat,MxPDof,MxPDof);  % Global mass assembly                                         +
-                                                                                                                                                                                                                               +
-                                                                                                              [CDat,RRndx] = PCDat(Xe,Elcon(:,ne),nn2nftp,Vdof);   % Convective data term                                      +
-                                                                                                              Vdata([RRndx]) = Vdata([RRndx])-CDat(:);                                                                         +
-                                                                                                                                                                                                                               +
-                                                                                                              [DDat,RRndx] = PDDatL(Xe,Elcon(:,ne),nn2nftp,Vdof);   % Diffusive data term                                      +
-                                                                                                              Vdata([RRndx]) = Vdata([RRndx]) + nu*DDat(:); % +nu*DDat;                                                        +
-                                                                                                              end;   % Loop ne                                                                                                 +
-                                                                                                              % END GLOBAL MATRIX ASSEMBLY                                                                                     +
-                                                                                                                                                                                                                               +
-                                                                                                              Vdatap=Vdata(EBCp.ndro)-pMat(EBCp.ndro,EBCp.dof)*EBCp.val;                                                       +
-                                                                                                              pMatr=pMat(EBCp.ndro,EBCp.ndro);                                                                                 +
-                                                                                                              Qs=Qp(EBCp.ndro);            % Non-Dirichlet (active) dofs                                                       +
-                                                                                                                                                                                                                               +
-                                                                                                              Pr=ilu_gmres_with_EBC(pMatr,Vdatap,[],GMRES,Qs,DropTol);                                                         +
-                                                                                                                                                                                                                               +
-                                                                                                              Qp=[Pr;EBCp.val];         % Augment active dofs with esential (Dirichlet) dofs                                   +
-                                                                                                              Qp=Qp(EBCp.p_vec_undo);      % Restore natural order                                                             +
-                                                                                                              Qp=reshape(Qp,np,NumNod);                                                                                        +
-                                                                                                              P= reshape(Qp(1,:),NumNy,NumNx);                                                                                 +
-                                                                                                              Px=reshape(Qp(2,:),NumNy,NumNx);                                                                                 +
-                                                                                                              Py=reshape(Qp(3,:),NumNy,NumNx);                                                                                 +
-                                                                                                              return;                                                                                                          +
-                                                                                                              % >>>>>>>>>>>>> End pressure recovery <<<<<<<<<<<<<                                                              +
-                                                                                                                                                                                                                               +
-                                                                                                              % -------------------- function pMassMat ----------------------------                                            +
-                                                                                                                                                                                                                               +
-                                                                                                              function [MM,Rndx,Cndx]=pMassMat(Xe,Elcon,nn2nftp)                                                               +
-                                                                                                              % Simple cubic gradient element "mass" matrix                                                                    +
-                                                                                                              % -------------- Constants and fixed data -----------------                                                      +
-                                                                                                              global GQ4;                                                                                                      +
-                                                                                                              xa=GQ4.xa; ya=GQ4.ya; wt=GQ4.wt; Nq=GQ4.size;                                                                    +
-                                                                                                              nn=[-1 -1; 1 -1; 1 1; -1 1]; % defines local nodal order                                                         +
-                                                                                                              nnd=4;                                                                                                           +
-                                                                                                              np=3; np4=nnd*np; NP=1:np;                                                                                       +
-                                                                                                              %                                                                                                                +
-                                                                                                              global ZG3412pm;                                                                                                 +
-                                                                                                              if (isempty(ZG3412pm)|size(ZG3412pm,2)~=Nq)                                                                      +
-                                                                                                              % Evaluate and save/cache the set of shape functions at quadrature pts.                                          +
-                                                                                                              ZG3412pm=cell(nnd,Nq);                                                                                           +
-                                                                                                              for k=1:Nq                                                                                                       +
-                                                                                                              for m=1:nnd                                                                                                      +
-                                                                                                              ZG3412pm{m,k}=Gr(nn(m,:),xa(k),ya(k));                                                                           +
-                                                                                                              end                                                                                                              +
-                                                                                                              end                                                                                                              +
-                                                                                                              end  % if(isempty(*))                                                                                            +
-                                                                                                              % --------------------- end fixed data -----------------                                                         +
-                                                                                                                                                                                                                               +
-                                                                                                              TGi=cell(nnd);                                                                                                   +
-                                                                                                              for m=1:nnd   % Loop over corner nodes                                                                           +
-                                                                                                              J=(Xe*GBL(nn(:,:),nn(m,1),nn(m,2)))'; % GBL is gradient of bilinear function                                     +
-                                                                                                              TGi{m} = blkdiag(1,J);                                                                                           +
-                                                                                                              end  % Loop m                                                                                                    +
-                                                                                                                                                                                                                               +
-                                                                                                              MM=zeros(np4,np4);  G=zeros(2,np4);   % Preallocate arrays                                                       +
-                                                                                                              for k=1:Nq                                                                                                       +
-                                                                                                              % Initialize functions and derivatives at the quadrature point (xa,ya).                                          +
-                                                                                                              J=(Xe*GBL(nn(:,:),xa(k),ya(k)))';         % transpose of Jacobian J at (xa,ya)                                   +
-                                                                                                              Det=J(1,1)*J(2,2)-J(1,2)*J(2,1);          % Determinant of J                                                     +
-                                                                                                              Ji=[J(2,2),-J(1,2); -J(2,1),J(1,1)]/Det;  % inverse of J                                                         +
-                                                                                                                                                                                                                               +
-                                                                                                              mm = 0;                                                                                                          +
-                                                                                                              for m=1:nnd                                                                                                      +
-                                                                                                              G(:,mm+NP) = Ji*ZG3412pm{m,k}*TGi{m};                                                                            +
-                                                                                                              mm = mm+np;                                                                                                      +
-                                                                                                              end  % loop m                                                                                                    +
-                                                                                                              MM = MM + G'*G*(wt(k)*Det);                                                                                      +
-                                                                                                              end        % end loop k (quadrature pts)                                                                         +
-                                                                                                                                                                                                                               +
-                                                                                                              gf=zeros(np4,1);          % array of global freedoms                                                             +
-                                                                                                              m=0;                                                                                                             +
-                                                                                                              for n=1:4                 % Loop over element nodes                                                              +
-                                                                                                              gf(m+NP)=(nn2nftp(1,Elcon(n))-1)+NP; % Get global freedoms                                                       +
-                                                                                                              m=m+np;                                                                                                          +
-                                                                                                              end                                                                                                              +
-                                                                                                                                                                                                                               +
-                                                                                                              Rndx=repmat(gf,1,np4);     % Row indices                                                                         +
-                                                                                                              Cndx=Rndx';                % Column indices                                                                      +
-                                                                                                                                                                                                                               +
-                                                                                                              MM = reshape(MM,1,np4*np4);                                                                                      +
-                                                                                                              Rndx=reshape(Rndx,1,np4*np4);                                                                                    +
-                                                                                                              Cndx=reshape(Cndx,1,np4*np4);                                                                                    +
-                                                                                                              return;                                                                                                          +
-                                                                                                                                                                                                                               +
-                                                                                                              % --------------------- funnction PCDat -----------------------                                                  +
-                                                                                                                                                                                                                               +
-                                                                                                              function [PC,Rndx]=PCDat(Xe,Elcon,nn2nftp,Vdof)                                                                  +
-                                                                                                              % Evaluate convective force data                                                                                 +
-                                                                                                              % ----------- Constants and fixed data ---------------                                                           +
-                                                                                                              global GQ5;                                                                                                      +
-                                                                                                              xa=GQ5.xa; ya=GQ5.ya; wt=GQ5.wt; Nq=GQ5.size;                                                                    +
-                                                                                                              nn=[-1 -1; 1 -1; 1 1; -1 1]; % defines local nodal order                                                         +
-                                                                                                              nnd=4;    % number of nodes                                                                                      +
-                                                                                                              np = 3;  np4=4*np;  NP=1:np;                                                                                     +
-                                                                                                              nd = 3;  nd4=4*nd;  ND=1:nd;                                                                                     +
-                                                                                                              %                                                                                                                +
-                                                                                                              global ZS3412pc; global ZSX3412pc; global ZSY3412pc; global ZG3412pc;                                            +
-                                                                                                              if (isempty(ZS3412pc)|size(ZS3412pc,2)~=Nq)                                                                      +
-                                                                                                              % Evaluate and save/cache the set of shape functions at quadrature pts.                                          +
-                                                                                                              ZS3412pc=cell(nnd,Nq); ZSX3412pc=cell(nnd,Nq);                                                                   +
-                                                                                                              ZSY3412pc=cell(nnd,Nq); ZG3412pc=cell(nnd,Nq);                                                                   +
-                                                                                                              for k=1:Nq                                                                                                       +
-                                                                                                              for m=1:nnd                                                                                                      +
-                                                                                                              ZS3412pc{m,k} =Sr(nn(m,:),xa(k),ya(k));                                                                          +
-                                                                                                              ZSX3412pc{m,k}=Sxr(nn(m,:),xa(k),ya(k));                                                                         +
-                                                                                                              ZSY3412pc{m,k}=Syr(nn(m,:),xa(k),ya(k));                                                                         +
-                                                                                                              ZG3412pc{m,k}=Gr(nn(m,:),xa(k),ya(k));                                                                           +
-                                                                                                              end  % loop m over nodes                                                                                         +
-                                                                                                              end  % loop k over Nq                                                                                            +
-                                                                                                              end  % if(isempty(*))                                                                                            +
-                                                                                                              % ----------------- end fixed data ------------------                                                            +
-                                                                                                                                                                                                                               +
-                                                                                                              Ti=cell(nnd);  TGi=cell(nnd);                                                                                    +
-                                                                                                              for m=1:nnd   % Loop over corner nodes                                                                           +
-                                                                                                              J=(Xe*GBL(nn(:,:),nn(m,1),nn(m,2)))';   % Jacobian at (xa,ya)                                                    +
-                                                                                                              Det=J(1,1)*J(2,2)-J(1,2)*J(2,1);      % Determinant of J & Jt                                                    +
-                                                                                                              JiD=[J(2,2),-J(1,2); -J(2,1),J(1,1)]; % inv(J)*det(J)                                                            +
-                                                                                                              Ti{m} = blkdiag(1,JiD');                                                                                         +
-                                                                                                              TGi{m} = blkdiag(1,J);                                                                                           +
-                                                                                                              end  % Loop m over corner nodes                                                                                  +
-                                                                                                                                                                                                                               +
-                                                                                                              PC=zeros(np4,1);                                                                                                 +
-                                                                                                              S=zeros(2,nd4);  Sx=zeros(2,nd4);  Sy=zeros(2,nd4);  G=zeros(2,np4);                                             +
-                                                                                                                                                                                                                               +
-                                                                                                              for k=1:Nq      % Loop over quadrature points                                                                    +
-                                                                                                              J=(Xe*GBL(nn(:,:),xa(k),ya(k)))';       % Jacobian at (xa,ya)                                                    +
-                                                                                                              Det=J(1,1)*J(2,2)-J(1,2)*J(2,1);      % Determinant of J & Jt                                                    +
-                                                                                                              Jtbd=(J/Det)';                        % transpose(J/D)                                                           +
-                                                                                                              JiD=[J(2,2),-J(1,2); -J(2,1),J(1,1)]; % inv(J)*det(J)                                                            +
-                                                                                                              Ji=JiD/Det;                           % inverse(J)                                                               +
-                                                                                                              for m=1:4         % Loop over element nodes                                                                      +
-                                                                                                              mm=nd*(m-1);                                                                                                     +
-                                                                                                              S(:,mm+ND) =Jtbd*ZS3412pc{m,k}*Ti{m};                                                                            +
-                                                                                                              Sx(:,mm+ND)=Jtbd*(Ji(1,1)*ZSX3412pc{m,k}+Ji(1,2)*ZSY3412pc{m,k})*Ti{m};                                          +
-                                                                                                              Sy(:,mm+ND)=Jtbd*(Ji(2,1)*ZSX3412pc{m,k}+Ji(2,2)*ZSY3412pc{m,k})*Ti{m};                                          +
-                                                                                                              mm=np*(m-1);                                                                                                     +
-                                                                                                              G(:,mm+NP)=Ji*ZG3412pc{m,k}*TGi{m};                                                                              +
-                                                                                                              end    % end loop over element nodes                                                                             +
-                                                                                                                                                                                                                               +
-                                                                                                              % Compute the fluid velocities at the quadrature point.                                                          +
-                                                                                                              U = S*Vdof(:);                                                                                                   +
-                                                                                                              Ux = Sx*Vdof(:);                                                                                                 +
-                                                                                                              Uy = Sy*Vdof(:);                                                                                                 +
-                                                                                                              UgU = U(1)*Ux+U(2)*Uy;   % U dot grad U                                                                          +
-                                                                                                              PC = PC + G'*UgU*(wt(k)*Det);                                                                                    +
-                                                                                                              end    % end loop over Nq quadrature points                                                                      +
-                                                                                                                                                                                                                               +
-                                                                                                              gf=zeros(1,np4);          % array of global freedoms                                                             +
-                                                                                                              m=0;                                                                                                             +
-                                                                                                              for n=1:4                 % Loop over element nodes                                                              +
-                                                                                                              gf(m+NP)=(nn2nftp(1,Elcon(n))-1)+NP; % Get global freedoms                                                       +
-                                                                                                              m=m+np;                                                                                                          +
-                                                                                                              end                                                                                                              +
-                                                                                                              Rndx=gf;                                                                                                         +
-                                                                                                              PC = reshape(PC,1,np4);                                                                                          +
-                                                                                                              return;                                                                                                          +
-                                                                                                                                                                                                                               +
-                                                                                                              % ----------------- function PDDatL -------------------------                                                    +
-                                                                                                                                                                                                                               +
-                                                                                                              function [PD,Rndx]=PDDatL(Xe,Elcon,nn2nftp,Vdof)                                                                 +
-                                                                                                              % Evaluate diffusive force data (Laplacian form)                                                                 +
-                                                                                                              % --------------- Constants and fixed data -------------------                                                   +
-                                                                                                              global GQ3;                                                                                                      +
-                                                                                                              xa=GQ3.xa; ya=GQ3.ya; wt=GQ3.wt; Nq=GQ3.size;                                                                    +
-                                                                                                              nnd=4;    % number of nodes                                                                                      +
-                                                                                                              nn=[-1 -1; 1 -1; 1 1; -1 1]; % defines local nodal order                                                         +
-                                                                                                              np = 3;  npdf=nnd*np;  NP=1:np;                                                                                  +
-                                                                                                              nd = 3;  nd4=nnd*nd;  ND=1:nd;                                                                                   +
-                                                                                                              global ZSXX3412pd; global ZSXY3412pd; global ZSYY3412pd; global ZG3412pd;                                        +
-                                                                                                              if (isempty(ZSXX3412pd)|size(ZSXX3412pd,2)~=Nq)                                                                  +
-                                                                                                              % Evaluate and save/cache the set of shape functions at quadrature pts.                                          +
-                                                                                                              ZSXX3412pd=cell(nnd,Nq); ZSXY3412pd=cell(nnd,Nq);                                                                +
-                                                                                                              ZSYY3412pd=cell(nnd,Nq);  ZG3412pd=cell(nnd,Nq);                                                                 +
-                                                                                                              global ZSYY3412pd;                                                                                               +
-                                                                                                              for k=1:Nq                                                                                                       +
-                                                                                                              for m=1:nnd                                                                                                      +
-                                                                                                              ZSXX3412pd{m,k}=Sxxr(nn(m,:),xa(k),ya(k));                                                                       +
-                                                                                                              ZSXY3412pd{m,k}=Sxyr(nn(m,:),xa(k),ya(k));                                                                       +
-                                                                                                              ZSYY3412pd{m,k}=Syyr(nn(m,:),xa(k),ya(k));                                                                       +
-                                                                                                              ZG3412pd{m,k}=Gr(nn(m,:),xa(k),ya(k));                                                                           +
-                                                                                                              end  % loop m over nodes                                                                                         +
-                                                                                                              end  % loop k over Nq                                                                                            +
-                                                                                                              end  % if(isempty(*))                                                                                            +
-                                                                                                              % ------------ end fixed data -------------------                                                                +
-                                                                                                              %                                                                                                                +
-                                                                                                              Ti=cell(nnd);  TGi=cell(nnd);                                                                                    +
-                                                                                                              for m=1:nnd   % Loop over corner nodes                                                                           +
-                                                                                                              J=(Xe*GBL(nn(:,:),nn(m,1),nn(m,2)))';   % Jacobian at (xa,ya)                                                    +
-                                                                                                              Det=J(1,1)*J(2,2)-J(1,2)*J(2,1);      % Determinant of J & Jt                                                    +
-                                                                                                              JiD=[J(2,2),-J(1,2); -J(2,1),J(1,1)]; % inv(J)*det(J)                                                            +
-                                                                                                              Ti{m} = blkdiag(1,JiD');                                                                                         +
-                                                                                                              TGi{m} = blkdiag(1,J);                                                                                           +
-                                                                                                              end  % Loop m over corner nodes                                                                                  +
-                                                                                                                                                                                                                               +
-                                                                                                              PD=zeros(npdf,1);                                                                                                +
-                                                                                                              Sxx=zeros(2,nd4);  Syy=zeros(2,nd4);  G=zeros(2,npdf);                                                           +
-                                                                                                              for k=1:Nq          % Loop over quadrature points                                                                +
-                                                                                                              Jt=(Xe*GBL(nn(:,:),xa(k),ya(k)));       % transpose of Jacobian at (xa,ya)                                       +
-                                                                                                              Det=Jt(1,1)*Jt(2,2)-Jt(1,2)*Jt(2,1);  % Determinant of Jt & J                                                    +
-                                                                                                              Jtd=Jt/Det;                                                                                                      +
-                                                                                                              JiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)];                                                                        +
-                                                                                                              Ji=JiD/Det;                                                                                                      +
-                                                                                                              for m=1:nnd        % Loop over element nodes                                                                     +
-                                                                                                              mm=nd*(m-1);    % This transform is approximate !!                                                               +
-                                                                                                              Sxx(:,mm+ND)=Jtd*(Ji(1,1)^2*ZSXX3412pd{m,k}+2*Ji(1,1)*Ji(1,2)*ZSXY3412pd{m,k}+Ji(1,2)^2*ZSYY3412pd{m,k})*Ti{m};  +
-                                                                                                              Syy(:,mm+ND)=Jtd*(Ji(2,1)^2*ZSXX3412pd{m,k}+2*Ji(2,1)*Ji(2,2)*ZSXY3412pd{m,k}+Ji(2,2)^2*ZSYY3412pd{m,k})*Ti{m};  +
-                                                                                                              mm=np*(m-1);                                                                                                     +
-                                                                                                              G(:,mm+NP) =Ji*ZG3412pd{m,k}*TGi{m};                                                                             +
-                                                                                                              end    % end loop over element nodes                                                                             +
-                                                                                                                                                                                                                               +
-                                                                                                              LapU = (Sxx+Syy)*Vdof(:);                                                                                        +
-                                                                                                              PD = PD+G'*LapU*(wt(k)*Det);                                                                                     +
-                                                                                                              end    % end loop over quadrature points                                                                         +
-                                                                                                                                                                                                                               +
-                                                                                                              gf=zeros(1,npdf);          % array of global freedoms                                                            +
-                                                                                                              m=0;  K=1:np;                                                                                                    +
-                                                                                                              for n=1:nnd                 % Loop over element nodes                                                            +
-                                                                                                              nfn=nn2nftp(1,Elcon(n))-1;  % Get global freedom                                                                 +
-                                                                                                              gf(m+NP)=(nn2nftp(1,Elcon(n))-1)+NP;                                                                             +
-                                                                                                              m=m+np;                                                                                                          +
-                                                                                                              end                                                                                                              +
-                                                                                                              Rndx=gf;                                                                                                         +
-                                                                                                              PD = reshape(PD,1,npdf);                                                                                         +
-                                                                                                              return;                                                                                                          +
-                                                                                                                                                                                                                               +
-                                                                                                              % ------------------------------------------------------------------------------                                 +
-                                                                                                              function gv=Gr(ni,q,r)                                                                                           +
-                                                                                                              qi=ni(1); q0=q*ni(1);                                                                                            +
-                                                                                                              ri=ni(2); r0=r*ni(2);                                                                                            +
-                                                                                                              % Cubic Hermite gradient                                                                                         +
-                                                                                                              gv=[1/8*qi*(1+r0)*(r0*(1-r0)+3*(1-q^2)), -1/8*(1+r0)*(1+q0)*(1-3*q0), ...                                        +
-                                                                                                              -1/8*qi/ri*(1-r^2)*(1+r0); ...                                                                                   +
-                                                                                                              1/8*ri*(1+q0)*(q0*(1-q0)+3*(1-r^2)),  -1/8/qi*ri*(1-q^2)*(1+q0), ...                                             +
-                                                                                                              -1/8*(1+q0)*(1+r0)*(1-3*r0)];                                                                                    +
-                                                                                                              return;                                                                                                          +
-                                                                                                                                                                                                                               +
-                                                                                                              function gx=Gxr(ni,q,r)                                                                                          +
-                                                                                                              %GRX - Cubic Hermite gradient basis functions for pressure gradient.                                             +
-                                                                                                              qi=ni(1); q0=q*ni(1);                                                                                            +
-                                                                                                              ri=ni(2); r0=r*ni(2);                                                                                            +
-                                                                                                              % x-derivative of irrotational vector                                                                            +
-                                                                                                              gx=[-3/4*qi^2*q0*(1+r0), 1/4*qi*(1+r0)*(1+3*q0), 0; ...                                                          +
-                                                                                                              1/8*qi*ri*(4-3*q0^2-3*r0^2), -1/8*ri*(1+q0)*(1-3*q0), -1/8*qi*(1+r0)*(1-3*r0)];                                  +
-                                                                                                              return;                                                                                                          +
-                                                                                                                                                                                                                               +
-                                                                                                              function gy=Gyr(ni,q,r)                                                                                          +
-                                                                                                              %GRY - Cubic Hermite gradient basis functions for pressure gradient.                                             +
-                                                                                                              qi=ni(1); q0=q*ni(1);                                                                                            +
-                                                                                                              ri=ni(2); r0=r*ni(2);                                                                                            +
-                                                                                                              % y-derivative of irrotational vector                                                                            +
-                                                                                                              gy=[1/8*qi*ri*(4-3*q0^2-3*r0^2), -1/8*ri*(1+q0)*(1-3*q0), -1/8*qi*(1+r0)*(1-3*r0); ...                           +
-                                                                                                              -3/4*ri^2*r0*(1+q0), 0, 1/4*ri*(1+q0)*(1+3*r0)];                                                                 +
-                                                                                                              return;                                                                                                          +
-                                                                                                                                                                                                                               +
-                                                                                                              % ------------------------------------------------------------------------------                                 +
-                                                                                                              function S=Sr(ni,q,r)                                                                                            +
-                                                                                                              %SR  Array of solenoidal basis functions.                                                                        +
-                                                                                                              qi=ni(1); q0=q*ni(1); q1=1+q0;                                                                                   +
-                                                                                                              ri=ni(2); r0=r*ni(2); r1=1+r0;                                                                                   +
-                                                                                                              % array of solenoidal vectors                                                                                    +
-                                                                                                              S=[ .125*ri*q1*(q0*(1-q0)+3*(1-r^2)), .125*q1*r1*(3*r0-1),    .125*ri/qi*q1^2*(1-q0); ...                        +
-                                                                                                              -.125*qi*r1*(r0*(1-r0)+3*(1-q^2)), .125*qi/ri*r1^2*(1-r0), .125*q1*r1*(3*q0-1)];                                 +
-                                                                                                              return;                                                                                                          +
-                                                                                                                                                                                                                               +
-                                                                                                              function S=Sxr(ni,q,r)                                                                                           +
-                                                                                                              %SXR  Array of x-derivatives of solenoidal basis functions.                                                      +
-                                                                                                              qi=ni(1); q0=q*ni(1); q1=1+q0;                                                                                   +
-                                                                                                              ri=ni(2); r0=r*ni(2); r1=1+r0;                                                                                   +
-                                                                                                              % array of solenoidal vectors                                                                                    +
-                                                                                                              S=[.125*qi*ri*(4-3*q^2-3*r^2), .125*qi*r1*(3*r0-1), -.125*ri*q1*(3*q0-1); ...                                    +
-                                                                                                              .75*qi^2*r1*q0,              0,                   .25*qi*r1*(3*q0+1)];                                           +
-                                                                                                              return;                                                                                                          +
-                                                                                                                                                                                                                               +
-                                                                                                              function s=Syr(ni,q,r)                                                                                           +
-                                                                                                              %SYR  Array of y-derivatives of solenoidal basis functions.                                                      +
-                                                                                                              qi=ni(1); q0=q*ni(1); q1=1+q0;                                                                                   +
-                                                                                                              ri=ni(2); r0=r*ni(2); r1=1+r0;                                                                                   +
-                                                                                                              % array of solenoidal vectors                                                                                    +
-                                                                                                              s=[-.75*ri^2*q1*r0,             .25*ri*q1*(3*r0+1),   0 ; ...                                                    +
-                                                                                                              -.125*qi*ri*(4-3*q^2-3*r^2), .125*qi*r1*(1-3*r0), .125*ri*q1*(3*q0-1)];                                          +
-                                                                                                              return;                                                                                                          +
-                                                                                                                                                                                                                               +
-                                                                                                              function s=Sxxr(ni,q,r)                                                                                          +
-                                                                                                              %SXXR  Array of second x-derivatives of solenoidal basis functions.                                              +
-                                                                                                              qi=ni(1); q0=q*ni(1); q1=1+q0;                                                                                   +
-                                                                                                              ri=ni(2); r0=r*ni(2); r1=1+r0;                                                                                   +
-                                                                                                              % array of solenoidal vectors                                                                                    +
-                                                                                                              s=[-3/4*qi^2*ri*q0, 0, -1/4*ri*qi*(1+3*q0); 3/4*qi^3*r1, 0, 3/4*qi^2*r1 ];                                       +
-                                                                                                              return;                                                                                                          +
-                                                                                                                                                                                                                               +
-                                                                                                              function s=Syyr(ni,q,r)                                                                                          +
-                                                                                                              %SYYR  Array of second y-derivatives of solenoidal basis functions.                                              +
-                                                                                                              qi=ni(1); q0=q*ni(1); q1=1+q0;                                                                                   +
-                                                                                                              ri=ni(2); r0=r*ni(2); r1=1+r0;                                                                                   +
-                                                                                                              % array of solenoidal vectors                                                                                    +
-                                                                                                              s=[-3/4*ri^3*q1, 3/4*ri^2*q1, 0; 3/4*qi*ri^2*r0, -1/4*qi*ri*(1+3*r0), 0 ];                                       +
-                                                                                                              return;                                                                                                          +
-                                                                                                                                                                                                                               +
-                                                                                                              function s=Sxyr(ni,q,r)                                                                                          +
-                                                                                                              %SXYR  Array of second (cross) xy-derivatives of solenoidal basis functions.                                     +
-                                                                                                              qi=ni(1); q0=q*ni(1); q1=1+q0;                                                                                   +
-                                                                                                              ri=ni(2); r0=r*ni(2); r1=1+r0;                                                                                   +
-                                                                                                              % array of solenoidal vectors                                                                                    +
-                                                                                                              s=[-3/4*qi*ri^2*r0, 1/4*qi*ri*(1+3*r0), 0; 3/4*qi^2*ri*q0, 0, 1/4*qi*ri*(1+3*q0)];                               +
-                                                                                                              return;                                                                                                          +
-                                                                                                                                                                                                                               +
-                                                                                                              function G=GBL(ni,q,r)                                                                                           +
-                                                                                                              % Transposed gradient (derivatives) of scalar bilinear mapping function.                                         +
-                                                                                                              % The parameter ni can be a vector of coordinate pairs.                                                          +
-                                                                                                              G=[.25*ni(:,1).*(1+ni(:,2).*r), .25*ni(:,2).*(1+ni(:,1).*q)];                                                    +
-                                                                                                              return;                                                                                                          +
-
+ -
+
-                                                                                                              function Q = ilu_gmres_with_EBC(mat,rhs,EBC,GMRES,Q0,DropTol)                                                    +
-                                                                                                              % ILU_GMRES_WITH_EBCx  Solves matrix equation mat*Q = rhs.                                                       +
-                                                                                                              %                                                                                                                +
-                                                                                                              % Solves the matrix equation mat*Q = rhs, optionally constrained                                                 +
-                                                                                                              % by Dirichlet boundary conditions described in diri_list, using                                                 +
-                                                                                                              % Matlab's preconditioned gmres sparse solver.  When Dirichlet                                                   +
-                                                                                                              % boundary conditions are provided, the routine enforces them by                                                 +
-                                                                                                              % reordering to partition out Dirichlet degrees of freedom.                                                      +
-                                                                                                              % usage:  Q = ilu_gmres_with_EBC(mat,rhs,EBC,GMRES,Q0,DropTol)                                                   +
-                                                                                                              %    or:  Q = ilu_gmres_with_EBC(mat,rhs,EBC,GMRES,Q0)                                                           +
-                                                                                                              %    or:  Q = ilu_gmres_with_EBC(mat,rhs,EBC,GMRES)                                                              +
-                                                                                                              %    or:  Q = ilu_gmres_with_EBC(mat,rhs,EBC)                                                                    +
-                                                                                                              %    or:  Q = ilu_gmres_with_EBC(mat,rhs)                                                                        +
-                                                                                                              %                                                                                                                +
-                                                                                                              %  mat  - matrix of linear system to be solved.                                                                  +
-                                                                                                              %  rhs  - right hand side of linear system.                                                                      +
-                                                                                                              %  EBC.dof, EBC.val - (optional) list of Dirichlet boundary                                                      +
-                                                                                                              %          conditions (constraints). May be empty ([]).                                                          +
-                                                                                                              %  GMRES - Structure specifying tolerance, max iterations and                                                    +
-                                                                                                              %          restarts. Use [] for default values.                                                                  +
-                                                                                                              %  Q0   - (optional) initial approximation to solution for restart,                                              +
-                                                                                                              %          may be empty ([]).                                                                                    +
-                                                                                                              %  DropTol - drop tolerance for luinc preconditioner (default=1e-6).                                             +
-                                                                                                              %                                                                                                                +
-                                                                                                              % The solution Q is reported with the original ordering restored.                                                +
-                                                                                                              %                                                                                                                +
-                                                                                                              % If specified and not empty, diri_list should have two columns                                                  +
-                                                                                                              %   total and one row for each diri degree of freedom.  The first                                                +
-                                                                                                              %   column of each row must contain global index of the degree of                                                +
-                                                                                                              %   freedom.  The second column contains the actual Dirichlet value.                                             +
-                                                                                                              %   If there are no Dirichlet nodes, omit this parameter if not using                                            +
-                                                                                                              %   the restart capabilities, or supply empty matrix ([]) as diri_list.                                          +
-                                                                                                              %                                                                                                                +
-                                                                                                              % Nominal values for GMRES for a large difficult problem might be:                                               +
-                                                                                                              %  GMRES.Tolerance  = 1.e-12,                                                                                    +
-                                                                                                              %  GMRES.MaxIterates = 75,                                                                                       +
-                                                                                                              %  GMRES.MaxRestarts = 14.                                                                                       +
-                                                                                                              %                                                                                                                +
-                                                                                                              % Jonas Holdeman,   revised February, 2009.                                                                      +
-                                                                                                                                                                                                                               +
-                                                                                                              % Drop tolerance for luinc preconditioner, nominal value - 1.e-6                                                 +
-                                                                                                              if (nargin<6 | isempty(DropTol) | DropTol<=0)                                                                    +
-                                                                                                              droptol = 1.e-6;        % default                                                                                +
-                                                                                                              else droptol = DropTol;   % assigned                                                                             +
-                                                                                                              end                                                                                                              +
-                                                                                                                                                                                                                               +
-                                                                                                              if nargin<=3 | isempty(GMRES)                                                                                    +
-                                                                                                              Tol = 1.e-12;  MaxIter = 75;  MaxRstrt = 14;                                                                     +
-                                                                                                              else                                                                                                             +
-                                                                                                              % Tol=tolerance for residual, increase it if solution takes too long.                                            +
-                                                                                                              Tol=GMRES.Tolerance;                                                                                             +
-                                                                                                              % MaxIter = maximum number of iterations before restart (MT used 5)                                              +
-                                                                                                              MaxIter=GMRES.MaxIterates;                                                                                       +
-                                                                                                              % MaxRstrt = maximum number of restarts before giving up (MT used 10)                                            +
-                                                                                                              MaxRstrt=GMRES.MaxRestarts;                                                                                      +
-                                                                                                              end                                                                                                              +
-                                                                                                                                                                                                                               +
-                                                                                                              % check arguments for reasonableness                                                                             +
-                                                                                                              tdof = size(mat,1);	 % total number of degrees of freedom                                                        +
-                                                                                                                                                                                                                               +
-                                                                                                              % good mat is a square tdof by tdof matrix                                                                       +
-                                                                                                              if (size(mat,2)~=size(mat,1) | size(size(mat),2)~=2)                                                             +
-                                                                                                              error('mat must be a square matrix')                                                                             +
-                                                                                                              end                                                                                                              +
-                                                                                                                                                                                                                               +
-                                                                                                              % valid rhs has the dimensions [tdof, 1]                                                                         +
-                                                                                                              if (size(rhs,1)~=tdof | size(rhs,2)~=1)                                                                          +
-                                                                                                              error('rhs must be a column matrix with the same number of rows as mat')                                         +
-                                                                                                              end                                                                                                              +
-                                                                                                                                                                                                                               +
-                                                                                                              % valid dimensions for optional diri_list                                                                        +
-                                                                                                              if nargin<=2                                                                                                     +
-                                                                                                              EBC=[];                                                                                                          +
-                                                                                                              elseif ~isempty(EBC) & (size(EBC.val,1)>=tdof ...                                                                +
-                                                                                                              | size(EBC.dof,1)>=tdof)                                                                                         +
-                                                                                                              error('check dimensions of EBC')                                                                                 +
-                                                                                                              end                                                                                                              +
-                                                                                                              % (optional) valid Q0 is empty or has the dimensions [tdof, 1]                                                   +
-                                                                                                              if nargin<=4                                                                                                     +
-                                                                                                              Q0=[];                                                                                                           +
-                                                                                                              elseif ~isempty(Q0) & (size(Q0,1)~=tdof | size(rhs,2)~=1)                                                        +
-                                                                                                              error('Q0 must be a column matrix with the same number of rows as mat')                                          +
-                                                                                                              end                                                                                                              +
-                                                                                                                                                                                                                               +
-                                                                                                              % handle the case of no Dirichlet dofs separately                                                                +
-                                                                                                              if isempty(EBC)                                                                                                  +
-                                                                                                              % skip diri partitioning, solve the system                                                                       +
-                                                                                                              [L,U] = luinc(mat,droptol);      % incomplete LU                                                                 +
-                                                                                                              Q = gmres(mat,rhs,MaxIter,Tol,MaxRstrt,L,U,Q0);	% GMRES                                                          +
-                                                                                                                                                                                                                               +
-                                                                                                              else                                                                                                             +
-                                                                                                              % Form list of all DOFs                                                                                          +
-                                                                                                              p_vec = [1:tdof]';                                                                                               +
-                                                                                                              % partion out diri dofs                                                                                          +
-                                                                                                              EBCdofs = EBC.dof(:,1);	 % list of dofs which are Dirichlet                                                      +
-                                                                                                              EBCvals = EBC.val(:,1);  % Dirichlet dof values                                                                  +
-                                                                                                                                                                                                                               +
-                                                                                                              % form a list of non-diri dofs                                                                                   +
-                                                                                                              ndro = p_vec(~ismember(p_vec, EBCdofs));	% list of non-diri dofs                                                 +
-                                                                                                                                                                                                                               +
-                                                                                                              % Move Dirichlet DOFs to right side                                                                              +
-                                                                                                              rhs_reduced = rhs(ndro) - mat(ndro, EBCdofs) * EBCvals;                                                          +
-                                                                                                                                                                                                                               +
-                                                                                                              % solve the reduced system (preconditioned gmres)                                                                +
-                                                                                                              A = mat(ndro,ndro);                                                                                              +
-                                                                                                                                                                                                                               +
-                                                                                                              % Compute incomplete LU preconditioner                                                                           +
-                                                                                                              [L,U] = luinc(A,droptol);      % incomplete LU                                                                   +
-                                                                                                                                                                                                                               +
-                                                                                                              % Remove Dirichlet DOFs from initial estimate                                                                    +
-                                                                                                              if ~isempty(Q0)  Q0=Q0(ndro);  end                                                                               +
-                                                                                                                                                                                                                               +
-                                                                                                              % solve the reduced system (preconditioned gmres)                                                                +
-                                                                                                              Q_reduced = gmres(A,rhs_reduced,MaxIter,Tol,MaxRstrt,L,U,Q0);                                                    +
-                                                                                                                                                                                                                               +
-                                                                                                              % insert the Dirichlet values into the solution                                                                  +
-                                                                                                              Q = [Q_reduced; EBCvals];                                                                                        +
-                                                                                                                                                                                                                               +
-                                                                                                              % calculate p_vec_undo to restore Q to the original dof ordering                                                 +
-                                                                                                              p_vec_undo = zeros(1,tdof);                                                                                      +
-                                                                                                              p_vec_undo([ndro;EBCdofs]) = [1:tdof];                                                                           +
-                                                                                                                                                                                                                               +
-                                                                                                              % restore the original ordering                                                                                  +
-                                                                                                              Q = Q(p_vec_undo);                                                                                               +
-                                                                                                              end                                                                                                              +
-
+ - + - File '''regrade.m''' for generating non-uniform mesh spacing. + - + -
+
-                                                                                                              %REGRADE  grade nodal points in array towards edge or center                                                     +
-                                                                                                              %                                                                                                                +
-                                                                                                              % Regrades array of points                                                                                       +
-                                                                                                              %                                                                                                                +
-                                                                                                              % x is array of nodal point coordinates in increasing order.                                                     +
-                                                                                                              % a is parameter which controls grading.                                                                         +
-                                                                                                              % e selects side or sides for refinement.                                                                        +
-                                                                                                              %                                                                                                                +
-                                                                                                              % if e=0: refine both sides, 1: refine upper, 2: refine lower.                                                   +
-                                                                                                              %                                                                                                                +
-                                                                                                              % if a=1 then return xarray unaltered.                                                                           +
-                                                                                                              % if a<1 then grade towards the edge(s)                                                                          +
-                                                                                                              % if a>1 then grade away from edge.                                                                              +
-                                                                                                                                                                                                                               +
-                                                                                                              ae=abs(a);                                                                                                       +
-                                                                                                              n=length(x);                                                                                                     +
-                                                                                                              y=x;                                                                                                             +
-                                                                                                              if ae==1 | n<3 | (e~=0 & e~=1 & e~=2) return;                                                                    +
-                                                                                                              end;                                                                                                             +
-                                                                                                              if e==0                                                                                                          +
-                                                                                                              Xmx=max(x);                                                                                                      +
-                                                                                                              Xmn=min(x);                                                                                                      +
-                                                                                                              Xc=(Xmx+Xmn)/2;                                                                                                  +
-                                                                                                              Xl=(Xmx-Xmn)/2;                                                                                                  +
-                                                                                                              for k=2:(n-1)                                                                                                    +
-                                                                                                              xk=x(k)-Xc;                                                                                                      +
-                                                                                                              y(k)=Xc+Xl*sign(xk)*(abs(xk/Xl))^ae;                                                                             +
-                                                                                                              end                                                                                                              +
-                                                                                                              elseif (e==1 & x(1)x(n))                                                                   +
-                                                                                                              for k=2:n-1                                                                                                      +
-                                                                                                              xk=x(k)-x(1);                                                                                                    +
-                                                                                                              y(k)=x(1)+(x(n)-x(1))*(abs(xk/(x(n)-x(1))))^ae;                                                                  +
-                                                                                                              end                                                                                                              +
-                                                                                                              else   % (e==2 & x(1)x(n))                                                                 +
-                                                                                                              for k=2:n-1                                                                                                      +
-                                                                                                              xk=x(k)-x(n);                                                                                                    +
-                                                                                                              y(k)=x(n)+(x(1)-x(n))*(abs(xk/(x(1)-x(n))))^ae;                                                                  +
-                                                                                                              end                                                                                                              +
-                                                                                                              end                                                                                                              +
-                                                                                                              return;                                                                                                          +
-
+

## Revision as of 04:19, 29 June 2011

Lid-driven cavity using pressure-free velocity formulation

```%LDCW            LID-DRIVEN CAVITY
% Finite element solution of the 2D Navier-Stokes equation using 4-node, 12 DOF,
%  (3-DOF/node), simple-cubic-derived rectangular Hermite basis for
%   the Lid-Driven Cavity problem.
%
% This could also be characterized as a VELOCITY-STREAM FUNCTION or
%   STREAM FUNCTION-VELOCITY method.
%
% Reference:  "A Hermite finite element method for incompressible fluid flow",
%    Int. J. Numer. Meth. Fluids, 64, P376-408 (2010).
%
% Simplified Wiki version
% The rectangular problem domain is defined between Cartesian
%   coordinates Xmin & Xmax and Ymin & Ymax.
% The computational grid has NumEx elements in the x-direction
%   and NumEy elements in the y-direction.
% The nodes and elements are numbered column-wise from the
%   upper left corner to the lower right corner.
%
%This script calls the user-defined functions:
% DMatW        - to evaluate element diffusion matrix
% CMatW        - to evaluate element convection matrix
% GetPresW     - to evaluate the pressure
% ilu_gmres_with_EBC - to solve the system with essential/Dirichlet BCs
%
% Jonas Holdeman   August 2007, revised June 2011

clear all;
disp('Lid-driven cavity');
disp(' Four-node, 12 DOF, simple-cubic stream function basis.');

% -------------------------------------------------------------
nd = 3; nd2=nd*nd;  % Number of DOF per node - do not change!!
% -------------------------------------------------------------
ETstart=clock;

% Parameters for GMRES solver
GMRES.Tolerance=1.e-14;
GMRES.MaxIterates=20;
GMRES.MaxRestarts=6;

% Optimal relaxation parameters for given Reynolds number
% (see IJNMF reference)
% Re          100   1000   3200   5000   7500  10000  12500
% RelxFac:  1.04    1.11   .860   .830   .780   .778   .730
% ExpCR1    1.488   .524   .192   .0378   --     --     --
% ExpCRO    1.624   .596   .390   .331   .243   .163   .133
% CritFac:  1.82    1.49   1.14  1.027   .942   .877   .804

% Define the problem geometry, set mesh bounds:
Xmin = 0.0; Xmax = 1.0; Ymin = 0.0; Ymax = 1.0;

% See below for explanation of use of parameters.
xgrd = .75; ygrd=.75;   % (xgrd = 1, ygrd=1 for uniform mesh)

% Set " RefineBoundary=1 " for additional refinement at boundary,
%  i.e., split first element along boundary into two.
RefineBoundary=1;

%     DEFINE THE MESH
% Set number of elements in each direction
NumEx = 16;   NumEy = NumEx;

% PLEASE CHANGE OR SET NUMBER OF ELEMENTS TO CHANGE/SET NUMBER OF NODES!
NumNx=NumEx+1;  NumNy=NumEy+1;

%   Define problem parameters:
% Lid velocity
Vlid=1.;

% Reynolds number
Re=1000.;

% factor for under/over-relaxation starting at iteration RelxStrt
RelxFac = 1.;  %

% Number of nonlinear iterations
MaxNLit=10; %

%--------------------------------------------------------

% Viscosity for specified Reynolds number
nu=Vlid*(Xmax-Xmin)/Re;

% if e=0: refine both sides, 1: refine upper, 2: refine lower
% if agrd=xgrd|ygrd is the parameter which controls grading, then
%   if agrd=1 then leave array unaltered.
%   if agrd<1 then refine (make finer) towards the ends
%   if agrd>1 then refine (make finer) towards the center.
%
%  Generate equally-spaced nodal coordinates and refine if desired.
if (RefineBoundary==1)
XNc=linspace(Xmin,Xmax,NumNx-2);
XNc=[XNc(1),(.62*XNc(1)+.38*XNc(2)),XNc(2:end-1),(.38*XNc(end-1)+.62*XNc(end)),XNc(end)];
YNc=linspace(Ymax,Ymin,NumNy-2);
YNc=[YNc(1),(.62*YNc(1)+.38*YNc(2)),YNc(2:end-1),(.38*YNc(end-1)+.62*YNc(end)),YNc(end)];
else
XNc=linspace(Xmin,Xmax,NumNx);
YNc=linspace(Ymax,Ymin,NumNy);
end
if xgrd ~= 1 XNc=regrade(XNc,xgrd,0); end;  % Refine mesh if desired
if ygrd ~= 1 YNc=regrade(YNc,ygrd,0); end;
[Xgrid,Ygrid]=meshgrid(XNc,YNc);% Generate the x- and y-coordinate meshes.

% Allocate storage for fields
psi0=zeros(NumNy,NumNx);
u0=zeros(NumNy,NumNx);
v0=zeros(NumNy,NumNx);

%--------------------Begin grid plot-----------------------
% ********************** FIGURE 1 *************************
% Plot the grid
figure(1);
clf;
orient portrait;  orient tall;
subplot(2,2,1);
hold on;
plot([Xmax;Xmin],[YNc;YNc],'k');
plot([XNc;XNc],[Ymax;Ymin],'k');
hold off;
axis([Xmin,Xmax,Ymin,Ymax]);
axis equal;
axis image;
title([num2str(NumNx) 'x' num2str(NumNy) ...
' node mesh for Lid-driven cavity']);
pause(.1);
%-------------- End plotting Figure 1 ----------------------

%Contour levels, Ghia, Ghia & Shin, Re=100, 400, 1000, 3200, ...
clGGS=[-.1175,-.1150,-.11,-.1,-.09,-.07,-.05,-.03,-.01,-1.e-4,-1.e-5,-1.e-7,-1.e-10,...
1.e-8,1.e-7,1.e-6,1.e-5,5.e-5,1.e-4,2.5e-4,5.e-4,1.e-3,1.5e-3,3.e-3];
CL=clGGS;   % Select contour level option
if (Vlid<0) CL=-CL; end

NumNod=NumNx*NumNy;     % total number of nodes
MaxDof=nd*NumNod;        % maximum number of degrees of freedom
EBC.Mxdof=nd*NumNod;        % maximum number of degrees of freedom

nn2nft=zeros(2,NumNod); % node number -> nf & nt
NodNdx=zeros(2,NumNod);
% Generate lists of active nodal indices, freedom number & type
ni=0;  nf=-nd+1;  nt=1;          %   ________
for nx=1:NumNx                   %  |        |
for ny=1:NumNy                %  |        |
ni=ni+1;                   %  |________|
NodNdx(:,ni)=[nx;ny];
nf=nf+nd;               % all nodes have 4 dofs
nn2nft(:,ni)=[nf;nt];   % dof number & type (all nodes type 1)
end;
end;
%NumNod=ni;     % total number of nodes
nf2nnt=zeros(2,MaxDof);  % (node, type) associated with dof
ndof=0; dd=[1:nd];
for n=1:NumNod
for k=1:nd
nf2nnt(:,ndof+k)=[n;k];
end
ndof=ndof+nd;
end

NumEl=NumEx*NumEy;

% Generate element connectivity, from upper left to lower right.
Elcon=zeros(4,NumEl);
ne=0;  LY=NumNy;
for nx=1:NumEx
for ny=1:NumEy
ne=ne+1;
Elcon(1,ne)=1+ny+(nx-1)*LY;
Elcon(2,ne)=1+ny+nx*LY;
Elcon(3,ne)=1+(ny-1)+nx*LY;
Elcon(4,ne)=1+(ny-1)+(nx-1)*LY;
end  % loop on ny
end  % loop on nx

% Begin essential boundary conditions, allocate space
MaxEBC = nd*2*(NumNx+NumNy-2);
EBC.dof=zeros(MaxEBC,1);  % Degree-of-freedom index
EBC.typ=zeros(MaxEBC,1);  % Dof type (1,2,3)
EBC.val=zeros(MaxEBC,1);  % Dof value

X1=XNc(2);  X2=XNc(NumNx-1);
nc=0;
for nf=1:MaxDof
ni=nf2nnt(1,nf);
nx=NodNdx(1,ni);
ny=NodNdx(2,ni);
x=XNc(nx);
y=YNc(ny);
if(x==Xmin | x==Xmax | y==Ymin)
nt=nf2nnt(2,nf);
switch nt;
case {1, 2, 3}
nc=nc+1; EBC.typ(nc)=nt; EBC.dof(nc)=nf; EBC.val(nc)=0;  % psi, u, v
end  % switch (type)
elseif (y==Ymax)
nt=nf2nnt(2,nf);
switch nt;
case {1, 3}
nc=nc+1; EBC.typ(nc)=nt; EBC.dof(nc)=nf; EBC.val(nc)=0;   % psi, v
case 2
nc=nc+1; EBC.typ(nc)=nt; EBC.dof(nc)=nf; EBC.val(nc)=Vlid;   % u
end  % switch (type)
end  % if (boundary)
end  % for nf
EBC.num=nc;

if (size(EBC.typ,1)>nc)   % Truncate arrays if necessary
EBC.typ=EBC.typ(1:nc);
EBC.dof=EBC.dof(1:nc);
EBC.val=EBC.val(1:nc);
end     % End ESSENTIAL (Dirichlet) boundary conditions

% partion out essential (Dirichlet) dofs
p_vec = [1:EBC.Mxdof]';         % List of all dofs
EBC.p_vec_undo = zeros(1,EBC.Mxdof);
% form a list of non-diri dofs
EBC.ndro = p_vec(~ismember(p_vec, EBC.dof));	% list of non-diri dofs
% calculate p_vec_undo to restore Q to the original dof ordering
EBC.p_vec_undo([EBC.ndro;EBC.dof]) = [1:EBC.Mxdof]; %p_vec';

Q=zeros(MaxDof,1); % Allocate space for solution (dof) vector

% Initialize fields to boundary conditions
for k=1:EBC.num
Q(EBC.dof(k))=EBC.val(k);
end;

errpsi=zeros(NumNy,NumNx);  % error correct for iteration

MxNL=max(1,MaxNLit);
np0=zeros(1,MxNL);     % Arrays for convergence info
nv0=zeros(1,MxNL);

Qs=[];

Dmat = spalloc(MaxDof,MaxDof,36*MaxDof);   % to save the diffusion matrix
Vdof=zeros(nd,4);
Xe=zeros(2,4);      % coordinates of element corners

NLitr=0; ND=1:nd;
while (NLitr<MaxNLit), NLitr=NLitr+1;   % <<< BEGIN NONLINEAR ITERATION

tclock=clock;   % Start assembly time <<<<<<<<<
% Generate and assemble element matrices
Mat=spalloc(MaxDof,MaxDof,36*MaxDof);
RHS=spalloc(MaxDof,1,MaxDof);
%RHS = zeros(MaxDof,1);
Emat=zeros(1,16*nd2);         % Values 144=4*4*(nd*nd)

% BEGIN GLOBAL MATRIX ASSEMBLY
for ne=1:NumEl

for k=1:4
ki=NodNdx(:,Elcon(k,ne));
Xe(:,k)=[XNc(ki(1));YNc(ki(2))];
end

if NLitr == 1
%     Fluid element diffusion matrix, save on first iteration
[DEmat,Rndx,Cndx] = DMatW(Xe,Elcon(:,ne),nn2nft);
Dmat=Dmat+sparse(Rndx,Cndx,DEmat,MaxDof,MaxDof);  % Global diffusion mat
end

if (NLitr>1)
%    Get stream function and velocities
for n=1:4
Vdof(ND,n)=Q((nn2nft(1,Elcon(n,ne))-1)+ND); % Loop over local element nodes
end
%     Fluid element convection matrix, first iteration uses Stokes equation.
[Emat,Rndx,Cndx] = CMatW(Xe,Elcon(:,ne),nn2nft,Vdof);
Mat=Mat+sparse(Rndx,Cndx,-Emat,MaxDof,MaxDof);  % Global convection assembly
end

end;  % loop ne over elements
% END GLOBAL MATRIX ASSEMBLY

Mat = Mat -nu*Dmat;    % Add in cached/saved global diffusion matrix

disp(['(' num2str(NLitr) ') Matrix assembly complete, elapsed time = '...
num2str(etime(clock,tclock)) ' sec']);  % Assembly time <<<<<<<<<<<
pause(1);

Q0 = Q;  % Save dof values

% Solve system
tclock=clock; %disp('start solution'); % Start solution time  <<<<<<<<<<<<<<

RHSr=RHS(EBC.ndro)-Mat(EBC.ndro,EBC.dof)*EBC.val;
Matr=Mat(EBC.ndro,EBC.ndro);
Qs=Q(EBC.ndro);

Qr=ilu_gmres_with_EBC(Matr,RHSr,[],GMRES,Qs);

Q=[Qr;EBC.val];        % Augment active dofs with esential (Dirichlet) dofs
Q=Q(EBC.p_vec_undo);   % Restore natural order

stime=etime(clock,tclock); % Solution time <<<<<<<<<<<<<<

% ****** APPLY RELAXATION FACTOR *********************
if(NLitr>1) Q=RelxFac*Q+(1-RelxFac)*Q0; end
% ****************************************************

% Compute change and copy dofs to field arrays
dsqp=0; dsqv=0;
for k=1:MaxDof
ni=nf2nnt(1,k); nx=NodNdx(1,ni); ny=NodNdx(2,ni);
switch nf2nnt(2,k) % switch on dof type
case 1
dsqp=dsqp+(Q(k)-Q0(k))^2; psi0(ny,nx)=Q(k);
errpsi(ny,nx)=Q0(k)-Q(k);
case 2
dsqv=dsqv+(Q(k)-Q0(k))^2; u0(ny,nx)=Q(k);
case 3
dsqv=dsqv+(Q(k)-Q0(k))^2; v0(ny,nx)=Q(k);
end  % switch on dof type
end  % for
np0(NLitr)=sqrt(dsqp);
nv0(NLitr)=sqrt(dsqv);

if (np0(NLitr)<=1e-15|nv0(NLitr)<=1e-15)
MaxNLit=NLitr; np0=np0(1:MaxNLit); nv0=nv0(1:MaxNLit);   end;
disp(['(' num2str(NLitr) ') Solution time for linear system = '...
num2str(etime(clock,tclock)) ' sec']); % Solution time <<<<<<<<<<<<

%---------- Begin plot of intermediate results ----------
% ********************** FIGURE 2 *************************
figure(1);

% Stream function (intermediate)
subplot(2,2,3);
contour(Xgrid,Ygrid,psi0,8,'k');  % Plot contours (trajectories)
axis([Xmin,Xmax,Ymin,Ymax]);
title(['Lid-driven cavity,  Re=' num2str(Re)]);
axis equal; axis image;

% Plot convergence info
subplot(2,2,2);
semilogy(1:NLitr,nv0(1:NLitr),'k-+',1:NLitr,np0(1:NLitr),'k-o');
xlabel('Nonlinear iteration number');
ylabel('Nonlinear correction');
axis square;
title(['Iteration conv.,  Re=' num2str(Re)]);
legend('U','Psi');

% Plot nonlinear iteration correction contours
subplot(2,2,4);
contour(Xgrid,Ygrid,errpsi,8,'k');  % Plot contours (trajectories)
axis([Xmin,Xmax,Ymin,Ymax]);
axis equal; axis image;
title(['Iteration correction']);
pause(1);
% ********************** END FIGURE 2 *************************
%----------  End plot of intermediate results  ---------

if (nv0(NLitr)<1e-15) break; end  % Terminate iteration if non-significant

end;   % <<< (while) END NONLINEAR ITERATION

format short g;
disp('Convergence results by iteration: velocity, stream function');
disp(['nv0:  ' num2str(nv0)]); disp(['np0:  ' num2str(np0)]);

% >>>>>>>>>>>>>> BEGIN PRESSURE RECOVERY <<<<<<<<<<<<<<<<<<
% Essential pressure boundary condition
% Index of node to apply pressure BC, value at node
PBCnx=fix((NumNx+1)/2);   % Apply at center of mesh
PBCny=fix((NumNy+1)/2);
PBCnod=0;
for k=1:NumNod
if (NodNdx(1,k)==PBCnx & NodNdx(2,k)==PBCny) PBCnod=k; break; end
end
else
EBCp.nodn = [PBCnod];  % Pressure BC node number
EBCp.val = [0];  % set P = 0.
end
% Cubic pressure
[P,Px,Py] = GetPresW(NumNod,NodNdx,Elcon,nn2nft,Xgrid,Ygrid,Q,EBCp,nu);
% ******************** END PRESSURE RECOVERY *********************

% ********************** CONTINUE FIGURE 1 *************************
figure(1);

% Stream function    (final)
subplot(2,2,3);
[CT,hn]=contour(Xgrid,Ygrid,psi0,CL,'k');  % Plot contours (trajectories)
clabel(CT,hn,CL([1,3,5,7,9,10,11,19,23]));
hold on;
plot([Xmin,Xmin,Xmax,Xmax,Xmin],[Ymax,Ymin,Ymin,Ymax,Ymax],'k');
hold off;
axis([Xmin,Xmax,Ymin,Ymax]);
axis equal;  axis image;
title(['Stream lines, ' num2str(NumNx) 'x' num2str(NumNy) ...
' mesh, Re=' num2str(Re)]);

% Plot pressure contours   (final)
subplot(2,2,4);
CPL=[-.002,0,.02,.05,.07,.09,.11,.12,.17,.3];
[CT,hn]=contour(Xgrid,Ygrid,P,CPL,'k');  % Plot pressure contours
clabel(CT,hn,CPL([3,5,7,10]));
hold on;
plot([Xmin,Xmin,Xmax,Xmax,Xmin],[Ymax,Ymin,Ymin,Ymax,Ymax],'k');
hold off;
axis([Xmin,Xmax,Ymin,Ymax]);
axis equal;  axis image;
title(['Simple cubic pressure contours, Re=' num2str(Re)]);
% ********************* END FIGURE 1 *************************

disp(['Total elapsed time = '...
num2str(etime(clock,ETstart)/60) ' min']); % Elapsed time from start <<<
```