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, Fold, Jg, vol] = hexa(A,X,u);
  31 F = -vol*Jg*u;
  32 err_F=big(F-Fold)
  33 eig_K=eig(K)  % one eigenvalue zero
  34
  35 % test F
  36 % for linear field mu with values V at nodes V*F should be -integral (grad mu) * u
  37 gradmu=[4,5,6]; c=1;
  38 V = gradmu*X + c;
  39 vol=hexa_volume(X)
  40 exact = - gradmu*u*vol;
  41 err_F = V*F - exact
  42
  43 if exist('vrrotvec2mat')
  44 % invariariant to random isometry
  45     S=vrrotvec2mat(vrrotvec(randn(3,1),randn(3,1)))*diag(sign(randn(3,1)));
  46     [K0, ~, ~] = hexa(A,X,u);
  47     [K1, ~, ~] = hexa(A,S*X,u);
  48     err_iso=norm(K0-K1,'fro')
  49     if err_iso > 1e-6
  50         warning(sprintf('err_iso %g too large',err_iso))
  51     end
  52
  53 else
  54     warning('No vrrotvec2mat, cannot test rotation invariance. Install Simulink 3D Animation.')
  55 end
  56
  57 % test stretch
  58 X3 = diag([1 1 2])*X0;
  59 [K3, ~, ~] = hexa(A,X3,u);
  60
  61 % test same results from matlab and fortran
  62 if exist('fortran/hexa_test.exe')
  63     disp('testing if same result in fortran')
  64     err=hexa_fortran(A,X3,u)
  65     if abs(err)<1e-6
  66         fprintf('error %g OK\n',err)
  67     else
  68         error(sprintf('error %g too large',err))
  69     end
  70 else
  71     warning('fortran/hexa_test.exe not available')
  72 end