# PFVT diffusion matrix

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search
```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;
```