PFVT Buoyancy matrix

From CFD-Wiki

function [Bm,RowNdx,ColNdx]=BMatSW(eu,Xe, Elcon, nn2vft, nn2tft)
% BMatSW - Rectangular(Hermite simple-cubic)element thermal bouyancy matrix
%  for segregated solution.
%
% Quadratic-complete, normal-conforming, solenoidal, Hermite basis for
%   velocity on 4-node rectangular elements with 3 DOF per node
%   and cubic Hermite temperature elements 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:
%   [Bm,RowNdx,ColNdx]=BMatSW(Xe, Elcon, nn2nft,eu)
%   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. 
%   nn2vft(n,1) - global freedom number for velocity at node n.
%   nn2vft(n,2) - global freedom type for velocity for node n.
%   nn2tft(n,1) - global freedom number for temperature at node n.
%   nn2tft(n,2) - global freedom type for temperature for node n.
%   eu          - class of shape function definitions (ELS4424r)
%
% Jonas Holdeman, December 2008     Revised March 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);  
%nedofs = nnodes*nndofs;  
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_TBr; 
if isempty(QQ_TBr) 
     QRorder = eu.mxpowr+et.mxpowr-2; % -2, 
    [QQ_TBr.xa, QQ_TBr.ya, QQ_TBr.wt, QQ_TBr.nq] = eu.hQuad(QRorder);
end  % if isempty...
xa = QQ_TBr.xa; ya = QQ_TBr.ya; wt = QQ_TBr.wt; Nq = QQ_TBr.nq;
% ------------------------------------------------------------------------
 
persistent ZZ_Stb; persistent ZZ_Gtb;  
if (isempty(ZZ_Stb)||size(ZZ_Stb,2)~=Nq)
% Evaluate and save/cache the set of shape functions at quadrature pts. 
ZZ_Stb=cell(nnd,Nq); ZZ_Gtb=cell(nnd,Nq);  
  for k=1:Nq
      for m=1:nnd
      ZZ_Stb{m,k}= eu.S(nn(m,:),xa(k),ya(k));
      ZZ_Gtb{m,k}= et.g(nn(m,:),xa(k),ya(k));
      end
  end
end
% --------------- End fixed data ----------------
 

affine = eu.isaffine(Xe);       % affine?
Ti=cell(nnd);  Tgi=cell(nnd);
if (affine==1)
  Jt=Xe*eu.Gm(nn(:,:),0,0);        % (J constant)
  JtiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)];   % det(J)*inv(J)
  if nV<4, 
    TI = blkdiag(1,JtiD); 
  else
    T2=[Jt(1,1)^2, 2*Jt(1,1)*Jt(1,2), Jt(1,2)^2; ... 
    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];
    TI=blkdiag(1,JtiD,T2); 
    Bxy=Xe*eu.DGm(nn(:,:),0,0); % Second cross derivatives
    TI(5,2:3)=Bxy([2,1])';
  end  % nV...
  for m=1:nnd, Ti{m}=TI; Tgi{m}=blkdiag(1,Jt'); end
  Det=Jt(1,1)*Jt(2,2)-Jt(1,2)*Jt(2,1);  % Determinant of Jt & J 
  Jtd=Jt/Det;
else       % ~affine
   
  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];
      TI=blkdiag(1,JtiD,T2); 
      Bxy=Xe*es.DGm(nn(:,:),nn(m,1),nn(m,2)); % Second cross derivatives
      TI(5,2:3)=Bxy([2,1])';
      Ti{m}=TI;
    end  % nV...
    Tgi{m}=blkdiag(1,Jt');
  end  % for m=... 
end  % if

Bm=zeros(nfV,nfT);  S=zeros(2,nfV); g=zeros(1,nfT);   % Allocate arrays
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;
  end
  % Initialize functions and derivatives at the quadrature point (xa,ya).
  mt = 0;  mv = 0;
for m=1:nnd  
  g(1,mt+NDT)=    ZZ_Gtb{m,k}*Tgi{m}; mt = mt + nT;
  S(:,mv+NDV)=Jtd*ZZ_Stb{m,k}*Ti{m}; mv = mv + nV;
end  % loop m
% Label rows by the test (weight) function index, columns by trial function index? 
% Submatrix ordered by psi,u,v 
  Bm = Bm + S(2,:)'*g*wt(k);   % Sy'*g

end  % loop k 
Bm = Bm*Det;  

gfr=zeros(nfV,1);  gfc=zeros(1,nfT);
mv = 0; mt=0;
for m=1:nnd
%  m, mv
  gfr(mv+NDV)=(nn2vft(Elcon(m),1)-1)+NDV;   mv = mv + nV; % get row dofs (V) 
  gfc(mt+NDT)=(nn2tft(Elcon(m),1)-1)+NDT;   mt = mt + nT; % get col dofs (T) 
end %  loop on k
%gfc=gfr'+3;
RowNdx=repmat(gfr,1,nfT);
ColNdx=repmat(gfc,nfV,1);
RowNdx=reshape(RowNdx,nfV*nfT,1);
ColNdx=reshape(ColNdx,nfV*nfT,1); 
Bm=reshape(Bm,nfV*nfT,1);
return;