CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Wiki > PFV convection matrix

PFV convection matrix

From CFD-Wiki

(Difference between revisions)
Jump to: navigation, search
(Created page with "Matlab function '''CMatW.m''' for pressure-free velocity convection matrix. <pre> function [Cm,RowNdx,ColNdx]=CMatW(Xe,Elcon,nn2nft,Vdof) %CMATW - Returns the affine-mapped ele...")
(Update to MatLab version R2012b and generalize to handle rectangular and triangular elements. Major revision.)
 
Line 2: Line 2:
<pre>
<pre>
-
function [Cm,RowNdx,ColNdx]=CMatW(Xe,Elcon,nn2nft,Vdof)
+
function [Cm,RowNdx,ColNdx,Rcm,RcNdx]=CMatW(eu,Xe,Elcon,nn2nft,Vdof)
-
%CMATW - Returns the affine-mapped element convection matrix for the simple cubic Hermite  
+
%CMatW - Returns the element convection matrix for the Hermite basis
-
basis functions on 4-node straight-sided quadrilateral elements with 3 DOF
+
%  functions with 3, 4, or 6 degrees-of-freedom and defined on a
-
per node using Gauss quadrature on the reference square and row/col indices.  
+
%  3-node (triangle) or 4-node (quadrilateral) element by the class
-
% The columns of the array Vdof must contain the 3 nodal degree-of-freedom  
+
instance es using Gauss quadrature on the reference element.  
-
vectors in the proper nodal order.  
+
% The columns of the array Vdof must contain the 3, 4, or 6 nodal  
-
% The degrees of freedom in Vdof are the stream function and the two components
+
degree-of-freedom vectors in the proper nodal order.  
-
%  u and v of the solenoidal velocity vector.
+
% The degrees of freedom in Vdof are the stream function, the two components
-
% The assumed nodal numbering starts with 1 at the lower left corner of the  
+
%  u and v of the solenoidal velocity vector, and possibly the second
-
element and proceeds counter-clockwise around the element.  
+
derivatives Pxx, Pxy, Pyy of the stream function.
%
%
% Usage:
% Usage:
-
%  [CM,Rndx,Cndx] = CMatW(Xe,Elcon,nn2nft,Vdof)
+
%  [CM,Rndx,Cndx] = CMatW(es,Xe,Elcon,nn2nft,Vdof)
 +
%  [CM,Rndx,Cndx,Rcm,RcNdx] = CMatW(es,Xe,Elcon,nn2nft,Vdof)
 +
%  eu      - handle for basis function definitions
%  Xe(1,:) -  x-coordinates of corner nodes of element.   
%  Xe(1,:) -  x-coordinates of corner nodes of element.   
%  Xe(2,:) -  y-coordinates of corner nodes of element.   
%  Xe(2,:) -  y-coordinates of corner nodes of element.   
%  Elcon - this element connectivity matrix  
%  Elcon - this element connectivity matrix  
%  nn2nft - global number and type of DOF at each node  
%  nn2nft - global number and type of DOF at each node  
-
%  Vdof - (3x4) array of stream function, velocity components, and second
+
%  Vdof   - (ndfxnnd) array of stream function, velocity components, and  
-
%    stream function derivatives to specify the velocity over the element.
+
%    perhaps second stream function derivatives to specify the velocity  
 +
%    over the element.
%
%
-
% Jonas Holdeman, August 2007, revised June 2011
+
% Indirectly may use (handle passed by eu):
 +
%  GQuad2  - function providing 2D rectangle quadrature rules.
 +
%  TQuad2  - function providing 2D triangle quadrature rules.
 +
%
 +
% Jonas Holdeman, August 2007, revised March 2013
-
% Constants and fixed data
+
% ----------------- Constants and fixed data -------------------------
-
nd = 3;  nd4=4*nd;  ND=1:nd;    % nd = number of dofs per node,
+
-
nn=[-1 -1; 1 -1; 1 1; -1 1];    % defines local nodal order
+
-
     
+
-
% Define 5-point quadrature data once, on first call.
+
-
% Gaussian weights and absissas to integrate 9th degree polynomials exactly.
+
-
global GQ5;
+
-
if (isempty(GQ5))  % 5-point quadrature data defined? If not, define arguments & weights.
+
-
  Aq=[-.906179845938664,-.538469310105683, .0,              .538469310105683, .906179845938664];
