# PFV Buoyancy matrix 2

Function BMat4424SW.m for quartic pressure-free buoyancy matrix

function [Bm,RowNdx,ColNdx]=BMat4424SW(Xe, Elcon, nn2vft, nn2tft)
% BMAT4424SW - Rectangular(Hermite simple-cubic)element thermal bouyancy matrix
%  for segregated solution.
%
% Quartic-complete, conforming, solenoidal, Hermite velocity basis
%   weights on 4-node rectangular elements with 6 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]=BMat4424SA(Xe, Elcon, nn2nft)
%   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(1,n) - global freedom number for velocity at node n.
%   nn2vft(2,n) - global freedom type for node n.
%   nn2tft(1,n) - global freedom number for temperature at node n.
%   nn2tft(2,n) - global freedom type for node n.
%
% 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.
nV = 6;   nfV = nnd*nV;    % nV = number of V dofs per node, nfV = number V dofs.
NDT = 1:nT;  NDV = 1:nV;

% Define 5-point quadrature data once, on first call.
% Gaussian weights and absissas to integrate 9th degree polynomials exactly.
global GQ5;
if (isempty(GQ5))   % Has 5-point quadrature data been defined? If not, define arguments & weights.
Aq=[-.906179845938664,-.538469310105683, .0,               .538469310105683, .906179845938664];
Hq=[ .236926885056189, .478628670499366, .568888888888889, .478628670499366, .236926885056189];
GQ5.size=25; GQ5.xa=[Aq;Aq;Aq;Aq;Aq]; GQ5.ya=GQ5.xa';
wt=[Hq;Hq;Hq;Hq;Hq]; GQ5.wt=wt.*wt';
end

xa=GQ5.xa; ya=GQ5.ya; wt=GQ5.wt; Nq=GQ5.size;

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

Ti=cell(nnd);  Tgi=cell(nnd);
for m=1:nnd   % Loop over corner nodes
Jt=Xe*GBL(nn(:,:),nn(m,1),nn(m,2));
JtiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)];   % det(J)*inv(J)
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*BLxy(nn(:,:),nn(m,1),nn(m,2));  % Second cross derivatives
T6(5,2:3)=Bxy([2,1])';
Ti{m}=T6;
Tgi{m} = blkdiag(1,Jt');
end  % Loop m

% Jt=[x_q, x_r; y_q, y_r];
Bm=zeros(nfV,nfT);  S=zeros(2,nfV); g=zeros(1,nfT);   % Preallocate arrays
for k=1:Nq
Jt=Xe*GBL(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;
% Initialize functions and derivatives at the quadrature point (xa,ya).
mt = 0;  mv = 0;
for m=1:nnd
g(1,mt+NDT)=    ZG4424b{m,k}*Tgi{m}; mt = mt + nT;
S(:,mv+NDV)=Jtd*ZS4424b{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(1,Elcon(m))-1)+NDV;   mv = mv + nV; % get row dofs (V)
gfc(mt+NDT)=(nn2tft(1,Elcon(m))-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;

% ------------------------------------------------------------------------------
function g=gr(ni,q,r)
%G  Cubic Hermite basis function for simple-cubic temperature.
qi=ni(1); q0=q*ni(1);
ri=ni(2); r0=r*ni(2);
% Scalar simple-cubic Hermite
g=[1/8*(1+q0)*(1+r0)*(2+q0*(1-q0)+r0*(1-r0)), -1/8/qi*(1+q0)*(1+r0)*(1-q^2),...
-1/8/ri*(1+q0)*(1+r0)*(1-r^2)];
return;

function S=Sr(ni,q,r)
%SR  Array of quartic solenoidal basis functions.
qi=ni(1); q0=q*ni(1); q1=1+q0;
ri=ni(2); r0=r*ni(2); r1=1+r0;
% array of solenoidal vectors
S=[3/64*ri*(1+q0)^2*(1-r^2)*((1+q0)*(8-9*q0+3*q^2)+(2-q0)*(1-5*r^2)), ...
-1/64*(1+q0)^2*(1+r0)^2*(2-q0)*(7-26*r0+15*r^2), ...
3/64*ri/qi*(1+q0)^3*(1-r^2)*(1-q0)*(5-3*q0), ...
3/64*ri/qi^2*(1+q0)^3*(1-r^2)*(1-q0)^2, ...
1/16/qi*(1+q0)^2*(1+r0)*(1-q0)*(1-3*r0), ...
1/64/ri*(1+q0)^2*(1+r0)^2*(1-r0)*(1-5*r0)*(2-q0); ...
-3/64*qi*(1+r0)^2*(1-q^2)*((1+r0)*(8-9*r0+3*r^2)+(2-r0)*(1-5*q^2)), ...
3/64*qi/ri*(1+r0)^3*(1-q^2)*(1-r0)*(5-3*r0), ...
-1/64*(1+r0)^2*(1+q0)^2*(2-r0)*(7-26*q0+15*q^2), ...
-1/64/qi*(1+r0)^2*(1+q0)^2*(1-q0)*(1-5*q0)*(2-r0), ...
-1/16/ri*(1+r0)^2*(1+q0)*(1-r0)*(1-3*q0), ...
-3/64*qi/ri^2*(1+r0)^3*(1-q^2)*(1-r0)^2];
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;

function D=BLxy(ni,q,r)
% Transposed second cross-derivative of scalar bilinear mapping function.
% The parameter ni can be a vector of coordinate pairs.
D=[.25*ni(:,1).*ni(:,2)];
return;