PFVT Tconvection matrix

From CFD-Wiki

function [TCm,RowNdx,ColNdx]=TCMatSW(eu,Xe, Elcon, nn2tft,Vdof)
% TCMatSW - Returns the element thermal convection matrix for
%   segregated solution using the simple-cubic Hermite basis Temperature 
%   functions on 4-node rectangular elements with 3 DOF per node using 
%   Gauss quadrature on the reference square. 
% Uses 4-node quartic velocity functions with 6 dof per node. 
% The columns of the array Vdof must contain the three degree-of-freedom 
%   vectors in the nodal order (psi,u,v). 
% The assumed nodal numbering starts with 1 at the lower left corner 
%   of the element and proceeds counter-clockwise around the element. 
%
% Usage:
%   [TCm,RowNdx,ColNdx]=TCMatSW(Xe, Elcon, nn2nft,Vdof)
%   Xe(1,:)     -  x-coordinates of 4 corner nodes of element.  
%   Xe(2,:)     -  y-coordinates of 4 corner nodes of element.  
%   Elcon(4)    - connectivity matrix for this element, list of nodes. 
%   nn2tft(n,1) - global freedom number for node n.
%   nn2tft(n,2) - global freedom type for node n.
%   Vdof        - (nx4) array of stream function and velocity components 
%                  (psi,u,v) to specify the velocity over the element.
%   eu          - class of shape function definitions 
%
% Jonas Holdeman, January 2009     Revised 2013

% ------------------ Constants and fixed data ----------------------------
%eu = ELS4424r;
et = ELG3412r;
nnd = eu.nnodes;     % number of nodes per element (4);
nV = eu.nndofs;      % nndofs = number of dofs per node, (3|6);  
nT = et.nndofs;      % nndofs = number of dofs per node, (3);  
nn = eu.nn;          % defines local nodal order, [-1 -1; 1 -1; 1 1; -1 1]
nfT = nnd*nT;        % nT = number of T dofs per node, nfT = number T dofs. 
nfV = nnd*nV;        % nV = number of V dofs per node, nfV = number V dofs.
NDT = 1:nT;  NDV = 1:nV; 

% ------------------------------------------------------------------------
% Define Gaussian quadrature data once, on first call.
persistent QQ_TCr; 
if isempty(QQ_TCr) 
    QRorder = eu.mxpowr+2*et.mxpowr-4; % -2, 
    [QQ_TCr.xa, QQ_TCr.ya, QQ_TCr.wt, QQ_TCr.nq] = eu.hQuad(QRorder);
end  % if isempty...
xa = QQ_TCr.xa; ya = QQ_TCr.ya; wt = QQ_TCr.wt; Nq = QQ_TCr.nq;
% ------------------------------------------------------------------------

persistent ZZ_Stc; persistent ZZ_gtc; persistent ZZ_Gtc; 
if (isempty(ZZ_Stc)||size(ZZ_Stc,2)~=Nq)
% Evaluate and save/cache the set of shape functions at quadrature pts. 
  ZZ_Stc=cell(nnd,Nq); ZZ_gtc=cell(nnd,Nq); ZZ_Gtc=cell(nnd,Nq);  
   for k=1:Nq
     for m=1:nnd
       ZZ_Stc{m,k} =eu.S(nn(m,:),xa(k),ya(k));
       ZZ_gtc{m,k}=et.g(nn(m,:),xa(k),ya(k));
       ZZ_Gtc{m,k}=et.G(nn(m,:),xa(k),ya(k));
     end
   end
end  % if(isempty(*))
% -------------------------- End fixed data ------------------------------
 
affine = eu.isaffine(Xe);       % affine?
Ti=cell(nnd);  TBi=cell(nnd);

if affine   % (J constant)
% Jt=[x_q, x_r; y_q, y_r];
  Jt=Xe*eu.Gm(nn(:,:),0,0); 
  JtiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)];   % det(J)*inv(J)
  if nV<4, 
     for m=1:nnd, Ti{m}=blkdiag(1,JtiD); end
  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];
    T6=blkdiag(1,JtiD,T2); 
    Bxy=Xe*eu.DGm(nn(:,:),0,0); % Second cross derivatives
    T6(5,2:3)=Bxy([2,1])';
    for m=1:nnd, Ti{m}=T6; end
  end  % nV...
  TB = blkdiag(1,Jt');
  for m=1:nnd, TBi{m}=TB;  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:nnd
    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 nV<4, 
       Ti{m}=blkdiag(1,JtiD);
    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];
      T6=blkdiag(1,JtiD,T2); 
      Bxy=Xe*eu.DGm(nn(:,:),nn(m,1),nn(m,2)); % Second cross derivatives
      T6(5,2:3)=Bxy([2,1])';
      Ti{m}=T6;
    end  % nV...
    TBi{m} = blkdiag(1,Jt');
  end  % for m=...
end  % if affine...

% Preallocate arrays
TCm=zeros(nfT,nfT);  S=zeros(2,nfV);  g=zeros(1,nfT);  G=zeros(2,nfT); 

% Begin loop over Gauss-Legendre quadrature points. 
for k=1:Nq  
  if ~affine
    Jt=Xe*es.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 arrays of functions and derivatives at the quadrature point (xa,ya).
mv = 0;  mt = 0;
for m=1:nnd  
  S(:,mv+NDV)= Jtd*ZZ_Stc{m,k}*Ti{m};  mv = mv+nV;
  g(1,mt+NDT)= ZZ_gtc{m,k}*TBi{m};
  G(:,mt+NDT)= Ji*ZZ_Gtc{m,k}*TBi{m};  mt = mt+nT;
end  % loop m
   
% Compute the fluid velocity at the quadrature point.
  U = S*Vdof(:);
  
  % Label rows by the test (weight) function index, columns by trial function index? 
% Submatrix ordered by T, Tx, Ty  
   TCm=TCm+(g'*(U(1)*G(1,:)+U(2)*G(2,:)))*(Det*wt(k));
end    % loop k over quadrature points 

gf=zeros(nfT,1);
m=0;
for k=1:nnd      % Loop over element nodes 
  gf(m+NDT)=(nn2tft(Elcon(k),1)-1)+NDT;  % Global freedoms for this node 
  m=m+nT;
end %  loop on k
RowNdx=repmat(gf,1,nfT);
ColNdx=RowNdx';
RowNdx=reshape(RowNdx,nfT*nfT,1);
ColNdx=reshape(ColNdx,nfT*nfT,1); 
TCm=reshape(TCm,nfT*nfT,1);
 
return;