← Back to branch summary

~logan/ubuntu/trusty/suitesparse/4.2.1-3ubuntu1

« back to all changes in this revision

Viewing changes to MESHND/meshsparse.m

  • Committer: Bazaar Package Importer
  • Author(s): Rafael Laboissiere, Rafael Laboissiere, Ondrej Certik
  • Date: 2008-02-21 14:46:50 UTC
  • mfrom: (1.1.2 upstream) (5.1.1 hardy)
  • Revision ID: james.westby@ubuntu.com-20080221144650-tgeppgj0t7s759i8
Tags: 3.1.0-3
[ Rafael Laboissiere ]
* Upload to unstable

[ Ondrej Certik ]
* XS-DM-Upload-Allowed: yes field added
* Ondrej Certik added to uploaders

Show diffs side-by-side

added added

removed removed

Lines of Context:
MESHND/meshsparse.m
 
1
function 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, Timothy A. Davis, Univ. of Florida
 
28
 
 
29
if (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
 
36
end
 
37
 
 
38
if (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
 
91
end
 
92
 
 
93
stencil = fix (stencil) ;
 
94
[npoints d] = size (stencil) ;
 
95
if (d == 2)
 
96
    % append zeros onto a 2D stencil to make it "3D"
 
97
    stencil = [stencil zeros(npoints,1)] ;
 
98
end
 
99
[npoints d] = size (stencil) ;
 
100
if (d ~= 3)
 
101
    error ('invalid stencil') ;
 
102
end
 
103
 
 
104
[m n k] = size (G) ;
 
105
i1 = 1:m ;
 
106
j1 = 1:n ;
 
107
k1 = 1:k ;
 
108
 
 
109
Ti = zeros (npoints*m*n*k, 1) ;
 
110
Tj = zeros (npoints*m*n*k, 1) ;
 
111
nz = 0 ;
 
112
 
 
113
for 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 ;
 
139
end
 
140
 
 
141
% convert the triplets into a sparse matrix
 
142
Ti = Ti (1:nz) ;
 
143
Tj = Tj (1:nz) ;
 
144
A = npoints * speye (m*n*k) - sparse (Ti, Tj, 1, m*n*k, m*n*k) ;