# PFV cubic velocity class

```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';
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)
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
```