quicwind/saddlepoint/saddle_test_3d.m

   1 % number of cells and spacing
   2 n = [5,5,5];
   3 h = [1,1,1];
   4 fprintf('linear array of %ix%ix%i cells\n',n(1),n(2),n(3))
   5 % dimensions
   6 d = size(n,2);
   7 % number of dof per cell
   8 factor = 2*d;
   9
  10 % creating grid
  11 X = regular_mesh(n,h,1);
  12 X = add_terrain_to_mesh(X,'hill','squash',0.1);
  13 CX = center_mesh(X);
  14 xx = CX{1}; yy = CX{2}; zz = CX{3};
  15
  16 % initial wind
  17 v0 = zeros(factor*prod(n),1);
  18 for e=1:prod(n)
  19     s=(e-1)*factor+1:e*factor; % span of local dofs
  20     % initial wind
  21     v0(s)=[1,1,0,0,0,0]';
  22 end
  23
  24 % create continuity and ground flux conditions
  25 C = sparse_C(n);
  26
  27 % static cell matrices
  28 E = [.5,.5,0,0,0,0;
  29     0,0,.5,.5,0,0;
  30     0,0,0,0,.5,.5]; % matrix of interpolated winds at the center of the cells
  31 D = ones(1,factor); % diverge operator
  32
  33 % moduli for penalization
  34 moduli = [1e3,1e3,1];
  35 % initialize resulting wind
  36 Pf = sparse(factor*prod(n),factor*prod(n));
  37 v = zeros(factor*prod(n),1);
  38 % solve mass-consistent using 3D indexing
  39 for k=1:n(3)
  40     for j=1:n(2)
  41         for i=1:n(1)
  42             ind = sub2ind(n,i,j,k);
  43             s=(ind-1)*factor+1:ind*factor;
  44             [A,B] = cell_A_B(h(1),h(2),X{3},i,j,k,moduli);
  45             Pf(s,s) = cell_P(A,B,D);
  46         end
  47     end
  48 end