PFV get pressure
From CFD-Wiki
(Difference between revisions)
(Created page with "Matlab function '''GetPresW.m''' to retrieve consistent pressure from velocity. <pre> function [P,Px,Py] = GetPresW(NumNod,NodNdx,Elcon,nn2nft,Xgrid,Ygrid,Q,EBCPr,nu) %GETPRESW...") |
(Update to MatLab version R2012b and generalize to handle rectangular and triangular elements. This is a major revision.) |
||
Line 2: | Line 2: | ||
<pre> | <pre> | ||
- | function [P,Px,Py] = | + | 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 | % 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 | % Outputs | ||
- | % | + | % P,Px,Py - pressure degrees of freedom |
- | + | ||
- | + | ||
% Uses | % 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 | 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, | + | 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 | for n=1:NumNod | ||
- | + | nn2pft(n,:)=[(n-1)*npd+1,1]; | |
- | + | end | |
- | + | % ---------------------- end fixed data ---------------------------------- | |
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
% Begin essential boundary conditions, allocate space | % Begin essential boundary conditions, allocate space | ||
- | + | EBCp.Mxdof=NumPdof; | |
- | EBCp.Mxdof= | + | |
% Essential boundary condition for pressure | % Essential boundary condition for pressure | ||
- | EBCp.dof = | + | EBCp.dof = nn2pft(EBCPr.nodn(1),1); % Degree-of-freedom index |
- | EBCp.val = EBCPr.val; | + | EBCp.val = EBCPr.val; % Pressure Dof value |
% partion out essential (Dirichlet) dofs | % partion out essential (Dirichlet) dofs | ||
- | p_vec = | + | p_vec = (1:EBCp.Mxdof)'; % List of all dofs |
EBCp.p_vec_undo = zeros(1,EBCp.Mxdof); | EBCp.p_vec_undo = zeros(1,EBCp.Mxdof); | ||
% form a list of non-diri dofs | % form a list of non-diri dofs | ||
EBCp.ndro = p_vec(~ismember(p_vec, EBCp.dof)); % 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 | % calculate p_vec_undo to restore Q to the original dof ordering | ||
- | EBCp.p_vec_undo([EBCp.ndro;EBCp.dof]) = | + | 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 | |
- | for | + | Vdof(ND,n)=Q((nn2nft(Elcon(ne,n),1)-1)+ND); % Loop over local nodes |
- | + | ||
- | + | ||
end | 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; % Loop ne | ||
% END GLOBAL MATRIX ASSEMBLY | % END GLOBAL MATRIX ASSEMBLY | ||
Line 131: | Line 108: | ||
Vdatap=Vdata(EBCp.ndro)-pMat(EBCp.ndro,EBCp.dof)*EBCp.val; | Vdatap=Vdata(EBCp.ndro)-pMat(EBCp.ndro,EBCp.dof)*EBCp.val; | ||
pMatr=pMat(EBCp.ndro,EBCp.ndro); | pMatr=pMat(EBCp.ndro,EBCp.ndro); | ||
- | Qs=Qp(EBCp.ndro); % Non-Dirichlet (active) dofs | + | %Qs=Qp(EBCp.ndro); % Non-Dirichlet (active) dofs |
- | Pr= | + | [Lm,Um] = ilu(pMatr,su); % incomplete LU |
+ | Pr = gmres(pMatr,Vdatap,GM.MIter,GM.Tol,GM.MRstrt,Lm,Um,[]); % GMRES | ||
- | Qp=[Pr;EBCp.val]; | + | Qp=[Pr;EBCp.val]; % Augment active dofs with esential (Dirichlet) dofs |
- | Qp=Qp(EBCp.p_vec_undo); | + | Qp=Qp(EBCp.p_vec_undo); % Restore natural order |
- | Qp=reshape(Qp, | + | if (nargout==3) |
- | P= | + | Qp=reshape(Qp,npd,NumNod); |
- | Px= | + | P = Qp(1,:); |
- | Py= | + | Px = Qp(2,:); |
+ | Py = Qp(3,:); | ||
+ | else | ||
+ | P = Qp; | ||
+ | end | ||
return; | return; | ||
% >>>>>>>>>>>>> End pressure recovery <<<<<<<<<<<<< | % >>>>>>>>>>>>> End pressure recovery <<<<<<<<<<<<< | ||
- | % | + | % ********************* function pMassMat ******************************** |
- | function [MM,Rndx,Cndx]=pMassMat(Xe,Elcon,nn2nftp) | + | function [MM,Rndx,Cndx]=pMassMat(Xe,Elcon,nn2nftp,ep) |
% Simple cubic gradient element "mass" matrix | % 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 -------------------------- | |
- | if (isempty( | + | 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. | % Evaluate and save/cache the set of shape functions at quadrature pts. | ||
- | + | ZZ_Gpmm=cell(npn,Nq); | |
for k=1:Nq | for k=1:Nq | ||
- | for m=1: | + | for m=1:npn |
- | + | ZZ_Gpmm{m,k}=ep.G(nn(m,:),xa(k),ya(k)); | |
end | end | ||
end | end | ||
end % if(isempty(*)) | end % if(isempty(*)) | ||
- | % --------------------- end fixed data ----------------- | + | % ------------------------ end fixed data -------------------------------- |
- | TGi=cell( | + | 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( | + | MM=zeros(npdf,npdf); G=zeros(2,npdf); % Preallocate arrays |
for k=1:Nq | for k=1:Nq | ||
+ | |||
% Initialize functions and derivatives at the quadrature point (xa,ya). | % Initialize functions and derivatives at the quadrature point (xa,ya). | ||
- | J=(Xe* | + | 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); | + | 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 | + | Ji=[J(2,2),-J(1,2); -J(2,1),J(1,1)]/Det; |
- | + | for m=1:npn | |
- | + | mm=(m-1)*npd; | |
- | for m=1: | + | G(:,mm+NP)=Ji*ZZ_Gpmm{m,k}*TGi{m}; |
- | G(:,mm+NP) = Ji* | + | |
- | + | ||
end % loop m | end % loop m | ||
MM = MM + G'*G*(wt(k)*Det); | MM = MM + G'*G*(wt(k)*Det); | ||
end % end loop k (quadrature pts) | end % end loop k (quadrature pts) | ||
- | gf=zeros( | + | gf=zeros(npdf,1); % array of global freedoms |
- | + | for n=1:npn % Loop over element nodes | |
- | for n=1: | + | m=(n-1)*npd; |
- | + | gf(m+NP)=(nn2nftp(Elcon(n),1)-1)+NP; % Get global freedoms | |
- | + | ||
end | end | ||
- | Rndx=repmat(gf,1, | + | Rndx=repmat(gf,1,npdf); % Row indices |
Cndx=Rndx'; % Column indices | Cndx=Rndx'; % Column indices | ||
- | MM = reshape(MM,1, | + | MM = reshape(MM,1,npdf*npdf); |
- | Rndx=reshape(Rndx,1, | + | Rndx=reshape(Rndx,1,npdf*npdf); |
- | Cndx=reshape(Cndx,1, | + | Cndx=reshape(Cndx,1,npdf*npdf); |
+ | |||
return; | return; | ||
- | % | + | % *********************** funnction PCDat ****************************** |
- | function [PC,Rndx]=PCDat(Xe,Elcon,nn2nftp,Vdof) | + | function [PC,Rndx]=PCDat(Xe,Elcon,nn2nftp,Vdof,ep,eu) |
% Evaluate convective force data | % Evaluate convective force data | ||
- | % | + | % ep = handle to class of definitions for cubic pressure functions: |
- | + | % ELG3412r (rectangle) or ELG2309t (triangle) | |
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
% | % | ||
- | + | % ----------- Constants and fixed data --------------- | |
- | if (isempty( | + | |
+ | 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. | % 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 k=1:Nq | ||
- | for m=1: | + | 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 m over nodes | ||
end % loop k over Nq | end % loop k over Nq | ||
end % if(isempty(*)) | end % if(isempty(*)) | ||
- | % ----------------- end fixed data ------------------ | + | % ----------------------- end fixed data --------------------------------- |
- | Ti=cell( | + | Ti=cell(nvn); TGi=cell(nvn); |
- | for m=1: | + | 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, | |
- | Ti{m} = | + | TT=blkdiag(1,JtiD); |
- | TGi{m} = blkdiag(1, | + | 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 | end % Loop m over corner nodes | ||
- | PC=zeros( | + | PC=zeros(npdf,1); |
- | S=zeros(2, | + | S=zeros(2,nvdf); Sx=zeros(2,nvdf); Sy=zeros(2,nvdf); G=zeros(2,npdf); |
for k=1:Nq % Loop over quadrature points | for k=1:Nq % Loop over quadrature points | ||
- | + | Jt=Xe*eu.Gm(nn(:,:),xa(k),ya(k)); % transpose of Jacobian at (xa,ya) | |
- | Det= | + | 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; | |
- | for m=1: | + | Ji=[Jti(1,1),Jti(2,1); Jti(1,2),Jti(2,2)]; |
- | mm= | + | for m=1:nvn % Loop over element nodes |
- | S(:,mm+ | + | mm=nvd*(m-1); |
- | Sx(:,mm+ | + | S(:,mm+NV) =Jtd*ZZ_Spcd{m,k}*Ti{m}; |
- | Sy(:,mm+ | + | Sx(:,mm+NV)=Jtd*(Jti(1,1)*ZZ_SXpcd{m,k}+Jti(2,1)*ZZ_SYpcd{m,k})*Ti{m}; |
- | mm= | + | Sy(:,mm+NV)=Jtd*(Jti(1,2)*ZZ_SXpcd{m,k}+Jti(2,2)*ZZ_SYpcd{m,k})*Ti{m}; |
- | G(:,mm+NP)=Ji* | + | mm=npd*(m-1); |
+ | G(:,mm+NP)=Ji*ZZ_Gpcd{m,k}*TGi{m}; | ||
end % end loop over element nodes | end % end loop over element nodes | ||
Line 262: | Line 282: | ||
Ux = Sx*Vdof(:); | Ux = Sx*Vdof(:); | ||
Uy = Sy*Vdof(:); | Uy = Sy*Vdof(:); | ||
- | UgU = U(1)*Ux+U(2)*Uy; | + | UgU = U(1)*Ux+U(2)*Uy; % U dot grad U |
+ | |||
PC = PC + G'*UgU*(wt(k)*Det); | PC = PC + G'*UgU*(wt(k)*Det); | ||
end % end loop over Nq quadrature points | end % end loop over Nq quadrature points | ||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
- | |||
gf=zeros(1,npdf); % array of global freedoms | gf=zeros(1,npdf); % array of global freedoms | ||
- | + | for n=1:nvn % Loop over element nodes | |
- | for n=1: | + | m=(n-1)*npd; |
- | + | gf(m+NP)=(nn2nftp(Elcon(n),1)-1)+NP; % Get global freedoms | |
- | + | ||
- | + | ||
end | end | ||
Rndx=gf; | Rndx=gf; | ||
- | + | PC = reshape(PC,1,npdf); | |
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
return; | return; | ||
</pre> | </pre> |
Latest revision as of 16:10, 15 March 2013
Matlab function GetPresW.m to retrieve consistent pressure from velocity.
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;