Mercurial > repos > public > sbplib

comparison +scheme/Staggered1DAcousticsVariable.m @ 653:2351a7690e8a feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Clean up and improve comments in StaggeredVariable
author Martin Almquist <malmquist@stanford.edu>
date Mon, 04 Dec 2017 11:18:13 -0800
parents be941bb0a11a
children 82024d32e333
comparison
equal deleted inserted replaced
652:be941bb0a11a 653:2351a7690e8a
26 26
27 % Pick out primal or dual component 27 % Pick out primal or dual component
28 e_primal 28 e_primal
29 e_dual 29 e_dual
30 30
31 % System matrices 31 % System matrices, cell matrices of function handles
32 A 32 A
33 B 33 B
34 34
35 % 35 % Variable coefficient matrices evaluated at boundaries.
36 A_l, A_r 36 A_l, A_r
37 B_l, B_r 37 B_l, B_r
38 end 38 end
39 39
40 40
47 % 47 %
48 % B = B^T and with diagonal entries = 0. 48 % B = B^T and with diagonal entries = 0.
49 % B = [0 b 49 % B = [0 b
50 % b 0] 50 % b 0]
51 % Here we store p on the primal grid and v on the dual 51 % Here we store p on the primal grid and v on the dual
52 % A and B are cell matrices of function handles
52 function obj = Staggered1DAcousticsVariable(g, order, A, B) 53 function obj = Staggered1DAcousticsVariable(g, order, A, B)
53 default_arg('B',{@(x)0*x, @(x)0*x+1; @(x)0*x+1, @(x)0*x}); 54 default_arg('B',{@(x)0*x, @(x)0*x+1; @(x)0*x+1, @(x)0*x});
54 default_arg('A',{@(x)0*x+1, @(x)0*x; @(x)0*x, @(x)0*x+1}); 55 default_arg('A',{@(x)0*x+1, @(x)0*x; @(x)0*x, @(x)0*x+1});
55 56
56 obj.order = order; 57 obj.order = order;
60 % Grids 61 % Grids
61 obj.m = g.size(); 62 obj.m = g.size();
62 xl = g.getBoundary('l'); 63 xl = g.getBoundary('l');
63 xr = g.getBoundary('r'); 64 xr = g.getBoundary('r');
64 xlim = {xl, xr}; 65 xlim = {xl, xr};
66 obj.grid = g;
67 obj.grid_primal = g.grids{1};
68 obj.grid_dual = g.grids{2};
69 x_primal = obj.grid_primal.points();
70 x_dual = obj.grid_dual.points();
65 71
66 % Boundary matrices 72 % Boundary matrices
67 obj.A_l = [A{1,1}(xl), A{1,2}(xl);.... 73 obj.A_l = [A{1,1}(xl), A{1,2}(xl);....
68 A{2,1}(xl), A{2,2}(xl)]; 74 A{2,1}(xl), A{2,2}(xl)];
69 obj.A_r = [A{1,1}(xr), A{1,2}(xr);.... 75 obj.A_r = [A{1,1}(xr), A{1,2}(xr);....
70 A{2,1}(xr), A{2,2}(xr)]; 76 A{2,1}(xr), A{2,2}(xr)];
71 obj.B_l = [B{1,1}(xl), B{1,2}(xl);.... 77 obj.B_l = [B{1,1}(xl), B{1,2}(xl);....
72 B{2,1}(xl), B{2,2}(xl)]; 78 B{2,1}(xl), B{2,2}(xl)];
73 obj.B_r = [B{1,1}(xr), B{1,2}(xr);.... 79 obj.B_r = [B{1,1}(xr), B{1,2}(xr);....
74 B{2,1}(xr), B{2,2}(xr)]; 80 B{2,1}(xr), B{2,2}(xr)];
75
76 obj.grid = g;
77 obj.grid_primal = g.grids{1};
78 obj.grid_dual = g.grids{2};
79 x_primal = obj.grid_primal.points();
80 x_dual = obj.grid_dual.points();
81 81
82 % Get operators 82 % Get operators
83 ops = sbp.D1StaggeredUpwind(obj.m, xlim, order); 83 ops = sbp.D1StaggeredUpwind(obj.m, xlim, order);
84 obj.h = ops.h; 84 obj.h = ops.h;
85 85
136 % BC on the form Lu - g = 0; 136 % BC on the form Lu - g = 0;
137 137
138 % Get boundary matrices 138 % Get boundary matrices
139 switch boundary 139 switch boundary
140 case 'l' 140 case 'l'
141 A = full(obj.A_l); 141 A = obj.A_l;
142 B = full(obj.B_l); 142 B = obj.B_l;
143 case 'r' 143 case 'r'
144 A = full(obj.A_r); 144 A = obj.A_r;
145 B = full(obj.B_r); 145 B = obj.B_r;
146 end 146 end
147 147
148 % Diagonalize B 148 % Diagonalize B
149 [T, Lambda] = eig(B); 149 [T, Lambda] = eig(B);
150 lambda = diag(Lambda); 150 lambda = diag(Lambda);