PFV diffusion matrix

From CFD-Wiki

Function DMatW.m for pressure-free velocity diffusion matrix

function [Dm,RowNdx,ColNdx]=DMatW(Xe,Elcon,nn2nft)
%DMATW - Returns the affine-mapped element diffusion 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. 

%
% Cubic-complete, fully-conforming, divergence-free, Hermite basis 
%   functions on 4-node rectangular elements with 3 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. 
% Uses second derivatives of stream function.
%
% Usage:
%   [Dm,Rndx,Cndx] = DMatW(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, 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 4-point quadrature data once, on first call. 
% Gaussian weights and absissas to integrate 7th degree polynomials exactly. 
global GQ4;
if (isempty(GQ4))       % Define 4-point quadrature data once, on first call. 
  Aq=[-.861136311594053,-.339981043584856,.339981043584856, .861136311594053]; %Abs
  Hq=[ .347854845137454, .652145154862546,.652145154862546, .347854845137454]; %Wts
  GQ4.size=16; GQ4.xa=[Aq;Aq;Aq;Aq]; GQ4.ya=GQ4.xa'; 
  wt=[Hq;Hq;Hq;Hq]; GQ4.wt=wt.*wt';
end

xa=GQ4.xa; ya=GQ4.ya; wt=GQ4.wt; Nq=GQ4.size;

% -----------------------------------------------
global Zs3412d2;  
if (isempty(Zs3412d2)|size(Zs3412d2,2)~=Nq)
% Evaluate and save/cache the set of shape functions at quadrature pts. 
  Zs3412d2=cell(4,Nq);  
    for k=1:Nq
      for m=1:4
       Zs3412d2{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));   % transpose of Jacobian at node m
  JtiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)];   % det(J)*inv(Jt)
  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];

Dm=zeros(nd4,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*Zs3412d2{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 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 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;