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
