PFV4 quartic velocity class

From CFD-Wiki

classdef ELS4424r < handle
    % ELS4424r
    %   Container class for 4-node modified quartic Hermite finite 
    %   elements on rectangle/quadrilateral (designated by 'S4424').
    %   The element has 24 degrees-of-freedom and is C1 continuous. 
    %   The vector element is the curl of the scalar element and is 
    %   C0 continuous. It is used for divergence-free vector fields 
    %   such as incompressible velocity and magnetic fields. 
    %
    %  Jonas Holdeman     January 2013
    
    properties (Constant)
        name = 'modified C1 quartic Hermite';
        designation = 'S4424r';
        shape = 'quadrilateral';
        nsides = 4;
        order = 4;       % order of completeness
        nnodes = 4;      % number of nodes
        nndofs = 6;      % max number of nodal dofs
        tndofe = 24;     % total number of dofs for element
        mxpowr = 5;;     % highest power/degree in s
        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 quartic-complete Hermite scalar stream function element
    % on the reference square [-1,1]x[-1,1]. The designated (4424), the
    % scalar function (4424) is C1 continuous with C0 curl and with 
    % 24 degrees-of-freedom, 6 dofs per vertex 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/64*(1+q0)^2*(1+r0)^2*(3*q0*(1-q0)^2*(2-r0)...
          +3*r0*(1-r0)^2*(2-q0)+4*(2-q0)*(2-r0)), ...
     -1/64/ri*(1+q0)^2*(1+r0)^3*(1-r0)*(2-q0)*(5-3*r0), ...
      1/64/qi*(1+q0)^3*(1-q0)*(1+r0)^2*(2-r0)*(5-3*q0),  ...
      1/64/qi^2*(1+q0)^3*(1-q0)^2*(1+r0)^2*(2-r0), ...
      1/16/qi/ri*(1+q0)^2*(1-q0)*(1+r0)^2*(1-r0), ...
      1/64/ri^2*(1+r0)^3*(1-r0)^2*(1+q0)^2*(2-q0)];
     end % s
     
    % Four-node cubic-complete Hermite solenoidal vector function element
    % on the reference square. The vector function is the curl of the 
    % quartic-complete element (4424) with 6 dofs per vertex 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=[3/64*ri*(1+q0)^2*(1-r^2)*((1+q0)*(8-9*q0+3*q^2)+(2-q0)*(1-5*r^2)), ...
    -1/64*(1+q0)^2*(1+r0)^2*(2-q0)*(7-26*r0+15*r^2), ...
     3/64*ri/qi*(1+q0)^3*(1-r^2)*(1-q0)*(5-3*q0), ...
     3/64*ri/qi^2*(1+q0)^3*(1-r^2)*(1-q0)^2, ...
     1/16/qi*(1+q0)^2*(1+r0)*(1-q0)*(1-3*r0), ... 
     1/64/ri*(1+q0)^2*(1+r0)^2*(1-r0)*(1-5*r0)*(2-q0); ... 
    -3/64*qi*(1+r0)^2*(1-q^2)*((1+r0)*(8-9*r0+3*r^2)+(2-r0)*(1-5*q^2)), ...
     3/64*qi/ri*(1+r0)^3*(1-q^2)*(1-r0)*(5-3*r0), ... 
    -1/64*(1+r0)^2*(1+q0)^2*(2-r0)*(7-26*q0+15*q^2), ... 
    -1/64/qi*(1+r0)^2*(1+q0)^2*(1-q0)*(1-5*q0)*(2-r0), ... 
    -1/16/ri*(1+r0)^2*(1+q0)*(1-r0)*(1-3*q0), ... 
    -3/64*qi/ri^2*(1+r0)^3*(1-q^2)*(1-r0)^2];
     end % S
        
   % First derivatives wrt q & r of four-node cubic-complete Hermite 
   % solenoidal vector function element on the reference square. 
   % The vector function is the curl of the quartic-complete element 
   % (4424) with 6 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=[9/64*qi*ri*(1-q^2)*(1-r^2)*(6-5*q^2-5*r^2), ...
      -3/64*qi*(1-q^2)*(1+r0)^2*(1-3*r0)*(7-5*r0), ...
       3/64*ri*(1+q0)^2*(1-r^2)*(1-3*q0)*(7-5*q0), ...
       3/64*ri/qi*(1+q0)^2*(1-r^2)*(1-q0)*(1-5*q0), ...
       1/16*(1+q0)*(1+r0)*(1-3*q0)*(1-3*r0), ...
       3/64*qi/ri*(1-q^2)*(1+r0)^2*(1-r0)*(1-5*r0); ...
       3/32*qi^2*q0*(1+r0)^2*(r0*(10*q^2+3*r^2-7)+20-20*q^2-6*r^2), ...
      -3/32*qi^2/ri*q0*(1+r0)^3*(1-r0)*(5-3*r0),...
       3/16*qi*(1-q^2)*(1+r0)^2*(2-r0)*(1+5*q0), ...
       1/16*(1+q0)*(1+r0)^2*(2-r0)*(1+2*q0-5*q^2), ...
       1/8*qi/ri*(1+r0)^2*(1-r0)*(1+3*q0), ...
       3/32*qi^2/ri^2*q0*(1+r0)^3*(1-r0)^2];

   SR=[-3/32*ri^2*r0*(1+q0)^2*(q0*(10*r^2+3*q^2-7)+20-20*r^2-6*q^2), ...
        3/16*ri*(1+q0)^2*(1-r^2)*(2-q0)*(1+5*r0),...
       -3/32*ri^2*r0/qi*(1+q0)^2*(1-q^2)*(5-3*q0), ...
       -3/32*ri^2/qi^2*r0*(1+q0)*(1-q^2)^2, ...
       -1/8*ri/qi*(1+q0)*(1-q^2)*(1+3*r0), ...
       -1/16*(1+q0)^2*(1+r0)*(2-q0)*(1+2*r0-5*r^2); ...
       -9/64*qi*ri*(1-q^2)*(1-r^2)*(6-5*q^2-5*r^2), ...
        3/64*qi*(1+r0)^2*(1-q^2)*(1-3*r0)*(7-5*r0), ...
       -3/64*ri*(1+q0)^2*(1-r^2)*(1-3*q0)*(7-5*q0), ...
       -3/64*ri/qi*(1+q0)*(1-q^2)*(1-r^2)*(1-5*q0), ...
       -1/16*(1+q0)*(1+r0)*(1-3*q0)*(1-3*r0), ...
       -3/64*qi/ri*(1-q^2)*(1-r^2)*(1+r0)*(1-5*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=[-9/32*qi^2*ri*q0*(1-r^2)*(11-10*q^2-5*r^2), ...
           3/32*qi^2*q0*(1+r0)^2*(1-3*r0)*(7-5*r0) , ...
          -9/16*qi*ri*(1-r^2)*(1-q^2)*(1+5*q0), ...
          -3/16*ri*(1+q0)*(1-r^2)*(1+2*q0-5*q^2), ...
          -1/8*qi*(1+r0)*(1+3*q0)*(1-3*r0), ...
          -3/32*qi^2/ri*q0*(1+r0)^2*(1-r0)*(1-5*r0); ...
           3/32*qi^3*(1+r0)^2*(6+(r0-2)*(30*q^2+3*r^2-7)), ...
          -3/32*qi^3/ri*(1+r0)^3*(1-r0)*(5-3*r0), ...
           3/16*qi^2*(1+r0)^2*(2-r0)*(5-2*q0-15*q^2), ...
           3/16*qi*(1+r0)^2*(2-r0)*(1-2*q0-5*q^2), ...
           3/8*qi^2/ri*(1+r0)^2*(1-r0), ...
           3/32*qi^3/ri^2*(1+r0)^3*(1-r0)^2];
       Sqr=[9/32*qi*ri^2*r0*(1-q^2)*(5*q^2-11+10*r^2), ...
            9/16*qi*ri*(1-q^2)*(1-r^2)*(5*r0+1), ...
           -3/32*ri^2*(q0+1)^2*r0*(-1+3*q0)*(-7+5*q0), ...
            3/32*ri^2/qi*(1-q^2)*(1+q0)*r0*(-1+5*q0), ...
            1/8*ri*(q0+1)*(-1+3*q0)*(1+3*r0), ...
            3/16*qi*(1-q^2)*(1+r0)*(5*r^2-2*r0-1); ...
           -9/32*qi^2*ri*q0*(1-r^2)*(5*r^2-11+10*q^2), ...
           -3/32*qi^2*q0*(3*r0-1)*(5*r0-7)*(1+r0)^2, ...
            9/16*qi*ri*(1-q^2)*(1-r^2)*(1+5*q0), ...
           -3/16*ri*(q0+1)*(1-r^2)*(5*q^2-2*q0-1), ...
           -1/8*qi*(1+r0)*(3*r0-1)*(1+3*q0), ...
           -3/32*qi^2/ri*q0*(5*r0-1)*(1-r^2)*(1+r0)];
       Srr=[-3/32*ri^3*(1+q0)^2*(6+(q0-2)*(30*r^2+3*q^2-7)), ...
             3/16*ri^2*(1+q0)^2*(2-q0)*(5-2*r0-15*r^2), ...
            -3/32*ri^3/qi*(1+q0)^2*(1-q^2)*(5-3*q0), ...
            -3/32*ri^3/qi^2*(1+q0)*(1-q^2)^2, ...
            -3/8*ri^2/qi*(1+q0)*(1-q^2), ...
            -3/16*ri*(1+q0)^2*(2-q0)*(1-2*r0-5*r^2); ...
             9/32*qi*ri^2*r0*(1-q^2)*(11-5*q^2-10*r^2), ...
            -9/16*qi*ri*(1-q^2)*(1-r^2)*(1+5*r0), ...
             3/32*ri^2*r0*(1+q0)^2*(1-3*q0)*(7-5*q0), ...
             3/32*ri^2/qi*(1+q0)*(1-q^2)*r0*(1-5*q0), ...
             1/8*ri*(1+q0)*(1-3*q0)*(1+3*r0), ...
             3/16*qi*(1-q^2)*(1+r0)*(1+2*r0-5*r^2)];
    end  % D2S  
    
    % Transpose of the Jacobian matrix at (q,r)
   function Jt=JacT(Xe,q,r)
      Jt=Xe*Gm(ELS4424r.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*ELS4424r.Gm(ELS4424r.nn(:,:),ELS4424r.nn(m,1),ELS4424r.nn(m,2)); 
     JtiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)];   % det(J)*inv(J)
     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];
     ti=blkdiag(1,JtiD,T2); 
     Bxy=Xe*ELS4424r.DGm(ELS4424r.nn(:,:),ELS4424r.nn(m,1), ...
           ELS4424r.nn(m,2));  % 2nd cross-derivative
     ti(5,2:3)=Bxy([2,1])';
   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

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