PFVT get pressure

From CFD-Wiki

function [P,Px,Py] = GetPres3W(eu,NodXY,Elcon,nn2nft,Q,EBCPr,~) 
% GetPres3W - Compute continuous simple-cubic pressure from velocity
%   field on general quadrilateral grid (bilinear geometric mapping) or
%   quadratic pressure for triangular grid (linear mapping).
% This version is restricted to 3-node triangles and 4-node quadrilaterals
% P,Px,Py must be reshaped or restructured for use in calling program with
%   P=reshape(P,NumNx,NumNy), etc. because it assumes that the grid may be
%   unstructured. 
%
% Usage
%   P = GetPres3W(eu,NodXY,Elcon,nn2nft,Q,EBCPr);
%   [P,Px,Py] = GetPres3W(eu,NodXY,Elcon,nn2nft,Q,EBCPr);
%
% Inputs
%   eu     - velocity class, (eg. ELS3412r, ELS4424r, ELS5424r, ELS2309t )
%   NodXY  - coordinates of nodes 
%   Elcon  - element connectivity, nodes in element 
%   nn2nft - global number and type of (non-pressure) DOF at each node 
%   Q      - array of DOFs for velocity elements 
%   EBCp   - essential pressure boundary condition structure
%     EBCp.nodn - node number of fixed pressure node 
%     EBCp.val  - value of pressure 
% Outputs 
%   P,Px,Py  - pressure degrees of freedom
% Uses 
%   ilu      - ilu preconditioner
%   gmres    - to solve the system 
% Indirectly may use (handle passed by eu):
%   GQuad2   - function providing 2D rectangle quadrature rules.
%   TQuad2   - function providing 2D triangle quadrature rules.
%   ELG3412r - basis function class defining the cubic/pressure elements(Q)
%   ELG2309t - basis function class defining the quadratic/pressure elements(T)
%
% Jonas Holdeman,   January 2007, revised March 2013  

% ------------------- Constants and fixed data ---------------------------

nvn = eu.nnodes;        % Number of nodes in elements (4)
nvd = eu.nndofs;        % number of velocity DOFs at nodes (3|4|6); 
if nvn==4, ep = ELG3412r;    % simple cubic voricity class definition (Q)
elseif nvn==3, ep=ELG2309t;  % quadratic voricity class definition (T)
else error(['pressure: ' num2str(nvn) ' nodes not supported']);
end
npd = ep.nndofs;        % Number DOFs in pressure fns (3, simple cubic)
ND=1:nvd;               % Number DOFs in velocity fns (bicubic-derived)

NumEl=size(Elcon,1);    % Number of elements 
NumNod=size(NodXY,1);   % Number of global nodes 
NumPdof=npd*NumNod;     % Max number pressure DOFs 

% Parameters for GMRES solver 
GM.Tol=1.e-11;
GM.MIter=30; 
GM.MRstrt=8;
% parameters for ilu preconditioner
% Decrease su.droptol if ilu preconditioner fails
su.type='ilutp';
su.droptol=1.e-5;

nn2pft = zeros(NumNod,2);
for n=1:NumNod
    nn2pft(n,:)=[(n-1)*npd+1,1];
end
% ---------------------- end fixed data ----------------------------------

