PFV V8xyzcW 3D linear FE derivatives

Function V8xyzcW.m for first derivatives of 3D linear divergence-free finite element. 

<pre>
function [Vx,Vy,Vz]=V8xyzcW(ni,q,r,s)
%V8XYZc -  3D, 8-node, 48 DOF linear solenoidal element on cube 
%Arrays of first-derivatives of solenoidal (velocity) vector basis functions.
% Usage
%   [Vx,Vy,Vz]=V8xyzcW(ni,q,r,s)
%   ni - coordinates of i-th node on reference cube 
%   q,r,s - point in reference cube for evaluation 
%   Vx,Vy,Vz - first derivatives of velocity vector 
%
% 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 vector first derivatives wrt x
Vx=[-3/64*ri*si*(r0^2-s0^2)*(3*q0^2+4*q0-1), ...
    1/128*qi*si*(8*r0+75*q0^2+63*r0^2+69*s0^2-24*r0*q0^2-63*s0^2*r0^2-61+81*s0^2*r0^2*q0^2-99*s0^2*q0^2-81*r0^2*q0^2), ...
    1/128*qi*ri*(-8*s0-75*q0^2-69*r0^2-63*s0^2+63*s0^2*r0^2+61-81*s0^2*r0^2*q0^2+24*s0*q0^2+81*s0^2*q0^2+99*r0^2*q0^2), ...
    3/64*qi*(q0^2-1)*(1+s0)*(1+r0)*(9*s0*r0-9*r0-9*s0+5), -3/128*ri*(r0^2-1)*(-15*q0^2+9*s0^2*q0^2-4*s0*q0-8*q0+5-3*s0^2), ...
   -3/128*si*(s0^2-1)*(-15*q0^2+9*r0^2*q0^2-4*r0*q0-8*q0+5-3*r0^2);...
   -1/64*qi*si*(1+r0)*(27*r0*q0-27*r0*s0^2*q0+4*r0-27*r0^2*q0+27*s0^2*r0^2*q0-4*r0^2+36*q0-36*s0^2*q0), ...
    3/32*qi^2*si/ri*(r0^2-1)*(r0+2)*q0, 3/64*qi^2*(1+r0)*q0*(11*r0-9*r0*s0^2-11*r0^2+9*s0^2*r0^2+12-12*s0^2), ...
   -3/64*qi^2/ri*(r0^2-1)*q0*(-5*r0+3*r0*s0^2-2*s0-4), 1/32*qi*(1+r0)^2*(1+s0)*(r0-2)*(9*s0*q0-9*q0-2), ...
   -3/64*qi*si/ri*(r0^2-1)*(s0^2-1)*(3*r0*q0-1); ...
    1/64*qi*ri*(1+s0)*(27*s0*q0-27*s0*r0^2*q0+4*s0-27*s0^2*q0+27*s0^2*r0^2*q0-4*s0^2+36*q0-36*r0^2*q0), ...
   -3/64*qi^2*(1+s0)*q0*(11*s0-9*s0*r0^2-11*s0^2+9*s0^2*r0^2+12-12*r0^2), -3/32*qi^2*ri/si*(s0^2-1)*(s0+2)*q0, ...
   -3/64*qi^2/si*(s0^2-1)*q0*(-5*s0+3*s0*r0^2-2*r0-4), -3/64*qi*ri/si*(s0^2-1)*(r0^2-1)*(3*s0*q0-1), ...
   1/32*qi*(1+s0)^2*(1+r0)*(s0-2)*(9*r0*q0-9*q0-2)];
