femwind/nd_mult.m

   1 function y=nd_mult(K,x)
   2 %y=nd_mult(K,x)
   3 % multiply vector x by matrix K from nd_assembly
   4 [n1,n2,n3,m1,m2,m3]=size(K);
   5 if any([m1,m2,m3]~=3), error('K must be 3D stencil'),end
   6 u=reshape(x,n1,n2,n3);  % make x into grid vector
   7 y=zeros(n1,n2,n3);
   8 for j3=-1:1
   9     for j2=-1:1
  10         for j1=-1:1
  11             for i3=max(1,1-j3):min(n3, n3-j3)
  12                 k3 = i3 + j3;
  13                 for i2=max(1,1-j2):min(n2, n2-j2)
  14                     k2 = i2 + j2;
  15                     for i1=max(1,1-j1):min(n1, n1-j1)
  16                         k1 = i1 + j1;
  17                         % contribution of K(i,j)*x(j)
  18                         %i1 and i2  and i3
  19                         y(i1,i2,i3)=y(i1,i2,i3)+K(i1,i2,i3,2 + j1,2 + j2,2 + j3)*x(k1,k2,k3);
  20                     end
  21                  end
  22             end
  23         end
  24     end
  25 end