CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Wiki > PFV get pressure 2

PFV get pressure 2

From CFD-Wiki

Revision as of 18:42, 4 July 2011 by Jonas Holdeman (Talk | contribs)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search

Matlab function GetPres4424W.m to retrieve consistent pressure from velocity.

function [P,Px,Py] = GetPres44243W(NumNod,NodNdx,Elcon,nn2nft,Xgrid,Ygrid,Q,EBCPr,nu)
%GETPRES44243G - Compute continuous simple cubic pressure from (quartic-based) 
%  velocity field on general quadrilateral grid (bilinear geometric mapping). 
%
% Usage
%   [P,PX,PY]=Getpres44243W(NumNod,NodNdx,Elcon,nn2nft,Xgrid,Ygrid,Q,EBDp,nu);
% 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 July 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 = 6;  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
% Set IsAffine=1 for affine elements, IsAffine=0 for general quad
  IsAffine=0;
% -------------- 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,36*MxPDof);   % allocate P mass matrix
Vdata = zeros(MxPDof,1);        % allocate velocity data 
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:nnd
     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)))';    % Jacobian at (xa,ya)
  Det=J(1,1)*J(2,2)-J(1,2)*J(2,1);  % Determinant of Jt & J 
  Ji=[J(2,2),-J(1,2); -J(2,1),J(1,1)]/Det;
  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 = 6;  nd4=4*nd;  ND=1:nd;