%
% array of solenoidal vector first derivatives wrt y 
Vy=[-3/32*ri^2*si/qi*(q0^2-1)*(q0+2)*r0, 1/64*ri*si*(1+q0)*(27*r0*q0*(1-s0^2)*(1-q0)+4*q0*(1-q0)+36*r0*(1-s0^2)), ...
    -3/64*ri^2*(1+q0)*r0*(11*q0*(1-q0)-9*s0^2*q0*(1-q0)+12*(1-s0^2)), 1/32*ri*(-2+q0)*(1+q0)^2*(1+s0)*(9*s0*r0-9*r0-2), ...
    -3/64*ri^2/qi*(q0^2-1)*r0*(-5*q0+3*s0^2*q0-2*s0-4), -3/64*ri*si/qi*(q0^2-1)*(s0^2-1)*(3*r0*q0-1); ...
     1/128*ri*si*(-8*q0-75*r0^2-69*s0^2+24*r0^2*q0+99*s0^2*r0^2+61+81*r0^2*q0^2*(1-s0^2)-63*q0^2*(1-s0^2)), ...
     3/64*qi*si*(q0^2-s0^2)*(3*r0^2+4*r0-1), ...
     1/128*qi*ri*(8*s0+69*q0^2+75*r0^2-24*s0*r0^2-81*s0^2*r0^2*(1-q0^2)-61+63*s0^2*(1-q0^2)-99*r0^2*q0^2), ...
    -1/128*qi*(q0^2-1)*(-45*r0^2+27*s0^2*r0^2-12*s0*r0-24*r0+15-9*s0^2), 3/64*ri*(r0^2-1)*(1+s0)*(1+q0)*(9*s0*q0-9*q0-9*s0+5), ...
    -3/128*si*(s0-1)*(1+s0)*(-15*r0^2+9*r0^2*q0^2-4*r0*q0-8*r0+5-3*q0^2); ...
     3/64*ri^2*(1+s0)*r0*(11*s0*(1-s0)-9*s0*q0^2*(1-s0)+12-12*q0^2), ...
    -1/64*qi*ri*(1+s0)*(27*s0*r0*(1-q0^2)*(1-s0)+4*s0-4*s0^2+36*r0*(1-q0^2)), 3/32*qi*ri^2/si*(s0^2-1)*(s0+2)*r0, ...
    -3/64*ri*qi/si*(s0^2-1)*(q0^2-1)*(3*s0*r0-1), -3/64*ri^2/si*(s0^2-1)*r0*(-5*s0+3*s0*q0^2-2*q0-4), ...
    1/32*ri*(1+s0)^2*(1+q0)*(s0-2)*(9*r0*q0-9*r0-2)];
%
% array of solenoidal vector first derivatives wrt z
Vz=[ 3/32*ri*si^2/qi*(q0^2-1)*(q0+2)*s0, 3/64*si^2*(1+q0)*s0*(11*q0*(1-q0)-9*r0^2*q0*(1-q0)+12*(1-r0^2)), ...
    -1/64*ri*si*(1+q0)*(27*s0*q0*(1-r0^2)-27*s0*q0^2*(1-r0^2)+4*q0*(1-q0)+36*s0*(1-r0^2)), ...
     1/32*si*(-2+q0)*(1+q0)^2*(1+r0)*(9*s0*r0-9*s0-2), -3/64*ri*si/qi*(q0^2-1)*(r0^2-1)*(3*s0*q0-1), ...
    -3/64*si^2/qi*(q0^2-1)*s0*(-5*q0+3*r0^2*q0-2*r0-4); ...
    -3/64*si^2*(1+r0)*s0*(11*r0-9*r0*q0^2-11*r0^2+9*r0^2*q0^2+12-12*q0^2), -3/32*qi*si^2/ri*(r0^2-1)*(r0+2)*s0, ...
       1/64*qi*si*(1+r0)*(27*s0*r0*(1-q0^2)+4*r0-27*s0*r0^2+27*s0*r0^2*q0^2-4*r0^2+36*s0-36*s0*q0^2), ...
    -3/64*qi*si/ri*(r0^2-1)*(q0-1)*(1+q0)*(3*s0*r0-1), 1/32*si*(r0-2)*(1+r0)^2*(1+q0)*(9*s0*q0-2-9*s0), ...
    -3/64*si^2/ri*(r0^2-1)*s0*(-5*r0+3*r0*q0^2-2*q0-4); ...
     1/128*ri*si*(8*q0+69*r0^2+75*s0^2-24*s0^2*q0-99*s0^2*r0^2-61-81*s0^2*q0^2*(1-r0^2)+63*q0^2*(1-r0^2)), ...
     1/128*qi*si*(-8*r0-69*q0^2-75*s0^2+61+24*r0*s0^2+81*s0^2*r0^2*(1-q0^2)+99*s0^2*q0^2-63*r0^2*(1-q0^2)), ...
    -1/64*qi*ri*(q0^2-r0^2)*(9*s0^2+12*s0-3), -1/128*qi*(q0^2-1)*(-45*s0^2+27*s0^2*r0^2-12*s0*r0-24*s0+15-9*r0^2), ...
    -1/128*ri*(3*r0-3)*(1+r0)*(-15*s0^2+9*s0^2*q0^2-4*s0*q0-8*s0+5-3*q0^2), ...
     3/64*si*(s0^2-1)*(1+r0)*(1+q0)*(-9*q0+9*r0*q0-9*r0+5)];
return;
</pre>