CFD Online Logo CFD Online URL
Home > Wiki > Code: Lid driven cavity using pressure free velocity form

Code: Lid driven cavity using pressure free velocity form

From CFD-Wiki

Revision as of 16:36, 15 March 2013 by Jonas Holdeman (Talk | contribs)
Jump to: navigation, search


Lid-driven cavity using pressure-free velocity formulation

This sample code uses four-node simple-cubic finite elements and simple iteration.


The incompressible Navier-Stokes equation is a differential algebraic equation, having the inconvenient feature that there is no explicit mechanism for advancing the pressure in time. Consequently, much effort has been expended to eliminate the pressure from all or part of the computational process. We show a simple, natural way of doing this.

The incompressible Navier-Stokes equation is composite, the sum of two orthogonal equations,

\frac{\partial\mathbf{v}}{\partial t}=\Pi^S(-\mathbf{v}\cdot\nabla\mathbf{v}+\nu\nabla^2\mathbf{v})+\mathbf{f}^S ,
\rho^{-1}\nabla p=\Pi^I(-\mathbf{v}\cdot\nabla\mathbf{v}+\nu\nabla^2\mathbf{v})+\mathbf{f}^I ,

where \Pi^S and \Pi^I are solenoidal and irrotational projection operators satisfying \Pi^S+\Pi^I=1 and \mathbf{f}^S and \mathbf{f}^I are the nonconservative and conservative parts of the body force. This result follows from the Helmholtz Theorem . The first equation is a pressureless governing equation for the velocity, while the second equation for the pressure is a functional of the velocity and is related to the pressure Poisson equation. The explicit functional forms of the projection operator in 2D and 3D are found from the Helmholtz Theorem, showing that these are integro-differential equations, and not particularly convenient for numerical computation.

Equivalent weak or variational forms of the equations, proved to produce the same velocity solution as the Navier-Stokes equation are

(\mathbf{w},\frac{\partial\mathbf{v}}{\partial t})=-(\mathbf{w},\mathbf{v}\cdot\nabla\mathbf{v})-\nu(\nabla\mathbf{w}: \nabla\mathbf{v})+(\mathbf{w},\mathbf{f}^S),
(\mathbf{g}_i,\nabla p)=-(\mathbf{g}_i,\mathbf{v}\cdot\nabla\mathbf{v}_j)-\nu(\nabla\mathbf{g}_i: \nabla\mathbf{v}_j)+(\mathbf{g}_i,\mathbf{f}^I)\,,

for divergence-free test functions \mathbf{w} and irrotational test functions \mathbf{g} satisfying appropriate boundary conditions. Here, the projections are accomplished by the orthogonality of the solenoidal and irrotational function spaces. The discrete form of this is emminently suited to finite element computation of divergence-free flow.

In the discrete case, it is desirable to choose basis functions for the velocity which reflect the essential feature of incompressible flow — the velocity elements must be divergence-free. While the velocity is the variable of interest, the existence of the stream function or vector potential is necessary by the Helmholtz Theorem. Further, to determine fluid flow in the absence of a pressure gradient, one can specify the difference of stream function values across a 2D channel, or the line integral of the tangential component of the vector potential around the channel in 3D, the flow being given by Stokes' Theorem. This leads naturally to the use of Hermite stream function (in 2D) or velocity potential elements (in 3D).

Involving, as it does, both stream function and velocity degrees-of-freedom, the method might be called a velocity-stream function or stream function-velocity method.

We now restrict discussion to 2D continuous Hermite finite elements which have at least first-derivative degrees-of-freedom. With this, one can draw a large number of candidate triangular and rectangular elements from the plate-bending literature. These elements have derivatives as components of the gradient. In 2D, the gradient and curl of a scalar are clearly orthogonal, given by the expressions,

\nabla\phi = \left[\frac{\partial \phi}{\partial x},\,\frac{\partial \phi}{\partial y}\right]^T, \quad
\nabla\times\phi = \left[\frac{\partial \phi}{\partial y},\,-\frac{\partial \phi}{\partial x}\right]^T.