+
-
  Hq=[ .236926885056189, .478628670499366, .568888888888889, .478628670499366, .236926885056189];
+
-
  GQ5.size=25; GQ5.xa=[Aq;Aq;Aq;Aq;Aq]; GQ5.ya=GQ5.xa';
+
-
  wt=[Hq;Hq;Hq;Hq;Hq]; GQ5.wt=wt.*wt';
+
-
end
+
-
  xa=GQ5.xa; ya=GQ5.ya; wt=GQ5.wt; Nq=GQ5.size; % Use GQ5 (5x5) for exact integration
+
nnodes = eu.nnodes;          % number of nodes per element (4);
 +
nndofs = eu.nndofs;          % nndofs = number of dofs per node, (3|6); 
 +
nedofs = nnodes*nndofs; 
 +
nn = eu.nn;          % defines local nodal order, [-1 -1; 1 -1; 1 1; -1 1]
 +
 +
% ------------------------------------------------------------------------
 +
persistent QQCM4;
 +
if isempty(QQCM4)
 +
    QRorder = 3*eu.mxpowr-1; % =9
 +
    [QQCM4.xa, QQCM4.ya, QQCM4.wt, QQCM4.nq] = eu.hQuad(QRorder);
 +
end  % if isempty...
 +
xa = QQCM4.xa; ya = QQCM4.ya; wt = QQCM4.wt; Nq = QQCM4.nq;
 +
% ------------------------------------------------------------------------
-
% -----------------------------------------------
+
persistent ZZ_Sc; persistent ZZ_SXc; persistent ZZ_SYc;  
-
global Zs3412D2c; global ZS3412c;
+
if (isempty(ZZ_Sc)||size(ZZ_Sc,2)~=Nq)
-
 
+
-
if (isempty(ZS3412c)|isempty(Zs3412D2c)|size(Zs3412D2c,2)~=Nq)
+
% Evaluate and save/cache the set of shape functions at quadrature pts.  
% Evaluate and save/cache the set of shape functions at quadrature pts.  
-
  Zs3412D2c=cell(4,Nq); ZS3412c=cell(4,Nq);
+
ZZ_Sc=cell(nnodes,Nq); ZZ_SXc=cell(nnodes,Nq); ZZ_SYc=cell(nnodes,Nq);  
   for k=1:Nq
   for k=1:Nq
-
       for m=1:4
+
       for m=1:nnodes
-
       ZS3412c{m,k}= Sr(nn(m,:),xa(k),ya(k));
+
       ZZ_Sc{m,k}= eu.S(nn(m,:),xa(k),ya(k));
-
       Zs3412D2c{m,k}=D3s(nn(m,:),xa(k),ya(k));
+
       [ZZ_SXc{m,k},ZZ_SYc{m,k}]=eu.DS(nn(m,:),xa(k),ya(k));
       end
       end
   end
   end
-
end  % if(isempty(*))
+
end  % if isempty...
-
% --------------- End fixed data ----------------
+
% ----------------------- End fixed data ---------------------------------
 +
 +
if (nargout<4), ItrType=0; else ItrType=1; end  % ItrType=1 for Newton iter
 +
affine = eu.isaffine(Xe);        % affine?
 +
%affine = (sum(abs(Xe(:,1)-Xe(:,2)+Xe(:,3)-Xe(:,4)))<4*eps);    % affine?
-
Ti=cell(4);
+
Ti=cell(nnodes);
-
%  
+
% Jt=[x_q, x_r; y_q, y_r];
-
for m=1:4
+
if affine  % (J constant)
-
   Jt=Xe*GBL(nn(:,:),nn(m,1),nn(m,2));  
+
   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)
   JtiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)];  % det(J)*inv(J)
-
   Ti{m}=blkdiag(1,JtiD);   
+
   if nndofs==3
-
end
+
    TT=blkdiag(1,JtiD);  
-
 
+
   elseif nndofs==4
-
% Move Jacobian evaluation inside k-loop for general convex quadrilateral.  
+
      TT=blkdiag(1,JtiD,Jt(1,1)*Jt(2,2)+Jt(1,2)*Jt(2,1));
-
% Jt=[x_q, x_r; y_q, y_r];
+
  else
