Mercurial > repos > public > sbplib
changeset 1068:e0ecce90f8cf feature/laplace_curvilinear_test
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 | 7a55a72729e6 |
| files | +scheme/LaplaceCurvilinearNewCorner.m |
| diffstat | 1 files changed, 33 insertions(+), 21 deletions(-) [+] |
line wrap: on
line diff
diff -r 9a858436f8fa -r e0ecce90f8cf +scheme/LaplaceCurvilinearNewCorner.m --- a/+scheme/LaplaceCurvilinearNewCorner.m Tue Jan 22 18:17:01 2019 -0800 +++ b/+scheme/LaplaceCurvilinearNewCorner.m Tue Jan 22 19:18:24 2019 -0800 @@ -58,11 +58,7 @@ function obj = LaplaceCurvilinearNewCorner(g, order, a, b, opSet) default_arg('opSet',@sbp.D2Variable); default_arg('a', 1); - default_arg('b', 1); - - if b ~=1 - error('Not implemented yet') - end + default_arg('b', @(x,y) 0*x + 1); % assert(isa(g, 'grid.Curvilinear')) if isa(a, 'function_handle') @@ -70,6 +66,16 @@ a = spdiag(a); end + if isa(b, 'function_handle') + b = grid.evalOn(g, b); + b = spdiag(b); + end + + % If b is scalar + if length(b) == 1 + b = b*speye(g.N(), g.N()); + end + dim = 2; m = g.size(); m_u = m(1); @@ -159,21 +165,23 @@ % Assemble full operators L_12 = spdiag(a12); - Duv = Du*L_12*Dv; - Dvu = Dv*L_12*Du; + Duv = Du*b*L_12*Dv; + Dvu = Dv*b*L_12*Du; Duu = sparse(m_tot); Dvv = sparse(m_tot); ind = grid.funcToMatrix(g, 1:m_tot); for i = 1:m_v - D = D2_u(a11(ind(:,i))); + b_a11 = b*a11; + D = D2_u(b_a11(ind(:,i))); p = ind(:,i); Duu(p,p) = D; end for i = 1:m_u - D = D2_v(a22(ind(i,:))); + b_a22 = b*a22; + D = D2_v(b_a22(ind(i,:))); p = ind(i,:); Dvv(p,p) = D; end @@ -281,6 +289,7 @@ J = obj.J; Hi = obj.Hi; a = obj.a; + b_b = e'*obj.b*e; switch type % Dirichlet boundary condition @@ -302,21 +311,21 @@ tau = tuning*(tau_R + tau_H); - closure = a*Hi*d*H_b*e' ... - -a*Hi*e*tau*H_b*e'; + closure = a*Hi*d*b_b*H_b*e' ... + -a*Hi*e*tau*b_b*H_b*e'; - penalty = -a*Hi*d*H_b ... - +a*Hi*e*tau*H_b; + penalty = -a*Hi*d*b_b*H_b ... + +a*Hi*e*tau*b_b*H_b; - % Neumann boundary condition + % Neumann boundary condition. Note that the penalty is for du/dn and not b*du/dn. case {'N','n','neumann'} tau1 = -1; tau2 = 0; tau = (tau1*e + tau2*d)*H_b; - closure = a*Hi*tau*d'; - penalty = -a*Hi*tau; + closure = a*Hi*tau*b_b*d'; + penalty = -a*Hi*tau*b_b; % Unknown, boundary condition @@ -367,6 +376,7 @@ % Coefficients, u K_u = u.K; J_u = u.J; + b_b_u = e_u'*u.b*e_u; % Boundary operators, v e_v = v.getBoundaryOperator('e', neighbour_boundary); @@ -379,6 +389,7 @@ % Coefficients, v K_v = v.K; J_v = v.J; + b_b_v = e_v'*v.b*e_v; %--- Penalty strength tau ------------- sigma_u = 0; @@ -401,7 +412,7 @@ tau_H_v(1,1) = dim*tau_H_v(1,1); tau_H_v(end,end) = dim*tau_H_v(end,end); - tau = 1/4*tuning*(tau_R_u + tau_H_u + tau_R_v + tau_H_v); + tau = 1/4*tuning*(b_b_u*(tau_R_u + tau_H_u) + b_b_v*(tau_R_v + tau_H_v)); %-------------------------------------- % Operators/coefficients that are only required from this side @@ -409,12 +420,12 @@ H_b = u.getBoundaryQuadrature(boundary); a = u.a; - closure = 1/2*a*Hi*d_u*H_b*e_u' ... - -1/2*a*Hi*e_u*H_b*d_u' ... + closure = 1/2*a*Hi*d_u*b_b_u*H_b*e_u' ... + -1/2*a*Hi*e_u*H_b*b_b_u*d_u' ... -a*Hi*e_u*tau*H_b*e_u'; - penalty = -1/2*a*Hi*d_u*H_b*e_v' ... - -1/2*a*Hi*e_u*H_b*d_v' ... + penalty = -1/2*a*Hi*d_u*b_b_u*H_b*e_v' ... + -1/2*a*Hi*e_u*H_b*b_b_v*d_v' ... +a*Hi*e_u*tau*H_b*e_v'; end @@ -423,6 +434,7 @@ % TODO: Make this work for curvilinear grids warning('LaplaceCurvilinear: Non-conforming grid interpolation only works for Cartesian grids.'); warning('LaplaceCurvilinear: Non-conforming interface uses Virtas penalty strength'); + warning('LaplaceCurvilinear: Non-conforming interface assumes that b is constant'); % User can request special interpolation operators by specifying type.interpOpSet default_field(type, 'interpOpSet', @sbp.InterpOpsOP);