Adopting continuous plate-bending elements, interchanging the derivative degrees-of-freedom and changing the sign of the appropriate one gives many families of stream function elements.

Taking the curl of the scalar stream function elements gives divergence-free velocity elements [1][2]. The requirement that the stream function elements be continuous assures that the normal component of the velocity is continuous across element interfaces, all that is necessary for vanishing divergence on these interfaces.

Boundary conditions are simple to apply. The stream function is constant on no-flow surfaces, with no-slip velocity conditions on surfaces. Stream function differences across open channels determine the flow. No boundary conditions are necessary on open boundaries [1], though consistent values may be used with some problems. These are all Dirichlet conditions.

The algebraic equations to be solved are simple to set up, but of course are non-linear, requiring iteration of the linearized equations.

The finite elements we will use here are apparently due to Melosh [3], but can also be found in Zienkiewitz [4]. These simple cubic-complete elements have three degrees-of-freedom at each of the four nodes. In the sample code we use this Hermite element for the pressure, and the modified form obtained by interchanging derivatives and the sign of one of them (though a simple bilinear element could be used for the pressure as well). The degrees-of-freedom are the pressure and pressure gardient, and the stream function and components of the solenoidal velocity for the modified element. The normal component of the velocity is continuous at element interfaces as is required, but the tangential velocity component may not be continuous.

The code implementing the lid-driven cavity problem is written for Matlab. The script below is problem-specific, and calls problem-independent functions to evaluate the element diffusion and convection matricies and evaluate the pressure from the resulting velocity field. These three functions accept general quadrilateral elements with straight sides as well as the rectangular elements used here. Other functions are a GMRES iterative solver using ILU preconditioning and incorporating the essential boundary conditions, and a function to produce non-uniform nodal spacing for the problem mesh.

This "educational code" is a simplified version of the code used in [1]. The user interface is the code itself. The user can experiment with changing the mesh, the Reynolds number, and the number of nonlinear iterations performed, as well as the relaxation factor. There are suggestions in the code regarding near-optimum choices for this factor as a function of Reynolds number. These values are given in the paper as well. For larger Reynolds numbers, a smaller relaxation factor speeds up convergence by smoothing the velocity factor (\mathbf{v}\cdot\nabla) in the convection term, but will impede convergence if made too small.

The output consists of graphic plots of contour levels of the stream function and the pressure levels.

This simplified version for this Wiki resulted from removal of computation of the vorticity, a restart capability, area weighting for the error, and production of publication-quality plots from one of the research codes used with the paper.

Lid-driven cavity Matlab script

% Finite element solution of the 2D Navier-Stokes equation using 4-node, 
%  12 DOF, (3-DOF/node), simple-cubic-derived rectangular Hermite basis for 
%   the Lid-Driven Cavity problem.
% This could also be characterized as a VELOCITY-STREAM FUNCTION or 
% Reference:  "A Hermite finite element method for incompressible fluid flow", 
%    Int. J. Numer. Meth. Fluids, 64, P376-408 (2010). 
% Simplified Wiki version 
% The rectangular problem domain is defined between Cartesian 
%   coordinates Xmin & Xmax and Ymin & Ymax.
% The computational grid has NumEx elements in the x-direction 
%   and NumEy elements in the y-direction. 
% The nodes and elements are numbered column-wise from the  
%   upper left corner to the lower right corner. 
% This script calls the user-defined functions:
%   ELS3412r     - class of velocity basis functions
%   DMatW        - to evaluate element diffusion matrix 
%   CMatW        - to evaluate element convection matrix
%   GetPresW     - to evaluate the pressure 
%   regrade      - to regrade the mesh 
% Uses
%   ilu          - incomplete LU preconditioner
%   gmres        - iterative solver
% Indirectly uses:
%    Gquad2      - Gauss integraion rules for rectangle
%    ELG3412r    - class of pressure basis functions
% Jonas Holdeman   August 2007, revised March 2013

clear all;
eu = ELS3412r;       % class of functions for velocity
disp('Lid-driven cavity');
disp([' Four-node, 12 DOF, ' ' basis.']);

