femwind/ndt_mult.m

   1 function y=ndt_mult(K,x,doprint)
   2 %y=nd_mult(K,x,doprint)
   3 % multiply vector x by matrix K from nd_assembly
   4 if ~exist('doprint','var')
   5     doprint=0;
   6 end
   7 [n1,n2,n3,m]=size(K);
   8 t = ndt_storage_table(m);
   9 u=reshape(x,n1,n2,n3);  % make x into grid vector if needed
  10 y=zeros(n1,n2,n3);
  11 % loop in WRF ordering (j,k,i)
  12 for i2=1:n2
  13     for i3=1:n3
  14         for i1=1:n1
  15             % global index of row = node i
  16             for j3=-1:1   % relative index j of neighbor = nonzero in the row i
  17                 for j2=-1:1
  18                     for j1=-1:1
  19                         % global index triple of neighbor node i+j out of bounds?
  20                         if ~(i1+j1<1 || i1+j1>n1 || i2+j2<1 || i2+j2>n2 || i3+j3<1 || i3+j3>n3 )
  21                             % offset m of the m where the entry (i,i+j) is stored, 0 if this row
  22                             % in fortran we won't have the 2+ because
  23                             % the array t will be indexed -1:1
  24                             m3 = t(3,2+j1,2+j2,2+j3);
  25                             m2 = t(2,2+j1,2+j2,2+j3);
  26                             m1 = t(1,2+j1,2+j2,2+j3);
  27                             % index of the matrix entry (i,i+j) in K(i+m,:)
  28                             jx = t(4,2+j1,2+j2,2+j3);
  29                             % contribution of K(i,i+j)*x(i+j) if index not out of bounds
  30                             y(i1,i2,i3)=y(i1,i2,i3)+K(i1+m1,i2+m2,i3+m3,jx)*u(i1+j1,i2+j2,i3+j3);
  31                             switch doprint
  32                                 case 1   % matlab ordering
  33                                     fprintf('i=(%g %g %g) j=(%g %g %g) m=(%g %g %g) jx=%g i+m=(%g %g %g) i+j=(%g %g %g)\n',...
  34                                            [ i1,i2,i3,    j1,j2,j3,    m1,m2,m3,  jx, i1+m1,i2+m2,i3+m3, i1+j1,i2+j2,i3+j3 ])
  35                                 case 2   % WRF storage ordering
  36                                     fprintf('ikj: row=(%g %g %g) rel.column=(%g %g %g) rel.row.stored=(%g %g %g) at %g\n',...
  37                                                     [ i1,i3,i2,             j1,j3,j2,             m1,m3,m2,    jx])
  38                                 case 3   % print code
  39                                     if i1 == 2 && i2 == 2 && i3 == 2
  40                                         fprintf('        kmat(i%s,k%s,j%s,%2i)*u(i%s,k%s,j%s) +  &\n',...
  41                                             pm(m1),pm(m3),pm(m2),jx,pm(j1),pm(j3),pm(j2))
  42                                     end
  43                             end
  44                         end
  45                     end
  46                  end
  47             end
  48         end
  49     end
  50 end
  51 end
  52
  53 function s=pm(j)
  54 switch j
  55     case 1
  56         s='+1';
  57     case -1
  58         s='-1';
  59     case 0
  60         s='  ';
  61     otherwise
  62         error('pm: argument must be in -1:1')
  63 end
  64 end