PFV convection matrix
From CFD-Wiki
(Difference between revisions)
(Created page with "Matlab function '''CMatW.m''' for pressure-free velocity convection matrix. <pre> function [Cm,RowNdx,ColNdx]=CMatW(Xe,Elcon,nn2nft,Vdof) %CMATW - Returns the affine-mapped ele...") |
(Update to MatLab version R2012b and generalize to handle rectangular and triangular elements. Major revision.) |
||
Line 2: | Line 2: | ||
<pre> | <pre> | ||
- | function [Cm,RowNdx,ColNdx]=CMatW(Xe,Elcon,nn2nft,Vdof) | + | function [Cm,RowNdx,ColNdx,Rcm,RcNdx]=CMatW(eu,Xe,Elcon,nn2nft,Vdof) |
- | % | + | %CMatW - Returns the element convection matrix for the Hermite basis |
- | % | + | % functions with 3, 4, or 6 degrees-of-freedom and defined on a |
- | % | + | % 3-node (triangle) or 4-node (quadrilateral) element by the class |
- | % The columns of the array Vdof must contain the 3 nodal degree-of-freedom | + | % instance es using Gauss quadrature on the reference element. |
- | + | % The columns of the array Vdof must contain the 3, 4, or 6 nodal | |
- | % The degrees of freedom in Vdof are the stream function | + | % degree-of-freedom vectors in the proper nodal order. |
- | % u and v of the solenoidal velocity vector | + | % The degrees of freedom in Vdof are the stream function, the two components |
- | + | % u and v of the solenoidal velocity vector, and possibly the second | |
- | % | + | % derivatives Pxx, Pxy, Pyy of the stream function. |
% | % | ||
% Usage: | % Usage: | ||
- | % [CM,Rndx,Cndx] = CMatW(Xe,Elcon,nn2nft,Vdof) | + | % [CM,Rndx,Cndx] = CMatW(es,Xe,Elcon,nn2nft,Vdof) |
+ | % [CM,Rndx,Cndx,Rcm,RcNdx] = CMatW(es,Xe,Elcon,nn2nft,Vdof) | ||
+ | % eu - handle for basis function definitions | ||
% Xe(1,:) - x-coordinates of corner nodes of element. | % Xe(1,:) - x-coordinates of corner nodes of element. | ||
% Xe(2,:) - y-coordinates of corner nodes of element. | % Xe(2,:) - y-coordinates of corner nodes of element. | ||
% Elcon - this element connectivity matrix | % Elcon - this element connectivity matrix | ||
% nn2nft - global number and type of DOF at each node | % nn2nft - global number and type of DOF at each node | ||
- | % Vdof | + | % Vdof - (ndfxnnd) array of stream function, velocity components, and |
- | % stream function derivatives to specify the velocity over the element. | + | % perhaps second stream function derivatives to specify the velocity |
+ | % over the element. | ||
% | % | ||
- | % Jonas Holdeman, August 2007, revised | + | % Indirectly may use (handle passed by eu): |
+ | % GQuad2 - function providing 2D rectangle quadrature rules. | ||
+ | % TQuad2 - function providing 2D triangle quadrature rules. | ||
+ | % | ||
+ | % Jonas Holdeman, August 2007, revised March 2013 | ||
- | % Constants and fixed data | + | % ----------------- Constants and fixed data ------------------------- |
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | nnodes = eu.nnodes; % number of nodes per element (4); | |
+ | nndofs = eu.nndofs; % nndofs = number of dofs per node, (3|6); | ||
+ | nedofs = nnodes*nndofs; | ||
+ | nn = eu.nn; % defines local nodal order, [-1 -1; 1 -1; 1 1; -1 1] | ||
+ | |||
+ | % ------------------------------------------------------------------------ | ||
+ | persistent QQCM4; | ||
+ | if isempty(QQCM4) | ||
+ | QRorder = 3*eu.mxpowr-1; % =9 | ||
+ | [QQCM4.xa, QQCM4.ya, QQCM4.wt, QQCM4.nq] = eu.hQuad(QRorder); | ||
+ | end % if isempty... | ||
+ | xa = QQCM4.xa; ya = QQCM4.ya; wt = QQCM4.wt; Nq = QQCM4.nq; | ||
+ | % ------------------------------------------------------------------------ | ||
- | + | persistent ZZ_Sc; persistent ZZ_SXc; persistent ZZ_SYc; | |
- | + | if (isempty(ZZ_Sc)||size(ZZ_Sc,2)~=Nq) | |
- | + | ||
- | if (isempty( | + | |
% 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_Sc=cell(nnodes,Nq); ZZ_SXc=cell(nnodes,Nq); ZZ_SYc=cell(nnodes,Nq); | |
for k=1:Nq | for k=1:Nq | ||
- | for m=1: | + | for m=1:nnodes |
- | + | ZZ_Sc{m,k}= eu.S(nn(m,:),xa(k),ya(k)); | |
- | + | [ZZ_SXc{m,k},ZZ_SYc{m,k}]=eu.DS(nn(m,:),xa(k),ya(k)); | |
end | end | ||
end | end | ||
- | end % if | + | end % if isempty... |
- | % --------------- End fixed data ---------------- | + | % ----------------------- End fixed data --------------------------------- |
+ | |||
+ | if (nargout<4), ItrType=0; else ItrType=1; end % ItrType=1 for Newton iter | ||
+ | affine = eu.isaffine(Xe); % affine? | ||
+ | %affine = (sum(abs(Xe(:,1)-Xe(:,2)+Xe(:,3)-Xe(:,4)))<4*eps); % affine? | ||
- | Ti=cell( | + | Ti=cell(nnodes); |
- | % | + | % Jt=[x_q, x_r; y_q, y_r]; |
- | + | if affine % (J constant) | |
- | Jt=Xe* | + | Jt=Xe*eu.Gm(nn(:,:),eu.cntr(1),eu.cntr(2)); |
JtiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)]; % det(J)*inv(J) | JtiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)]; % det(J)*inv(J) | ||
- | + | if nndofs==3 | |
- | end | + | TT=blkdiag(1,JtiD); |
- | + | elseif nndofs==4 | |
- | % | + | TT=blkdiag(1,JtiD,Jt(1,1)*Jt(2,2)+Jt(1,2)*Jt(2,1)); |
- | % Jt=[ | + | else |
- | + | T2=[Jt(1,1)^2, 2*Jt(1,1)*Jt(1,2), Jt(1,2)^2; ... % alt | |
- | Cm=zeros( | + | Jt(1,1)*Jt(2,1), Jt(1,1)*Jt(2,2)+Jt(1,2)*Jt(2,1), Jt(1,2)*Jt(2,2); ... |
- | S=zeros(2, | + | Jt(2,1)^2, 2*Jt(2,1)*Jt(2,2), Jt(2,2)^2]; |
- | + | TT=blkdiag(1,JtiD,T2); | |
+ | Bxy=Xe*eu.DGm(nn(:,:),0,0); % Second cross derivatives | ||
+ | TT(5,2)= Bxy(2); | ||
+ | TT(5,3)=-Bxy(1); | ||
+ | end % nndofs... | ||
+ | for m=1:nnodes, Ti{m}=TT; end | ||
+ | Det=Jt(1,1)*Jt(2,2)-Jt(1,2)*Jt(2,1); % Determinant of Jt & J | ||
+ | Jtd=Jt/Det; | ||
+ | Ji=[Jt(2,2),-Jt(2,1); -Jt(1,2),Jt(1,1)]/Det; | ||
+ | else | ||
+ | for m=1:nnodes | ||
+ | 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 nndofs==3 | ||
+ | TT=blkdiag(1,JtiD); | ||
+ | elseif nndofs==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(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]; | ||
+ | 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 % nndofs... | ||
+ | Ti{m}=TT; | ||
+ | end % for m=... | ||
+ | end % affine | ||
+ | % Allocate arrays | ||
+ | Cm=zeros(nedofs,nedofs); Rcm=zeros(nedofs,1); | ||
+ | S=zeros(2,nedofs); Sx=zeros(2,nedofs); Sy=zeros(2,nedofs); | ||
+ | ND=1:nndofs; | ||
% Begin loop over Gauss-Legendre quadrature points. | % Begin loop over Gauss-Legendre quadrature points. | ||
for k=1:Nq | for k=1:Nq | ||
- | + | if ~affine | |
- | + | 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; | |
- | + | Ji=[Jt(2,2),-Jt(2,1); -Jt(1,2),Jt(1,1)]/Det; | |
- | + | end | |
- | + | ||
% Initialize functions and derivatives at the quadrature point (xa,ya). | % Initialize functions and derivatives at the quadrature point (xa,ya). | ||
- | for m=1: | + | for m=1:nnodes |
- | mm= | + | mm=nndofs*(m-1); |
- | S(:,mm+ND) = Jtd* | + | S(:,mm+ND)=Jtd*ZZ_Sc{m,k}*Ti{m}; |
- | + | Sx(:,mm+ND)=Jtd*(Ji(1,1)*ZZ_SXc{m,k}+Ji(1,2)*ZZ_SYc{m,k})*Ti{m}; | |
- | Sx(:,mm+ND) = | + | Sy(:,mm+ND)=Jtd*(Ji(2,1)*ZZ_SXc{m,k}+Ji(2,2)*ZZ_SYc{m,k})*Ti{m}; |
- | Sy(:,mm+ND) = | + | |
end % loop m | end % loop m | ||
- | % Compute the fluid velocity at the quadrature point. | + | % Compute the fluid velocity at the quadrature point. |
U = S*Vdof(:); | U = S*Vdof(:); | ||
+ | |||
% Submatrix ordered by psi,u,v | % Submatrix ordered by psi,u,v | ||
- | Cm = Cm + S'*(U(1)*Sx+U(2)*Sy)*(wt(k)*Det); | + | Cm = Cm + S'*(U(1)*Sx+U(2)*Sy)*(wt(k)*Det); |
+ | |||
+ | if (ItrType~=0) % iteration type is Newton | ||
+ | Ux=Sx*Vdof(:); | ||
+ | Uy=Sy*Vdof(:); | ||
+ | Cm = Cm + S'*(Ux*S(1,:)+Uy*S(2,:))*(wt(k)*Det); | ||
+ | Rcm=Rcm + S'*(U(1)*Ux+U(2)*Uy)*(wt(k)*Det); | ||
+ | end % Cm & Rcm for Newton iteration | ||
end % end loop k over quadrature points | end % end loop k over quadrature points | ||
- | gf=zeros( | + | gf=zeros(nedofs,1); |
m=0; | m=0; | ||
- | for n=1: | + | for n=1:nnodes % Loop over element nodes |
- | gf(m+ND)=(nn2nft( | + | gf(m+ND)=(nn2nft(Elcon(n),1)-1)+ND; % Get global freedoms |
- | m=m+ | + | m=m+nndofs; |
end | end | ||
- | RowNdx=repmat(gf,1, | + | RowNdx=repmat(gf,1,nedofs); % Row indices |
ColNdx=RowNdx'; % Col indices | ColNdx=RowNdx'; % Col indices | ||
- | Cm = reshape(Cm, | + | Cm = reshape(Cm,nedofs*nedofs,1); |
- | RowNdx=reshape(RowNdx, | + | RowNdx=reshape(RowNdx,nedofs*nedofs,1); |
- | ColNdx=reshape(ColNdx, | + | ColNdx=reshape(ColNdx,nedofs*nedofs,1); |
- | + | if(ItrType==0), Rcm=[]; RcNdx=[]; else RcNdx=gf; end | |
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | + | ||
- | |||
- | |||
- | |||
- | |||
return; | return; | ||
</pre> | </pre> |
Latest revision as of 16:06, 15 March 2013
Matlab function CMatW.m for pressure-free velocity convection matrix.
function [Cm,RowNdx,ColNdx,Rcm,RcNdx]=CMatW(eu,Xe,Elcon,nn2nft,Vdof) %CMatW - Returns the element convection matrix for the Hermite basis % functions with 3, 4, or 6 degrees-of-freedom and defined on a % 3-node (triangle) or 4-node (quadrilateral) element by the class % instance es using Gauss quadrature on the reference element. % The columns of the array Vdof must contain the 3, 4, or 6 nodal % degree-of-freedom vectors in the proper nodal order. % The degrees of freedom in Vdof are the stream function, the two components % u and v of the solenoidal velocity vector, and possibly the second % derivatives Pxx, Pxy, Pyy of the stream function. % % Usage: % [CM,Rndx,Cndx] = CMatW(es,Xe,Elcon,nn2nft,Vdof) % [CM,Rndx,Cndx,Rcm,RcNdx] = CMatW(es,Xe,Elcon,nn2nft,Vdof) % eu - handle for basis function definitions % 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 - (ndfxnnd) array of stream function, velocity components, and % perhaps second stream function derivatives to specify the velocity % over the element. % % Indirectly may use (handle passed by eu): % GQuad2 - function providing 2D rectangle quadrature rules. % TQuad2 - function providing 2D triangle quadrature rules. % % Jonas Holdeman, August 2007, revised March 2013 % ----------------- Constants and fixed data ------------------------- nnodes = eu.nnodes; % number of nodes per element (4); nndofs = eu.nndofs; % nndofs = number of dofs per node, (3|6); nedofs = nnodes*nndofs; nn = eu.nn; % defines local nodal order, [-1 -1; 1 -1; 1 1; -1 1] % ------------------------------------------------------------------------ persistent QQCM4; if isempty(QQCM4) QRorder = 3*eu.mxpowr-1; % =9 [QQCM4.xa, QQCM4.ya, QQCM4.wt, QQCM4.nq] = eu.hQuad(QRorder); end % if isempty... xa = QQCM4.xa; ya = QQCM4.ya; wt = QQCM4.wt; Nq = QQCM4.nq; % ------------------------------------------------------------------------ persistent ZZ_Sc; persistent ZZ_SXc; persistent ZZ_SYc; if (isempty(ZZ_Sc)||size(ZZ_Sc,2)~=Nq) % Evaluate and save/cache the set of shape functions at quadrature pts. ZZ_Sc=cell(nnodes,Nq); ZZ_SXc=cell(nnodes,Nq); ZZ_SYc=cell(nnodes,Nq); for k=1:Nq for m=1:nnodes ZZ_Sc{m,k}= eu.S(nn(m,:),xa(k),ya(k)); [ZZ_SXc{m,k},ZZ_SYc{m,k}]=eu.DS(nn(m,:),xa(k),ya(k)); end end end % if isempty... % ----------------------- End fixed data --------------------------------- if (nargout<4), ItrType=0; else ItrType=1; end % ItrType=1 for Newton iter affine = eu.isaffine(Xe); % affine? %affine = (sum(abs(Xe(:,1)-Xe(:,2)+Xe(:,3)-Xe(:,4)))<4*eps); % affine? Ti=cell(nnodes); % Jt=[x_q, x_r; y_q, y_r]; if affine % (J constant) Jt=Xe*eu.Gm(nn(:,:),eu.cntr(1),eu.cntr(2)); JtiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)]; % det(J)*inv(J) if nndofs==3 TT=blkdiag(1,JtiD); elseif nndofs==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(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]; TT=blkdiag(1,JtiD,T2); Bxy=Xe*eu.DGm(nn(:,:),0,0); % Second cross derivatives TT(5,2)= Bxy(2); TT(5,3)=-Bxy(1); end % nndofs... for m=1:nnodes, Ti{m}=TT; end Det=Jt(1,1)*Jt(2,2)-Jt(1,2)*Jt(2,1); % Determinant of Jt & J Jtd=Jt/Det; Ji=[Jt(2,2),-Jt(2,1); -Jt(1,2),Jt(1,1)]/Det; else for m=1:nnodes 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 nndofs==3 TT=blkdiag(1,JtiD); elseif nndofs==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(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]; 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 % nndofs... Ti{m}=TT; end % for m=... end % affine % Allocate arrays Cm=zeros(nedofs,nedofs); Rcm=zeros(nedofs,1); S=zeros(2,nedofs); Sx=zeros(2,nedofs); Sy=zeros(2,nedofs); ND=1:nndofs; % Begin loop over Gauss-Legendre quadrature points. for k=1:Nq if ~affine 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; Ji=[Jt(2,2),-Jt(2,1); -Jt(1,2),Jt(1,1)]/Det; end % Initialize functions and derivatives at the quadrature point (xa,ya). for m=1:nnodes mm=nndofs*(m-1); S(:,mm+ND)=Jtd*ZZ_Sc{m,k}*Ti{m}; Sx(:,mm+ND)=Jtd*(Ji(1,1)*ZZ_SXc{m,k}+Ji(1,2)*ZZ_SYc{m,k})*Ti{m}; Sy(:,mm+ND)=Jtd*(Ji(2,1)*ZZ_SXc{m,k}+Ji(2,2)*ZZ_SYc{m,k})*Ti{m}; 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); if (ItrType~=0) % iteration type is Newton Ux=Sx*Vdof(:); Uy=Sy*Vdof(:); Cm = Cm + S'*(Ux*S(1,:)+Uy*S(2,:))*(wt(k)*Det); Rcm=Rcm + S'*(U(1)*Ux+U(2)*Uy)*(wt(k)*Det); end % Cm & Rcm for Newton iteration end % end loop k over quadrature points gf=zeros(nedofs,1); m=0; for n=1:nnodes % Loop over element nodes gf(m+ND)=(nn2nft(Elcon(n),1)-1)+ND; % Get global freedoms m=m+nndofs; end RowNdx=repmat(gf,1,nedofs); % Row indices ColNdx=RowNdx'; % Col indices Cm = reshape(Cm,nedofs*nedofs,1); RowNdx=reshape(RowNdx,nedofs*nedofs,1); ColNdx=reshape(ColNdx,nedofs*nedofs,1); if(ItrType==0), Rcm=[]; RcNdx=[]; else RcNdx=gf; end return;