# PFV diffusion matrix 2

Function DMat4424W.m for quartic pressure-free velocity diffusion matrix

```function [Dm,RowNdx,ColNdx]=DMat4424W(Xe,Elcon,nn2nft)
%DMAT4424W - General quadrilateral Hermite (C1 quartic-derived)element convection matrix.
%
% Quartic-complete, fully-conforming, divergence-free, Hermite basis
%   functions on 4-node rectangular elements with 6 DOF per node using
%   Gauss quadrature on the 2x2 reference square.
% The assumed nodal numbering starts with 1 at the lower left corner
%   of the element and proceeds counter-clockwise around the element.
%
% Usage:
%   [Dm,Rndx,Cndx] = DMat4424W(Xe,Elcon,nn2nft)
%   Xe(1,:) -  x-coordinates of corner nodes of element.
%   Xe(2,:) -  y-coordinates of corner nodes of element.
%      and shape of the element. It is constant for affine elements.
%   Elcon - connectivity matrix for this element.
%   nn2nft - global number and type of DOF at each node
%
% Jonas Holdeman, January 2007, revised July 2011

% Constants and fixed data
nd = 6;  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))   % Has 5-point quadrature data been 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;

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

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

Ti=cell(4);
% Jt=[x_q, x_r; y_q, y_r];
for m=1:4   % Loop over corner nodes

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)
T2=[Jt(1,1)^2, 2*Jt(1,1)*Jt(1,2), Jt(1,2)^2; ...
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];
T6=blkdiag(1,JtiD,T2);
Bxy=Xe*BLxy(nn(:,:),nn(m,1),nn(m,2));  % Second cross derivatives
T6(5,2:3)=Bxy([2,1])';
Ti{m}=T6;
end  % Loop m

Dm=zeros(nd4,nd4); Fx=zeros(2,nd4); Fy=zeros(2,nd4);  Sx=zeros(2,nd4); Sy=zeros(2,nd4);  % Pre-allocate arrays

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
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);
Ds = TL*Zs4424D2d{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

Dm = Dm+(Sx'*Sx+Sy'*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

Dm = reshape(Dm,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 quartic stream function.
qi=ni(1); q0=q*ni(1); q1=1+q0;
ri=ni(2); r0=r*ni(2); r1=1+r0;
p2=[-3/32*qi^2*q0*(1+r0)^2*(r0*(10*q^2+3*r^2-7)+20-20*q^2-6*r^2), 3/32*qi^2/ri*q0*(1+r0)^3*(1-r0)*(5-3*r0), ...
-3/16*qi*(1-q^2)*(1+r0)^2*(2-r0)*(1+5*q0), -1/16*(1+q0)*(1+r0)^2*(2-r0)*(1+2*q0-5*q^2), ...
-1/8*qi/ri*(1+r0)^2*(1-r0)*(1+3*q0), -3/32*qi^2/ri^2*q0*(1+r0)^3*(1-r0)^2; ...
9/64*qi*ri*(1-q^2)*(1-r^2)*(6-5*q^2-5*r^2), -3/64*qi*(1-q^2)*(1+r0)^2*(1-3*r0)*(7-5*r0), ...
3/64*ri*(1+q0)^2*(1-r^2)*(1-3*q0)*(7-5*q0), 3/64*ri/qi*(1+q0)^2*(1-r^2)*(1-q0)*(1-5*q0), ...
1/16*(1+q0)*(1+r0)*(1-3*q0)*(1-3*r0), 3/64*qi/ri*(1-q^2)*(1+r0)^2*(1-r0)*(1-5*r0); ...
-3/32*ri^2*r0*(1+q0)^2*(q0*(10*r^2+3*q^2-7)+20-20*r^2-6*q^2), 3/16*ri*(1+q0)^2*(1-r^2)*(2-q0)*(1+5*r0), ...
-3/32*ri^2*r0/qi*(1+q0)^2*(1-q^2)*(5-3*q0), -3/32*ri^2/qi^2*r0*(1+q0)*(1-q^2)^2, ...
-1/8*ri/qi*(1+q0)*(1-q^2)*(1+3*r0), -1/16*(1+q0)^2*(1+r0)*(2-q0)*(1+2*r0-5*r^2)] ;
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;

function D=BLxy(ni,q,r)
% Transposed second cross-derivative of scalar bilinear mapping function.
% The parameter ni can be a vector of coordinate pairs.
D=[.25*ni(:,1).*ni(:,2)];
return;

```