1function A = meshsparse (G, stencil) 2%MESHSPARSE convert a 2D or 3D mesh into a sparse matrix matrix. 3% 4% Example: 5% A = meshsparse (G) 6% A = meshsparse (G,5) % 2D 5-point stencil (default for 2D case) 7% A = meshsparse (G,9) % 2D 9-point stencil 8% A = meshsparse (G,7) % 3D 7-point stencil (default for 3D case) 9% A = meshsparse (G,27) % 3D 27-point stencil 10% A = meshsparse (G,stencil) % user-provided stencil 11% 12% To create a sparse matrix for an m-by-n 2D mesh or m-by-n-by-k 3D mesh, use 13% 14% A = meshsparse (meshnd (m,n)) ; 15% A = meshsparse (meshnd (m,n,k)) ; 16% 17% G is an m-by-n-by-k matrix, with entries numbered 1 to m*n*k (with k=1 for 18% the 2D case). The entries in G can appear in any order, but no duplicate 19% entries can appear. That is sort(G(:))' must equal 1:m*n*k. A is returned as 20% a sparse matrix with m*n*k rows and columns whose pattern depends on the 21% stencil. The number of nonzeros in most rows/columns of A is equal to the 22% number of points in the stencil. For examples on how to specify your own 23% stencil, see the contents of meshsparse.m. 24% 25% See also meshnd. 26 27% Copyright 2007-2009, Timothy A. Davis, http://www.suitesparse.com 28 29if (nargin < 2) 30 [m n k] = size (G) ; 31 if (k == 1) 32 stencil = 5 ; % 2D default is a 5-point stencil 33 else 34 stencil = 7 ; % 3D default is a 7-point stencil 35 end 36end 37 38if (numel (stencil) == 1) 39 40 % create the stencil 41 42 if (stencil == 5) 43 44 % 5-point stencil (2D) 45 stencil = [ 46 -1 0 % north 47 1 0 % south 48 0 1 % east 49 0 -1 ] ; % west 50 51 elseif (stencil == 9) 52 53 % 9-point stencil (2D) 54 stencil = [ 55 -1 0 % north 56 1 0 % south 57 0 1 % east 58 0 -1 % west 59 -1 -1 % north-west 60 -1 1 % north-east 61 1 -1 % south-west 62 1 1 ] ; % south-east 63 64 elseif (stencil == 7) 65 66 % 7-point stencil (3D) 67 stencil = [ 68 -1 0 0 % north 69 1 0 0 % south 70 0 1 0 % east 71 0 -1 0 % west 72 0 0 -1 % up 73 0 0 1] ; % down 74 75 elseif (stencil == 27) 76 77 % 27-point stencil (3D) 78 stencil = zeros (26, 3) ; 79 t = 0 ; 80 for i = -1:1 81 for j = -1:1 82 for k = -1:1 83 if (~(i == 0 & j == 0 & k == 0)) %#ok 84 t = t + 1 ; 85 stencil (t,:) = [i j k] ; 86 end 87 end 88 end 89 end 90 end 91end 92 93stencil = fix (stencil) ; 94[npoints d] = size (stencil) ; 95if (d == 2) 96 % append zeros onto a 2D stencil to make it "3D" 97 stencil = [stencil zeros(npoints,1)] ; 98end 99[npoints d] = size (stencil) ; 100if (d ~= 3) 101 error ('invalid stencil') ; 102end 103 104[m n k] = size (G) ; 105i1 = 1:m ; 106j1 = 1:n ; 107k1 = 1:k ; 108 109Ti = zeros (npoints*m*n*k, 1) ; 110Tj = zeros (npoints*m*n*k, 1) ; 111nz = 0 ; 112 113for point = 1:npoints 114 115 % find the overlapping rows of G 116 idelta = stencil (point,1) ; 117 i2 = i1 + idelta ; 118 ki = find (i2 >= 1 & i2 <= m) ; 119 120 % find the overlapping columns of G 121 jdelta = stencil (point,2) ; 122 j2 = j1 + jdelta ; 123 kj = find (j2 >= 1 & j2 <= n) ; 124 125 % find the overlapping slices of G 126 kdelta = stencil (point,3) ; 127 k2 = k1 + kdelta ; 128 kk = find (k2 >= 1 & k2 <= k) ; 129 130 % find the nodes in G the shifted G that touch 131 g2 = G (i2 (ki), j2 (kj), k2 (kk)) ; % shifted mesh 132 g1 = G (i1 (ki), j1 (kj), k1 (kk)) ; % unshifted mesh 133 134 % place the edges in the triplet list 135 e = numel (g1) ; 136 Ti ((nz+1):(nz+e)) = g1 (:) ; 137 Tj ((nz+1):(nz+e)) = g2 (:) ; 138 nz = nz + e ; 139end 140 141% convert the triplets into a sparse matrix 142Ti = Ti (1:nz) ; 143Tj = Tj (1:nz) ; 144A = npoints * speye (m*n*k) - sparse (Ti, Tj, 1, m*n*k, m*n*k) ; 145