Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinear.m @ 560:6132c52bf923 feature/grids/laplace_refactor

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 29 Aug 2017 13:04:07 +0200
parents e3f587bfd7f8 54c775c3348a
children 3a13916f8ff0
comparison
equal deleted inserted replaced
559:e3f587bfd7f8 560:6132c52bf923
233 [e, d_n, d_t, a_n, a_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t] = obj.get_boundary_ops(boundary); 233 [e, d_n, d_t, a_n, a_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t] = obj.get_boundary_ops(boundary);
234 234
235 A_n = spdiag(a_n); 235 A_n = spdiag(a_n);
236 A_t = spdiag(a_t); 236 A_t = spdiag(a_t);
237 237
238 F = (s * A_n * d_n' + s * A_t*d_t')'; 238 F = s*(A_n*d_n' + A_t*d_t')';
239 239
240 u = obj; 240 u = obj;
241 241
242 b1 = gamm*u.lambda./u.a11.^2; 242 b1 = gamm*u.lambda./u.a11.^2;
243 b2 = gamm*u.lambda./u.a22.^2; 243 b2 = gamm*u.lambda./u.a22.^2;
254 254
255 % Neumann boundary condition 255 % Neumann boundary condition
256 case {'N','n','neumann'} 256 case {'N','n','neumann'}
257 A_n = spdiag(a_n); 257 A_n = spdiag(a_n);
258 A_t = spdiag(a_t); 258 A_t = spdiag(a_t);
259 d = (A_n * d_n' + A_t*d_t')'; 259 d = s*(A_n * d_n' + A_t*d_t')';
260 260
261 tau1 = -s; 261 tau1 = -1;
262 tau2 = 0; 262 tau2 = 0;
263 tau = obj.a * obj.Ji*(tau1*e + tau2*d); 263 tau = s*obj.a*obj.Ji*(tau1*e + tau2*d);
264 264
265 closure = halfnorm_inv*tau*d'; 265 closure = halfnorm_inv*tau*d';
266 penalty = -halfnorm_inv*tau; 266 penalty = -halfnorm_inv*tau;
267 267
268 % Characteristic boundary condition 268 % Characteristic boundary condition
297 A_n_u = spdiag(a_n_u); 297 A_n_u = spdiag(a_n_u);
298 A_t_u = spdiag(a_t_u); 298 A_t_u = spdiag(a_t_u);
299 A_n_v = spdiag(a_n_v); 299 A_n_v = spdiag(a_n_v);
300 A_t_v = spdiag(a_t_v); 300 A_t_v = spdiag(a_t_v);
301 301
302 F_u = (s_u * A_n_u * d_n_u' + s_u * A_t_u*d_t_u')'; 302 F_u = s_u*(A_n_u * d_n_u' + A_t_u*d_t_u')';
303 F_v = (s_v * A_n_v * d_n_v' + s_v * A_t_v*d_t_v')'; 303 F_v = s_v*(A_n_v * d_n_v' + A_t_v*d_t_v')';
304 304
305 u = obj; 305 u = obj;
306 v = neighbour_scheme; 306 v = neighbour_scheme;
307 307
308 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2; 308 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2;