Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinear.m @ 554:6e71d4f8b155 feature/grids/laplace_refactor

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Better names for the metric coefficients
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 29 Aug 2017 11:21:02 +0200
parents 63d5c07943dd
children 2a856a589510
comparison
equal deleted inserted replaced
553:63d5c07943dd 554:6e71d4f8b155
168 obj.e_w = e_w; 168 obj.e_w = e_w;
169 obj.e_e = e_e; 169 obj.e_e = e_e;
170 obj.e_s = e_s; 170 obj.e_s = e_s;
171 obj.e_n = e_n; 171 obj.e_n = e_n;
172 172
173
173 obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv; 174 obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv;
174 obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv; 175 obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv;
175 176
176 177
177 % Misc. 178 % Misc.
201 % neighbour_boundary is a string specifying which boundary to interface to. 202 % neighbour_boundary is a string specifying which boundary to interface to.
202 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) 203 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
203 default_arg('type','neumann'); 204 default_arg('type','neumann');
204 default_arg('parameter', []); 205 default_arg('parameter', []);
205 206
206 [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv , ~, ~, ~, scale_factor] = obj.get_boundary_ops(boundary); 207 [e, d_n, d_t, a_n, a_t, s, gamm, halfnorm_inv , ~, ~, ~, scale_factor] = obj.get_boundary_ops(boundary);
207 switch type 208 switch type
208 % Dirichlet boundary condition 209 % Dirichlet boundary condition
209 case {'D','d','dirichlet'} 210 case {'D','d','dirichlet'}
210 % v denotes the solution in the neighbour domain 211 % v denotes the solution in the neighbour domain
211 tuning = 1.2; 212 tuning = 1.2;
212 % tuning = 20.2; 213 % tuning = 20.2;
213 [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t] = obj.get_boundary_ops(boundary); 214 [e, d_n, d_t, a_n, a_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t] = obj.get_boundary_ops(boundary);
214 215
215 a_n = spdiag(coeff_n); 216 A_n = spdiag(a_n);
216 a_t = spdiag(coeff_t); 217 A_t = spdiag(a_t);
217 218
218 F = (s * a_n * d_n' + s * a_t*d_t')'; 219 F = (s * A_n * d_n' + s * A_t*d_t')';
219 220
220 u = obj; 221 u = obj;
221 222
222 b1 = gamm*u.lambda./u.a11.^2; 223 b1 = gamm*u.lambda./u.a11.^2;
223 b2 = gamm*u.lambda./u.a22.^2; 224 b2 = gamm*u.lambda./u.a22.^2;
232 penalty = -obj.Ji*obj.a * penalty_parameter_1; 233 penalty = -obj.Ji*obj.a * penalty_parameter_1;
233 234
234 235
235 % Neumann boundary condition 236 % Neumann boundary condition
236 case {'N','n','neumann'} 237 case {'N','n','neumann'}
237 a_n = spdiag(coeff_n); 238 A_n = spdiag(a_n);
238 a_t = spdiag(coeff_t); 239 A_t = spdiag(a_t);
239 d = (a_n * d_n' + a_t*d_t')'; 240 d = (A_n * d_n' + A_t*d_t')';
240 241
241 tau1 = -s; 242 tau1 = -s;
242 tau2 = 0; 243 tau2 = 0;
243 tau = obj.a * obj.Ji*(tau1*e + tau2*d); 244 tau = obj.a * obj.Ji*(tau1*e + tau2*d);
244 245
248 % Characteristic boundary condition 249 % Characteristic boundary condition
249 case {'characteristic', 'char', 'c'} 250 case {'characteristic', 'char', 'c'}
250 default_arg('parameter', 1); 251 default_arg('parameter', 1);
251 beta = parameter; 252 beta = parameter;
252 253
253 a_n = spdiag(coeff_n); 254 A_n = spdiag(a_n);
254 a_t = spdiag(coeff_t); 255 A_t = spdiag(a_t);
255 d = s*(a_n * d_n' + a_t*d_t')'; % outward facing normal derivative 256 d = s*(A_n * d_n' + A_t*d_t')'; % outward facing normal derivative
256 257
257 tau = -obj.a * 1/beta*obj.Ji*e; 258 tau = -obj.a * 1/beta*obj.Ji*e;
258 259
259 closure{1} = halfnorm_inv*tau*spdiag(scale_factor)*e'; 260 closure{1} = halfnorm_inv*tau*spdiag(scale_factor)*e';
260 closure{2} = halfnorm_inv*tau*beta*d'; 261 closure{2} = halfnorm_inv*tau*beta*d';
269 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 270 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
270 % u denotes the solution in the own domain 271 % u denotes the solution in the own domain
271 % v denotes the solution in the neighbour domain 272 % v denotes the solution in the neighbour domain
272 tuning = 1.2; 273 tuning = 1.2;
273 % tuning = 20.2; 274 % tuning = 20.2;
274 [e_u, d_n_u, d_t_u, coeff_n_u, coeff_t_u, s_u, gamm_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t, I_u] = obj.get_boundary_ops(boundary); 275 [e_u, d_n_u, d_t_u, a_n_u, a_t_u, s_u, gamm_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t, I_u] = obj.get_boundary_ops(boundary);
275 [e_v, d_n_v, d_t_v, coeff_n_v, coeff_t_v, s_v, gamm_v, halfnorm_inv_v_n, halfnorm_inv_v_t, halfnorm_v_t, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 276 [e_v, d_n_v, d_t_v, a_n_v, a_t_v, s_v, gamm_v, halfnorm_inv_v_n, halfnorm_inv_v_t, halfnorm_v_t, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
276 277
277 a_n_u = spdiag(coeff_n_u); 278 A_n_u = spdiag(a_n_u);
278 a_t_u = spdiag(coeff_t_u); 279 A_t_u = spdiag(a_t_u);
279 a_n_v = spdiag(coeff_n_v); 280 A_n_v = spdiag(a_n_v);
280 a_t_v = spdiag(coeff_t_v); 281 A_t_v = spdiag(a_t_v);
281 282
282 F_u = (s_u * a_n_u * d_n_u' + s_u * a_t_u*d_t_u')'; 283 F_u = (s_u * A_n_u * d_n_u' + s_u * A_t_u*d_t_u')';
283 F_v = (s_v * a_n_v * d_n_v' + s_v * a_t_v*d_t_v')'; 284 F_v = (s_v * A_n_v * d_n_v' + s_v * A_t_v*d_t_v')';
284 285
285 u = obj; 286 u = obj;
286 v = neighbour_scheme; 287 v = neighbour_scheme;
287 288
288 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2; 289 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2;
305 306
306 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 307 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
307 % The right boundary is considered the positive boundary 308 % The right boundary is considered the positive boundary
308 % 309 %
309 % I -- the indecies of the boundary points in the grid matrix 310 % I -- the indecies of the boundary points in the grid matrix
310 function [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t, I, scale_factor] = get_boundary_ops(obj, boundary) 311 function [e, d_n, d_t, a_n, a_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t, I, scale_factor] = get_boundary_ops(obj, boundary)
311 312
312 % gridMatrix = zeros(obj.m(2),obj.m(1)); 313 % gridMatrix = zeros(obj.m(2),obj.m(1));
313 % gridMatrix(:) = 1:numel(gridMatrix); 314 % gridMatrix(:) = 1:numel(gridMatrix);
314 315
315 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m)); 316 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
320 d_n = obj.du_w; 321 d_n = obj.du_w;
321 d_t = obj.dv_w; 322 d_t = obj.dv_w;
322 s = -1; 323 s = -1;
323 324
324 I = ind(1,:); 325 I = ind(1,:);
325 coeff_n = obj.a11(I); 326 a_n = obj.a11(I);
326 coeff_t = obj.a12(I); 327 a_t = obj.a12(I);
327 scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2); 328 scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2);
328 case 'e' 329 case 'e'
329 e = obj.e_e; 330 e = obj.e_e;
330 d_n = obj.du_e; 331 d_n = obj.du_e;
331 d_t = obj.dv_e; 332 d_t = obj.dv_e;
332 s = 1; 333 s = 1;
333 334
334 I = ind(end,:); 335 I = ind(end,:);
335 coeff_n = obj.a11(I); 336 a_n = obj.a11(I);
336 coeff_t = obj.a12(I); 337 a_t = obj.a12(I);
337 scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2); 338 scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2);
338 case 's' 339 case 's'
339 e = obj.e_s; 340 e = obj.e_s;
340 d_n = obj.dv_s; 341 d_n = obj.dv_s;
341 d_t = obj.du_s; 342 d_t = obj.du_s;
342 s = -1; 343 s = -1;
343 344
344 I = ind(:,1)'; 345 I = ind(:,1)';
345 coeff_n = obj.a22(I); 346 a_n = obj.a22(I);
346 coeff_t = obj.a12(I); 347 a_t = obj.a12(I);
347 scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2); 348 scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2);
348 case 'n' 349 case 'n'
349 e = obj.e_n; 350 e = obj.e_n;
350 d_n = obj.dv_n; 351 d_n = obj.dv_n;
351 d_t = obj.du_n; 352 d_t = obj.du_n;
352 s = 1; 353 s = 1;
353 354
354 I = ind(:,end)'; 355 I = ind(:,end)';
355 coeff_n = obj.a22(I); 356 a_n = obj.a22(I);
356 coeff_t = obj.a12(I); 357 a_t = obj.a12(I);
357 scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2); 358 scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2);
358 otherwise 359 otherwise
359 error('No such boundary: boundary = %s',boundary); 360 error('No such boundary: boundary = %s',boundary);
360 end 361 end
361 362