PFV4 convection matrix

From CFD-Wiki

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;