PFV V8cW 3D linear FE

From CFD-Wiki

Function V8cW.m for 3D linear divergence-free finite element.

function S=V8cW(ni,q,r,s)
%V8CW - 3D, 8-node, 48 DOF solenoidal velocity element
%Array of solenoidal velocity basis functions, corner nodes.
% Usage
%   V = V8cW(ni,q,r,s)
%   ni - coordinates of i-th node on reference cube 
%   q,r,s - point in reference cube for evaluation 
%   V - velocity vector at (q,r,s)
%
% Jonas Holdeman, January 2011
% --------------------------------------------------------------

qi=ni(1); q0=q*ni(1); 
ri=ni(2); r0=r*ni(2);
si=ni(3); s0=s*ni(3);
% array of solenoidal vectors 
S=[ 3/64*ri*si/qi*(1-q0^2)*(q0+2)*(r0^2-s0^2), 1/128*si*(1+q0)*(q0*(1-q0)*(33*s0^2+27*r0^2*(1-s0^2)+8*r0-25)-36*(1-r0^2)*(1-s0^2)), ...
   -1/128*ri*(1+q0)*(q0*(1-q0)*(33*r0^2+27*s0^2*(1-r0^2)+8*s0-25)-36*(1-r0^2)*(1-s0^2)), ...
    1/64*(1+q0)^2*(1+r0)*(1+s0)*(2-q0)*(4-(9-9*r0)*(1-s0)), 3/128*ri/qi*(1-q0^2)*(1-r0^2)*(3*q0*(1-s0^2)+2*q0+2*s0+4), ...
    3/128*si/qi*(1-q0^2)*(1-s0^2)*(3*q0*(1-r0^2)+2*q0+2*r0+4); ...
   -1/128*si*(1+r0)*(r0*(1-r0)*(33*s0^2+27*q0^2*(1-s0^2)+8*q0-25)-36*(1-q0^2)*(1-s0^2)), ...
    3/64*qi*si/ri*(1-r0^2)*(r0+2)*(s0^2-q0^2), 1/128*qi*(1+r0)*(r0*(1-r0)*(33*q0^2+27*s0^2*(1-q0^2)+8*s0-25)-36*(1-q0^2)*(1-s0^2)), ...
    3/128*qi/ri*(1-r0^2)*(1-q0^2)*(3*r0*(1-s0^2)+2*r0+2*s0+4), 1/64*(1+r0)^2*(1+s0)*(1+q0)*(2-r0)*(4-9*(1-s0)*(1-q0)), ...
    3/128*si/ri*(1-r0^2)*(1-s0^2)*(3*r0*(1-q0^2)+2*r0+2*q0+4); ...
    1/128*ri*(1+s0)*(s0*(1-s0)*(33*r0^2+27*q0^2*(1-r0^2)+8*q0-25)-36*(1-q0^2)*(1-r0^2)), ...
   -1/128* qi*(1+s0)*(s0*(1-s0)*(33*q0^2+27*r0^2*(1-q0^2)+8*r0-25)-36*(1-q0^2)*(1-r0^2)), 3/64*qi*ri/si*(1-s0^2)*(s0+2)*(q0^2-r0^2), ...
    3/128*qi/si*(1-s0^2)*(1-q0^2)*(3*s0*(1-r0^2)+2*s0+2*r0+4), 3/128*ri/si*(1-s0^2)*(1-r0^2)*(3*s0*(1-q0^2)+2*s0+2*q0+4), ...
    1/64*(1+s0)^2*(1+q0)*(1+r0)*(2-s0)*(4-9*(1-q0)*(1-r0))]; 
return;