PFVT diffusion matrix
From CFD-Wiki
Revision as of 20:55, 15 March 2013 by Jonas Holdeman (Talk | contribs)
function [Dm,RowNdx,ColNdx]=DMatW(eu,Xe,Elcon,nn2nft) %DMatW - Returns the element diffusion matrix for the Hermite basis % functions with 3, 4, or 6 degrees-of-freedom and defined on a % 3-node (triangle) or 4-node (quadrilateral) element by the class % instance es using Gauss quadrature on the reference element. % % Usage: % [Dm,Rndx,Cndx] = DMatW(Xe,Elcon,nn2nft,es) % es - reference for basis function definitions % Xe(1,:) - x-coordinates of corner nodes of element. % Xe(2,:) - y-coordinates of corner nodes of element. % Elcon - connectivity matrix for this element. % nn2nft - global DOF and type of DOF at each node % % Indirectly may use (handle passed by eu): % GQuad2 - function providing 2D rectangle quadrature rules. % TQuad2 - function providing 2D triangle quadrature rules. % % Jonas Holdeman, January 2007, revised March 2013 % ------------------- Constants and fixed data --------------------------- nnodes = eu.nnodes; % number of nodes per element (4); nndofs = eu.nndofs; % nndofs = number of dofs per node, (3|6); nedofs=nnodes*nndofs; % nndofs = number of dofs per node, nn = eu.nn; % defines local nodal order, [-1 -1; 1 -1; 1 1; -1 1] % ------------------------------------------------------------------------ persistent QQDM4; if isempty(QQDM4) QRorder = 2*(eu.mxpowr-1)+1; % =9 [QQDM4.xa, QQDM4.ya, QQDM4.wt, QQDM4.nq] = eu.hQuad(QRorder); end % if isempty... xa = QQDM4.xa; ya = QQDM4.ya; wt = QQDM4.wt; Nq = QQDM4.nq; % ------------------------------------------------------------------------ persistent ZZ_SXd; persistent ZZ_SYd; if (isempty(ZZ_SXd)||isempty(ZZ_SYd)||size(ZZ_SXd,2)~=Nq) % Evaluate and save/cache the set of shape functions at quadrature pts. ZZ_SXd=cell(nnodes,Nq); ZZ_SYd=cell(nnodes,Nq); for k=1:Nq for m=1:nnodes [ZZ_SXd{m,k},ZZ_SYd{m,k}]=eu.DS(nn(m,:),xa(k),ya(k)); end end end % if(isempty(*)) % -------------------------- End fixed data ------------------------------ affine = eu.isaffine(Xe); % affine? %affine = (sum(abs(Xe(:,1)-Xe(:,2)+Xe(:,3)-Xe(:,4)))<4*eps); % affine? Ti=cell(nnodes); % Jt=[x_q, x_r; y_q, y_r]; if affine % (J constant) Jt=Xe*eu.Gm(nn(:,:),eu.cntr(1),eu.cntr(2)); JtiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)]; % det(J)*inv(J) if nndofs==3 TT=blkdiag(1,JtiD); elseif nndofs==4 TT=blkdiag(1,JtiD,Jt(1,1)*Jt(2,2)+Jt(1,2)*Jt(2,1)); 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]; TT=blkdiag(1,JtiD,T2); Bxy=Xe*eu.DGm(nn(:,:),0,0); % Second cross derivatives TT(5,2)= Bxy(2); TT(5,3)=-Bxy(1); end % nndofs... for m=1:nnodes, Ti{m}=TT; 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:nnodes % Loop over corner nodes 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 nndofs==3, TT=blkdiag(1,JtiD); elseif nndofs==4 TT=blkdiag(1,JtiD,(Jt(1,1)*Jt(2,2)+Jt(1,2)*Jt(2,1))); 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]; TT=blkdiag(1,JtiD,T2); Bxy=Xe*eu.DGm(nn(:,:),nn(m,1),nn(m,2)); % 2nd cross derivatives TT(5,2)= Bxy(2); TT(5,3)=-Bxy(1); end Ti{m}=TT; end % Loop m end % Allocate arrays Dm=zeros(nedofs,nedofs); Sx=zeros(2,nedofs); Sy=zeros(2,nedofs); ND=1:nndofs; for k=1:Nq if ~affine Jt=Xe*eu.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 functions and derivatives at the quadrature point (xa,ya). for m=1:nnodes mm=nndofs*(m-1); Sx(:,mm+ND)=Jtd*(Ji(1,1)*ZZ_SXd{m,k}+Ji(1,2)*ZZ_SYd{m,k})*Ti{m}; Sy(:,mm+ND)=Jtd*(Ji(2,1)*ZZ_SXd{m,k}+Ji(2,2)*ZZ_SYd{m,k})*Ti{m}; end % loop m Dm = Dm+(Sx'*Sx+Sy'*Sy)*(wt(k)*Det); end % end loop k over quadrature points gf=zeros(nedofs,1); m=0; for n=1:nnodes % Loop over element nodes gf(m+ND)=(nn2nft(Elcon(n),1)-1)+ND; % Get global freedoms m=m+nndofs; end RowNdx=repmat(gf,1,nedofs); % Row indices ColNdx=RowNdx'; % Col indices Dm = reshape(Dm,nedofs*nedofs,1); RowNdx=reshape(RowNdx,nedofs*nedofs,1); ColNdx=reshape(ColNdx,nedofs*nedofs,1); return;