Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinear.m @ 1210:5e7692ed7c7c feature/laplace_curvilinear_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add CG interface coupling
author Martin Almquist <malmquist@stanford.edu>
date Sun, 22 Sep 2019 19:05:17 -0700
parents 7259b3f3e227
children
comparison
equal deleted inserted replaced
1141:91058813b6e7 1210:5e7692ed7c7c
7 grid 7 grid
8 8
9 order % Order of accuracy for the approximation 9 order % Order of accuracy for the approximation
10 10
11 a,b % Parameters of the operator 11 a,b % Parameters of the operator
12 weight % Parameter in front of time derivative (e.g. u_tt in wave equation) here: 1/a.
12 13
13 14
14 % Inner products and operators for physical coordinates 15 % Inner products and operators for physical coordinates
15 D % Laplace operator 16 D % Laplace operator
16 H, Hi % Inner product 17 H, Hi % Inner product
59 default_arg('b', @(x,y) 0*x + 1); 60 default_arg('b', @(x,y) 0*x + 1);
60 61
61 % assert(isa(g, 'grid.Curvilinear')) 62 % assert(isa(g, 'grid.Curvilinear'))
62 if isa(a, 'function_handle') 63 if isa(a, 'function_handle')
63 a = grid.evalOn(g, a); 64 a = grid.evalOn(g, a);
65 end
66
67 % If a is scalar
68 if length(a) == 1
69 a = a*ones(g.N(), 1);
64 end 70 end
65 a = spdiag(a); 71 a = spdiag(a);
66 72
67 if isa(b, 'function_handle') 73 if isa(b, 'function_handle')
68 b = grid.evalOn(g, b); 74 b = grid.evalOn(g, b);
237 obj.order = order; 243 obj.order = order;
238 obj.grid = g; 244 obj.grid = g;
239 obj.dim = dim; 245 obj.dim = dim;
240 246
241 obj.a = a; 247 obj.a = a;
248 obj.weight = inv(a);
242 obj.b = b; 249 obj.b = b;
243 obj.a11 = a11; 250 obj.a11 = a11;
244 obj.a12 = a12; 251 obj.a12 = a12;
245 obj.a22 = a22; 252 obj.a22 = a22;
246 obj.s_w = spdiag(s_w); 253 obj.s_w = spdiag(s_w);
325 % Fields: 332 % Fields:
326 % -- tuning: penalty strength, defaults to 1.2 333 % -- tuning: penalty strength, defaults to 1.2
327 % -- interpolation: type of interpolation, default 'none' 334 % -- interpolation: type of interpolation, default 'none'
328 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) 335 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
329 336
330 % error('Not implemented') 337 defaultType.coupling = 'sat';
331
332 defaultType.tuning = 1.0; 338 defaultType.tuning = 1.0;
333 defaultType.interpolation = 'none'; 339 defaultType.interpolation = 'none';
334 default_struct('type', defaultType); 340 default_struct('type', defaultType);
335 341
336 switch type.interpolation 342 switch type.coupling
337 case {'none', ''} 343 case {'cg', 'CG'}
338 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type); 344 [closure, penalty] = interfaceCG(obj,boundary,neighbour_scheme,neighbour_boundary,type);
339 case {'op','OP'} 345 case {'sat', 'SAT'}
340 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type); 346 switch type.interpolation
347 case {'none', ''}
348 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type);
349 case {'op','OP'}
350 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type);
351 otherwise
352 error('Unknown type of interpolation: %s ', type.interpolation);
353 end
341 otherwise 354 otherwise
342 error('Unknown type of interpolation: %s ', type.interpolation); 355 error('Unknown type of coupling: %s ', type.coupling);
343 end 356 end
344 end 357 end
345 358
346 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type) 359 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
347 tuning = type.tuning; 360 tuning = type.tuning;
412 -a*Hi*e_u*tau*H_b*e_u'; 425 -a*Hi*e_u*tau*H_b*e_u';
413 426
414 penalty = -1/2*a*Hi*d_u*b_b_u*H_b*e_v' ... 427 penalty = -1/2*a*Hi*d_u*b_b_u*H_b*e_v' ...
415 -1/2*a*Hi*e_u*H_b*b_b_v*d_v' ... 428 -1/2*a*Hi*e_u*H_b*b_b_v*d_v' ...
416 +a*Hi*e_u*tau*H_b*e_v'; 429 +a*Hi*e_u*tau*H_b*e_v';
430 end
431
432 function [closure, penalty] = interfaceCG(obj,boundary,neighbour_scheme,neighbour_boundary,type)
433
434 % There is no penalty, only a closure. And the closure is the same as for Neumann BC
435 e = obj.getBoundaryOperator('e', boundary);
436 d = obj.getBoundaryOperator('d', boundary);
437 H_b = obj.getBoundaryQuadrature(boundary);
438
439 e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
440
441 Hi = obj.Hi;
442 a = obj.a;
443 b_b = e'*obj.b*e;
444
445 tau1 = -1;
446 tau2 = 0;
447 tau = (tau1*e + tau2*d)*H_b;
448
449 closure = a*Hi*tau*b_b*d';
450
451 % Zero penalty of correct dimensions
452 penalty = 0*e*e_v';
417 end 453 end
418 454
419 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type) 455 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
420 456
421 % TODO: Make this work for curvilinear grids 457 % TODO: Make this work for curvilinear grids