quicwind/Amul_v.m

   1 function res = Amul_v(X,A,v)
   2 % multiplies the penalty weight diagonal matrix and a vector
   3 % the operation is res = A * v, where A is a diagonal matrix
   4 % set:
   5 %    X   mesh
   6 %    A   matrix to multiply, must be diagonal
   7 %    v   vector to multiply
   8
   9 % Error checking
  10 check_mesh(X);
  11
  12 % Initialize result
  13 res = zeros(size(v));
  14
  15 % Mesh sizes and vars
  16 nx = size(X{1},1)-1;
  17 ny = size(X{1},2)-1;
  18 nz = size(X{1},3)-1;
  19 x = X{1}; y = X{2}; z = X{3};
  20
  21 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% A * v %%%%%%%%%%%%%%%%%%%%%%%%%%%%
  22 for k=1:nz
  23     for j=1:ny
  24         for i=1:(nx+1)
  25             res(i + (j-1)*(nx+1) + (k-1)*(nx+1)*ny) = ...
  26                 A(i + (j-1)*(nx+1) + (k-1)*(nx+1)*ny,i + (j-1)*(nx+1) + (k-1)*(nx+1)*ny) * v(i + (j-1)*(nx+1) + (k-1)*(nx+1)*ny);
  27         end
  28     end
  29 end
  30 add_factor = (nx+1)*ny*nz;
  31 for k=1:nz
  32     for j=1:(ny+1)
  33         for i=1:nx
  34             res(add_factor + i + (j-1)*nx + (k-1)*nx*(ny+1)) = ...
  35                 A(add_factor + i + (j-1)*nx + (k-1)*nx*(ny+1),add_factor + i + (j-1)*nx + (k-1)*nx*(ny+1)) * v(add_factor + i + (j-1)*nx + (k-1)*nx*(ny+1));
  36         end
  37     end
  38 end
  39 add_factor = add_factor + nx*(ny+1)*nz;
  40 for k=1:(nz+1)
  41     for j=1:ny
  42         for i=1:nx
  43             res(add_factor + i + (j-1)*nx + (k-1)*nx*ny) = ...
  44                 A(add_factor + i + (j-1)*nx + (k-1)*nx*ny,add_factor + i + (j-1)*nx + (k-1)*nx*ny) * v(add_factor + i + (j-1)*nx + (k-1)*nx*ny);
  45         end
  46     end
  47 end