femwind/hexa.m

   1 function [Kloc,Floc,Jg]=hexa(A,X,u0)
   2 % create local stiffness matrix for hexa 3d element
   3 % in:
   4 %   A   coefficient matrix size 3x3, symmetric positive definite
   5 %   X   nodes coordinates size 3x8, one each column is one node
   6 %   u0   column vector of input wind at the center of the element
   7 % out:
   8 %   Kloc   local stiffness matrix
   9 %   Floc   local divergence load vector
  10 %   Jg     gradient at center of function with values V is V'*Jg
  11
  12 % basis functions on reference element [-1,1]^3
  13 Nb = 8;  % number of basis functions
  14 ib =[  % coordinates of basis functions
  15     -1    -1    -1
  16      1    -1    -1
  17     -1     1    -1
  18      1     1    -1
  19     -1    -1     1
  20      1    -1     1
  21     -1     1     1
  22      1     1     1];
  23 % the value of basis function k at x is
  24 % bf= @(k,x) (1+ib(k,1)*x(1))*(1+ib(k,2)*x(2))*(1+ib(k,3)*x(3))/8;
  25 %check_symmetry(A,'A',eps)
  26 % gaussian quadrature nodes
  27 g=0.5773502691896257;
  28 s = g*ib;
  29 Ng=Nb;  %  number of Gauss points
  30 s(Ng+1,:)=0; % extra point at center
  31
  32 Kloc = zeros(Nb);
  33 Floc = zeros(Nb,1);
  34 for j=1:Ng+1
  35     gradf = gradbfs(s(j,:));
  36     Jx = X*gradf; % Jacobian at s
  37     [q,r]=qr(Jx);
  38     Jg   = (gradf/r)*q'; % gradf*inv(Jx)
  39     adetJx = abs(prod(diag(r))); % det(Jx)
  40     if j<=Ng % contribution to stiffness
  41         K_at_s = Jg * A * Jg' * adetJx;
  42         Kloc = Kloc + K_at_s;
  43     elseif ~isempty(u0)   % contribution to divergence load
  44         % volume = determinant at midpoint times the volume of reference element
  45         % exact if the element is linear image of reference
  46         % or all sides planar, and sides in two directions are parallel
  47         vol = adetJx*8;
  48         Floc = Floc - Jg * u0 * vol;
  49     end
  50 end
  51 %check_symmetry(Kloc,'Kloc',eps)
  52 end
  53