% ------------------------------------------------------------------------
  ndf = eu.nndofs;       % nndofs = number of velocity dofs per node, (3); 
  nnd = eu.nnodes;       % nnodes = number of element nodes
  nd2=ndf*ndf;  % Number of DOF per node - do not change!!
% ------------------------------------------------------------------------

% Parameters for GMRES solver 
% parameters for ilu preconditioner
% Decrease su.droptol if ilu preconditioner fails

% Optimal relaxation parameters for given Reynolds number
% (see IJNMF reference)
% Re          100   1000   3200   5000   7500  10000  12500 
% RelxFac:  1.04    1.11   .860   .830   .780   .778   .730 
% ExpCR1    1.488   .524   .192   .0378   --     --     -- 
% ExpCRO    1.624   .596   .390   .331   .243   .163   .133
% CritFac:  1.82    1.49   1.14  1.027   .942   .877   .804 

% Define the problem geometry, set mesh bounds:
Xmin = 0.0; Xmax = 1.0; Ymin = 0.0; Ymax = 1.0; 

% Set mesh grading parameters (set to 1 if no grading).
% See below for explanation of use of parameters. 
xgrd = .75; ygrd=.75;   % (xgrd = 1, ygrd=1 for uniform mesh) 

% Set " RefineBoundary=1 " for additional refinement at boundary, 
%  i.e., split first element along boundary into two. 

% Set number of elements in each direction
NumEx = 18;   NumEy = NumEx;

NumNx=NumEx+1;  NumNy=NumEy+1;

%   Define problem parameters: 
 % Lid velocity

 % Reynolds number

% factor for under/over-relaxation starting at iteration RelxStrt 
RelxFac = 1.;  % 

% Start with simple non-linear iteration, then switch to Newton method 
%   when non-linear corrections are less than ItThreshold. 
%   CAUTION! Large Re may require very small threshold. 
% If ItThreshold = 0 (or sufficiently small) method will never switch.
ItThreshold = .1; %1.e-2; 

% Number of nonlinear iterations
MaxNLit=12; %20; %


 % Viscosity for specified Reynolds number
