# PFVT Buoyancy matrix

```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;
```