femwind/femwind_test.m

   1 disp('femwind_test')
   2
   3 do_plot=0
   4
   5 % dimensions in elements
   6 sc=1; % mesh scale
   7 n = sc*[20,10,5];
   8 sc2=10;
   9 n(1:2)=n(1:2)*sc2
  10 h = [1,1,1]/sc;
  11 fprintf('linear array of %ix%ix%i cells\n',n(1),n(2),n(3))
  12 da=[1 1 5]
  13 string_diag_A=sprintf('%g %g %g',da);
  14 A = diag(da);
  15 lambda = zeros(prod(n+1),1); % placeholder solution
  16
  17 % creating the grid
  18 X = regular_mesh(n,h,1.5^(1/sc));
  19 X = add_terrain_to_mesh(X,'hill','squash',0.3);
  20 CX = center_mesh(X);
  21
  22 % initial wind at the centers of the elements
  23 U0={ones(n),zeros(n),zeros(n)};
  24
  25 if do_plot
  26
  27     % show mesh
  28     hold off
  29     figure
  30     plot_mesh_3d(X), hold on,
  31     title('The wind mesh, wind vector in centers, lambda in corners')
  32
  33     % show initial wind
  34     figure
  35     plot_mesh_3d(X,[1,n(1),1,n(2),1,1]), hold on,
  36     plot_wind_3d(CX,U0)
  37     hold off
  38     axis equal
  39     title('Initial wind')
  40
  41     % show initial wind
  42     figure
  43     plot_mesh_3d(X,[1,n(1),1,n(2),1,1]), hold on,
  44     plot_wind_3d(CX,U0,1)
  45     hold off
  46     axis equal
  47     title('Initial wind lowest layer')
  48 end
  49
  50 % assemble sparse system matrix
  51 [K,F,~] = sparse_assembly(A,X,U0,lambda);
  52
  53 % dirichlet boundary conditions
  54 [K,F]=apply_boundary_conditions(K,F,X);
  55
  56 % solve the equations
  57 lambda = sparse_solve(K,F,X,'s');
  58
  59 % assemble final wind
  60 [~,~,W] = sparse_assembly(A,X,U0,lambda);
  61
  62 if do_plot
  63     % plot resulting wind
  64     figure
  65     plot_mesh_3d(X,[1,n(1),1,n(2),1,1]), hold on,
  66     plot_wind_3d(CX,W)
  67     hold off
  68     axis equal
  69     title(['Final wind a=',string_diag_A])
  70
  71     figure
  72     plot_mesh_3d(X,[1,n(1),1,n(2),1,1]), hold on,
  73     plot_wind_3d(CX,W,1)
  74     hold off
  75     axis equal
  76     title(['Final wind lowest layer a=',string_diag_A])
  77
  78     figure
  79     plot_mesh_3d(X,[1,n(1),1,n(2),1,1]), hold on,
  80     plot_wind_3d(CX,W,1:2)
  81     hold off
  82     axis equal
  83     title(['Final wind lowest layers a=',string_diag_A])
  84 end
  85 % condition_number=scond(K)
  86 n