% Begin essential boundary conditions, allocate space 
EBCp.Mxdof=NumPdof;
% Essential boundary condition for pressure 
EBCp.dof = nn2pft(EBCPr.nodn(1),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';

% Allocate space for pressure matrix, velocity data  
pMat = spalloc(NumPdof,NumPdof,36*NumPdof);   % allocate P mass matrix
Vdata = zeros(NumPdof,1);        % allocate velocity data 
Vdof = zeros(nvd,nvn);             % Nodal velocity DOFs 
Xe = zeros(2,nvn);
   
% BEGIN GLOBAL MATRIX ASSEMBLY
for ne=1:NumEl
    
     Xe(1:2,1:nvn)=NodXY(Elcon(ne,1:nvn),1:2)'; 

% Get stream function and velocities
   for n=1:nvn  
     Vdof(ND,n)=Q((nn2nft(Elcon(ne,n),1)-1)+ND); % Loop over local nodes
   end
   
%       Pressure "mass" matrix
   [Emat,Rndx,Cndx] = pMassMat(Xe,Elcon(ne,:),nn2pft,ep);  
   pMat=pMat+sparse(Rndx,Cndx,Emat,NumPdof,NumPdof);  % Global mass assembly 
%       Convective data term   
   [CDat,RRndx] = PCDat(Xe,Elcon(ne,:),nn2pft,Vdof,ep,eu); 
   Vdata(RRndx) = Vdata(RRndx)-CDat(:);
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 

[Lm,Um] = ilu(pMatr,su);      % incomplete LU
Pr = gmres(pMatr,Vdatap,GM.MIter,GM.Tol,GM.MRstrt,Lm,Um,[]);	% GMRES

Qp=[Pr;EBCp.val];      % Augment active dofs with esential (Dirichlet) dofs
Qp=Qp(EBCp.p_vec_undo);    % Restore natural order
if (nargout==3)
   Qp=reshape(Qp,npd,NumNod);
   P =  Qp(1,:);
   Px = Qp(2,:); 
   Py = Qp(3,:); 
else
    P = Qp;
end
return;
% >>>>>>>>>>>>> End pressure recovery <<<<<<<<<<<<<

% ********************* function pMassMat ********************************

function [MM,Rndx,Cndx]=pMassMat(Xe,Elcon,nn2nftp,ep)
% Simple cubic gradient element "mass" matrix 
% ep = handle to class of definitions for cubic pressure functions: 
% ELG3412r (rectangle) or ELG2309t (triangle)
%
% -------------------- Constants and fixed data --------------------------
npn = ep.nnodes;        % number of velocity/vorticity element nodes (4)
npd = ep.nndofs;        % number of vorticity DOFs at nodes (3); 
nn = ep.nn;           % defines local nodal order [-1 -1; 1 -1; 1 1; -1 1]
npdf=npn*npd; 
NP=1:npd; 
% ------------------------------------------------------------------------

persistent QQ_prMMp;       % quadrature rules
if isempty(QQ_prMMp) 
   QRord = (2*ep.mxpowr+1); % quadtature rule order 
  [QQ_prMMp.xa, QQ_prMMp.ya, QQ_prMMp.wt, QQ_prMMp.nq]=ep.hQuad(QRord);
end  % if isempty...
xa = QQ_prMMp.xa; ya = QQ_prMMp.ya; wt = QQ_prMMp.wt; Nq = QQ_prMMp.nq;
% ------------------------------------------------------------------------
persistent ZZ_Gpmm;        % pressure functions
if (isempty(ZZ_Gpmm)||size(ZZ_Gpmm,2)~=Nq)
% Evaluate and save/cache the set of shape functions at quadrature pts. 
  ZZ_Gpmm=cell(npn,Nq); 
  for k=1:Nq
    for m=1:npn
      ZZ_Gpmm{m,k}=ep.G(nn(m,:),xa(k),ya(k));
    end
  end
end  % if(isempty(*))
% ------------------------ end fixed data --------------------------------

TGi=cell(npn);
for m=1:npn   % Loop over corner nodes, GBL is gradient of bilinear fn
  J=(Xe*ep.Gm(nn(:,:),nn(m,1),nn(m,2)))'; %
  TGi{m} = blkdiag(1,J);
end  % Loop m 

MM=zeros(npdf,npdf);  G=zeros(2,npdf);   % Preallocate arrays
for k=1:Nq  
    
% Initialize functions and derivatives at the quadrature point (xa,ya).
  J=(Xe*ep.Gm(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;
  for m=1:npn  
     mm=(m-1)*npd;
    G(:,mm+NP)=Ji*ZZ_Gpmm{m,k}*TGi{m};
  end  % loop m
  MM = MM + G'*G*(wt(k)*Det);
end        % end loop k (quadrature pts)

gf=zeros(npdf,1);           % array of global freedoms 
for n=1:npn             % Loop over element nodes 
   m=(n-1)*npd;
   gf(m+NP)=(nn2nftp(Elcon(n),1)-1)+NP; % Get global freedoms
end

Rndx=repmat(gf,1,npdf);     % Row indices
Cndx=Rndx';                % Column indices
 
MM = reshape(MM,1,npdf*npdf);
Rndx=reshape(Rndx,1,npdf*npdf);
Cndx=reshape(Cndx,1,npdf*npdf);   

return;

% *********************** funnction PCDat ******************************

function [PC,Rndx]=PCDat(Xe,Elcon,nn2nftp,Vdof,ep,eu)
% Evaluate convective force data 
% ep = handle to class of definitions for cubic pressure functions: 
%  ELG3412r (rectangle) or ELG2309t (triangle)
%
% ----------- Constants and fixed data ---------------

nvn = eu.nnodes;        % number of velocity element nodes   (4)
nvd = eu.nndofs;        % number of velocity DOFs at nodes (3|4|6); 
npd = ep.nndofs;        % number of vorticity DOFs at nodes (3);   
nn = eu.nn;            % defines local nodal order [-1 -1; 1 -1; 1 1; -1 1]
npdf=nvn*npd;  
nvdf=nvn*nvd;  
NP=1:npd;
NV=1:nvd;
% ------------------------------------------------------------------------
persistent QQ_prPCp;       % quadrature rules
if isempty(QQ_prPCp) 
   QRord = (eu.mxpowr+ep.mxpowr-1); % quadtature rule order 
  [QQ_prPCp.xa, QQ_prPCp.ya, QQ_prPCp.wt, QQ_prPCp.nq]=eu.hQuad(QRord);
end  % if isempty...
xa = QQ_prPCp.xa; ya = QQ_prPCp.ya; wt = QQ_prPCp.wt; Nq = QQ_prPCp.nq;
% ------------------------------------------------------------------------
persistent ZZ_Spcd; persistent ZZ_SXpcd; persistent ZZ_SYpcd; 
persistent ZZ_Gpcd;  
if (isempty(ZZ_Spcd)||isempty(ZZ_Gpcd)||size(ZZ_Spcd,2)~=Nq)
% Evaluate and save/cache the set of shape functions at quadrature pts. 
  ZZ_Spcd=cell(nvn,Nq); ZZ_SXpcd=cell(nvn,Nq); 
  ZZ_SYpcd=cell(nvn,Nq); ZZ_Gpcd=cell(nvn,Nq);  
  for k=1:Nq
    for m=1:nvn
      ZZ_Spcd{m,k} =eu.S(nn(m,:),xa(k),ya(k));
      [ZZ_SXpcd{m,k},ZZ_SYpcd{m,k}]=eu.DS(nn(m,:),xa(k),ya(k));
      ZZ_Gpcd{m,k}=ep.G(nn(m,:),xa(k),ya(k));
    end  % loop m over nodes  
  end  % loop k over Nq
end  % if(isempty(*))
% ----------------------- end fixed data ---------------------------------

Ti=cell(nvn);  TGi=cell(nvn);
for m=1:nvn   % Loop over corner nodes 
    Jt=Xe*eu.Gm(nn(:,:),nn(m,1),nn(m,2)); 
    JtiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)];   % det(J)*inv(J)
  if nvd==3, 
    TT=blkdiag(1,JtiD);
  elseif nvd==4
    TT=blkdiag(1,JtiD,(Jt(1,1)*Jt(2,2)+Jt(1,2)*Jt(2,1)) );
  else
    T2=[Jt(1,1)^2, 2*Jt(1,1)*Jt(2,1), Jt(2,1)^2; ...  % alt
    Jt(1,1)*Jt(1,2), Jt(1,1)*Jt(2,2)+Jt(1,2)*Jt(2,1), Jt(2,1)*Jt(2,2); ...
    Jt(1,2)^2, 2*Jt(1,2)*Jt(2,2), Jt(2,2)^2];
    TT=blkdiag(1,JtiD,T2); 
    Bxy=Xe*eu.DGm(nn(:,:),nn(m,1),nn(m,2));  % Second cross derivatives
    TT(5,2)= Bxy(2);
    TT(5,3)=-Bxy(1);
  end
  Ti{m}=TT;
% 
%   J=[Jt(1,1),Jt(2,1); Jt(1,2),Jt(2,2)];   % evaluate J from transpose
   TGi{m} = blkdiag(1,Jt');
end  % Loop m over corner nodes

PC=zeros(npdf,1);
S=zeros(2,nvdf);  Sx=zeros(2,nvdf);  Sy=zeros(2,nvdf);  G=zeros(2,npdf);

for k=1:Nq      % Loop over quadrature points 
  Jt=Xe*eu.Gm(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;
  JtiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)];
  Jti=JtiD/Det;
  Ji=[Jti(1,1),Jti(2,1); Jti(1,2),Jti(2,2)]; 
  for m=1:nvn         % Loop over element nodes 
    mm=nvd*(m-1);
    S(:,mm+NV) =Jtd*ZZ_Spcd{m,k}*Ti{m};
    Sx(:,mm+NV)=Jtd*(Jti(1,1)*ZZ_SXpcd{m,k}+Jti(2,1)*ZZ_SYpcd{m,k})*Ti{m};
    Sy(:,mm+NV)=Jtd*(Jti(1,2)*ZZ_SXpcd{m,k}+Jti(2,2)*ZZ_SYpcd{m,k})*Ti{m};
     mm=npd*(m-1);
    G(:,mm+NP)=Ji*ZZ_Gpcd{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,npdf);          % array of global freedoms 
for n=1:nvn            % Loop over element nodes 
   m=(n-1)*npd;
   gf(m+NP)=(nn2nftp(Elcon(n),1)-1)+NP; % Get global freedoms
end
Rndx=gf;
PC = reshape(PC,1,npdf);
return;