quicwind/mat_div3.m

   1 function Dmat=mat_div3(X,h)
   2 % Dmat=mat_div3(X)
   3 % Compute divergence matrix for a staggered mesh
   4 %    Note: this will work only for meshes on a standard Cartesian
   5 %    coordinate system
   6 % arguments:
   7 %    X      staggered mesh (midpoints of sides)
   8 %    h(1:3) mesh step
   9 % output:
  10 %    Dmat   divergence matrix (sparse) on mesh
  11 %
  12
  13 % Input parameter checking
  14 if (size(X) ~= 3)
  15     error('mat_div3 only works in 3 dimensions')
  16 end
  17 if ~(isequal(size(X{1}),size(X{2}),size(X{3})))
  18     error('Size of mesh must be the same in for u,w,z for mat_div3')
  19 end
  20 if ~exist('h','var')
  21     h=[1,1,1];
  22 end
  23
  24 % Mesh computations
  25 x = X{1}; y = X{2}; z = X{3};
  26 [nx1, ny1, nz1] = size(x);
  27 nx = nx1 - 1; ny = ny1 - 1; nz = nz1 - 1;
  28 n = nx * ny * nz1 + nx * ny1 * nz + nx1 * ny * nz;
  29
  30 % Compute number of nonzeros
  31 %   Each boundary cell accounts for 1 nonzero, each non-boundary cell adds
  32 %   2 (in each dimension)
  33 nnz = (nx1 - 2) * ny * nz + ny * nz + ...
  34     (ny1 - 2) * nx * nz + nx * nz + ...
  35     (nz1 - 2) * nx * ny + nx * ny;
  36 nnz = 2 * nnz;
  37
  38 % For sparse matrix
  39 I = zeros(nnz,1);
  40 J = zeros(nnz,1);
  41 D = zeros(nnz,1); % Value of Dmat at index (I,J)
  42
  43 D_col = 1;
  44 D_entry = 1; % Current entry of Dmat
  45 for ux=1:numel(X)
  46
  47     if (ux == 1)
  48         surf_area = h(2) * h(3);
  49     elseif (ux == 2)
  50         surf_area = h(3) * h(1);
  51     else
  52         surf_area = h(1) * h(2);
  53     end
  54
  55     % Size of the grid for the current direction
  56     grid_size = [nx,ny,nz];
  57     grid_size(ux) = grid_size(ux) + 1;
  58     X_len = grid_size(1); Y_len = grid_size(2); Z_len = grid_size(3);
  59
  60     % Main loop, iterates through all columns of Dmat
  61     for kx=1:Z_len
  62         for jx=1:Y_len
  63             for ix=1:X_len
  64
  65                 % Edge cases and relevant indices
  66                 if (ux == 1)
  67                     low_row_ind = ix - 1 + (jx - 1) * nx + (kx - 1) * ny * nx;
  68                     row_ind = ix + (jx - 1) * nx + (kx - 1) * ny * nx;
  69                     is_lower_edge = ix == 1;
  70                     is_upper_edge = ix == X_len;
  71                 elseif (ux == 2)
  72                     low_row_ind = ix + (jx - 2) * nx + (kx - 1) * ny * nx;
  73                     row_ind = ix + (jx - 1) * nx + (kx - 1) * ny * nx;
  74                     is_lower_edge = jx == 1;
  75                     is_upper_edge = jx == Y_len;
  76                 else
  77                     low_row_ind = ix + (jx - 1) * nx + (kx - 2) * ny * nx;
  78                     row_ind = ix + (jx - 1) * nx + (kx - 1) * ny * nx;
  79                     is_lower_edge = kx == 1;
  80                     is_upper_edge = kx == Z_len;
  81                 end
  82
  83                 if is_lower_edge
  84                     % Special case for top boundary
  85                     D(D_entry) = -1 * surf_area;
  86                     I(D_entry) = row_ind;
  87                     J(D_entry) = D_col;
  88                     D_entry = D_entry + 1;
  89                 elseif is_upper_edge
  90                     % Special case for lower boundary
  91                     D(D_entry) = 1 * surf_area;
  92                     I(D_entry) = low_row_ind;
  93                     J(D_entry) = D_col;
  94                     D_entry = D_entry + 1;
  95                 else
  96                     % Other elements
  97                     D(D_entry) = 1 * surf_area;
  98                     I(D_entry) = low_row_ind;
  99                     J(D_entry) = D_col;
 100                     D_entry = D_entry + 1;
 101                     D(D_entry) = -1 / h(ux);
 102                     I(D_entry) = row_ind;
 103                     J(D_entry) = D_col;
 104                     D_entry = D_entry + 1;
 105                 end
 106             D_col = D_col + 1;
 107             end
 108         end
 109     end
 110 end
 111 Dmat = sparse(I, J, D, nx * nz *ny, n);
 112