Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinearNewCorner.m @ 1068:e0ecce90f8cf feature/laplace_curvilinear_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add variable b. Tested for interface and Dirichlet, but not Neumann yet.
author Martin Almquist <malmquist@stanford.edu>
date Tue, 22 Jan 2019 19:18:24 -0800
parents 9a858436f8fa
children d5290a056049
comparison
equal deleted inserted replaced
1067:9a858436f8fa 1068:e0ecce90f8cf
56 % TODO: Implement proper H, it should be the real physical quadrature, the logic quadrature may be but in a separate variable (H_logic?) 56 % TODO: Implement proper H, it should be the real physical quadrature, the logic quadrature may be but in a separate variable (H_logic?)
57 57
58 function obj = LaplaceCurvilinearNewCorner(g, order, a, b, opSet) 58 function obj = LaplaceCurvilinearNewCorner(g, order, a, b, opSet)
59 default_arg('opSet',@sbp.D2Variable); 59 default_arg('opSet',@sbp.D2Variable);
60 default_arg('a', 1); 60 default_arg('a', 1);
61 default_arg('b', 1); 61 default_arg('b', @(x,y) 0*x + 1);
62
63 if b ~=1
64 error('Not implemented yet')
65 end
66 62
67 % assert(isa(g, 'grid.Curvilinear')) 63 % assert(isa(g, 'grid.Curvilinear'))
68 if isa(a, 'function_handle') 64 if isa(a, 'function_handle')
69 a = grid.evalOn(g, a); 65 a = grid.evalOn(g, a);
70 a = spdiag(a); 66 a = spdiag(a);
67 end
68
69 if isa(b, 'function_handle')
70 b = grid.evalOn(g, b);
71 b = spdiag(b);
72 end
73
74 % If b is scalar
75 if length(b) == 1
76 b = b*speye(g.N(), g.N());
71 end 77 end
72 78
73 dim = 2; 79 dim = 2;
74 m = g.size(); 80 m = g.size();
75 m_u = m(1); 81 m_u = m(1);
157 obj.y_u = y_u; 163 obj.y_u = y_u;
158 obj.y_v = y_v; 164 obj.y_v = y_v;
159 165
160 % Assemble full operators 166 % Assemble full operators
161 L_12 = spdiag(a12); 167 L_12 = spdiag(a12);
162 Duv = Du*L_12*Dv; 168 Duv = Du*b*L_12*Dv;
163 Dvu = Dv*L_12*Du; 169 Dvu = Dv*b*L_12*Du;
164 170
165 Duu = sparse(m_tot); 171 Duu = sparse(m_tot);
166 Dvv = sparse(m_tot); 172 Dvv = sparse(m_tot);
167 ind = grid.funcToMatrix(g, 1:m_tot); 173 ind = grid.funcToMatrix(g, 1:m_tot);
168 174
169 for i = 1:m_v 175 for i = 1:m_v
170 D = D2_u(a11(ind(:,i))); 176 b_a11 = b*a11;
177 D = D2_u(b_a11(ind(:,i)));
171 p = ind(:,i); 178 p = ind(:,i);
172 Duu(p,p) = D; 179 Duu(p,p) = D;
173 end 180 end
174 181
175 for i = 1:m_u 182 for i = 1:m_u
176 D = D2_v(a22(ind(i,:))); 183 b_a22 = b*a22;
184 D = D2_v(b_a22(ind(i,:)));
177 p = ind(i,:); 185 p = ind(i,:);
178 Dvv(p,p) = D; 186 Dvv(p,p) = D;
179 end 187 end
180 188
181 189
279 287
280 K = obj.K; 288 K = obj.K;
281 J = obj.J; 289 J = obj.J;
282 Hi = obj.Hi; 290 Hi = obj.Hi;
283 a = obj.a; 291 a = obj.a;
292 b_b = e'*obj.b*e;
284 293
285 switch type 294 switch type
286 % Dirichlet boundary condition 295 % Dirichlet boundary condition
287 case {'D','d','dirichlet'} 296 case {'D','d','dirichlet'}
288 tuning = 1.0; 297 tuning = 1.0;
300 tau_H(1,1) = obj.dim*tau_H(1,1); 309 tau_H(1,1) = obj.dim*tau_H(1,1);
301 tau_H(end,end) = obj.dim*tau_H(end,end); 310 tau_H(end,end) = obj.dim*tau_H(end,end);
302 311
303 tau = tuning*(tau_R + tau_H); 312 tau = tuning*(tau_R + tau_H);
304 313
305 closure = a*Hi*d*H_b*e' ... 314 closure = a*Hi*d*b_b*H_b*e' ...
306 -a*Hi*e*tau*H_b*e'; 315 -a*Hi*e*tau*b_b*H_b*e';
307 316
308 penalty = -a*Hi*d*H_b ... 317 penalty = -a*Hi*d*b_b*H_b ...
309 +a*Hi*e*tau*H_b; 318 +a*Hi*e*tau*b_b*H_b;
310 319
311 320
312 % Neumann boundary condition 321 % Neumann boundary condition. Note that the penalty is for du/dn and not b*du/dn.
313 case {'N','n','neumann'} 322 case {'N','n','neumann'}
314 tau1 = -1; 323 tau1 = -1;
315 tau2 = 0; 324 tau2 = 0;
316 tau = (tau1*e + tau2*d)*H_b; 325 tau = (tau1*e + tau2*d)*H_b;
317 326
318 closure = a*Hi*tau*d'; 327 closure = a*Hi*tau*b_b*d';
319 penalty = -a*Hi*tau; 328 penalty = -a*Hi*tau*b_b;
320 329
321 330
322 % Unknown, boundary condition 331 % Unknown, boundary condition
323 otherwise 332 otherwise
324 error('No such boundary condition: type = %s',type); 333 error('No such boundary condition: type = %s',type);
365 m_u = u.getBoundaryNumber(boundary); 374 m_u = u.getBoundaryNumber(boundary);
366 375
367 % Coefficients, u 376 % Coefficients, u
368 K_u = u.K; 377 K_u = u.K;
369 J_u = u.J; 378 J_u = u.J;
379 b_b_u = e_u'*u.b*e_u;
370 380
371 % Boundary operators, v 381 % Boundary operators, v
372 e_v = v.getBoundaryOperator('e', neighbour_boundary); 382 e_v = v.getBoundaryOperator('e', neighbour_boundary);
373 d_v = v.getBoundaryOperator('d', neighbour_boundary); 383 d_v = v.getBoundaryOperator('d', neighbour_boundary);
374 gamm_v = v.getBoundaryBorrowing(neighbour_boundary); 384 gamm_v = v.getBoundaryBorrowing(neighbour_boundary);
377 m_v = v.getBoundaryNumber(neighbour_boundary); 387 m_v = v.getBoundaryNumber(neighbour_boundary);
378 388
379 % Coefficients, v 389 % Coefficients, v
380 K_v = v.K; 390 K_v = v.K;
381 J_v = v.J; 391 J_v = v.J;
392 b_b_v = e_v'*v.b*e_v;
382 393
383 %--- Penalty strength tau ------------- 394 %--- Penalty strength tau -------------
384 sigma_u = 0; 395 sigma_u = 0;
385 sigma_v = 0; 396 sigma_v = 0;
386 for i = 1:obj.dim 397 for i = 1:obj.dim
399 410
400 tau_H_v = 1/th_H_v*sigma_v; 411 tau_H_v = 1/th_H_v*sigma_v;
401 tau_H_v(1,1) = dim*tau_H_v(1,1); 412 tau_H_v(1,1) = dim*tau_H_v(1,1);
402 tau_H_v(end,end) = dim*tau_H_v(end,end); 413 tau_H_v(end,end) = dim*tau_H_v(end,end);
403 414
404 tau = 1/4*tuning*(tau_R_u + tau_H_u + tau_R_v + tau_H_v); 415 tau = 1/4*tuning*(b_b_u*(tau_R_u + tau_H_u) + b_b_v*(tau_R_v + tau_H_v));
405 %-------------------------------------- 416 %--------------------------------------
406 417
407 % Operators/coefficients that are only required from this side 418 % Operators/coefficients that are only required from this side
408 Hi = u.Hi; 419 Hi = u.Hi;
409 H_b = u.getBoundaryQuadrature(boundary); 420 H_b = u.getBoundaryQuadrature(boundary);
410 a = u.a; 421 a = u.a;
411 422
412 closure = 1/2*a*Hi*d_u*H_b*e_u' ... 423 closure = 1/2*a*Hi*d_u*b_b_u*H_b*e_u' ...
413 -1/2*a*Hi*e_u*H_b*d_u' ... 424 -1/2*a*Hi*e_u*H_b*b_b_u*d_u' ...
414 -a*Hi*e_u*tau*H_b*e_u'; 425 -a*Hi*e_u*tau*H_b*e_u';
415 426
416 penalty = -1/2*a*Hi*d_u*H_b*e_v' ... 427 penalty = -1/2*a*Hi*d_u*b_b_u*H_b*e_v' ...
417 -1/2*a*Hi*e_u*H_b*d_v' ... 428 -1/2*a*Hi*e_u*H_b*b_b_v*d_v' ...
418 +a*Hi*e_u*tau*H_b*e_v'; 429 +a*Hi*e_u*tau*H_b*e_v';
419 end 430 end
420 431
421 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type) 432 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
422 433
423 % TODO: Make this work for curvilinear grids 434 % TODO: Make this work for curvilinear grids
424 warning('LaplaceCurvilinear: Non-conforming grid interpolation only works for Cartesian grids.'); 435 warning('LaplaceCurvilinear: Non-conforming grid interpolation only works for Cartesian grids.');
425 warning('LaplaceCurvilinear: Non-conforming interface uses Virtas penalty strength'); 436 warning('LaplaceCurvilinear: Non-conforming interface uses Virtas penalty strength');
437 warning('LaplaceCurvilinear: Non-conforming interface assumes that b is constant');
426 438
427 % User can request special interpolation operators by specifying type.interpOpSet 439 % User can request special interpolation operators by specifying type.interpOpSet
428 default_field(type, 'interpOpSet', @sbp.InterpOpsOP); 440 default_field(type, 'interpOpSet', @sbp.InterpOpsOP);
429 interpOpSet = type.interpOpSet; 441 interpOpSet = type.interpOpSet;
430 tuning = type.tuning; 442 tuning = type.tuning;