https://cfd-online.com/W/index.php?title=PFV4_get_pressure&feed=atom&action=history
PFV4 get pressure - Revision history
2024-03-28T14:48:45Z
Revision history for this page on the wiki
MediaWiki 1.16.5
https://cfd-online.com/W/index.php?title=PFV4_get_pressure&diff=17567&oldid=prev
Jonas Holdeman: pressure recovery
2013-03-15T19:29:51Z
<p>pressure recovery</p>
<p><b>New page</b></p><div><pre><br />
function [P,Px,Py] = GetPres3W(eu,NodXY,Elcon,nn2nft,Q,EBCPr,~) <br />
% GetPres3W - Compute continuous simple-cubic pressure from velocity<br />
% field on general quadrilateral grid (bilinear geometric mapping) or<br />
% quadratic pressure for triangular grid (linear mapping).<br />
% This version is restricted to 3-node triangles and 4-node quadrilaterals<br />
% P,Px,Py must be reshaped or restructured for use in calling program with<br />
% P=reshape(P,NumNx,NumNy), etc. because it assumes that the grid may be<br />
% unstructured. <br />
%<br />
% Usage<br />
% P = GetPres3W(eu,NodXY,Elcon,nn2nft,Q,EBCPr);<br />
% [P,Px,Py] = GetPres3W(eu,NodXY,Elcon,nn2nft,Q,EBCPr);<br />
%<br />
% Inputs<br />
% eu - velocity class, (eg. ELS3412r, ELS4424r, ELS5424r, ELS2309t )<br />
% NodXY - coordinates of nodes <br />
% Elcon - element connectivity, nodes in element <br />
% nn2nft - global number and type of (non-pressure) DOF at each node <br />
% Q - array of DOFs for velocity elements <br />
% EBCp - essential pressure boundary condition structure<br />
% EBCp.nodn - node number of fixed pressure node <br />
% EBCp.val - value of pressure <br />
% Outputs <br />
% P,Px,Py - pressure degrees of freedom<br />
% Uses <br />
% ilu - ilu preconditioner<br />
% gmres - to solve the system <br />
% Indirectly may use (handle passed by eu):<br />
% GQuad2 - function providing 2D rectangle quadrature rules.<br />
% TQuad2 - function providing 2D triangle quadrature rules.<br />
% ELG3412r - basis function class defining the cubic/pressure elements(Q)<br />
% ELG2309t - basis function class defining the quadratic/pressure elements(T)<br />
%<br />
% Jonas Holdeman, January 2007, revised March 2013 <br />
<br />
% ------------------- Constants and fixed data ---------------------------<br />
<br />
nvn = eu.nnodes; % Number of nodes in elements (4)<br />
nvd = eu.nndofs; % number of velocity DOFs at nodes (3|4|6); <br />
if nvn==4, ep = ELG3412r; % simple cubic voricity class definition (Q)<br />
elseif nvn==3, ep=ELG2309t; % quadratic voricity class definition (T)<br />
else error(['pressure: ' num2str(nvn) ' nodes not supported']);<br />
end<br />
npd = ep.nndofs; % Number DOFs in pressure fns (3, simple cubic)<br />
ND=1:nvd; % Number DOFs in velocity fns (bicubic-derived)<br />
<br />
NumEl=size(Elcon,1); % Number of elements <br />
NumNod=size(NodXY,1); % Number of global nodes <br />
NumPdof=npd*NumNod; % Max number pressure DOFs <br />
<br />
% Parameters for GMRES solver <br />
GM.Tol=1.e-11;<br />
GM.MIter=30; <br />
GM.MRstrt=8;<br />
% parameters for ilu preconditioner<br />
% Decrease su.droptol if ilu preconditioner fails<br />
su.type='ilutp';<br />
su.droptol=1.e-5;<br />
<br />
nn2pft = zeros(NumNod,2);<br />
for n=1:NumNod<br />
nn2pft(n,:)=[(n-1)*npd+1,1];<br />
end<br />
% ---------------------- end fixed data ----------------------------------<br />
<br />
% Begin essential boundary conditions, allocate space <br />
EBCp.Mxdof=NumPdof;<br />
% Essential boundary condition for pressure <br />
EBCp.dof = nn2pft(EBCPr.nodn(1),1); % Degree-of-freedom index<br />
EBCp.val = EBCPr.val; % Pressure Dof value <br />
<br />
% partion out essential (Dirichlet) dofs<br />
p_vec = (1:EBCp.Mxdof)'; % List of all dofs<br />
EBCp.p_vec_undo = zeros(1,EBCp.Mxdof);<br />
% form a list of non-diri dofs<br />
EBCp.ndro = p_vec(~ismember(p_vec, EBCp.dof)); % list of non-diri dofs<br />
% calculate p_vec_undo to restore Q to the original dof ordering<br />
EBCp.p_vec_undo([EBCp.ndro;EBCp.dof]) = (1:EBCp.Mxdof); %p_vec';<br />
<br />
% Allocate space for pressure matrix, velocity data <br />
pMat = spalloc(NumPdof,NumPdof,36*NumPdof); % allocate P mass matrix<br />
Vdata = zeros(NumPdof,1); % allocate velocity data <br />
Vdof = zeros(nvd,nvn); % Nodal velocity DOFs <br />
Xe = zeros(2,nvn);<br />
<br />
% BEGIN GLOBAL MATRIX ASSEMBLY<br />
for ne=1:NumEl<br />
<br />
Xe(1:2,1:nvn)=NodXY(Elcon(ne,1:nvn),1:2)'; <br />
<br />
% Get stream function and velocities<br />
for n=1:nvn <br />
Vdof(ND,n)=Q((nn2nft(Elcon(ne,n),1)-1)+ND); % Loop over local nodes<br />
end<br />
<br />
% Pressure "mass" matrix<br />
[Emat,Rndx,Cndx] = pMassMat(Xe,Elcon(ne,:),nn2pft,ep); <br />
pMat=pMat+sparse(Rndx,Cndx,Emat,NumPdof,NumPdof); % Global mass assembly <br />
% Convective data term <br />
[CDat,RRndx] = PCDat(Xe,Elcon(ne,:),nn2pft,Vdof,ep,eu); <br />
Vdata(RRndx) = Vdata(RRndx)-CDat(:);<br />
end; % Loop ne <br />
% END GLOBAL MATRIX ASSEMBLY<br />
<br />
Vdatap=Vdata(EBCp.ndro)-pMat(EBCp.ndro,EBCp.dof)*EBCp.val;<br />
pMatr=pMat(EBCp.ndro,EBCp.ndro);<br />
%Qs=Qp(EBCp.ndro); % Non-Dirichlet (active) dofs <br />
<br />
[Lm,Um] = ilu(pMatr,su); % incomplete LU<br />
Pr = gmres(pMatr,Vdatap,GM.MIter,GM.Tol,GM.MRstrt,Lm,Um,[]); % GMRES<br />
<br />
Qp=[Pr;EBCp.val]; % Augment active dofs with esential (Dirichlet) dofs<br />
Qp=Qp(EBCp.p_vec_undo); % Restore natural order<br />
if (nargout==3)<br />
Qp=reshape(Qp,npd,NumNod);<br />
P = Qp(1,:);<br />
Px = Qp(2,:); <br />
Py = Qp(3,:); <br />
else<br />
P = Qp;<br />
end<br />
return;<br />
% >>>>>>>>>>>>> End pressure recovery <<<<<<<<<<<<<<br />
<br />
% ********************* function pMassMat ********************************<br />
<br />
function [MM,Rndx,Cndx]=pMassMat(Xe,Elcon,nn2nftp,ep)<br />
% Simple cubic gradient element "mass" matrix <br />
% ep = handle to class of definitions for cubic pressure functions: <br />
% ELG3412r (rectangle) or ELG2309t (triangle)<br />
%<br />
% -------------------- Constants and fixed data --------------------------<br />
npn = ep.nnodes; % number of velocity/vorticity element nodes (4)<br />
npd = ep.nndofs; % number of vorticity DOFs at nodes (3); <br />
nn = ep.nn; % defines local nodal order [-1 -1; 1 -1; 1 1; -1 1]<br />
npdf=npn*npd; <br />
NP=1:npd; <br />
% ------------------------------------------------------------------------<br />
<br />
persistent QQ_prMMp; % quadrature rules<br />
if isempty(QQ_prMMp) <br />
QRord = (2*ep.mxpowr+1); % quadtature rule order <br />
[QQ_prMMp.xa, QQ_prMMp.ya, QQ_prMMp.wt, QQ_prMMp.nq]=ep.hQuad(QRord);<br />
end % if isempty...<br />
xa = QQ_prMMp.xa; ya = QQ_prMMp.ya; wt = QQ_prMMp.wt; Nq = QQ_prMMp.nq;<br />
% ------------------------------------------------------------------------<br />
persistent ZZ_Gpmm; % pressure functions<br />
if (isempty(ZZ_Gpmm)||size(ZZ_Gpmm,2)~=Nq)<br />
% Evaluate and save/cache the set of shape functions at quadrature pts. <br />
ZZ_Gpmm=cell(npn,Nq); <br />
for k=1:Nq<br />
for m=1:npn<br />
ZZ_Gpmm{m,k}=ep.G(nn(m,:),xa(k),ya(k));<br />
end<br />
end<br />
end % if(isempty(*))<br />
% ------------------------ end fixed data --------------------------------<br />
<br />
TGi=cell(npn);<br />
for m=1:npn % Loop over corner nodes, GBL is gradient of bilinear fn<br />
J=(Xe*ep.Gm(nn(:,:),nn(m,1),nn(m,2)))'; %<br />
TGi{m} = blkdiag(1,J);<br />
end % Loop m <br />
<br />
MM=zeros(npdf,npdf); G=zeros(2,npdf); % Preallocate arrays<br />
for k=1:Nq <br />
<br />
% Initialize functions and derivatives at the quadrature point (xa,ya).<br />
J=(Xe*ep.Gm(nn(:,:),xa(k),ya(k)))'; % Jacobian at (xa,ya)<br />
Det=J(1,1)*J(2,2)-J(1,2)*J(2,1); % Determinant of Jt & J <br />
Ji=[J(2,2),-J(1,2); -J(2,1),J(1,1)]/Det;<br />
for m=1:npn <br />
mm=(m-1)*npd;<br />
G(:,mm+NP)=Ji*ZZ_Gpmm{m,k}*TGi{m};<br />
end % loop m<br />
MM = MM + G'*G*(wt(k)*Det);<br />
end % end loop k (quadrature pts)<br />
<br />
gf=zeros(npdf,1); % array of global freedoms <br />
for n=1:npn % Loop over element nodes <br />
m=(n-1)*npd;<br />
gf(m+NP)=(nn2nftp(Elcon(n),1)-1)+NP; % Get global freedoms<br />
end<br />
<br />
Rndx=repmat(gf,1,npdf); % Row indices<br />
Cndx=Rndx'; % Column indices<br />
<br />
MM = reshape(MM,1,npdf*npdf);<br />
Rndx=reshape(Rndx,1,npdf*npdf);<br />
Cndx=reshape(Cndx,1,npdf*npdf); <br />
<br />
return;<br />
<br />
% *********************** funnction PCDat ******************************<br />
<br />
function [PC,Rndx]=PCDat(Xe,Elcon,nn2nftp,Vdof,ep,eu)<br />
% Evaluate convective force data <br />
% ep = handle to class of definitions for cubic pressure functions: <br />
% ELG3412r (rectangle) or ELG2309t (triangle)<br />
%<br />
% ----------- Constants and fixed data ---------------<br />
<br />
nvn = eu.nnodes; % number of velocity element nodes (4)<br />
nvd = eu.nndofs; % number of velocity DOFs at nodes (3|4|6); <br />
npd = ep.nndofs; % number of vorticity DOFs at nodes (3); <br />
nn = eu.nn; % defines local nodal order [-1 -1; 1 -1; 1 1; -1 1]<br />
npdf=nvn*npd; <br />
nvdf=nvn*nvd; <br />
NP=1:npd;<br />
NV=1:nvd;<br />
% ------------------------------------------------------------------------<br />
persistent QQ_prPCp; % quadrature rules<br />
if isempty(QQ_prPCp) <br />
QRord = (eu.mxpowr+ep.mxpowr-1); % quadtature rule order <br />
[QQ_prPCp.xa, QQ_prPCp.ya, QQ_prPCp.wt, QQ_prPCp.nq]=eu.hQuad(QRord);<br />
end % if isempty...<br />
xa = QQ_prPCp.xa; ya = QQ_prPCp.ya; wt = QQ_prPCp.wt; Nq = QQ_prPCp.nq;<br />
% ------------------------------------------------------------------------<br />
persistent ZZ_Spcd; persistent ZZ_SXpcd; persistent ZZ_SYpcd; <br />
persistent ZZ_Gpcd; <br />
if (isempty(ZZ_Spcd)||isempty(ZZ_Gpcd)||size(ZZ_Spcd,2)~=Nq)<br />
% Evaluate and save/cache the set of shape functions at quadrature pts. <br />
ZZ_Spcd=cell(nvn,Nq); ZZ_SXpcd=cell(nvn,Nq); <br />
ZZ_SYpcd=cell(nvn,Nq); ZZ_Gpcd=cell(nvn,Nq); <br />
for k=1:Nq<br />
for m=1:nvn<br />
ZZ_Spcd{m,k} =eu.S(nn(m,:),xa(k),ya(k));<br />
[ZZ_SXpcd{m,k},ZZ_SYpcd{m,k}]=eu.DS(nn(m,:),xa(k),ya(k));<br />
ZZ_Gpcd{m,k}=ep.G(nn(m,:),xa(k),ya(k));<br />
end % loop m over nodes <br />
end % loop k over Nq<br />
end % if(isempty(*))<br />
% ----------------------- end fixed data ---------------------------------<br />
<br />
Ti=cell(nvn); TGi=cell(nvn);<br />
for m=1:nvn % Loop over corner nodes <br />
Jt=Xe*eu.Gm(nn(:,:),nn(m,1),nn(m,2)); <br />
JtiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)]; % det(J)*inv(J)<br />
if nvd==3, <br />
TT=blkdiag(1,JtiD);<br />
elseif nvd==4<br />
TT=blkdiag(1,JtiD,(Jt(1,1)*Jt(2,2)+Jt(1,2)*Jt(2,1)) );<br />
else<br />
T2=[Jt(1,1)^2, 2*Jt(1,1)*Jt(2,1), Jt(2,1)^2; ... % alt<br />
Jt(1,1)*Jt(1,2), Jt(1,1)*Jt(2,2)+Jt(1,2)*Jt(2,1), Jt(2,1)*Jt(2,2); ...<br />
Jt(1,2)^2, 2*Jt(1,2)*Jt(2,2), Jt(2,2)^2];<br />
TT=blkdiag(1,JtiD,T2); <br />
Bxy=Xe*eu.DGm(nn(:,:),nn(m,1),nn(m,2)); % Second cross derivatives<br />
TT(5,2)= Bxy(2);<br />
TT(5,3)=-Bxy(1);<br />
end<br />
Ti{m}=TT;<br />
% <br />
% J=[Jt(1,1),Jt(2,1); Jt(1,2),Jt(2,2)]; % evaluate J from transpose<br />
TGi{m} = blkdiag(1,Jt');<br />
end % Loop m over corner nodes<br />
<br />
PC=zeros(npdf,1);<br />
S=zeros(2,nvdf); Sx=zeros(2,nvdf); Sy=zeros(2,nvdf); G=zeros(2,npdf);<br />
<br />
for k=1:Nq % Loop over quadrature points <br />
Jt=Xe*eu.Gm(nn(:,:),xa(k),ya(k)); % transpose of Jacobian at (xa,ya)<br />
Det=Jt(1,1)*Jt(2,2)-Jt(1,2)*Jt(2,1); % Determinant of Jt & J <br />
Jtd=Jt/Det;<br />
JtiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)];<br />
Jti=JtiD/Det;<br />
Ji=[Jti(1,1),Jti(2,1); Jti(1,2),Jti(2,2)]; <br />
for m=1:nvn % Loop over element nodes <br />
mm=nvd*(m-1);<br />
S(:,mm+NV) =Jtd*ZZ_Spcd{m,k}*Ti{m};<br />
Sx(:,mm+NV)=Jtd*(Jti(1,1)*ZZ_SXpcd{m,k}+Jti(2,1)*ZZ_SYpcd{m,k})*Ti{m};<br />
Sy(:,mm+NV)=Jtd*(Jti(1,2)*ZZ_SXpcd{m,k}+Jti(2,2)*ZZ_SYpcd{m,k})*Ti{m};<br />
mm=npd*(m-1);<br />
G(:,mm+NP)=Ji*ZZ_Gpcd{m,k}*TGi{m}; <br />
end % end loop over element nodes<br />
<br />
% Compute the fluid velocities at the quadrature point.<br />
U = S*Vdof(:);<br />
Ux = Sx*Vdof(:);<br />
Uy = Sy*Vdof(:);<br />
UgU = U(1)*Ux+U(2)*Uy; % U dot grad U <br />
<br />
PC = PC + G'*UgU*(wt(k)*Det);<br />
end % end loop over Nq quadrature points<br />
<br />
gf=zeros(1,npdf); % array of global freedoms <br />
for n=1:nvn % Loop over element nodes <br />
m=(n-1)*npd;<br />
gf(m+NP)=(nn2nftp(Elcon(n),1)-1)+NP; % Get global freedoms<br />
end<br />
Rndx=gf;<br />
PC = reshape(PC,1,npdf);<br />
return;<br />
</pre></div>
Jonas Holdeman