PFVT Tdiffusion matrix
From CFD-Wiki
Revision as of 21:02, 15 March 2013 by Jonas Holdeman (Talk | contribs)
function [TDm,RowNdx,ColNdx]=TDMatSW(et, Xe, Elcon, nn2tft) % TDMatSW - Affine (Hermite simple-cubic)element thermal diffusion % matrix for SEGREGATED solutions. % % Quadratic-complete, tangential-conforming, irrotational, Hermite basis % functions on 4-node rectangular elements with 3 DOF per node 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: % [TDm,RowNdx,ColNdx]=TDMatSW(Xe, Elcon, nn2tft) % 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. % nn2tft(n,1) - global temperature freedom number for node n. % nn2tft(n,2) - global temperature freedom type for node n (not used here). % % Jonas Holdeman, December 2008 Revised March 2013 % Constants and fixed data %et = ELG3412r; nnd = et.nnodes; % number of nodes per element (4); nT = et.nndofs; % nndofs = number of dofs per node, (3); nn = et.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. NDT = 1:nT; % ------------------------------------------------------------------------ % Define Gaussian quadrature data once, on first call. persistent QQ_TDr; if isempty(QQ_TDr) QRorder = 2*et.mxpowr-2; % -2, [QQ_TDr.xa, QQ_TDr.ya, QQ_TDr.wt, QQ_TDr.nq] = et.hQuad(QRorder); end % if isempty... xa = QQ_TDr.xa; ya = QQ_TDr.ya; wt = QQ_TDr.wt; Nq = QQ_TDr.nq; % ------------------------------------------------------------------------ persistent ZZ_Gtd; if (isempty(ZZ_Gtd)||size(ZZ_Gtd,2)~=Nq) % Evaluate and save/cache the set of shape functions at quadrature pts. ZZ_Gtd=cell(nnd,Nq); for k=1:Nq for m=1:nnd ZZ_Gtd{m,k}=et.G(nn(m,:),xa(k),ya(k)); end end end % if(isempty(*)) % -------------------------- End fixed data ------------------------------ affine = et.isaffine(Xe); % affine? TBi=cell(nnd); if (affine) % (J constant) % Jt=[x_q, x_r; y_q, y_r]; J = (Xe*et.Gm(nn(:,:),0,0))'; Det = J(1,1)*J(2,2)-J(1,2)*J(2,1); % Determinant of Jt & J Tb = blkdiag(1,J); for m=1:nnd, TBi{m}=Tb; end Ji = [J(2,2),-J(1,2); -J(2,1),J(1,1)]/Det; % det(J)*inv(J) else for m=1:nnd % Loop over corner nodes J=(Xe*et.Gm(nn(:,:),nn(m,1),nn(m,2)))'; % Gradient of bilinear function/element TBi{m} = blkdiag(1,J); end % Loop m end % if TDm=zeros(nfT,nfT); G=zeros(2,nfT); % Preallocate arrays for k=1:Nq % Initialize functions and derivatives at the quadrature point (xa,ya). if ~affine J=(Xe*et.Gm(nn(:,:),nn(m,1),nn(m,2)))'; % Gradient of bilinear fn Det=J(1,1)*J(2,2)-J(1,2)*J(2,1); % Determinant of Jt & J Ji=[J(2,2),-J(1,2); -J(2,1),J(1,1)]/Det; end mm = 0; for m=1:nnd G(:,mm+NDT)=Ji*ZZ_Gtd{m,k}*TBi{m}; mm = mm+nT; end % loop m % Submatrix ordered by T,Tx,Ty TDm = TDm + (G(1,:)'*G(1,:) + G(2,:)'*G(2,:))*(Det*wt(k)); end % loop k gf=zeros(nfT,1); mm=0; for m=1:nnd gf(mm+NDT)=(nn2tft(Elcon(m),1)-1)+NDT; % get first global freedom (T) mm = mm+nT; end % loop on m RowNdx=repmat(gf,1,nfT); ColNdx=RowNdx'; RowNdx=reshape(RowNdx,nfT*nfT,1); ColNdx=reshape(ColNdx,nfT*nfT,1); TDm=reshape(TDm,nfT*nfT,1); return;