PFVT Buoyancy matrix
From CFD-Wiki
Revision as of 21:07, 15 March 2013 by Jonas Holdeman (Talk | contribs)
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;