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