femwind/ndt_convert.m

   1 function K2=ndt_convert(K,st2)
   2 % K2=K(K,st2)
   3 % convert ndt matrix from one storage format to another
   4 % st = storage type of output, 27 or 14 or 0=sparse
   5 switch st2
   6     case {27,'nons'}
   7         st2=27;
   8     case {14,'symm'}
   9         st2=14
  10     case {0,'sparse'}
  11         st2=0;
  12     otherwise
  13         st2,error('allowable values of st: 27 or nons, 14 or symm, 0 or sparse')
  14 end
  15 [n1,n2,n3,st1]=size(K);
  16 t1=ndt_storage_table(st1);
  17 if st2, % ndt storage
  18     t2=ndt_storage_table(st2);
  19     K2=zeros(n1,n2,n3,st2);
  20 else  % sparse storage
  21     n=n1*n2*n3;
  22     mnz=n*27; % max nonzeros
  23     ii=zeros(mnz,1);
  24     jj=zeros(mnz,1);
  25     vv=zeros(mnz,1);
  26     nn=0;
  27     % K2=sparse([],[],[],n,n,mnz);
  28 end
  29 g=@(i1,i2,i3)i1+n1*(i2-1+n2*(i3-1));  % flattened index in 1:n by 1:n matrix
  30 for i3=1:n3
  31     for i2=1:n2
  32         for i1=1:n1
  33             for j3=-1:1
  34                 for j2=-1:1
  35                     for j1=-1:1
  36                         % global index of neighbor node k
  37                         k3 = i3+j3;
  38                         k2 = i2+j2;
  39                         k1 = i1+j1;
  40                         % relative index of the row where the entry (i,k)
  41                         % is stored, 0 0 0 = this row
  42                         % in fortran we won't have the 2+ because
  43                         % the array t will be indexed -1:1
  44                         % row m where the entry (i,k) is stored in K
  45                         m3 = i3+t1(3,2+j1,2+j2,2+j3);
  46                         m2 = i2+t1(2,2+j1,2+j2,2+j3);
  47                         m1 = i1+t1(1,2+j1,2+j2,2+j3);
  48                         % index of the matrix entry (i,k) in K(m,:)
  49                         jx = t1(4,2+j1,2+j2,2+j3);
  50                         if ~(m1<1 || m1>n1 || m2<1 || m2>n2 || m3<1 || m3>n3 || ...
  51                              k1<1 || k1>n1 || k2<1 || k2>n2 || k3<1 || k3>n3 )
  52                            % j is within bounds
  53                             v = K(m1,m2,m3,jx); % value of the entry
  54                             if st2  % ndt storage scheme
  55                                 % entry (i,k) is stored in K2(l,p)
  56                                 l3=i3+t2(3,2+j1,2+j2,2+j3); % row + offset
  57                                 l2=i2+t2(2,2+j1,2+j2,2+j3); % row + offset
  58                                 l1=i1+t2(1,2+j1,2+j2,2+j3); % row + offset
  59                                 if ~(l1<1 || l1>n1 || l2<1 || l2>n2 || l3<1 || l3>n3)
  60                                     px = t2(4,2+j1,2+j2,2+j3); % index in the row
  61                                     K2(l1,l2,l3,px)=v;
  62                                 end
  63                             else % sparse
  64                                 i=g(i1,i2,i3);
  65                                 k=g(k1,k2,k3);
  66                                 % K2(i,k)=v;
  67                                 nn=nn+1;
  68                                 ii(nn)=i;
  69                                 jj(nn)=k;
  70                                 vv(nn)=v;
  71                             end
  72                         end
  73                     end
  74                 end
  75             end
  76         end
  77     end
  78 end
  79 if st2==0
  80     K2 = sparse(ii(1:nn),jj(1:nn),vv(1:nn),n,n);
  81 end
  82 end
  83