-
 
+
    T2=[Jt(1,1)^2, 2*Jt(1,1)*Jt(1,2), Jt(1,2)^2; ...  % alt
-
Cm=zeros(nd4,nd4); Rcm=zeros(nd4,1);
+
    Jt(1,1)*Jt(2,1), Jt(1,1)*Jt(2,2)+Jt(1,2)*Jt(2,1), Jt(1,2)*Jt(2,2); ...
-
S=zeros(2,nd4); Sx=zeros(2,nd4); Sy=zeros(2,nd4);  % Pre-allocate arrays
+
    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
 +
    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; ...  % 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(:,:),nn(m,1),nn(m,2)); % Second cross derivatives
 +
      TT(5,2)= Bxy(2);
 +
      TT(5,3)=-Bxy(1);
 +
    end  % nndofs...
 +
    Ti{m}=TT;
 +
  end  % for m=...
 +
end  % affine
 +
% Allocate arrays
 +
Cm=zeros(nedofs,nedofs); Rcm=zeros(nedofs,1);  
 +
S=zeros(2,nedofs); Sx=zeros(2,nedofs); Sy=zeros(2,nedofs);   
 +
ND=1:nndofs;
% Begin loop over Gauss-Legendre quadrature points.  
% Begin loop over Gauss-Legendre quadrature points.  
for k=1:Nq   
for k=1:Nq   
-
    
+
   if ~affine
-
  Jt=Xe*GBL(nn(:,:),xa(k),ya(k));         % transpose of Jacobian at (xa,ya)
+
    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  
+
    Det=Jt(1,1)*Jt(2,2)-Jt(1,2)*Jt(2,1); % Determinant of Jt & J  
-
  Jtd=Jt/Det;
+
    Jtd=Jt/Det;
-
  TL=[Jt(2,2)^2, -2*Jt(2,1)*Jt(2,2),  Jt(2,1)^2; ...
+
    Ji=[Jt(2,2),-Jt(2,1); -Jt(1,2),Jt(1,1)]/Det;
-
    -Jt(1,2)*Jt(2,2), Jt(1,1)*Jt(2,2)+Jt(2,1)*Jt(1,2), -Jt(1,1)*Jt(2,1); ...
+
  end
-
      Jt(1,2)^2, -2*Jt(1,1)*Jt(1,2),  Jt(1,1)^2 ]/Det^2;
+
    
    
% Initialize functions and derivatives at the quadrature point (xa,ya).
% Initialize functions and derivatives at the quadrature point (xa,ya).
-
   for m=1:4
+
   for m=1:nnodes
-
     mm=nd*(m-1);
+
     mm=nndofs*(m-1);
-
     S(:,mm+ND) = Jtd*ZS3412c{m,k}*Ti{m};
+
     S(:,mm+ND)=Jtd*ZZ_Sc{m,k}*Ti{m};
-
    Ds = TL*Zs3412D2c{m,k}*Ti{m};  
+
     Sx(:,mm+ND)=Jtd*(Ji(1,1)*ZZ_SXc{m,k}+Ji(1,2)*ZZ_SYc{m,k})*Ti{m};
-
     Sx(:,mm+ND) = [Ds(2,:); -Ds(1,:)];    % [Pyx, -Pxx]
+
     Sy(:,mm+ND)=Jtd*(Ji(2,1)*ZZ_SXc{m,k}+Ji(2,2)*ZZ_SYc{m,k})*Ti{m};
-
     Sy(:,mm+ND) = [Ds(3,:); -Ds(2,:)];    % [Pyy, -Pxy]
+
   end  % loop m
   end  % loop m
    
    
-
% Compute the fluid velocity at the quadrature point.
+
  % Compute the fluid velocity at the quadrature point.
   U = S*Vdof(:);
   U = S*Vdof(:);
 +
 
% Submatrix ordered by psi,u,v  
% Submatrix ordered by psi,u,v  
-
   Cm = Cm + S'*(U(1)*Sx+U(2)*Sy)*(wt(k)*Det);
