femwind/hexa_sub.m

   1 function p=hexa_sub(X,K,i)
   2 % p=hexa_sub(X,K,i)
   3 % extract one submatrix from global matrix from i to i+1
   4 % in:
   5 %    X mesh coordinates, cell vector 3
   6 %    K global matrix
   7 %    i index of lower left front corner
   8 % out:
   9 %    p.S2 2D submatrix of K
  10 %    p.X  3 by 8 matrix of coordinates, same as in hexa
  11 %    p.S6 submatrix of K as 6D
  12 n= size(X{1});
  13 nn = prod(n);
  14 if any(size(K)~=nn),
  15     error('wrong size K')
  16 end
  17 p.S6=zeros(2,2,2,2,2,2);
  18 p.S2=zeros(8,8);
  19 p.X=zeros(3,8);
  20 for i1=0:1
  21     for i2=0:1
  22         for i3=0:1
  23             ix=sub2ind(n,i(1)+i1,i(2)+i2,i(3)+i3); % global index
  24             ik=sub2ind([2,2,2],1+i1,1+i2,1+i3);    % local index
  25             p.X(1,ik)=X{1}(i(1)+i1,i(2)+i2,i(3)+i3);
  26             p.X(2,ik)=X{2}(i(1)+i1,i(2)+i2,i(3)+i3);
  27             p.X(3,ik)=X{3}(i(1)+i1,i(2)+i2,i(3)+i3);
  28             for j1=0:1
  29                 for j2=0:1
  30                     for j3=0:1
  31                         jx=sub2ind(n,i(1)+j1,i(2)+j2,i(3)+j3);
  32                         jk=sub2ind([2,2,2],1+j1,1+j2,1+j3);
  33                         p.S2(ik,jk)=K(ix,jx);
  34                         p.S6(i1+1,i2+1,i3+1,j1+1,j2+1,j3+1)=K(ix,jx);
  35                     end
  36                 end
  37             end
  38         end
  39     end
  40 end
  41 end