quicwind/Dmul_v.m

   1 function res = Dmul_v(X,trans,v)
   2 % multiplies the divergence matrix (or its transpose) and a vector
   3 % the operation is res = op(D) * v, where op(D) = D or D^T
   4 % set:
   5 %    X       mesh
   6 %    trans   indicates whether the matrix is transposed or not:
   7 %               trans = 'T' or trans = 't': op(D) = D^T
   8 %               trans = 'N' or trans = 'n': op(D) = D
   9 %            default is 'N'
  10 %    v       vector to multiply
  11
  12 % Error checking
  13 check_mesh(X);
  14 if ~exist('trans','var')
  15     trans = 'N';
  16 end
  17 if (trans ~= 'T' && trans ~= 't' && trans ~= 'N' && trans ~= 'n')
  18     error('trans must be N, T, n, or t')
  19 end
  20
  21 % Mesh sizes and vars
  22 nx = size(X{1},1)-1;
  23 ny = size(X{1},2)-1;
  24 nz = size(X{1},3)-1;
  25 x = X{1}; y = X{2}; z = X{3};
  26
  27 % Initialize result
  28 if (trans == 'T' || trans == 't')
  29     res = zeros((nx+1)*ny*nz+nx*(ny+1)*nz+nx*ny*(nz+1),1);
  30 else
  31     res = zeros(nx*ny*nz,1);
  32 end
  33
  34 if (trans == 't' || trans == 'T')
  35
  36 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% D^T * v %%%%%%%%%%%%%%%%%%%%%%%%%%%%
  37
  38 helper_x = 1;
  39 helper_y = 1;
  40 helper_z = 1;
  41 for k=1:nz
  42     for j=1:ny
  43         for i=1:(nx+1)
  44             if (i == 1)
  45                 res(ny*(nx+1)*(k-1)+(nx+1)*(j-1)+1) = ...
  46                     -v(nx * (helper_x - 1) + 1);
  47             elseif (i == (nx+1))
  48                 res(ny*(nx+1)*(k-1)+(nx+1)*(j-1)+nx+1) = ...
  49                     v(nx * helper_x);
  50             else
  51                 res(ny*(nx+1)*(k-1)+(nx+1)*(j-1)+i) = ...
  52                     v(nx * (helper_x-1) + i-1)-v(nx * (helper_x-1)+i);
  53             end
  54         end
  55         helper_x = helper_x + 1;
  56     end
  57 end
  58 add_factor = (nx+1)*ny*nz;
  59 for k=1:nz
  60     for j=1:(ny+1)
  61         for i=1:nx
  62             if (j == 1)
  63                 res(add_factor + (ny+1)*nx*(k-1)+nx*(j-1)+i) = ...
  64                     -v(helper_y - 1 + i);
  65             elseif (j == (ny+1))
  66                 res(add_factor + (ny+1) * nx * k - nx + i) = ...
  67                     v(helper_y);
  68                 helper_y = helper_y + 1;
  69             else
  70                 res(add_factor + (ny + 1) * nx * (k-1) + nx * (j - 1) + i) = ...
  71                     v(helper_y)-v(helper_y + nx);
  72                 helper_y = helper_y + 1;
  73             end
  74         end
  75     end
  76 end
  77 add_factor = add_factor + (ny+1)*nx*nz;
  78 for k=1:(nz+1)
  79     for j=1:ny
  80         for i=1:nx
  81             if (k == 1)
  82                 res(add_factor + (ny+1)*nx*(k-1)+nx*(j-1)+i) = ...
  83                     -v(helper_z - 1 + (j - 1) * nx + i);
  84             elseif (k == (nz+1))
  85                 res(add_factor + nx * ny * nz + (j-1)*nx + i) = ...
  86                     v(helper_z);
  87                 helper_z = helper_z + 1;
  88             else
  89                 res(add_factor + ny * nx * (k-1) + nx * (j-1) + i) = ...
  90                     v(helper_z)-v(helper_z + nx*ny);
  91                 helper_z = helper_z + 1;
  92             end
  93         end
  94     end
  95 end
  96
  97 else % trans = 'N' || trans = 'n'
  98
  99 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% D * v %%%%%%%%%%%%%%%%%%%%%%%%%%%%
 100 for k=1:nz
 101     for j=1:ny
 102         for i=1:(nx+1)
 103             if (i ~= 1)
 104                 res((i-1) + (j-1)*nx + (k-1)*nx*(ny)) = ...
 105                     res((i-1) + (j-1)*nx + (k-1)*nx*(ny)) ...
 106                     - v((i-1) + (j-1)*(nx+1) + (k-1)*(nx+1)*(ny));
 107             end
 108             if (i ~= (nx+1))
 109                 res(i + (j-1)*nx + (k-1)*nx*(ny)) = ...
 110                     res(i + (j-1)*nx + (k-1)*nx*(ny)) ...
 111                     + v(i + 1 + (j-1)*(nx+1) + (k-1)*(nx+1)*(ny));
 112             end
 113         end
 114     end
 115 end
 116 add_factor = (nx+1)*ny*nz;
 117 for k=1:nz
 118     for j=1:(ny+1)
 119         for i=1:nx
 120             if (j ~= 1)
 121                 res(i + (j-2)*nx + (k-1)*nx*ny) = ...
 122                     res(i + (j-2)*nx + (k-1)*nx*ny) ...
 123                     - v(add_factor + i + (j-2)*nx + (k-1)*nx*(ny+1));
 124             end
 125             if (j ~= (ny+1))
 126                 res(i + (j-1)*nx + (k-1)*nx*ny) = ...
 127                     res(i + (j-1)*nx + (k-1)*nx*ny) ...
 128                     + v(add_factor + i + j*nx + (k-1)*nx*(ny+1));
 129             end
 130         end
 131     end
 132 end
 133 add_factor = add_factor + nx*(ny+1)*nz;
 134 for k=1:(nz+1)
 135     for j=1:ny
 136         for i=1:nx
 137             if (k ~= 1)
 138                 res(i + (j-1)*nx + (k-2)*nx*ny) = ...
 139                     res(i + (j-1)*nx + (k-2)*nx*ny) ...
 140                     - v(add_factor + i + (j-1)*nx + (k-2)*nx*ny);
 141             end
 142             if (k ~= (nz+1))
 143                 res(i + (j-1)*nx + (k-1)*nx*ny) = ...
 144                     res(i + (j-1)*nx + (k-1)*nx*ny) ...
 145                     + v(add_factor + i + (j-1)*nx + (k)*nx*ny);
 146             end
 147         end
 148     end
 149 end
 150 end
 151
 152
 153