PFV cubic velocity class

From CFD-Wiki

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), ...
        -1/8*ri*(1+q0)*(1-3*q0)];
    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)
      Jt=Xe*Gm(ELS3412r.nn(:,:),q,r); 
   end  % JacT
   
   % Test to see if transformation is affine, returns True or False
   function isit=isaffine(Xe)
       isit=sum(abs(Xe(:,1)-Xe(:,2)+Xe(:,3)-Xe(:,4)))<4*eps;
   end  % isaffine
   
   % Post-multiplying matrix Ti^-1
   function ti=Ti(Xe,m)
      Jt=Xe*ELS3412r.Gm(ELS3412r.nn(:,:),ELS3412r.nn(m,1),ELS3412r.nn(m,2)); 
      JtiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)];   % det(J)*inv(J)
      ti=blkdiag(1,JtiD);
   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)
       g=.25*(1-ni(:,1).*q)*(1+ni(:,2).*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