femwind/ndt_storage_table.m

   1 function t=ndt_storage_table(m)
   2 % t=ndt_storage_table(m)
   3 % returns 6 storage indices for each of 27 neigbors as array size (6,3,3,3)
   4 % t(1:3,...) row offsets 0 = this this row
   5
   6 switch m
   7     case 27
   8         % no row offset
   9         t=zeros(4,3,3,3);
  10         % flat indexing within the row
  11         t(4,:,:,:)=reshape(1:27,3,3,3);
  12     case 14
  13         g=@(j1,j2,j3)1+(j1+1)+3*(j2+1+3*(j3+1)); % storage function
  14         g0 = g(0,0,0);
  15         n=0;
  16         for j3=-1:1
  17             for j2=-1:1
  18                 for j1=-1:1
  19                     gj=g(j1,j2,j3);
  20                     if gj >= g0
  21                         % we are in the upper triangle including diagonal
  22                         % connecting from 0 0 0 to j1 j2 j3
  23                         n=n+1;
  24                         t(1:3,2+j1,2+j2,2+j3)=0; % is centered at 2 not 0
  25                         t(4  ,2+j1,2+j2,2+j3)=n;
  26                         % disp([gj,j1,j2,j3,squeeze(t(:,2+j1,2+j2,2+j3))'])
  27                     end
  28                 end
  29             end
  30         end
  31         if n ~= 14  % just checking
  32             nu, error('wrong size of upper triangular + diagonal part')
  33         end
  34         % now fill the lower triangular entries by index
  35         % to the upper triangular entry on the other row
  36         for j3=-1:1
  37             for j2=-1:1
  38                 for j1=-1:1
  39                     gj=g(j1,j2,j3);
  40                     if gj < g0
  41                         % we are in lower triangle
  42                         % conecting from j1 j2 j3 to 0 0 0
  43                         t(1:3,2+j1,2+j2,2+j3)=[j1; j2; j3]; % from
  44                         t(4  ,2+j1,2+j2,2+j3)=t(4  ,2-j1,2-j2,2-j3); % to
  45                         % disp([gj,j1,j2,j3,squeeze(t(:,2+j1,2+j2,2+j3))'])
  46                         if t(4  ,2+j1,2+j2,2+j3) == 0
  47                             error('no target')
  48                         end
  49                     end
  50                 end
  51             end
  52         end
  53     otherwise
  54         m,error('unknown storage scheme, must be 14 or 27')
  55 end