femwind/ndt_test.m

   1 format compact
   2 % test for ndt_assembly and ndt_mult
   3 nel = [5,4,3]
   4 %nel=[1 1 1]
   5 n=nel+1
   6 h = [1,1,1]
   7 expand=1.3
   8 A = diag([1 1 1])
   9 lambda=[]
  10 params=[]
  11 u0=[];
  12
  13 X = regular_mesh(nel,h,expand);
  14 % X = add_terrain_to_mesh(X, 'hill', 'shift', 0.1)
  15 X = add_terrain_to_mesh(X, 'hill', 'squash', 0.1)  % more thorough testing
  16
  17 plot_mesh_3d(X)
  18
  19 K=nd_assembly(A,X,lambda,params);
  20 [n1,n2,n3,m1,m2,m3]=size(K);
  21 K27 = reshape(K,n1,n2,n3,m1*m2*m3);  % make to n x 27 nd format
  22
  23 disp('converting to reduced storage format with 14 numbers per row')
  24 K14 = ndt_convert(K27,14);
  25 disp('multiplying st 14 by all ones')
  26 x1=ones(n1,n2,n3);
  27 y1=ndt_mult(K14,x1,3);
  28 K14_err_zero=big(y1)  % should be zero
  29 if abs(K14_err_zero)>1e-10, error('should be zero'),end
  30
  31 disp('testing multiplication by ones should give zero')
  32 x1=ones(n1,n2,n3);
  33 y27=ndt_mult(K27,x1);
  34 K27_err_zero=big(y27)  % should be zero
  35
  36
  37 disp('compare matrix vector multiply nd and ndt 27 numbers per row')
  38 xr=rand(n1,n2,n3);
  39 y27=ndt_mult(K27,xr);
  40 yn=nd_mult(K,xr);
  41 Kn_27_err=big(yn-y27)
  42
  43 y14=ndt_mult(K14,xr);
  44 K14_err_rand=big(y14-y27)  % should be zero
  45 if abs(K14_err_rand)>1e-10, error('should be zero'),end
  46
  47 disp('convert to sparse compare matrix-vector multiply')
  48 Ks = ndt_convert(K,'sparse');
  49 ys=Ks*xr(:);
  50 K_Ks_rand_err=big(ys(:)-y27(:))
  51 if abs(K_Ks_rand_err)>1e-10, error('should be zero'),end
  52
  53 if exist('fortran/ndt_boundary_conditions_test.exe')
  54     disp('testing if ndt_boundary_conditions gives same result in fortran')
  55     ndt_boundary_conditions_fortran(K14);
  56 else
  57     warning('fortran/ndt_boundary_conditions_test.exe not available')
  58 end
  59
  60 disp('testing multiply vec by st 14 vs st 27')
  61
  62 % test same results for ndt_mult from matlab and fortran
  63 if exist('fortran/ndt_mult_test.exe')
  64     disp('testing if ndt_mult same result in fortran')
  65     err=ndt_mult_fortran(K14,xr);
  66     if abs(err)<1e-6
  67         fprintf('error %g OK\n',err)
  68     else
  69         error(sprintf('error %g too large',err))
  70     end
  71 else
  72     warning('fortran/ndt_mult_test.exe not available')
  73 end
  74
  75 disp('sparse assembly vs ndt_assembly test')
  76 K_sparse=sparse_assembly(A,X,u0,lambda,params);
  77 err_mat_sparse=big(Ks-K_sparse);
  78 if abs(err_mat_sparse)>1e-10, error('should be zero'),end
  79
  80 K27a = ndt_assembly(A,X,u0,lambda,params,27);
  81 K27a_err = big(K27 - K27a)
  82 if abs(K27a_err)>1e-10, error('should be zero'),end
  83
  84
  85 K14a = ndt_assembly(A,X,u0,lambda,params,14);
  86 disp('multiplying ndt_assembly st 14 by all ones')
  87 x=ones(n1,n2,n3);
  88 y1=ndt_mult(K14,x1);
  89 K14a_err_zero=big(y1)  % should be zero
  90 if K14a_err_zero > 1e-6
  91     error('residual should be zero')
  92 end
  93
  94 disp('comparing st 14 with converged st 27')
  95 K14a_err = big(K14 - K14a)
  96 if abs(K14a_err)>1e-10, error('should be zero'),end
  97
  98 disp('ndt_assembly_fortran')
  99 err = ndt_assembly_fortran(A,X,u0,lambda,params,14)
 100 if err > 1e-6
 101     error('ndt_assembly_fortran error too large')
 102 end