Mercurial > repos > public > sbplib

comparison spdiagsVariablePeriodic.m @ 799:8c65ef13df89 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Bort med martins hemska whitespace
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 26 Jul 2018 17:51:48 -0700
parents 50e77b15d841
children 87ea9cac3287
comparison
equal deleted inserted replaced
797:5cf9fdf4c98f 799:8c65ef13df89
1 function A = spdiagsVariablePeriodic(vals,diags) 1 function A = spdiagsVariablePeriodic(vals,diags)
2 % Creates an m x m periodic discretization matrix. 2 % Creates an m x m periodic discretization matrix.
3 % vals - m x ndiags matrix of values 3 % vals - m x ndiags matrix of values
4 % diags - 1 x ndiags vector of the 'center diagonals' that vals end up on 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 5 % vals that are not on main diagonal are going to spill over to
6 % off-diagonal corners. 6 % off-diagonal corners.
7 7
8 default_arg('diags',0); 8 default_arg('diags',0);
9 9
10 [m, ~] = size(vals); 10 [m, ~] = size(vals);
11 11
12 A = sparse(m,m); 12 A = sparse(m,m);
13 13
14 for i = 1:length(diags) 14 for i = 1:length(diags)
15 15
16 d = diags(i); 16 d = diags(i);
17 a = vals(:,i); 17 a = vals(:,i);
18 18
19 % Sub-diagonals 19 % Sub-diagonals
20 if d < 0 20 if d < 0
21 a_bulk = a(1+abs(d):end); 21 a_bulk = a(1+abs(d):end);
22 a_corner = a(1:1+abs(d)-1); 22 a_corner = a(1:1+abs(d)-1);
23 corner_diag = m-abs(d); 23 corner_diag = m-abs(d);
24 A = A + spdiagVariable(a_bulk, d); 24 A = A + spdiagVariable(a_bulk, d);
25 A = A + spdiagVariable(a_corner, corner_diag); 25 A = A + spdiagVariable(a_corner, corner_diag);
26 26
27 % Super-diagonals 27 % Super-diagonals
28 elseif d > 0 28 elseif d > 0
29 a_bulk = a(1:end-d); 29 a_bulk = a(1:end-d);
30 a_corner = a(end-d+1:end); 30 a_corner = a(end-d+1:end);
31 corner_diag = -m + d; 31 corner_diag = -m + d;
32 A = A + spdiagVariable(a_bulk, d); 32 A = A + spdiagVariable(a_bulk, d);
33 A = A + spdiagVariable(a_corner, corner_diag); 33 A = A + spdiagVariable(a_corner, corner_diag);
34 34
35 % Main diagonal 35 % Main diagonal
36 else 36 else
37 A = A + spdiagVariable(a, 0); 37 A = A + spdiagVariable(a, 0);