Mercurial > repos > public > sbplib

comparison +scheme/Staggered1DAcousticsVariable.m @ 672:82024d32e333 feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Correct char BC in staggered acoustics 1D variable coeff.
author Martin Almquist <malmquist@stanford.edu>
date Wed, 27 Dec 2017 21:16:06 +0100
parents 2351a7690e8a
children ec0ac76006e2
comparison
equal deleted inserted replaced
671:993aac771efd 672:82024d32e333
142 B = obj.B_l; 142 B = obj.B_l;
143 case 'r' 143 case 'r'
144 A = obj.A_r; 144 A = obj.A_r;
145 B = obj.B_r; 145 B = obj.B_r;
146 end 146 end
147 C = inv(A)*B;
147 148
148 % Diagonalize B 149 % Diagonalize B
149 [T, Lambda] = eig(B); 150 [T, Lambda] = eig(B);
150 lambda = diag(Lambda); 151 lambda = diag(Lambda);
151 152
153 Iplus = lambda > 0; 154 Iplus = lambda > 0;
154 Iminus = lambda < 0; 155 Iminus = lambda < 0;
155 156
156 switch boundary 157 switch boundary
157 case 'l' 158 case 'l'
158 Iout = Iminus;
159 Iin = Iplus; 159 Iin = Iplus;
160 case 'r' 160 case 'r'
161 Iout = Iplus;
162 Iin = Iminus; 161 Iin = Iminus;
163 end 162 end
164
165 Tin = T(:,Iin); 163 Tin = T(:,Iin);
166 Tout = T(:,Iout);
167 164
168 switch type 165 switch type
169 case 'p' 166 case 'p'
170 L = [1, 0]; 167 L = [1, 0];
171 case 'v' 168 case 'v'
172 L = [0, 1]; 169 L = [0, 1];
173 case 'characteristic' 170 case 'characteristic'
174 L = Tin'; 171 % Diagonalize C
172 [T_C, Lambda_C] = eig(C);
173 lambda_C = diag(Lambda_C);
174 % Identify in- and outgoing characteristic variables
175 Iplus = lambda_C > 0;
176 Iminus = lambda_C < 0;
177
178 switch boundary
179 case 'l'
180 Iin_C = Iplus;
181 case 'r'
182 Iin_C = Iminus;
183 end
184 T_C_inv = inv(T_C);
185 L = T_C_inv(Iin_C,:);
175 otherwise 186 otherwise
176 error('Boundary condition not implemented.'); 187 error('Boundary condition not implemented.');
177 end 188 end
178 189
179 % Penalty parameters 190 % Penalty parameters
180 sigma = [0; 0]; 191 sigma = [0; 0];
181 sigma(Iin) = lambda(Iin); 192 sigma(Iin) = lambda(Iin);
193
194 % Sparsify
195 A = sparse(A);
196 T = sparse(T);
197 sigma = sparse(sigma);
198 L = sparse(L);
199 Tin = sparse(Tin);
200
182 switch boundary 201 switch boundary
183 case 'l' 202 case 'l'
184 tau = -1*obj.e_l * inv(A) * T * sigma * inv(L*Tin); 203 tau = -1*obj.e_l * inv(A) * T * sigma * inv(L*Tin);
185 closure = obj.Hi*tau*L*obj.e_l'; 204 closure = obj.Hi*tau*L*obj.e_l';
186 205