PFV Tdiffusion matrix 2
From CFD-Wiki
Revision as of 16:08, 7 July 2011 by Jonas Holdeman (Talk | contribs)
Function TDMat4424SW.m for simple-cubic pressure-free thermal diffusion matrix
function [TDm,RowNdx,ColNdx]=TDMat4424SW(Xe, Elcon, nn2tft) % TDMAT4424SW - General (Hermite simple-cubic)element thermal diffusion % matrix for SEGREGATED solutions. % % Simple-cubic Hermite basis, gradient is tangential-conforming, irrotational, % function 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]=TDMat4424SW(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(1,n) - global temperature freedom number for node n. % nn2tft(2,n) - global temperature freedom type for node n (not used here). % % 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. NDT = 1:nT; % Define 4-point quadrature data once, on first call. % Gaussian weights and absissas to integrate 7th degree polynomials exactly. global GQ4; if (isempty(GQ4)) % Define 4-point quadrature data once, on first call. Aq=[-.861136311594053,-.339981043584856,.339981043584856, .861136311594053]; %Abs Hq=[ .347854845137454, .652145154862546,.652145154862546, .347854845137454]; %Wts GQ4.size=16; GQ4.xa=[Aq;Aq;Aq;Aq]; GQ4.ya=GQ4.xa'; W=[Hq;Hq;Hq;Hq]; GQ4.wt=W.*W'; % 4x4 end xa=GQ4.xa; ya=GQ4.ya; wt=GQ4.wt; Nq=GQ4.size; global ZG3412td; if (isempty(ZG3412td)|size(ZG3412td,2)~=Nq) % Evaluate and save/cache the set of shape functions at quadrature pts. ZG3412td=cell(nnd,Nq); for k=1:Nq for m=1:nnd ZG3412td{m,k}=Gr(nn(m,:),xa(k),ya(k)); end end end % if(isempty(*)) % --------------- End fixed data ---------------- TBi=cell(nnd); for m=1:nnd % Loop over corner nodes J=(Xe*GBL(nn(:,:),nn(m,1),nn(m,2)))'; % Gradient of bilinear function/element TBi{m} = blkdiag(1,J); end % Loop m % Jt=[x_q, x_r; y_q, y_r]; TDm=zeros(12,12); G=zeros(2,12); % Preallocate arrays for k=1:Nq % Initialize functions and derivatives at the quadrature point (xa,ya). J=(Xe*GBL(nn(:,:),xa(k),ya(k)))'; % Jacobian at (xa,ya) 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; mm = 0; for m=1:nnd G(:,mm+NDT)=Ji*ZG3412td{m,k}*TBi{m}; mm = mm+nT; end % loop m % Label rows by the test (weight) function index, columns by trial function index? % Submatrix ordered by T,Tx,Ty TDm = TDm + (G(1,:)'*G(1,:) + G(2,:)'*G(2,:))*wt(k); end % loop k TDm = TDm*Det; gf=zeros(nfT,1); mm=0; for m=1:nnd gf(mm+NDT)=(nn2tft(1,Elcon(m))-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; function G=Gr(ni,q,r) %GR gradient array of scalar simple-cubic basis functions. qi=ni(1); q0=q*ni(1); ri=ni(2); r0=r*ni(2); % array of irrotational vectors G=[ 1/8*qi*(1+r0)*(r0*(1-r0)+3*(1-q^2)), -1/8*(1+r0)*(1+q0)*(1-3*q0), ... -1/8*qi/ri*(1-r^2)*(1+r0); 1/8*ri*(1+q0)*(q0*(1-q0)+3*(1-r^2)), ... -1/8/qi*ri*(1-q^2)*(1+q0), -1/8*(1+q0)*(1+r0)*(1-3*r0)]; 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;