PFV 3D convection matrix

From CFD-Wiki

Function CMat3D8W.m for pressure-free velocity convection matrix

function [Cm,RowNdx,ColNdx]=CMat3D8W(Xe, Elcon, nn2nft, Vdof)
% CMat3D8W - Returns the element convection matrix for the  
%   3D linear-complete, normal-conforming, divergence-free, Hermite basis 
%   functions on 8-node rectangular hexahedral elements with 6 DOF per node 
%   using Gauss quadrature on the 2x2x2 reference cube. 
% The 8 columns of the array V3dof must contain the six degree-of-freedom 
%   vectors in the nodal order (A,B,C,u,v,w). 
% The assumed nodal numbering starts with 1 at the lower left corner (-1,-1,-1) 
%   of the element. 
%
% Usage:
%   [Cm,Rndx,Cndx] = CMat3D8W(Xe, Elcon, nn2nft,V3dof)
%   [Cm,Rndx,Cndx,Rcm,RcNdx] = CMat3D8W(Xe, Elcon, nn2nft,V3dof)
%   Xe(1,:) -  x-coordinates of 8 corner nodes of element.  
%   Xe(2,:) -  y-coordinates of 8 corner nodes of element.  
%   Xe(3,:) -  z-coordinates of 8 corner nodes of element.  
%   Elcon(8) - connectivity matrix for this element, list of nodes. 
%   nn2nft(1,n) - global freedom number for node n.
%   nn2nft(2,n) - global freedom type for node n.
%   Vdof(6,8)  - VP & velocity Dofs at 8 nodes.
%
% Calls:
%   V8cW(nc,x,y,z), V8xyzcW(nc,x,y,z)
%
% Jonas Holdeman, July 2011
%

% Constants and fixed data
nc=[-1,-1,-1; 1,-1,-1; -1,1,-1; 1,1,-1; -1,-1,1; 1,-1,1; -1,1,1; 1,1,1]; % defines corner nodal order
ndfn=6;        % number of degrees of freedom per node. 
nne=8;         % number of nodes per element
ndfe=ndfn*nne;  % number of degrees of freedom per element

% Define 5-point quadrature data once, on first call. 
% Gaussian weights and absissas to integrate 9th degree polynomials exactly. 
global GQC5;
if (isempty(GQC5))   % Has 5-point quadrature data been defined? If not, define arguments & weights. 
   Aq=[-.906179845938664,-.538469310105683, .0,               .538469310105683, .906179845938664];
   Hq=[ .236926885056189, .478628670499366, .568888888888889, .478628670499366, .236926885056189];
   GQC5.xa=zeros(125,1); GQC5.ya=zeros(125,1); GQC5.za=zeros(125,1); GQC5.wt=zeros(125,1); 
   nr=0; 
   for nz=1:5; for ny=1:5; for nx=1:5
       nr=nr+1; GQC5.xa(nr)=Aq(nx); GQC5.ya(nr)=Aq(ny); GQC5.za(nr)=Aq(nz); 
       GQC5.wt(nr)=Hq(nx)*Hq(ny)*Hq(nz);
   end; end; end
   GQC5.size=nr; 
end

xa=GQC5.xa; ya=GQC5.ya; za=GQC5.za; W=GQC5.wt; Nq=GQC5.size;

% ---------------------------------------------------
global Z3_V8c; global Z3_V8xc; global Z3_V8yc; global Z3_V8zc;
if (isempty(Z3_V8c) | isempty(Z3_V8xc) | size(Z3_V8xc,2)~=Nq)
  Z3_V8c=cell(nne,Nq); Z3_V8xc=cell(nne,Nq); 
  Z3_V8yc=cell(nne,Nq); Z3_V8zc=cell(nne,Nq);
  for k=1:Nq
    for m=1:nne
      ncm=nc(m,:);
      Z3_V8c{m,k}=V8cW(ncm,xa(k),ya(k),za(k));
      [Z3_V8xc{m,k},Z3_V8yc{m,k},Z3_V8zc{m,k}]=V8xyzcW(ncm,xa(k),ya(k),za(k));
    end
  end
end   % if (isempty(*))
% ----------------- End fixed data ------------------
 
Ti=cell(nne);
for m=1:nne
  Jt=Xe*GTL(nc(:,:),nc(m,1),nc(m,2),nc(m,3));
  Det=det(Jt);
  JtiD=inv(Jt)*Det;
  J=Jt';
  Ti{m}=blkdiag(J,JtiD);
end  % loop m

Cm=zeros(ndfe,ndfe); S=zeros(3,ndfe);   % Preallocate arrays

% Begin loop over Gauss-Legendre quadrature points. 
for k=1:Nq  
   
  Jt=Xe*GTL(nc(:,:),xa(k),ya(k),za(k));    % transpose of Jacobian at (xa,ya,za)
  Det=det(Jt);
  JtbD=Jt/Det;
  Jti=inv(Jt); 
  JtiD=Jti*Det; 
  Ji=Jti';
%
% Compute mapped element Si and the fluid velocity at the quadrature point (xa,ya,za).
  Ua=[0;0;0]; 
  for m=1:nne         % velocity & derivatives at 8 corner nodes
    mm=ndfn*(m-1);  mm3=mm+1:mm+ndfn;  ncm=nc(m,:);
    S(:,mm3)= JtbD*Z3_V8c{m,k}*Ti{m};  
    Ua = Ua + S(:,mm+1:mm+ndfn)*Vdof(:,m); % A,B,C,u,v,w
  end
  Ub=Jti*Ua;
  UgS=zeros(3,ndfe);
  for m=1:nne         % velocity & derivatives at 8 corner nodes
    mm=ndfn*(m-1);  mm3=mm+1:mm+ndfn; 
    UgS(:,mm3)=JtbD*(Ub(1)*Z3_V8xc{m,k}+Ub(2)*Z3_V8yc{m,k}+Ub(3)*Z3_V8zc{m,k})*Ti{m};
  end 
  Cm = Cm + S'*UgS*W(k)*Det;
  
end    % loop k over quadrature points 

gf=zeros(ndfe,1);
m=0;
for k=1:nne
  m=m+1; gf(m)=nn2nft(1,Elcon(k));  % get global freedom number 
  for k1=2:ndfn
    m=m+1; gf(m)=gf(m-1)+1;  % next 
  end  % if
end %  loop on k
RowNdx=repmat(gf,1,ndfe);
ColNdx=RowNdx';
RowNdx=reshape(RowNdx,ndfe*ndfe,1);
ColNdx=reshape(ColNdx,ndfe*ndfe,1); 
Cm=reshape(Cm,ndfe*ndfe,1);
 
return;

% -----------------------------------------------------------------------

function G=GTL(ni,q,r,s)
% Transposed gradient (derivatives) of scalar trilinear mapping function. 
% The parameter ni can be a vector of coordinate pairs. 
G=[.125*ni(:,1).*(1+ni(:,2).*r).*(1+ni(:,3).*s), .125*ni(:,2).*(1+ni(:,1).*q).*(1+ni(:,3).*s), ... 
    .125*ni(:,3).*(1+ni(:,1).*q).*(1+ni(:,2).*r)];
return;