PFVT Tconvection matrix
From CFD-Wiki
Revision as of 21:04, 15 March 2013 by Jonas Holdeman (Talk | contribs)
function [TCm,RowNdx,ColNdx]=TCMatSW(eu,Xe, Elcon, nn2tft,Vdof) % TCMatSW - Returns the element thermal convection matrix for % segregated solution using the simple-cubic Hermite basis Temperature % functions on 4-node rectangular elements with 3 DOF per node using % Gauss quadrature on the reference square. % Uses 4-node quartic velocity functions with 6 dof per node. % The columns of the array Vdof must contain the three degree-of-freedom % vectors in the nodal order (psi,u,v). % The assumed nodal numbering starts with 1 at the lower left corner % of the element and proceeds counter-clockwise around the element. % % Usage: % [TCm,RowNdx,ColNdx]=TCMatSW(Xe, Elcon, nn2nft,Vdof) % 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 freedom number for node n. % nn2tft(n,2) - global freedom type for node n. % Vdof - (nx4) array of stream function and velocity components % (psi,u,v) to specify the velocity over the element. % eu - class of shape function definitions % % Jonas Holdeman, January 2009 Revised 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); 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_TCr; if isempty(QQ_TCr) QRorder = eu.mxpowr+2*et.mxpowr-4; % -2, [QQ_TCr.xa, QQ_TCr.ya, QQ_TCr.wt, QQ_TCr.nq] = eu.hQuad(QRorder); end % if isempty... xa = QQ_TCr.xa; ya = QQ_TCr.ya; wt = QQ_TCr.wt; Nq = QQ_TCr.nq; % ------------------------------------------------------------------------ persistent ZZ_Stc; persistent ZZ_gtc; persistent ZZ_Gtc; if (isempty(ZZ_Stc)||size(ZZ_Stc,2)~=Nq) % Evaluate and save/cache the set of shape functions at quadrature pts. ZZ_Stc=cell(nnd,Nq); ZZ_gtc=cell(nnd,Nq); ZZ_Gtc=cell(nnd,Nq); for k=1:Nq for m=1:nnd ZZ_Stc{m,k} =eu.S(nn(m,:),xa(k),ya(k)); ZZ_gtc{m,k}=et.g(nn(m,:),xa(k),ya(k)); ZZ_Gtc{m,k}=et.G(nn(m,:),xa(k),ya(k)); end end end % if(isempty(*)) % -------------------------- End fixed data ------------------------------ affine = eu.isaffine(Xe); % affine? Ti=cell(nnd); TBi=cell(nnd); if affine % (J constant) % Jt=[x_q, x_r; y_q, y_r]; Jt=Xe*eu.Gm(nn(:,:),0,0); JtiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)]; % det(J)*inv(J) if nV<4, for m=1:nnd, Ti{m}=blkdiag(1,JtiD); end 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]; T6=blkdiag(1,JtiD,T2); Bxy=Xe*eu.DGm(nn(:,:),0,0); % Second cross derivatives T6(5,2:3)=Bxy([2,1])'; for m=1:nnd, Ti{m}=T6; end end % nV... TB = blkdiag(1,Jt'); for m=1:nnd, TBi{m}=TB; end Det=Jt(1,1)*Jt(2,2)-Jt(1,2)*Jt(2,1); % Determinant of Jt & J Jtd=Jt/Det; Ji=[Jt(2,2),-Jt(2,1); -Jt(1,2),Jt(1,1)]/Det; else 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]; T6=blkdiag(1,JtiD,T2); Bxy=Xe*eu.DGm(nn(:,:),nn(m,1),nn(m,2)); % Second cross derivatives T6(5,2:3)=Bxy([2,1])'; Ti{m}=T6; end % nV... TBi{m} = blkdiag(1,Jt'); end % for m=... end % if affine... % Preallocate arrays TCm=zeros(nfT,nfT); S=zeros(2,nfV); g=zeros(1,nfT); G=zeros(2,nfT); % Begin loop over Gauss-Legendre quadrature points. 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; Ji=[Jt(2,2),-Jt(2,1); -Jt(1,2),Jt(1,1)]/Det; end % Initialize arrays of functions and derivatives at the quadrature point (xa,ya). mv = 0; mt = 0; for m=1:nnd S(:,mv+NDV)= Jtd*ZZ_Stc{m,k}*Ti{m}; mv = mv+nV; g(1,mt+NDT)= ZZ_gtc{m,k}*TBi{m}; G(:,mt+NDT)= Ji*ZZ_Gtc{m,k}*TBi{m}; mt = mt+nT; end % loop m % Compute the fluid velocity at the quadrature point. U = S*Vdof(:); % Label rows by the test (weight) function index, columns by trial function index? % Submatrix ordered by T, Tx, Ty TCm=TCm+(g'*(U(1)*G(1,:)+U(2)*G(2,:)))*(Det*wt(k)); end % loop k over quadrature points gf=zeros(nfT,1); m=0; for k=1:nnd % Loop over element nodes gf(m+NDT)=(nn2tft(Elcon(k),1)-1)+NDT; % Global freedoms for this node m=m+nT; end % loop on k RowNdx=repmat(gf,1,nfT); ColNdx=RowNdx'; RowNdx=reshape(RowNdx,nfT*nfT,1); ColNdx=reshape(ColNdx,nfT*nfT,1); TCm=reshape(TCm,nfT*nfT,1); return;