femwind/hexa_test.m

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