femwind/ndt_convert.m

   1 function K2=ndt_convert(K,n1,n2,n3,st)
   2 % K2=ndt_convert(K,st)
   3 % K2=ndt_convert(K,n1,n2,n3,st)
   4 % convert ndt matrix from one storage format to another
   5 % st = storage type of output, 27 or 14 or 0=sparse
   6 % n1 n2 n3 = dimensions of output, if K is sparse
   7 if issparse(K) % storage type of input:
   8     st1 = 0;
   9     [m,n]=size(K);
  10     if n1*n2*n3 ~= n
  11         n1,n2,n3,n
  12         error('input must satisfy n1*n2*n3 == n')
  13     end
  14     if(big(K-K'))
  15         error('Sparse input matrix must be symmetric')
  16     end
  17     Kg=K;
  18     t1=ndt_storage_table(27);  % need to loop over all output entries
  19 else
  20     st = n1;
  21     [n1,n2,n3,st1]=size(K);
  22     t1=ndt_storage_table(st1);
  23 end
  24 switch st %  set storage type of output: 0 or 14 or 27
  25     case {27,'nons'}
  26         st2=27;
  27     case {14,'symm'}
  28         st2=14;
  29     case {0,'sparse'}
  30         st2=0;
  31     otherwise
  32         st,error('allowable values of st: 27 or nons, 14 or symm, 0 or sparse')
  33 end
  34 if st2, % ndt storage
  35     t2=ndt_storage_table(st2);
  36     K2=zeros(n1,n2,n3,st2);
  37 else  % sparse storage
  38     n=n1*n2*n3;
  39     mnz=n*27; % max nonzeros
  40     ii=zeros(mnz,1);
  41     jj=zeros(mnz,1);
  42     vv=zeros(mnz,1);
  43     nn=0;
  44     % K2=sparse([],[],[],n,n,mnz);
  45 end
  46 g=@(i1,i2,i3)i1+n1*(i2-1+n2*(i3-1));  % flattened index in 1:n by 1:n matrix
  47 for i3=1:n3
  48     for i2=1:n2
  49         for i1=1:n1
  50             for j3=-1:1
  51                 for j2=-1:1
  52                     for j1=-1:1
  53                         % global index of neighbor node k
  54                         k3 = i3+j3;
  55                         k2 = i2+j2;
  56                         k1 = i1+j1;
  57                         % relative index of the row where the entry (i,k)
  58                         % is stored, 0 0 0 = this row
  59                         % in fortran we won't have the 2+ because
  60                         % the array t will be indexed -1:1
  61                         % row m where the entry (i,k) is stored in K
  62                         m3 = i3+t1(3,2+j1,2+j2,2+j3);
  63                         m2 = i2+t1(2,2+j1,2+j2,2+j3);
  64                         m1 = i1+t1(1,2+j1,2+j2,2+j3);
  65                         % index of the matrix entry (i,k) in K(m,:)
  66                         jx = t1(4,2+j1,2+j2,2+j3);
  67                         if ~(m1<1 || m1>n1 || m2<1 || m2>n2 || m3<1 || m3>n3 || ...
  68                              k1<1 || k1>n1 || k2<1 || k2>n2 || k3<1 || k3>n3 )
  69                             % j is within bounds
  70                             if st1   % ndt scheme
  71                                 v = K(m1,m2,m3,jx); % value of the entry in ndt
  72                             else  % sparse
  73                                 ig=g(k1,k2,k3);
  74                                 jg=g(m1,m2,m3);
  75                                 v = K(ig,jg);
  76                                 Kg(ig,jg)=0;  % mark location used
  77                             end
  78                             if st2  % ndt storage scheme
  79                                 % entry (i,k) is stored in K2(l,p)
  80                                 l3=i3+t2(3,2+j1,2+j2,2+j3); % row + offset
  81                                 l2=i2+t2(2,2+j1,2+j2,2+j3); % row + offset
  82                                 l1=i1+t2(1,2+j1,2+j2,2+j3); % row + offset
  83                                 if ~(l1<1 || l1>n1 || l2<1 || l2>n2 || l3<1 || l3>n3)
  84                                     px = t2(4,2+j1,2+j2,2+j3); % index in the row
  85                                     K2(l1,l2,l3,px)=v;
  86                                 end
  87                             else % sparse
  88                                 i=g(i1,i2,i3);
  89                                 k=g(k1,k2,k3);
  90                                 % K2(i,k)=v;
  91                                 nn=nn+1;
  92                                 ii(nn)=i;
  93                                 jj(nn)=k;
  94                                 vv(nn)=v;
  95                             end
  96                         end
  97                     end
  98                 end
  99             end
 100         end
 101     end
 102 end
 103 if st1==0
 104     if(nnz(Kg))
 105         error('K is not in the 3D grid form, nonzeros left')
 106     end
 107 end
 108 if st2==0
 109     K2 = sparse(ii(1:nn),jj(1:nn),vv(1:nn),n,n);
 110 end
 111 end
 112