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;
