PFV Tdiffusion matrix 2

From CFD-Wiki

Function TDMat4424SW.m for simple-cubic pressure-free thermal diffusion matrix

function [TDm,RowNdx,ColNdx]=TDMat4424SW(Xe, Elcon, nn2tft)
% TDMAT4424SW - General (Hermite simple-cubic)element thermal diffusion 
%   matrix for SEGREGATED solutions.
%
% Simple-cubic Hermite basis, gradient is tangential-conforming, irrotational, 
%   function 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]=TDMat4424SW(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(1,n) - global temperature freedom number for node n.
%   nn2tft(2,n) - global temperature freedom type for node n (not used here).
%
% Jonas Holdeman, December 2008,  revised July 2011

% Constants and fixed data
nn=[-1 -1; 1 -1; 1 1; -1 1]; % defines local nodal order
nnd = 4;                   % nnd = number of nodes in element
nT = 3;   nfT = nnd*nT;    % nT = number of T dofs per node, nfT = number T dofs. 
NDT = 1:nT;  

% 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'; 
   W=[Hq;Hq;Hq;Hq]; GQ4.wt=W.*W';     % 4x4 
end

xa=GQ4.xa; ya=GQ4.ya; wt=GQ4.wt; Nq=GQ4.size;
 
global ZG3412td;  
if (isempty(ZG3412td)|size(ZG3412td,2)~=Nq)
% Evaluate and save/cache the set of shape functions at quadrature pts. 
  ZG3412td=cell(nnd,Nq);  
  for k=1:Nq
    for m=1:nnd
      ZG3412td{m,k}=Gr(nn(m,:),xa(k),ya(k));
    end
  end
end  % if(isempty(*))
% --------------- End fixed data ----------------

TBi=cell(nnd);
for m=1:nnd   % Loop over corner nodes 
  J=(Xe*GBL(nn(:,:),nn(m,1),nn(m,2)))';  % Gradient of bilinear function/element
  TBi{m} = blkdiag(1,J);
end  % Loop m 

% Jt=[x_q, x_r; y_q, y_r];
TDm=zeros(12,12);  G=zeros(2,12);   % Preallocate arrays
for k=1:Nq  
% Initialize functions and derivatives at the quadrature point (xa,ya).
  J=(Xe*GBL(nn(:,:),xa(k),ya(k)))';    % Jacobian at (xa,ya)
  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;
mm = 0;
for m=1:nnd  
  G(:,mm+NDT)=Ji*ZG3412td{m,k}*TBi{m};
  mm = mm+nT;
end  % loop m
% Label rows by the test (weight) function index, columns by trial function index? 
% Submatrix ordered by T,Tx,Ty 
  TDm = TDm + (G(1,:)'*G(1,:) + G(2,:)'*G(2,:))*wt(k);
end  % loop k 
TDm = TDm*Det;  

gf=zeros(nfT,1);
mm=0; 
for m=1:nnd
  gf(mm+NDT)=(nn2tft(1,Elcon(m))-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;

function G=Gr(ni,q,r)
%GR  gradient array of scalar simple-cubic basis functions.
   qi=ni(1); q0=q*ni(1);  
   ri=ni(2); r0=r*ni(2); 
   % array of irrotational vectors 
G=[ 1/8*qi*(1+r0)*(r0*(1-r0)+3*(1-q^2)), -1/8*(1+r0)*(1+q0)*(1-3*q0), ...
   -1/8*qi/ri*(1-r^2)*(1+r0); 1/8*ri*(1+q0)*(q0*(1-q0)+3*(1-r^2)), ...
   -1/8/qi*ri*(1-q^2)*(1+q0), -1/8*(1+q0)*(1+r0)*(1-3*r0)];
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;