+
   Cm = Cm + S'*(U(1)*Sx+U(2)*Sy)*(wt(k)*Det);  
 +
 
 +
  if (ItrType~=0)      % iteration type is Newton 
 +
    Ux=Sx*Vdof(:);
 +
    Uy=Sy*Vdof(:);
 +
    Cm = Cm + S'*(Ux*S(1,:)+Uy*S(2,:))*(wt(k)*Det);
 +
    Rcm=Rcm + S'*(U(1)*Ux+U(2)*Uy)*(wt(k)*Det);
 +
  end              % Cm & Rcm for Newton iteration
end    % end loop k over quadrature points  
end    % end loop k over quadrature points  
-
gf=zeros(nd4,1);
+
gf=zeros(nedofs,1);
m=0;  
m=0;  
-
for n=1:4                 % Loop over element nodes  
+
for n=1:nnodes                 % Loop over element nodes  
-
   gf(m+ND)=(nn2nft(1,Elcon(n))-1)+ND;  % Get global freedoms
+
   gf(m+ND)=(nn2nft(Elcon(n),1)-1)+ND;  % Get global freedoms
-
   m=m+nd;
+
   m=m+nndofs;
end
end
-
RowNdx=repmat(gf,1,nd4);     % Row indices
+
RowNdx=repmat(gf,1,nedofs);   % Row indices
ColNdx=RowNdx';              % Col indices
ColNdx=RowNdx';              % Col indices
   
   
-
Cm = reshape(Cm,nd4*nd4,1);
+
Cm = reshape(Cm,nedofs*nedofs,1);
-
RowNdx=reshape(RowNdx,nd4*nd4,1);
+
RowNdx=reshape(RowNdx,nedofs*nedofs,1);
-
ColNdx=reshape(ColNdx,nd4*nd4,1);  
+
ColNdx=reshape(ColNdx,nedofs*nedofs,1);  
-
return;
+
if(ItrType==0), Rcm=[]; RcNdx=[]; else RcNdx=gf; end
-
 
+
-
% ----------------------------------------------------------------------------
+
-
 
+
-
function P2=D3s(ni,q,r)
+
-
% Second derivatives [Pxx; Pxy; Pyy] of simple cubic stream function.
+
-
qi=ni(1); q0=q*ni(1); q1=1+q0;
+
-
ri=ni(2); r0=r*ni(2); r1=1+r0;
+
-
  P2=[-.75*qi^2*(r0+1)*q0, 0, -.25*qi*(r0+1)*(3*q0+1); ...
+
-
      .125*qi*ri*(4-3*(q^2+r^2)), .125*qi*(r0+1)*(3*r0-1), -.125*ri*(q0+1)*(3*q0-1); ...
+
-
      -.75*ri^2*(q0+1)*r0, .25*ri*(q0+1)*(3*r0+1), 0] ;
+
-
return;
+
-
 
+
-
function S=Sr(ni,q,r)
+
-
%S  Array of solenoidal basis functions on rectangle.
+
-
  qi=ni(1); q0=q*ni(1); q1=1+q0;
+
-
  ri=ni(2); r0=r*ni(2); r1=1+r0;
+
-
  % array of solenoidal vectors
+
-
S=[ .125*ri*q1*(q0*(1-q0)+3*(1-r^2)), .125*q1*r1*(3*r0-1),    .125*ri/qi*q1^2*(1-q0); ...
+
-
      -.125*qi*r1*(r0*(1-r0)+3*(1-q^2)), .125*qi/ri*r1^2*(1-r0), .125*q1*r1*(3*q0-1)];
+
-
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;
return;
</pre>
</pre>

Latest revision as of 16:06, 15 March 2013

Matlab function CMatW.m for pressure-free velocity convection matrix.