%
global ZS4424pc; global ZSX4424pc; global ZSY4424pc; global ZG4424pc;  
if (isempty(ZS4424pc)|size(ZS4424pc,2)~=Nq)
% Evaluate and save/cache the set of shape functions at quadrature pts. 
  ZS4424pc=cell(nnd,Nq); ZSX4424pc=cell(nnd,Nq); 
  ZSY4424pc=cell(nnd,Nq); ZG4424pc=cell(nnd,Nq);  
  for k=1:Nq
    for m=1:nnd
      ZS4424pc{m,k} =Sr(nn(m,:),xa(k),ya(k));
      ZSX4424pc{m,k}=Sxr(nn(m,:),xa(k),ya(k));
      ZSY4424pc{m,k}=Syr(nn(m,:),xa(k),ya(k));
      ZG4424pc{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 
  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)
  T2=[Jt(1,1)^2, 2*Jt(1,1)*Jt(1,2), Jt(1,2)^2; ...  % alt
    Jt(1,1)*Jt(2,1), Jt(1,1)*Jt(2,2)+Jt(1,2)*Jt(2,1), Jt(1,2)*Jt(2,2); ...
    Jt(2,1)^2, 2*Jt(2,1)*Jt(2,2), Jt(2,2)^2];
  T6=blkdiag(1,JtiD,T2); 
  Bxy=Xe*BLxy(nn(:,:),nn(m,1),nn(m,2));  % Second cross derivatives
  T6(5,2:3)=Bxy([2,1])';
  Ti{m}=T6;
  TGi{m} = blkdiag(1,Jt');
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 
  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);
    S(:,mm+ND) =Jtd*ZS4424pc{m,k}*Ti{m};
    Sx(:,mm+ND)=Jtd*(Ji(1,1)*ZSX4424pc{m,k}+Ji(2,1)*ZSY4424pc{m,k})*Ti{m};
    Sy(:,mm+ND)=Jtd*(Ji(1,2)*ZSX4424pc{m,k}+Ji(2,2)*ZSY4424pc{m,k})*Ti{m};
     mm=np*(m-1);
    G(:,mm+NP)=Ji*ZG4424pc{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 = 6;  nd4=nnd*nd;  ND=1:nd;
global ZSXX4424pd; global ZSXY4424pd; global ZSYY4424pd; global ZG4424pd;  
if (isempty(ZSXX4424pd)|size(ZSXX4424pd,2)~=Nq)
% Evaluate and save/cache the set of shape functions at quadrature pts. 
ZSXX4424pd=cell(nnd,Nq); ZSXY4424pd=cell(nnd,Nq); 
ZSYY4424pd=cell(nnd,Nq); ZG4424pd=cell(nnd,Nq);  
  for k=1:Nq
    for m=1:nnd
      ZSXX4424pd{m,k}=Sxxr(nn(m,:),xa(k),ya(k));
      ZSXY4424pd{m,k}=Sxyr(nn(m,:),xa(k),ya(k));
      ZSYY4424pd{m,k}=Syyr(nn(m,:),xa(k),ya(k));
      ZG4424pd{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 
  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)
  T2=[Jt(1,1)^2, 2*Jt(1,1)*Jt(1,2), Jt(1,2)^2; ...  % alt
    Jt(1,1)*Jt(2,1), Jt(1,1)*Jt(2,2)+Jt(1,2)*Jt(2,1), Jt(1,2)*Jt(2,2); ...
    Jt(2,1)^2, 2*Jt(2,1)*Jt(2,2), Jt(2,2)^2];
  T6=blkdiag(1,JtiD,T2); 
  Bxy=Xe*BLxy(nn(:,:),nn(m,1),nn(m,2));  % Second cross derivatives
  T6(5,2:3)=Bxy([2,1])';
  Ti{m}=T6;
  TGi{m} = blkdiag(1,Jt');
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)=Ji(1,1)^2*Jtd*ZSXX4424pd{m,k}*Ti{m};
%    Syy(:,mm+ND)=Ji(2,2)^2*Jtd*ZSYY4424pd{m,k}*Ti{m};
    Sxx(:,mm+ND)=Jtd*(Ji(1,1)^2*ZSXX4424pd{m,k}+2*Ji(1,1)*Ji(1,2)*ZSXY4424pd{m,k}+Ji(1,2)^2*ZSYY4424pd{m,k})*Ti{m};
    Syy(:,mm+ND)=Jtd*(Ji(2,1)^2*ZSXX4424pd{m,k}+2*Ji(2,1)*Ji(2,2)*ZSXY4424pd{m,k}+Ji(2,2)^2*ZSYY4424pd{m,k})*Ti{m};
    mm=np*(m-1);
    G(:,mm+NP) =Ji*ZG4424pd{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)
%GR  Gradient of cubic Hermite basis functions for pressuret.
   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 - x-derivative of cubic Hermite gradient basis functions for pressure.
   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 - y-derivative of cubic Hermite gradient basis functions for pressuret.
   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 (quartic).
   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/64*ri*(1+q0)^2*(1-r^2)*((1+q0)*(8-9*q0+3*q^2)+(2-q0)*(1-5*r^2)), ... 
   -1/64*(1+q0)^2*(1+r0)^2*(2-q0)*(7-26*r0+15*r^2), ...
    3/64*ri/qi*(1+q0)^3*(1-r^2)*(1-q0)*(5-3*q0), ...
    3/64*ri/qi^2*(1+q0)^3*(1-r^2)*(1-q0)^2, ...
    1/16/qi*(1+q0)^2*(1+r0)*(1-q0)*(1-3*r0), ... 
    1/64/ri*(1+q0)^2*(1+r0)^2*(1-r0)*(1-5*r0)*(2-q0); ... 
   -3/64*qi*(1+r0)^2*(1-q^2)*((1+r0)*(8-9*r0+3*r^2)+(2-r0)*(1-5*q^2)), ... 
    3/64*qi/ri*(1+r0)^3*(1-q^2)*(1-r0)*(5-3*r0), ... 
   -1/64*(1+r0)^2*(1+q0)^2*(2-r0)*(7-26*q0+15*q^2), ... 
   -1/64/qi*(1+r0)^2*(1+q0)^2*(1-q0)*(1-5*q0)*(2-r0), ... 
   -1/16/ri*(1+r0)^2*(1+q0)*(1-r0)*(1-3*q0), ... 
   -3/64*qi/ri^2*(1+r0)^3*(1-q^2)*(1-r0)^2];
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=[9/64*qi*ri*(1-q^2)*(1-r^2)*(6-5*q^2-5*r^2), -3/64*qi*(1-q^2)*(1+r0)^2*(1-3*r0)*(7-5*r0), ...
    3/64*ri*(1+q0)^2*(1-r^2)*(1-3*q0)*(7-5*q0), 3/64*ri/qi*(1+q0)^2*(1-r^2)*(1-q0)*(1-5*q0), ...
    1/16*(1+q0)*(1+r0)*(1-3*q0)*(1-3*r0), 3/64*qi/ri*(1-q^2)*(1+r0)^2*(1-r0)*(1-5*r0); ...
    3/32*qi^2*q0*(1+r0)^2*(r0*(10*q^2+3*r^2-7)+20-20*q^2-6*r^2), -3/32*qi^2/ri*q0*(1+r0)^3*(1-r0)*(5-3*r0),...
    3/16*qi*(1-q^2)*(1+r0)^2*(2-r0)*(1+5*q0), 1/16*(1+q0)*(1+r0)^2*(2-r0)*(1+2*q0-5*q^2), ...
    1/8*qi/ri*(1+r0)^2*(1-r0)*(1+3*q0), 3/32*qi^2/ri^2*q0*(1+r0)^3*(1-r0)^2];
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=[-3/32*ri^2*r0*(1+q0)^2*(q0*(10*r^2+3*q^2-7)+20-20*r^2-6*q^2), 3/16*ri*(1+q0)^2*(1-r^2)*(2-q0)*(1+5*r0),...
      -3/32*ri^2*r0/qi*(1+q0)^2*(1-q^2)*(5-3*q0), -3/32*ri^2/qi^2*r0*(1+q0)*(1-q^2)^2, ...
      -1/8*ri/qi*(1+q0)*(1-q^2)*(1+3*r0), -1/16*(1+q0)^2*(1+r0)*(2-q0)*(1+2*r0-5*r^2); ...
      -9/64*qi*ri*(1-q^2)*(1-r^2)*(6-5*q^2-5*r^2), 3/64*qi*(1+r0)^2*(1-q^2)*(1-3*r0)*(7-5*r0), ...
      -3/64*ri*(1+q0)^2*(1-r^2)*(1-3*q0)*(7-5*q0), -3/64*ri/qi*(1+q0)*(1-q^2)*(1-r^2)*(1-5*q0), ...
      -1/16*(1+q0)*(1+r0)*(1-3*q0)*(1-3*r0), -3/64*qi/ri*(1-q^2)*(1-r^2)*(1+r0)*(1-5*r0)];
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=[-9/32*qi^2*ri*q0*(1-r^2)*(11-10*q^2-5*r^2), 3/32*qi^2*q0*(1+r0)^2*(1-3*r0)*(7-5*r0) , ...
   -9/16*qi*ri*(1-r^2)*(1-q^2)*(1+5*q0), -3/16*ri*(1+q0)*(1-r^2)*(1+2*q0-5*q^2), ...
   -1/8*qi*(1+r0)*(1+3*q0)*(1-3*r0), -3/32*qi^2/ri*q0*(1+r0)^2*(1-r0)*(1-5*r0); ...
    3/32*qi^3*(1+r0)^2*(6+(r0-2)*(30*q^2+3*r^2-7)), -3/32*qi^3/ri*(1+r0)^3*(1-r0)*(5-3*r0), ...
    3/16*qi^2*(1+r0)^2*(2-r0)*(5-2*q0-15*q^2), 3/16*qi*(1+r0)^2*(2-r0)*(1-2*q0-5*q^2), ...
    3/8*qi^2/ri*(1+r0)^2*(1-r0), 3/32*qi^3/ri^2*(1+r0)^3*(1-r0)^2];
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/32*ri^3*(1+q0)^2*(6+(q0-2)*(30*r^2+3*q^2-7)), 3/16*ri^2*(1+q0)^2*(2-q0)*(5-2*r0-15*r^2), ...
   -3/32*ri^3/qi*(1+q0)^2*(1-q^2)*(5-3*q0), -3/32*ri^3/qi^2*(1+q0)*(1-q^2)^2, ...
   -3/8*ri^2/qi*(1+q0)*(1-q^2), -3/16*ri*(1+q0)^2*(2-q0)*(1-2*r0-5*r^2); ...
    9/32*qi*ri^2*r0*(1-q^2)*(11-5*q^2-10*r^2), -9/16*qi*ri*(1-q^2)*(1-r^2)*(1+5*r0), ...
    3/32*ri^2*r0*(1+q0)^2*(1-3*q0)*(7-5*q0), 3/32*ri^2/qi*(1+q0)*(1-q^2)*r0*(1-5*q0), ...
    1/8*ri*(1+q0)*(1-3*q0)*(1+3*r0), 3/16*qi*(1-q^2)*(1+r0)*(1+2*r0-5*r^2)];
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=[9/32*qi*ri^2*r0*(1-q^2)*(5*q^2-11+10*r^2), 9/16*qi*ri*(1-q^2)*(1-r^2)*(5*r0+1), ...
  -3/32*ri^2*(q0+1)^2*r0*(-1+3*q0)*(-7+5*q0), 3/32*ri^2/qi*(1-q^2)*(1+q0)*r0*(-1+5*q0), ...
   1/8*ri*(q0+1)*(-1+3*q0)*(1+3*r0), 3/16*qi*(1-q^2)*(1+r0)*(5*r^2-2*r0-1); ...
  -9/32*qi^2*ri*q0*(1-r^2)*(5*r^2-11+10*q^2), -3/32*qi^2*q0*(3*r0-1)*(5*r0-7)*(1+r0)^2, ...
   9/16*qi*ri*(1-q^2)*(1-r^2)*(1+5*q0), -3/16*ri*(q0+1)*(1-r^2)*(5*q^2-2*q0-1), ...
  -1/8*qi*(1+r0)*(3*r0-1)*(1+3*q0), -3/32*qi^2/ri*q0*(5*r0-1)*(1-r^2)*(1+r0)];
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 D=BLxy(ni,q,r)
% Transposed second cross-derivative of scalar bilinear mapping function. 
% The parameter ni can be a vector of coordinate pairs. 
  D=[.25*ni(:,1).*ni(:,2)]; 
return;
My wiki