Mercurial > repos > public > sbplib

comparison spdiagsPeriodic.m @ 812:6b83dcb46f54 feature/grids

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with feature/poroelastic
author Martin Almquist <malmquist@stanford.edu>
date Fri, 27 Jul 2018 10:31:51 -0700
parents bbf303c1f0cf
children
comparison
equal deleted inserted replaced
798:e76321b89c1e 812:6b83dcb46f54
1 function A = spdiagsPeriodic(vals,diags)
2 % Creates an m x m periodic discretization matrix.
3 % vals - m x ndiags matrix of values
4 % diags - 1 x ndiags vector of the 'center diagonals' that vals end up on
5 % vals that are not on main diagonal are going to spill over to
6 % off-diagonal corners.
7
8 default_arg('diags',0);
9
10 [m, ~] = size(vals);
11
12 A = sparse(m,m);
13
14 for i = 1:length(diags)
15
16 d = diags(i);
17 a = vals(:,i);
18
19 % Sub-diagonals
20 if d < 0
21 a_bulk = a(1+abs(d):end);
22 a_corner = a(1:1+abs(d)-1);
23 corner_diag = m-abs(d);
24 A = A + spdiagVariable(a_bulk, d);
25 A = A + spdiagVariable(a_corner, corner_diag);
26
27 % Super-diagonals
28 elseif d > 0
29 a_bulk = a(1:end-d);
30 a_corner = a(end-d+1:end);
31 corner_diag = -m + d;
32 A = A + spdiagVariable(a_bulk, d);
33 A = A + spdiagVariable(a_corner, corner_diag);
34
35 % Main diagonal
36 else
37 A = A + spdiagVariable(a, 0);
38 end
39
40 end
41
42 end
43
44 function A = spdiagVariable(a,i)
45 default_arg('i',0);
46
47 if isrow(a)
48 a = a';
49 end
50
51 n = length(a)+abs(i);
52
53 if i > 0
54 a = [sparse(i,1); a];
55 elseif i < 0
56 a = [a; sparse(abs(i),1)];
57 end
58
59 A = spdiags(a,i,n,n);
60 end