# PFV convection matrix

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)

Matlab function CMatW.m for pressure-free velocity convection matrix.

```function [Cm,RowNdx,ColNdx]=CMatW(Xe,Elcon,nn2nft,Vdof)
%CMATW - Returns the affine-mapped element convection matrix for the simple cubic Hermite
%   basis functions on 4-node straight-sided quadrilateral elements with 3 DOF
%   per node using Gauss quadrature on the reference square and row/col indices.
% The columns of the array Vdof must contain the 3 nodal degree-of-freedom
%   vectors in the proper nodal order.
% The degrees of freedom in Vdof are the stream function and the two components
%   u and v of the solenoidal velocity vector.
% The assumed nodal numbering starts with 1 at the lower left corner of the
%   element and proceeds counter-clockwise around the element.
%
% Usage:
%   [CM,Rndx,Cndx] = CMatW(Xe,Elcon,nn2nft,Vdof)
%   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  - (3x4) array of stream function, velocity components, and second
%     stream function derivatives to specify the velocity over the element.
%
% Jonas Holdeman, August 2007, revised  June 2011

% Constants and fixed data
nd = 3;  nd4=4*nd;  ND=1:nd;    % nd = number of dofs per node,
nn=[-1 -1; 1 -1; 1 1; -1 1];    % defines local nodal order

% Define 5-point quadrature data once, on first call.
% Gaussian weights and absissas to integrate 9th degree polynomials exactly.
global GQ5;
if (isempty(GQ5))   % 5-point quadrature data 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';
end

xa=GQ5.xa; ya=GQ5.ya; wt=GQ5.wt; Nq=GQ5.size;  % Use GQ5 (5x5) for exact integration

% -----------------------------------------------
global Zs3412D2c;  global ZS3412c;

if (isempty(ZS3412c)|isempty(Zs3412D2c)|size(Zs3412D2c,2)~=Nq)
% Evaluate and save/cache the set of shape functions at quadrature pts.
Zs3412D2c=cell(4,Nq);  ZS3412c=cell(4,Nq);
for k=1:Nq
for m=1:4
ZS3412c{m,k}= Sr(nn(m,:),xa(k),ya(k));
Zs3412D2c{m,k}=D3s(nn(m,:),xa(k),ya(k));
end
end
end  % if(isempty(*))

% --------------- End fixed data ----------------

Ti=cell(4);
%
for m=1:4
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)
Ti{m}=blkdiag(1,JtiD);
end

% Move Jacobian evaluation inside k-loop for general convex quadrilateral.
% Jt=[x_q, x_r; y_q, y_r];

Cm=zeros(nd4,nd4); Rcm=zeros(nd4,1);
S=zeros(2,nd4); Sx=zeros(2,nd4); Sy=zeros(2,nd4);  % Pre-allocate arrays

% Begin loop over Gauss-Legendre quadrature points.
for k=1:Nq

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;
TL=[Jt(2,2)^2, -2*Jt(2,1)*Jt(2,2),  Jt(2,1)^2; ...
-Jt(1,2)*Jt(2,2), Jt(1,1)*Jt(2,2)+Jt(2,1)*Jt(1,2), -Jt(1,1)*Jt(2,1); ...
Jt(1,2)^2, -2*Jt(1,1)*Jt(1,2),  Jt(1,1)^2 ]/Det^2;

% Initialize functions and derivatives at the quadrature point (xa,ya).
for m=1:4
mm=nd*(m-1);
S(:,mm+ND) = Jtd*ZS3412c{m,k}*Ti{m};
Ds = TL*Zs3412D2c{m,k}*Ti{m};
Sx(:,mm+ND) = [Ds(2,:); -Ds(1,:)];     % [Pyx, -Pxx]
Sy(:,mm+ND) = [Ds(3,:); -Ds(2,:)];     % [Pyy, -Pxy]
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);
end    % end loop k over quadrature points

gf=zeros(nd4,1);
m=0;
for n=1:4                 % Loop over element nodes
gf(m+ND)=(nn2nft(1,Elcon(n))-1)+ND;  % Get global freedoms
m=m+nd;
end

RowNdx=repmat(gf,1,nd4);      % Row indices
ColNdx=RowNdx';               % Col indices

Cm = reshape(Cm,nd4*nd4,1);
RowNdx=reshape(RowNdx,nd4*nd4,1);
ColNdx=reshape(ColNdx,nd4*nd4,1);
return;

% ----------------------------------------------------------------------------

function P2=D3s(ni,q,r)
% Second derivatives [Pxx; Pxy; Pyy] of simple cubic stream function.
qi=ni(1); q0=q*ni(1); q1=1+q0;
ri=ni(2); r0=r*ni(2); r1=1+r0;
P2=[-.75*qi^2*(r0+1)*q0, 0, -.25*qi*(r0+1)*(3*q0+1); ...
.125*qi*ri*(4-3*(q^2+r^2)), .125*qi*(r0+1)*(3*r0-1), -.125*ri*(q0+1)*(3*q0-1); ...
-.75*ri^2*(q0+1)*r0, .25*ri*(q0+1)*(3*r0+1), 0] ;
return;

function S=Sr(ni,q,r)
%S  Array of solenoidal basis functions on rectangle.
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=[ .125*ri*q1*(q0*(1-q0)+3*(1-r^2)), .125*q1*r1*(3*r0-1),    .125*ri/qi*q1^2*(1-q0); ...
-.125*qi*r1*(r0*(1-r0)+3*(1-q^2)), .125*qi/ri*r1^2*(1-r0), .125*q1*r1*(3*q0-1)];
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;
```