PFVT Tdiffusion matrix

From CFD-Wiki

function [TDm,RowNdx,ColNdx]=TDMatSW(et, Xe, Elcon, nn2tft)
% TDMatSW - Affine (Hermite simple-cubic)element thermal diffusion 
%   matrix for SEGREGATED solutions.
%
% Quadratic-complete, tangential-conforming, irrotational, 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. 
%
% Usage:
%   [TDm,RowNdx,ColNdx]=TDMatSW(Xe, Elcon, nn2tft)
%   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 temperature freedom number for node n.
%   nn2tft(n,2) - global temperature freedom type for node n (not used here).
%
% Jonas Holdeman, December 2008    Revised March 2013


% Constants and fixed data
%et = ELG3412r;
nnd = et.nnodes;         % number of nodes per element (4);
nT = et.nndofs;          % nndofs = number of dofs per node, (3);  
nn = et.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. 
NDT = 1:nT;  

% ------------------------------------------------------------------------
% Define Gaussian quadrature data once, on first call.
persistent QQ_TDr; 
if isempty(QQ_TDr) 
     QRorder = 2*et.mxpowr-2; % -2, 
    [QQ_TDr.xa, QQ_TDr.ya, QQ_TDr.wt, QQ_TDr.nq] = et.hQuad(QRorder);
end  % if isempty...
xa = QQ_TDr.xa; ya = QQ_TDr.ya; wt = QQ_TDr.wt; Nq = QQ_TDr.nq;
% ------------------------------------------------------------------------
 
persistent ZZ_Gtd;  
if (isempty(ZZ_Gtd)||size(ZZ_Gtd,2)~=Nq)
% Evaluate and save/cache the set of shape functions at quadrature pts. 
  ZZ_Gtd=cell(nnd,Nq);  
  for k=1:Nq
    for m=1:nnd
      ZZ_Gtd{m,k}=et.G(nn(m,:),xa(k),ya(k));
    end
  end
end  % if(isempty(*))
% -------------------------- End fixed data ------------------------------

affine = et.isaffine(Xe);       % affine?
TBi=cell(nnd);

if (affine)       % (J constant)
% Jt=[x_q, x_r; y_q, y_r];
  J = (Xe*et.Gm(nn(:,:),0,0))'; 
  Det = J(1,1)*J(2,2)-J(1,2)*J(2,1);  % Determinant of Jt & J 
  Tb = blkdiag(1,J);
  for m=1:nnd,  TBi{m}=Tb;  end 
  Ji = [J(2,2),-J(1,2); -J(2,1),J(1,1)]/Det;   % det(J)*inv(J)
else   
  for m=1:nnd   % Loop over corner nodes 
    J=(Xe*et.Gm(nn(:,:),nn(m,1),nn(m,2)))'; % Gradient of bilinear function/element
    TBi{m} = blkdiag(1,J);
  end  % Loop m 
end  % if

TDm=zeros(nfT,nfT);  G=zeros(2,nfT);   % Preallocate arrays
for k=1:Nq  
% Initialize functions and derivatives at the quadrature point (xa,ya).
  if ~affine
    J=(Xe*et.Gm(nn(:,:),nn(m,1),nn(m,2)))';   % Gradient of bilinear fn
    Det=J(1,1)*J(2,2)-J(1,2)*J(2,1);          % Determinant of Jt & J 
    Ji=[J(2,2),-J(1,2); -J(2,1),J(1,1)]/Det;
  end
  mm = 0;
  for m=1:nnd  
    G(:,mm+NDT)=Ji*ZZ_Gtd{m,k}*TBi{m};
    mm = mm+nT;
  end  % loop m

% Submatrix ordered by T,Tx,Ty 
  TDm = TDm + (G(1,:)'*G(1,:) + G(2,:)'*G(2,:))*(Det*wt(k));
end  % loop k 

gf=zeros(nfT,1);
mm=0; 
for m=1:nnd
  gf(mm+NDT)=(nn2tft(Elcon(m),1)-1)+NDT;  % get first global freedom (T) 
  mm = mm+nT;
end %  loop on m
RowNdx=repmat(gf,1,nfT);
ColNdx=RowNdx';
RowNdx=reshape(RowNdx,nfT*nfT,1);
ColNdx=reshape(ColNdx,nfT*nfT,1); 
TDm=reshape(TDm,nfT*nfT,1);
return;