function [Cm,RowNdx,ColNdx,Rcm,RcNdx]=CMatW(eu,Xe,Elcon,nn2nft,Vdof)
%CMatW - Returns the element convection 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. 
% The columns of the array Vdof must contain the 3, 4, or 6 nodal 
%   degree-of-freedom vectors in the proper nodal order. 
% The degrees of freedom in Vdof are the stream function, the two components
%   u and v of the solenoidal velocity vector, and possibly the second 
%   derivatives Pxx, Pxy, Pyy of the stream function.
%
% Usage:
%   [CM,Rndx,Cndx] = CMatW(es,Xe,Elcon,nn2nft,Vdof)
%   [CM,Rndx,Cndx,Rcm,RcNdx] = CMatW(es,Xe,Elcon,nn2nft,Vdof)
%   eu      - handle for basis function definitions
%   Xe(1,:) -  x-coordinates of corner nodes of element.  
%   Xe(2,:) -  y-coordinates of corner nodes of element.  
%   Elcon - this element connectivity matrix 
%   nn2nft - global number and type of DOF at each node 
%   Vdof    - (ndfxnnd) array of stream function, velocity components, and 
%     perhaps second stream function derivatives to specify the velocity 
%     over the element.
%
% Indirectly may use (handle passed by eu):
%   GQuad2   - function providing 2D rectangle quadrature rules.
%   TQuad2   - function providing 2D triangle quadrature rules.
%
% Jonas Holdeman, August 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;  
nn = eu.nn;           % defines local nodal order, [-1 -1; 1 -1; 1 1; -1 1]
 
% ------------------------------------------------------------------------
persistent QQCM4; 
if isempty(QQCM4) 
     QRorder = 3*eu.mxpowr-1; % =9
    [QQCM4.xa, QQCM4.ya, QQCM4.wt, QQCM4.nq] = eu.hQuad(QRorder);
end  % if isempty...
xa = QQCM4.xa; ya = QQCM4.ya; wt = QQCM4.wt; Nq = QQCM4.nq;
% ------------------------------------------------------------------------

persistent ZZ_Sc; persistent ZZ_SXc; persistent ZZ_SYc; 
if (isempty(ZZ_Sc)||size(ZZ_Sc,2)~=Nq)
% Evaluate and save/cache the set of shape functions at quadrature pts. 
ZZ_Sc=cell(nnodes,Nq); ZZ_SXc=cell(nnodes,Nq); ZZ_SYc=cell(nnodes,Nq); 
   for k=1:Nq
      for m=1:nnodes
      ZZ_Sc{m,k}= eu.S(nn(m,:),xa(k),ya(k));
      [ZZ_SXc{m,k},ZZ_SYc{m,k}]=eu.DS(nn(m,:),xa(k),ya(k));
      end
   end
end  % if isempty...

% ----------------------- End fixed data ---------------------------------
 
if (nargout<4), ItrType=0; else ItrType=1; end  % ItrType=1 for Newton iter
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
    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; ...  % 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(:,:),nn(m,1),nn(m,2)); % Second cross derivatives
      TT(5,2)= Bxy(2);
      TT(5,3)=-Bxy(1);
    end  % nndofs...
    Ti{m}=TT;
  end  % for m=...
end  % affine
% Allocate arrays
Cm=zeros(nedofs,nedofs); Rcm=zeros(nedofs,1); 
S=zeros(2,nedofs); Sx=zeros(2,nedofs); Sy=zeros(2,nedofs);  
ND=1:nndofs; 
% Begin loop over Gauss-Legendre quadrature points. 
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);
    S(:,mm+ND)=Jtd*ZZ_Sc{m,k}*Ti{m};
    Sx(:,mm+ND)=Jtd*(Ji(1,1)*ZZ_SXc{m,k}+Ji(1,2)*ZZ_SYc{m,k})*Ti{m};
    Sy(:,mm+ND)=Jtd*(Ji(2,1)*ZZ_SXc{m,k}+Ji(2,2)*ZZ_SYc{m,k})*Ti{m};
  end  % loop m
  
  % Compute the fluid velocity at the quadrature point.
  U = S*Vdof(:);
   
% Submatrix ordered by psi,u,v 
  Cm = Cm + S'*(U(1)*Sx+U(2)*Sy)*(wt(k)*Det); 

  if (ItrType~=0)       % iteration type is Newton  
    Ux=Sx*Vdof(:);
    Uy=Sy*Vdof(:);
    Cm = Cm + S'*(Ux*S(1,:)+Uy*S(2,:))*(wt(k)*Det);
    Rcm=Rcm + S'*(U(1)*Ux+U(2)*Uy)*(wt(k)*Det);
  end               % Cm & Rcm for Newton iteration 
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
 
Cm = reshape(Cm,nedofs*nedofs,1);
RowNdx=reshape(RowNdx,nedofs*nedofs,1);
ColNdx=reshape(ColNdx,nedofs*nedofs,1);   
if(ItrType==0), Rcm=[]; RcNdx=[]; else RcNdx=gf; end

return;
My wiki