Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinearNewCorner.m @ 1067:9a858436f8fa feature/laplace_curvilinear_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Implement new penalty strength for interface. Bugfix missing coeff a in Dirichlet penalty.
author Martin Almquist <malmquist@stanford.edu>
date Tue, 22 Jan 2019 18:17:01 -0800
parents d64062bed5fb
children e0ecce90f8cf
comparison
equal deleted inserted replaced
1066:d64062bed5fb 1067:9a858436f8fa
272 272
273 e = obj.getBoundaryOperator('e', boundary); 273 e = obj.getBoundaryOperator('e', boundary);
274 d = obj.getBoundaryOperator('d', boundary); 274 d = obj.getBoundaryOperator('d', boundary);
275 H_b = obj.getBoundaryQuadrature(boundary); 275 H_b = obj.getBoundaryQuadrature(boundary);
276 s_b = obj.getBoundaryScaling(boundary); 276 s_b = obj.getBoundaryScaling(boundary);
277 [th_H, th_M, th_R] = obj.getBoundaryBorrowing(boundary); 277 [th_H, ~, th_R] = obj.getBoundaryBorrowing(boundary);
278 m = obj.getBoundaryNumber(boundary); 278 m = obj.getBoundaryNumber(boundary);
279 279
280 K = obj.K; 280 K = obj.K;
281 J = obj.J; 281 J = obj.J;
282 Hi = obj.Hi; 282 Hi = obj.Hi;
283 a_b = e'*obj.a*e; 283 a = obj.a;
284 284
285 switch type 285 switch type
286 % Dirichlet boundary condition 286 % Dirichlet boundary condition
287 case {'D','d','dirichlet'} 287 case {'D','d','dirichlet'}
288 tuning = 1.0; 288 tuning = 1.0;
300 tau_H(1,1) = obj.dim*tau_H(1,1); 300 tau_H(1,1) = obj.dim*tau_H(1,1);
301 tau_H(end,end) = obj.dim*tau_H(end,end); 301 tau_H(end,end) = obj.dim*tau_H(end,end);
302 302
303 tau = tuning*(tau_R + tau_H); 303 tau = tuning*(tau_R + tau_H);
304 304
305 closure = Hi*d*a_b*H_b*e' ... 305 closure = a*Hi*d*H_b*e' ...
306 -Hi*e*tau*H_b*e'; 306 -a*Hi*e*tau*H_b*e';
307 307
308 penalty = -Hi*d*a_b*H_b ... 308 penalty = -a*Hi*d*H_b ...
309 +Hi*e*tau*H_b; 309 +a*Hi*e*tau*H_b;
310 310
311 311
312 % Neumann boundary condition 312 % Neumann boundary condition
313 case {'N','n','neumann'} 313 case {'N','n','neumann'}
314 tau1 = -1; 314 tau1 = -1;
315 tau2 = 0; 315 tau2 = 0;
316 tau = (tau1*e + tau2*d)*H_b; 316 tau = (tau1*e + tau2*d)*H_b;
317 317
318 closure = obj.a*Hi*tau*d'; 318 closure = a*Hi*tau*d';
319 penalty = -obj.a*Hi*tau; 319 penalty = -a*Hi*tau;
320 320
321 321
322 % Unknown, boundary condition 322 % Unknown, boundary condition
323 otherwise 323 otherwise
324 error('No such boundary condition: type = %s',type); 324 error('No such boundary condition: type = %s',type);
331 % -- interpolation: type of interpolation, default 'none' 331 % -- interpolation: type of interpolation, default 'none'
332 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) 332 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
333 333
334 % error('Not implemented') 334 % error('Not implemented')
335 335
336 defaultType.tuning = 1.2; 336 defaultType.tuning = 1.0;
337 defaultType.interpolation = 'none'; 337 defaultType.interpolation = 'none';
338 default_struct('type', defaultType); 338 default_struct('type', defaultType);
339 339
340 switch type.interpolation 340 switch type.interpolation
341 case {'none', ''} 341 case {'none', ''}
348 end 348 end
349 349
350 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type) 350 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
351 tuning = type.tuning; 351 tuning = type.tuning;
352 352
353 dim = obj.dim;
353 % u denotes the solution in the own domain 354 % u denotes the solution in the own domain
354 % v denotes the solution in the neighbour domain 355 % v denotes the solution in the neighbour domain
355 e_u = obj.getBoundaryOperator('e', boundary);
356 d_u = obj.getBoundaryOperator('d', boundary);
357 H_b_u = obj.getBoundaryQuadrature(boundary);
358 I_u = obj.getBoundaryIndices(boundary);
359 gamm_u = obj.getBoundaryBorrowing(boundary);
360
361 e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
362 d_v = neighbour_scheme.getBoundaryOperator('d', neighbour_boundary);
363 H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary);
364 I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary);
365 gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary);
366
367 u = obj; 356 u = obj;
368 v = neighbour_scheme; 357 v = neighbour_scheme;
369 358
370 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2; 359 % Boundary operators, u
371 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2; 360 e_u = u.getBoundaryOperator('e', boundary);
372 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2; 361 d_u = u.getBoundaryOperator('d', boundary);
373 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2; 362 gamm_u = u.getBoundaryBorrowing(boundary);
374 363 s_b_u = u.getBoundaryScaling(boundary);
375 tau1 = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v); 364 [th_H_u, ~, th_R_u] = u.getBoundaryBorrowing(boundary);
376 tau1 = tuning * spdiag(tau1); 365 m_u = u.getBoundaryNumber(boundary);
377 tau2 = 1/2; 366
378 367 % Coefficients, u
379 sig1 = -1/2; 368 K_u = u.K;
380 sig2 = 0; 369 J_u = u.J;
381 370
382 tau = (e_u*tau1 + tau2*d_u)*H_b_u; 371 % Boundary operators, v
383 sig = (sig1*e_u + sig2*d_u)*H_b_u; 372 e_v = v.getBoundaryOperator('e', neighbour_boundary);
384 373 d_v = v.getBoundaryOperator('d', neighbour_boundary);
385 closure = obj.a*obj.Hi*( tau*e_u' + sig*d_u'); 374 gamm_v = v.getBoundaryBorrowing(neighbour_boundary);
386 penalty = obj.a*obj.Hi*(-tau*e_v' + sig*d_v'); 375 s_b_v = v.getBoundaryScaling(neighbour_boundary);
376 [th_H_v, ~, th_R_v] = v.getBoundaryBorrowing(neighbour_boundary);
377 m_v = v.getBoundaryNumber(neighbour_boundary);
378
379 % Coefficients, v
380 K_v = v.K;
381 J_v = v.J;
382
383 %--- Penalty strength tau -------------
384 sigma_u = 0;
385 sigma_v = 0;
386 for i = 1:obj.dim
387 sigma_u = sigma_u + e_u'*J_u*K_u{i,m_u}*K_u{i,m_u}*e_u;
388 sigma_v = sigma_v + e_v'*J_v*K_v{i,m_v}*K_v{i,m_v}*e_v;
389 end
390 sigma_u = sigma_u/s_b_u;
391 sigma_v = sigma_v/s_b_v;
392
393 tau_R_u = 1/th_R_u*sigma_u;
394 tau_R_v = 1/th_R_v*sigma_v;
395
396 tau_H_u = 1/th_H_u*sigma_u;
397 tau_H_u(1,1) = dim*tau_H_u(1,1);
398 tau_H_u(end,end) = dim*tau_H_u(end,end);
399
400 tau_H_v = 1/th_H_v*sigma_v;
401 tau_H_v(1,1) = dim*tau_H_v(1,1);
402 tau_H_v(end,end) = dim*tau_H_v(end,end);
403
404 tau = 1/4*tuning*(tau_R_u + tau_H_u + tau_R_v + tau_H_v);
405 %--------------------------------------
406
407 % Operators/coefficients that are only required from this side
408 Hi = u.Hi;
409 H_b = u.getBoundaryQuadrature(boundary);
410 a = u.a;
411
412 closure = 1/2*a*Hi*d_u*H_b*e_u' ...
413 -1/2*a*Hi*e_u*H_b*d_u' ...
414 -a*Hi*e_u*tau*H_b*e_u';
415
416 penalty = -1/2*a*Hi*d_u*H_b*e_v' ...
417 -1/2*a*Hi*e_u*H_b*d_v' ...
418 +a*Hi*e_u*tau*H_b*e_v';
387 end 419 end
388 420
389 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type) 421 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
390 422
391 % TODO: Make this work for curvilinear grids 423 % TODO: Make this work for curvilinear grids
392 warning('LaplaceCurvilinear: Non-conforming grid interpolation only works for Cartesian grids.'); 424 warning('LaplaceCurvilinear: Non-conforming grid interpolation only works for Cartesian grids.');
425 warning('LaplaceCurvilinear: Non-conforming interface uses Virtas penalty strength');
393 426
394 % User can request special interpolation operators by specifying type.interpOpSet 427 % User can request special interpolation operators by specifying type.interpOpSet
395 default_field(type, 'interpOpSet', @sbp.InterpOpsOP); 428 default_field(type, 'interpOpSet', @sbp.InterpOpsOP);
396 interpOpSet = type.interpOpSet; 429 interpOpSet = type.interpOpSet;
397 tuning = type.tuning; 430 tuning = type.tuning;