% Grade the mesh spacing if desired, call regrade(x,agrd,e). 
% if e=0: refine both sides, 1: refine upper, 2: refine lower
% if agrd=xgrd|ygrd is the parameter which controls grading, then
%   if agrd=1 then leave array unaltered.
%   if agrd<1 then refine (make finer) towards the ends
%   if agrd>1 then refine (make finer) towards the center.
%  Generate equally-spaced nodal coordinates and refine if desired.
if (RefineBoundary==1)
  XNc=[XNc(1),(.62*XNc(1)+.38*XNc(2)),XNc(2:end-1),(.38*XNc(end-1) ...
  YNc=[YNc(1),(.62*YNc(1)+.38*YNc(2)),YNc(2:end-1),(.38*YNc(end-1) ...
if xgrd ~= 1, XNc=regrade(XNc,xgrd,0); end;  % Refine mesh if desired
if ygrd ~= 1, YNc=regrade(YNc,ygrd,0); end;
[Xgrid,Ygrid]=meshgrid(XNc,YNc);% Generate the x- and y-coordinate meshes.

% Allocate storage for fields 

%--------------------Begin grid plot-----------------------
% ********************** FIGURE 1 *************************
% Plot the grid 
orient portrait;  orient tall; 
hold on;
hold off;
axis equal;
axis image;
title([num2str(NumNx) 'x' num2str(NumNy) ...
      ' node mesh for Lid-driven cavity']);

%-------------- End plotting Figure 1 ----------------------

%Contour levels, Ghia, Ghia & Shin, Re=100, 400, 1000, 3200, ...
clGGS=[-.1175,-.1150,-.11,-.1,-.09,-.07,-.05,-.03,-.01,-1.e-4,-1.e-5, ...
       -1.e-7,-1.e-10,1.e-8,1.e-7,1.e-6,1.e-5,5.e-5,1.e-4,2.5e-4, ...
CL=clGGS;   % Select contour level option
if (Vlid<0), CL=-CL; end

NumNod=NumNx*NumNy;     % total number of nodes
MaxDof=ndf*NumNod;        % maximum number of degrees of freedom
EBC.Mxdof=ndf*NumNod;        % maximum number of degrees of freedom

nn2nft=zeros(NumNod,2); % node number -> nf & nt
% Generate lists of active nodal indices, freedom number & type 
ni=0;  nf=-ndf+1;  nt=1;          %   ________
for nx=1:NumNx                   %  |        |
   for ny=1:NumNy                %  |        |
      ni=ni+1;                   %  |________|
      nf=nf+ndf;               % all nodes have 4 dofs 
      nn2nft(ni,:)=[nf,nt];   % dof number & type (all nodes type 1)
%NumNod=ni;     % total number of nodes
nf2nnt=zeros(MaxDof,2);  % (node, type) associated with dof
ndof=0; dd = 1:ndf;
for n=1:NumNod
  for k=1:ndf


% Generate element connectivity, from upper left to lower right. 
ne=0;  LY=NumNy;
for nx=1:NumEx
  for ny=1:NumEy
  end  % loop on ny
end  % loop on nx

% Begin essential boundary conditions, allocate space 
MaxEBC = ndf*2*(NumNx+NumNy-2);
EBC.dof=zeros(MaxEBC,1);  % Degree-of-freedom index  
EBC.val=zeros(MaxEBC,1);  % Dof value 

 X1=XNc(2);  X2=XNc(NumNx-1);
for nf=1:MaxDof
  x=NodXY(ni,1); y=NodXY(ni,2);
  if(x==Xmin || x==Xmax || y==Ymin)
    switch nt;
    case {1, 2, 3}
      nc=nc+1; EBC.dof(nc)=nf; EBC.val(nc)=0;  % psi, u, v 
    end  % switch (type)
  elseif (y==Ymax)
    switch nt;
    case {1, 3}
      nc=nc+1; EBC.dof(nc)=nf; EBC.val(nc)=0;   % psi, v 
    case 2
      nc=nc+1; EBC.dof(nc)=nf; EBC.val(nc)=Vlid;   % u
    end  % switch (type) 
  end  % if (boundary)
end  % for nf 
if (size(EBC.dof,1)>nc)   % Truncate arrays if necessary 
end     % End ESSENTIAL (Dirichlet) boundary conditions

% partion out essential (Dirichlet) dofs
p_vec = (1:EBC.Mxdof)';         % List of all dofs
EBC.p_vec_undo = zeros(1,EBC.Mxdof);
% form a list of non-diri dofs
EBC.ndro = p_vec(~ismember(p_vec, EBC.dof));	% list of non-diri dofs
% calculate p_vec_undo to restore Q to the original dof ordering
EBC.p_vec_undo([EBC.ndro;EBC.dof]) = (1:EBC.Mxdof); %p_vec';

Q=zeros(MaxDof,1); % Allocate space for solution (dof) vector

% Initialize fields to boundary conditions
for k=1:EBC.num

errpsi=zeros(NumNy,NumNx);  % error correct for iteration

np0=zeros(1,MxNL);     % Arrays for convergence info

Dmat = spalloc(MaxDof,MaxDof,36*MaxDof);   % to save the diffusion matrix
Xe=zeros(2,nnd);      % coordinates of element corners 

NLitr=0; ND=1:ndf;
while (NLitr<MaxNLit), NLitr=NLitr+1;   % <<< BEGIN NONLINEAR ITERATION 

    % <<<<<<Newton?<<<<<<<<<<<<<<*******
if(ItType==0 && NLitr>1 && nv0(NLitr-1)<ItThreshold && ...
  ItType=1; disp(' >>> Begin Newton method >>>');
tclock=clock;   % Start assembly time <<<<<<<<<
% Generate and assemble element matrices

for ne=1:NumEl   
  if NLitr == 1    
%     Fluid element diffusion matrix, save on first iteration    
     [Emat,Rndx,Cndx] = DMatW(eu,Xe,Elcon(ne,:),nn2nft);
     Dmat=Dmat+sparse(Rndx,Cndx,Emat,MaxDof,MaxDof);  % Global diffusion mat 
   if (NLitr>1) 
%    Get stream function and velocities, loop over local element nodes
     for n=1:nnd  
%     Fluid element convection matrix, first iteration uses Stokes equation
   if ItType==0
%      Convection term for simple iteration
       [Emat,Rndx,Cndx] = CMatW(eu,Xe,Elcon(ne,:),nn2nft,Vdof);  
%      Convection term for Newton iteration
       [Emat,Rndx,Cndx,Rcm,RcNdx] = CMatW(eu,Xe,Elcon(ne,:),nn2nft,Vdof);
       RHS = RHS + sparse(RcNdx,1,Rcm,MaxDof,1); 
   end  % if ItType...
% Global convection assembly 
   end  % if NLitr...

end;  % loop ne over elements 

Mat = Mat + nu*Dmat;    % Add in cached/saved global diffusion matrix 

disp(['(' num2str(NLitr) ') Matrix assembly complete, elapsed time = '...
      num2str(etime(clock,tclock)) ' sec']);  % Assembly time <<<<<<<<<<<

Q0 = Q;  % Save dof values 

% Solve system
tclock=clock; %disp('start solution'); % Start solution time  <<<<<<<<<<<<


[Lm,Um] = ilu(Matr,su);      % incomplete LU
Qr = gmres(Matr,RHSr,GM.MIter,GM.Tol,GM.MRstrt,Lm,Um,Qs);	% GMRES

Q=[Qr;EBC.val];        % Augment active dofs with esential (Dirichlet) dofs
Q=Q(EBC.p_vec_undo);   % Restore natural order
stime=etime(clock,tclock); % Solution time <<<<<<<<<<<<<<

% ****** APPLY RELAXATION FACTOR *********************
if(NLitr>1), Q=RelxFac*Q+(1-RelxFac)*Q0; end
% ****************************************************

% Compute change and copy dofs to field arrays
dsqp=0; dsqv=0;
for k=1:MaxDof
  ni=nf2nnt(k,1); nx=NodNdx(ni,1); ny=NodNdx(ni,2);
  switch nf2nnt(k,2) % switch on dof type 
    case 1
      dsqp=dsqp+(Q(k)-Q0(k))^2; psi0(ny,nx)=Q(k);
    case 2
      dsqv=dsqv+(Q(k)-Q0(k))^2; u0(ny,nx)=Q(k);
    case 3
      dsqv=dsqv+(Q(k)-Q0(k))^2; v0(ny,nx)=Q(k);
  end  % switch on dof type 
end  % for 
np0(NLitr)=sqrt(dsqp); dP=np0(NLitr);

if (np0(NLitr)<=1e-15||nv0(NLitr)<=1e-15) 
  MaxNLit=NLitr; np0=np0(1:MaxNLit); nv0=nv0(1:MaxNLit); 
disp(['(' num2str(NLitr) ') Solution time for linear system = '...
     num2str(etime(clock,tclock)) ' sec,' ' dV = ' num2str(nv0(NLitr)) ...
     ', dP = ' num2str(np0(NLitr))]);         % Solution time <<<<<<<<<<<<
%---------- Begin plot of intermediate results ----------
% ********************** FIGURE 2 *************************

% Stream function (intermediate) 
contour(Xgrid,Ygrid,psi0,8,'k');  % Plot contours (trajectories)
title(['Lid-driven cavity,  Re=' num2str(Re)]);
axis equal; axis image;

% Plot convergence info 
xlabel('Nonlinear iteration number');
ylabel('Nonlinear correction');
axis square; 
title(['Iteration conv.,  Re=' num2str(Re)]);

% Plot nonlinear iteration correction contours 
if dP>0
contour(Xgrid,Ygrid,errpsi,8,'k');  % Plot contours (trajectories)
axis equal; axis image;
title('Iteration correction');
end  % if dP

% ********************** END FIGURE 2 *************************
%----------  End plot of intermediate results  ---------

if (nv0(NLitr)<1e-15), break; end  % Terminate iteration if non-significant 

end;   % <<< (while) END NONLINEAR ITERATION

format short g;
disp('Convergence results by iteration: velocity, stream function');
disp(['nv0:  ' num2str(nv0)]); disp(['np0:  ' num2str(np0)]); 

% >>>>>>>>>>>>>> BEGIN PRESSURE RECOVERY <<<<<<<<<<<<<<<<<<
% Essential pressure boundary condition 
% Index of node to apply pressure BC, value at node
PBCnx=fix((NumNx+1)/2);   % Apply at center of mesh
for k=1:NumNod
  if (NodNdx(k,1)==PBCnx && NodNdx(k,2)==PBCny), PBCnod=k; break; end
if (PBCnod==0), error('Pressure BC node not found'); 
  EBCp.nodn = PBCnod;  % Pressure BC node number
  EBCp.val = 0;  % set P = 0.
% Cubic pressure 
[P,Px,Py] = GetPres3W(eu,NodXY,Elcon,nn2nft,Q,EBCp,0); 
% Copy pressure to structured mesh for plotting
P =reshape(P,NumNy,NumNx);
% ******************** END PRESSURE RECOVERY *********************

% ********************** CONTINUE FIGURE 1 *************************

% Stream function    (final)
[CT,hn]=contour(Xgrid,Ygrid,psi0,CL,'k');  % Plot contours (trajectories)
hold on;
shading interp  %flat;
hold off;
axis equal;  axis image;
title(['Stream lines, ' num2str(NumNx) 'x' num2str(NumNy) ...
    ' mesh, Re=' num2str(Re)]);

% Plot pressure contours   (final)
[CT,hn]=contour(Xgrid,Ygrid,P,CPL,'k');  % Plot pressure contours
hold on;
%shading interp  %flat;
hold off;
axis equal;  axis image;
title(['Simple cubic pressure contours, Re=' num2str(Re)]);
% ********************* END FIGURE 1 *************************

disp(['Total elapsed time = '...
   num2str(etime(clock,ETstart)/60) ' min']); % Elapsed time from start <<<
beep; pause(.25); beep; pause(.25); beep;

Diffusion matrix for pressure-free velocity method (DMatW.m)

Convection matrix for pressure-free velocity method (CMatW.m)

Modified-cubic velocity element class for pressure-free velocity method (ELS3412f.m)

classdef ELS3412r < handle
    % ELS3412r
    %   Container class for 4-node modified simple-cubic Hermite finite 
    %   elements on rectangle/quadrilateral (designated by 'S3412').
    %   The vector element is used for divergence-free vector fields
    %   such as incompressible velocity and magnetic fields. 
    %  Jonas Holdeman     January 2013
    properties (Constant)
        name = 'modified simple-cubic Hermite';
        designation = 'S3412r';
        shape = 'quadrilateral';
        nsides = 4;
        order = 3;       % order of completeness
        nnodes = 4;      % number of nodes
        nndofs = 3;      % max number of nodal dofs
        tndofe = 12;     % total number of dofs for element
        mxpowr = 3;      % highest power/degree in s (for quadrature rule)
        hQuad = @GQuad2; % handle to quadrature rules on rectangles
        cntr = [0 0];    % reference element centroid
        nn = [-1 -1; 1 -1; 1 1; -1 1];  % standard nodal order of coords
    end  % properties (Constant)
    methods  (Static)
    % Four-node cubic-complete Hermite scalar stream function 
    % element on the reference square. The vector function is the curl 
    % of this simple-cubic element (S3412r) with 3 dofs per corner node.
    % Parameter ni is one of coordinate pairs (-1,-1),(1,-1),(1,1),(-1,1) 
    function s=s(ni,q,r)
       qi=ni(1); q0=q*ni(1); 
       ri=ni(2); r0=r*ni(2); 
       s=[1/8*(1+q0)*(1+r0)*(2+q0*(1-q0)+r0*(1-r0)), ...
         -1/8*ri*(1+q0)*(1+r0)*(1-r^2), 1/8*qi*(1+q0)*(1+r0)*(1-q^2)];
     end % s

    % Four-node quadratic-complete Hermite solenoidal vector function 
    % element on the reference square. The vector function is the curl 
    % of the cubic-complete element (3412) with 3 dofs per corner node.
    % Parameter ni is one of coordinate pairs (-1,-1),(1,-1),(1,1),(-1,1) 
     function S=S(ni,q,r)
        qi=ni(1); q0=q*ni(1); 
        ri=ni(2); r0=r*ni(2); 
        S=[1/8*ri*(1+q0)*(q0*(1-q0)+3*(1-r^2)), ...
          -1/8*(1+q0)*(1+r0)*(1-3*r0), 1/8*qi*ri*(1-q^2)*(1+q0); ...
          -1/8*qi*(1+r0)*(r0*(1-r0)+3*(1-q^2)), ...
           1/8*qi*ri*(1-r^2)*(1+r0), -1/8*(1+r0)*(1+q0)*(1-3*q0)];
     end % S
   % First derivatives wrt q & r of four-node quadratic-complete Hermite 
   % solenoidal vector function element on the reference square. 
   % The vector function is the curl of the cubic-complete element 
   % (3412) with 3 dofs at each corner node.
   % SQ = array of q-derivatives of solenoidal vectors 
   % SR = array of r-derivatives of solenoidal vectors 
   % Parameter ni is one of coordinate pairs (-1,-1),(1,-1),(1,1),(-1,1)
    function [SQ,SR]=DS(ni,q,r)
       qi=ni(1); q0=q*ni(1); 
       ri=ni(2); r0=r*ni(2); 
    SQ=[ 1/8*qi*ri*(4-3*q0^2-3*r0^2), -1/8*qi*(1+r0)*(1-3*r0), ...
         1/8*ri*(1+q0)*(1-3*q0); ...
         3/4*qi^2*q0*(1+r0), 0, 1/4*qi*(1+r0)*(1+3*q0)];
    SR=[-3/4*ri^2*r0*(1+q0), 1/4*ri*(1+q0)*(1+3*r0), 0; ...
        -1/8*qi*ri*(4-3*q0^2-3*r0^2), 1/8*qi*(1+r0)*(1-3*r0), ...
    end  % DS

   % Second derivatives wrt {qq,qr,rr} of four-node quadratic-complete 
   % Hermite solenoidal vector function element on the reference square. 
    function [Sqq,Sqr,Srr]=D2S(ni,q,r)
       qi=ni(1); q0=q*ni(1); 
       ri=ni(2); r0=r*ni(2); 
       Sqq=[-3/4*qi^2*ri*q0, 0, -1/4*ri*qi*(1+3*q0); ...
             3/4*qi^3*(1+r0), 0, 3/4*qi^2*(1+r0)];
       Sqr=[-3/4*qi*ri^2*r0,1/4*qi*ri*(1+3*r0), 0; ...
             3/4*qi^2*ri*q0, 0, 1/4*qi*ri*(1+3*q0)];
       Srr=[-3/4*ri^3*(1+q0), 3/4*ri^2*(1+q0), 0; ...
             3/4*qi*ri^2*r0,-1/4*qi*ri*(1+3*r0), 0];
    end  % D2S
   % Transpose of the Jacobian matrix at (q,r)
   function Jt=JacT(Xe,q,r)
   end  % JacT
   % Test to see if transformation is affine, returns True or False
   function isit=isaffine(Xe)
   end  % isaffine
   % Post-multiplying matrix Ti^-1
   function ti=Ti(Xe,m)
      JtiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)];   % det(J)*inv(J)
   end  % Ti

   % Bilinear mapping function from (q,r) in the reference square 
   % [-1.1]x[-1,1] to (x,y) in the straight-sided quadrilateral 
   % finite elements. 
   % The parameter ni can be a vector of coordinate pairs. 
   function g=gm(ni,q,r)
   end  % gm
   % Transposed gradient (derivatives) of scalar bilinear mapping function
   % The parameter ni can be a vector of coordinate pairs. 
   function G=Gm(ni,q,r)
      G=[.25*ni(:,1).*(1+ni(:,2).*r), .25*ni(:,2).*(1+ni(:,1).*q)];
   end  % Gm

   % Second (cross) derivative of scalar bilinear mapping function
   % The parameter ni can be a vector of coordinate pairs. 
   function D=DGm(ni,~,~)
      D = .25*ni(:,1).*ni(:,2); 
   end  % DGm
    end  % method (Static)
end  % classdef

Simple-cubic pressure element class for pressure-free velocity method (ELG3412f.m)

Consistent pressure for pressure-free velocity method (GetPresW.m)

Gauss quadrature rules on rectangle (GQuad2.m)

function [ xa,ya,wt,nq ] = GQuad2( ord )
% Returns arrays of integration points and weights for specified degree
% Usage:
%   [xa,ya,wt,nq] = GQuad2( ord )
% where
%   ord - highest degree of term integrated exactly on interval [-1,1]
  if  ord<1,     error('integration order < 1');
  elseif ord>15, error('max quad data order = 15');
  nt = fix(ord/2)+1;
    switch nt;       
       case {1, 2}  % degree 3
   Aq=[-0.5773502691896257, 0.5773502691896257];
   Hq=[ 1.0,                1.0];  
        case 3    % degree 5
   Aq=[-0.7745966692414833, 0.000000000000000, 0.7745966692414833];
   Hq=[ 0.5555555555555558, 0.8888888888888888,0.5555555555555558];
        case 4    % degree 7
   Aq=[-0.8611363115940526,-0.3399810435848563, 0.3399810435848563, ...
   Hq=[ 0.3478548451374538, 0.6521451548625461, 0.6521451548625461, ...
        case 5    % degree 9
   Aq=[-0.9061798459386640,-0.5384693101056831, 0.0000000000000000, ...
        0.5384693101056831, 0.9061798459386640];
   Hq=[ .23692688505618909, .47862867049936647, .56888888888888889, ...
        .47862867049936647, .23692688505618909];
        case 6    % degree 11
   Aq=[-.932469514203152,-.661209386466265,-.238619186083197, ...
        .238619186083197, .661209386466265, .932469514203152];
   Hq=[ .171324492379170, .360761573048139, .467913934572691, ...
        .467913934572691, .360761573048139, .171324492379170];
        case 7   % degree 13
   Aq=[-0.9491079123427585,-0.7415311855993945,-0.4058451513773972, ...
        0.0000000000000000, 0.4058451513773972, 0.7415311855993945, ...
   Hq=[ 0.1294849661688699, 0.2797053914892767, 0.3818300505051190, ...
        0.4179591836734694, 0.3818300505051190, 0.2797053914892767, ...
        case 8   % degree 15
   Aq=[-0.9602898564975363,-0.7966664774136267,-0.5255324099163290, ...
       -0.1834346424956498, 0.1834346424956498, 0.5255324099163290, ...
        0.7966664774136267, 0.9602898564975363];
   Hq=[ 0.1012285362903761, 0.2223810344533745, 0.3137066458778873, ...
        0.3626837833783620, 0.3626837833783620, 0.3137066458778873, ...
        0.2223810344533745, 0.1012285362903761];

    xa = repmat(Aq,nt,1);
    ya = xa';
    wt = repmat(Hq,nt,1);
    wt = wt.*wt';
    nq = nt*nt;


Grade node spacing (regrade.m)


[1] Holdeman, J. T. (2010), "A Hermite finite element method for incompressible fluid flow", Int. J. Numer. Meth. Fluids, 64: 376-408.

[2] Holdeman, J. T. and Kim, J.W. (2010), "Computation of incompressible thermal flows using Hermite finite elements", Comput. Methods Appl. Mech. Engr., 199: 3297-3304.

[3] Melosh, R. J. (1963), "Basis of derivation of matricies for the direct stifness method", J.A.I.A.A., 1: 1631-1637.

[4] Zienkiewicz, O. C. (1971), The Finite Element Method in Engineering Science, McGraw-Hill, London.

My wiki