femwind/hexa_test.m

   1 format compact
   2 A=eye(3);
   3 % lexicographic unit cube
   4 %    x y z
   5 X0 = [0 0 0  %1
   6      1 0 0  %2
   7      0 1 0  %3
   8      1 1 0  %4
   9      0 0 1  %5
  10      1 0 1  %6
  11      0 1 1  %7
  12      1 1 1  %8
  13      ]';
  14
  15 %   7-----8
  16 %  /|    /|
  17 % 5-----6 |
  18 % | |   | |
  19 % | 3---|-4
  20 % |/    |/
  21 % 1-----2
  22
  23 % linear transformation
  24 T = rand(3);
  25 % T = magic(3);
  26 X = T*X0+rand(3,1)*ones(1,8);
  27
  28 u = rand(3,1);
  29 [K, F, Jg] = hexa(A,X,u);
  30 eig(K)  % one eigenvalue zero
  31
  32 % test F
  33 % for linear field mu with values V at nodes V*F should be -integral (grad mu) * u
  34 gradmu=[4,5,6]; c=1;
  35 V = gradmu*X + c;
  36 vol=hexa_volume(X)
  37 exact = - gradmu*u*vol;
  38 err_F = V*F - exact
  39
  40 if exist('vrrotvec2mat')
  41 % invariariant to random isometry
  42     S=vrrotvec2mat(vrrotvec(randn(3,1),randn(3,1)))*diag(sign(randn(3,1)));
  43     [K0, ~, ~] = hexa(A,X,u);
  44     [K1, ~, ~] = hexa(A,S*X,u);
  45     err_iso=norm(K0-K1,'fro')
  46 else
  47     warning('No vrrotvec2mat, cannot test rotation invariance. Install Simulink 3D Animation.')
  48 end
  49
  50 % test stretch
  51 X3 = diag([1 1 2])*X0;
  52 [K3, ~, ~] = hexa(A,X3,u);
  53
  54 % test same results from matlab and fortran
  55 if exist('fortran/hexa_test.exe')
  56     disp('testing if same result in fortran')
  57     err=hexa_fortran(A,X3,u)
  58     if abs(err)<1e-6
  59         fprintf('error %g OK\n',err)
  60     else
  61         error(sprintf('error %g too large',err))
  62     end
  63 else
  64     warning('fortran/hexa_test.exe not available')
  65 end