PFV4 cubic pressure class

From CFD-Wiki

classdef ELG3412r < handle
    % ELG3412r
    %   Container class for 4-node simple cubic Hermite finite elements
    %   on rectangle/quadrilateral (designated by 'G3412r').
    %   The scalar element is used for scalar fields such as pressure
    %   or temperature. 
    %   The vector element is used for irrotational vector fields
    %   such as pressure or thermal gradients or irrotational fluid flow. 
    %
    %  Jonas Holdeman     January 2013
    
    properties  (Constant)
        name = 'Simple-cubic Hermite';
        designation = 'G3412r';
        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 g
        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
    
 
    methods  (Static)
    
    % Four-node cubic-complete Hermite scalar potential function 
    % element on the reference square. The vector function is the gradient 
    % of this simple-cubic element (3412) with 3 dofs per corner node.
    % Parameter ni is one of coordinate pairs (-1,-1),(1,-1),(1,1),(-1,1) 
    function g=g(ni,q,r)
       qi=ni(1); q0=q*ni(1); 
       ri=ni(2); r0=r*ni(2); 
       g=[1/8*(1+q0)*(1+r0)*(2+q0*(1-q0)+r0*(1-r0)), ...
         -1/8*qi*(1+q0)*(1+r0)*(1-q^2), -1/8*ri*(1+q0)*(1+r0)*(1-r^2)];
     end % g

    % Four-node quadratic-complete Hermite irrotational vector function 
    % element on the reference square. The vector function is the gradient 
    % 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 gv=G(ni,q,r)
        qi=ni(1); q0=q*ni(1); 
        ri=ni(2); r0=r*ni(2); 
        gv=[1/8*qi*(1+r0)*(r0*(1-r0)+3*(1-q^2)), ...
           -1/8*(1+r0)*(1+q0)*(1-3*q0), ...
           -1/8*qi/ri*(1-r^2)*(1+r0); ...
            1/8*ri*(1+q0)*(q0*(1-q0)+3*(1-r^2)), ...
           -1/8/qi*ri*(1-q^2)*(1+q0), ...
           -1/8*(1+q0)*(1+r0)*(1-3*r0)];
     end % G
        
   % First derivatives wrt q & r of four-node quadratic-complete Hermite 
   % gradient vector function element on the reference square. 
   % The vector function is the gradient of the cubic-complete element 
   % (3412) with 3 dofs at each corner node.
   % Gq = array of q-derivatives of irrotational vectors 
   % Gr = array of r-derivatives of irrotational vectors 
   % Parameter ni is one of coordinate pairs (-1,-1),(1,-1),(1,1),(-1,1)
    function [Gq,Gr]=DG(ni,q,r)
       qi=ni(1); q0=q*ni(1); 
       ri=ni(2); r0=r*ni(2); 
       Gq=[-3/4*qi^2*q0*(1+r0), 1/4*qi*(1+r0)*(1+3*q0), 0; ...
           1/8*qi*ri*(4-3*q0^2-3*r0^2), -1/8*ri*(1+q0)*(1-3*q0), ...
          -1/8*qi*(1+r0)*(1-3*r0)];
       Gr=[1/8*qi*ri*(4-3*q0^2-3*r0^2), -1/8*ri*(1+q0)*(1-3*q0), ...
          -1/8*qi*(1+r0)*(1-3*r0); ...
          -3/4*ri^2*r0*(1+q0), 0, 1/4*ri*(1+q0)*(1+3*r0)];
    end  % DG
    
   % Transpose of the Jacobian matrix at (q,r)
   function Jt=JacT(Xe,q,r)
      Jt=Xe*Gm(ELG3412rr.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*ELG3412r.Gm(ELS32r.nn(:,:),ELG3412r.nn(m,1),ELG3412r.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