Mercurial > repos > public > sbplib
changeset 491:26125cfefe11 feature/quantumTriangles
Schrodinger 2d runs without exploding with fixed coordinates
| author | Ylva Rydin <ylva.rydin@telia.com> |
|---|---|
| date | Fri, 10 Feb 2017 14:29:53 +0100 |
| parents | b13d44271ead |
| children | 6b95a894cbd7 |
| files | +grid/Schrodinger2dCurve.m +scheme/Schrodinger2dCurve.m +scheme/Schrodinger2dCurve.m~ |
| diffstat | 3 files changed, 552 insertions(+), 290 deletions(-) [+] |
line wrap: on
line diff
--- a/+grid/Schrodinger2dCurve.m Thu Feb 09 11:41:21 2017 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,290 +0,0 @@ -classdef Schrodinger2dCurve < scheme.Scheme - properties - m % Number of points in each direction, possibly a vector - h % Grid spacing - - grid - - order % Order accuracy for the approximation - - D % non-stabalized scheme operator - M % Derivative norm - H % Discrete norm - Hi - H_u, H_v % Norms in the x and y directions - Hu,Hv % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. - Hi_u, Hi_v - Hiu, Hiv - D1_v D1_u - D2_v D2_u - Du Dv - - - e_w, e_e, e_s, e_n - du_w, dv_w - du_e, dv_e - du_s, dv_s - du_n, dv_n - g_1 - g_2 - - p,p_tau - end - - methods - function obj = Schrodinger2dCurve(g ,order, opSet, p,p_tau) - default_arg('opSet',@sbp.D2Variable); - default_arg('c', 1); - - assert(isa(g, 'grid.Curvilinear')) - - obj.p=p; - obj.p_tau=p_tau; - obj.c=1; - - m = g.size(); - m_u = m(1); - m_v = m(2); - m_tot = g.N(); - - h = g.scaling(); - h_u = h(1); - h_v = h(2); - - % Operators - ops_u = opSet(m_u, {0, 1}, order); - ops_v = opSet(m_v, {0, 1}, order); - - I_u = speye(m_u); - I_v = speye(m_v); - - D1_u = ops_u.D1; - - H_u = ops_u.H; - Hi_u = ops_u.HI; - e_l_u = ops_u.e_l; - e_r_u = ops_u.e_r; - d1_l_u = ops_u.d1_l; - d1_r_u = ops_u.d1_r; - - obj.D1_v = ops_v.D1; - obj.D2_v = ops_v.D2; - H_v = ops_v.H; - Hi_v = ops_v.HI; - e_l_v = ops_v.e_l; - e_r_v = ops_v.e_r; - d1_l_v = ops_v.d1_l; - d1_r_v = ops_v.d1_r; - - obj.Du = kr(D1_u,I_v); - obj.Dv = kr(I_u,D1_v); - - obj.H = kr(H_u,H_v); - obj.Hi = kr(Hi_u,Hi_v); - obj.Hu = kr(H_u,I_v); - obj.Hv = kr(I_u,H_v); - obj.Hiu = kr(Hi_u,I_v); - obj.Hiv = kr(I_u,Hi_v); - - obj.e_w = kr(e_l_u,I_v); - obj.e_e = kr(e_r_u,I_v); - obj.e_s = kr(I_u,e_l_v); - obj.e_n = kr(I_u,e_r_v); - obj.du_w = kr(d1_l_u,I_v); - obj.dv_w = (obj.e_w'*Dv)'; - obj.du_e = kr(d1_r_u,I_v); - obj.dv_e = (obj.e_e'*Dv)'; - obj.du_s = (obj.e_s'*Du)'; - obj.dv_s = kr(I_u,d1_l_v); - obj.du_n = (obj.e_n'*Du)'; - obj.dv_n = kr(I_u,d1_r_v); - -% obj.x_u = x_u; -% obj.x_v = x_v; -% obj.y_u = y_u; -% obj.y_v = y_v; - - obj.m = m; - obj.h = [h_u h_v]; - obj.order = order; - obj.grid = g; - - - end - - - function [D ]= d_fun(obj,t) - % Metric derivatives - ti = parametrization.Ti.points(obj.p.p1(t),obj.p.p2(t),obj.p.p3,obj.p.p4); - ti_tau = parametrization.Ti.points(obj.p_tau.p1(t),obj.p_tau.p2(t),obj.p_tau.p3,obj.p_tau.p4); - - coords = parametrization.ti.map(); - coords_tau = parametrization.ti_tau.map(); - x = coords(:,1); - y = coords(:,2); - - x_tau = coords_tau(:,1); - y_tau = coords_tau(:,2); - - x_u = obj.Du*x; - x_v = obj.Dv*x; - y_u = obj.Du*y; - y_v = obj.Dv*y; - - J = x_u.*y_v - x_v.*y_u; - a11 = 1./J.* (x_v.^2 + y_v.^2); - a12 = -1./J .* (x_u.*x_v + y_u.*y_v); - a22 = 1./J .* (x_u.^2 + y_u.^2); - lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2)); - - % Assemble full operators - L_12 = spdiags(a12, 0, m_tot, m_tot); - Duv = obj.Du*L_12*obj.Dv; - Dvu = obj.Dv*L_12*obj.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))); - p = ind(:,i); - Duu(p,p) = D; - end - - for i = 1:m_u - D = D2_v(a22(ind(i,:))); - p = ind(i,:); - Dvv(p,p) = D; - end - - - J = spdiags(J, 0, m_tot, m_tot); - Ji = spdiags(1./J, 0, m_tot, m_tot); - obj.g_1 = x_tau.*y_v-y_tau.*x_v; - obj.g_2 = -x_tau.*y_u + y_tau.*x_u; - - %Add the flux splitting - D = Ji*(obj.g_1*obj.Du + obj.g_2*obj.Dv + 1i*obj.c^2*(Duu + Duv + Dvu + Dvv)); - -% obj.gamm_u = h_u*ops_u.borrowing.M.d1; -% obj.gamm_v = h_v*ops_v.borrowing.M.d1; - - end - - % Closure functions return the opertors applied to the own doamin to close the boundary - % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. - % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. - % type is a string specifying the type of boundary condition if there are several. - % data is a function returning the data that should be applied at the boundary. - % neighbour_scheme is an instance of Scheme that should be interfaced to. - % neighbour_boundary is a string specifying which boundary to interface to. - function [closure, penalty] = boundary_condition(obj, boundary, parameter) - default_arg('type','neumann'); - default_arg('parameter', []); - - % v denotes the solution in the neighbour domain - tuning = 1.2; - % tuning = 20.2; - [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g] = obj.get_boundary_ops(boundary); - - a_n = spdiag(coeff_n); - a_t = spdiag(coeff_t); - - F = (s * a_n * d_n' + s * a_t*d_t')'; - - u = obj; - - b1 = gamm*u.lambda./u.a11.^2; - b2 = gamm*u.lambda./u.a22.^2; - - tau1 = -1./b1 - 1./b2; - tau1 = tuning * spdiag(tau1); - sig1 = 1; - - a = e'*g; - tau2 = (-1*s*a - abs(a))/4; - - penalty_parameter_1 = halfnorm_inv_n*(tau1 + sig1*halfnorm_inv_t*F*e'*halfnorm_t)*e; - penalty_parameter_2 = halfnorm_inv_n(tau2)*e; - - closure = obj.Ji*obj.c^2 * penalty_parameter_1*e' +obj.Ji* penalty_parameter_2*e'; - penalty = -obj.Ji*obj.c^2 * penalty_parameter_1; - - end - - function [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t, I, scale_factor,g] = get_boundary_ops(obj, boundary) - - % gridMatrix = zeros(obj.m(2),obj.m(1)); - % gridMatrix(:) = 1:numel(gridMatrix); - - ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m)); - - switch boundary - case 'w' - e = obj.e_w; - d_n = obj.du_w; - d_t = obj.dv_w; - s = -1; - - I = ind(1,:); - coeff_n = obj.a11(I); - coeff_t = obj.a12(I); - scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2); - case 'e' - e = obj.e_e; - d_n = obj.du_e; - d_t = obj.dv_e; - s = 1; - - I = ind(end,:); - coeff_n = obj.a11(I); - coeff_t = obj.a12(I); - scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2); - case 's' - e = obj.e_s; - d_n = obj.dv_s; - d_t = obj.du_s; - s = -1; - - I = ind(:,1)'; - coeff_n = obj.a22(I); - coeff_t = obj.a12(I); - scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2); - case 'n' - e = obj.e_n; - d_n = obj.dv_n; - d_t = obj.du_n; - s = 1; - - I = ind(:,end)'; - coeff_n = obj.a22(I); - coeff_t = obj.a12(I); - scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2); - otherwise - error('No such boundary: boundary = %s',boundary); - end - - switch boundary - case {'w','e'} - halfnorm_inv_n = obj.Hiu; - halfnorm_inv_t = obj.Hiv; - halfnorm_t = obj.Hv; - gamm = obj.gamm_u; - g=obj.g_1; - case {'s','n'} - halfnorm_inv_n = obj.Hiv; - halfnorm_inv_t = obj.Hiu; - halfnorm_t = obj.Hu; - gamm = obj.gamm_v; - g=obj.g_2; - end - end - - function N = size(obj) - N = prod(obj.m); - end - - - end -end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Schrodinger2dCurve.m Fri Feb 10 14:29:53 2017 +0100 @@ -0,0 +1,276 @@ +classdef Schrodinger2dCurve < scheme.Scheme + properties + m % Number of points in each direction, possibly a vector + h % Grid spacing + + grid + xm, ym + + order % Order accuracy for the approximation + + D % non-stabalized scheme operator + M % Derivative norm + H % Discrete norm + Hi + H_u, H_v % Norms in the x and y directions + Hu,Hv % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. + Hi_u, Hi_v + Hiu, Hiv + D1_v, D1_u + D2_v, D2_u + Du, Dv + + + e_w, e_e, e_s, e_n + du_w, dv_w + du_e, dv_e + du_s, dv_s + du_n, dv_n + g_1, g_2 + c + a11, a12, a22 + m_tot, m_u, m_v + p,p_tau + Ji + end + + methods + function obj = Schrodinger2dCurve(g ,order, opSet,p,p_tau) + default_arg('opSet',@sbp.D2Variable); + default_arg('c', 1); + + obj.p=p; + obj.p_tau=p_tau; + obj.c=1; + + m = g.size(); + obj.m_u = m(1); + obj.m_v = m(2); + obj.m_tot = g.N(); + + h = g.scaling(); + h_u = h(1); + h_v = h(2); + + % Operators + ops_u = opSet(obj.m_u, {0, 1}, order); + ops_v = opSet(obj.m_v, {0, 1}, order); + + I_u = speye(obj.m_u); + I_v = speye(obj.m_v); + + obj.D1_u = ops_u.D1; + obj.D2_u = ops_u.D2; + + H_u = ops_u.H; + Hi_u = ops_u.HI; + e_l_u = ops_u.e_l; + e_r_u = ops_u.e_r; + d1_l_u = ops_u.d1_l; + d1_r_u = ops_u.d1_r; + + obj.D1_v = ops_v.D1; + obj.D2_v = ops_v.D2; + H_v = ops_v.H; + Hi_v = ops_v.HI; + e_l_v = ops_v.e_l; + e_r_v = ops_v.e_r; + d1_l_v = ops_v.d1_l; + d1_r_v = ops_v.d1_r; + + obj.Du = kr(obj.D1_u,I_v); + obj.Dv = kr(I_u,obj.D1_v); + + obj.H = kr(H_u,H_v); + obj.Hi = kr(Hi_u,Hi_v); + obj.Hu = kr(H_u,I_v); + obj.Hv = kr(I_u,H_v); + obj.Hiu = kr(Hi_u,I_v); + obj.Hiv = kr(I_u,Hi_v); + + obj.e_w = kr(e_l_u,I_v); + obj.e_e = kr(e_r_u,I_v); + obj.e_s = kr(I_u,e_l_v); + obj.e_n = kr(I_u,e_r_v); + obj.du_w = kr(d1_l_u,I_v); + obj.dv_w = (obj.e_w'*obj.Dv)'; + obj.du_e = kr(d1_r_u,I_v); + obj.dv_e = (obj.e_e'*obj.Dv)'; + obj.du_s = (obj.e_s'*obj.Du)'; + obj.dv_s = kr(I_u,d1_l_v); + obj.du_n = (obj.e_n'*obj.Du)'; + obj.dv_n = kr(I_u,d1_r_v); + +% obj.x_u = x_u; +% obj.x_v = x_v; +% obj.y_u = y_u; +% obj.y_v = y_v; + + obj.m = m; + obj.h = [h_u h_v]; + obj.order = order; + obj.grid = g; + + + end + + + function [D ]= d_fun(obj,t) + % Metric derivatives + ti = parametrization.Ti.points(obj.p.p1(t),obj.p.p2(t),obj.p.p3(t),obj.p.p4(t)); + ti_tau = parametrization.Ti.points(obj.p_tau.p1(t),obj.p_tau.p2(t),obj.p_tau.p3(t),obj.p_tau.p4(t)); + + lcoords=points(obj.grid); + [obj.xm,obj.ym]= ti.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2)); + [x_tau,y_tau]= ti_tau.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2)); + x = reshape(obj.xm,obj.m_tot,1); + y = reshape(obj.ym,obj.m_tot,1); + + x_tau = reshape(x_tau,obj.m_tot,1); + y_tau = reshape(y_tau,obj.m_tot,1); + + x_u = obj.Du*x; + x_v = obj.Dv*x; + y_u = obj.Du*y; + y_v = obj.Dv*y; + + J = x_u.*y_v - x_v.*y_u; + a11 = 1./J.* (x_v.^2 + y_v.^2); + a12 = -1./J .* (x_u.*x_v + y_u.*y_v); + a22 = 1./J .* (x_u.^2 + y_u.^2); + + obj.a11=a11; + obj.a12=a12; + obj.a22=a22; + + % Assemble full operators + L_12 = spdiags(a12, 0, obj.m_tot, obj.m_tot); + Duv = obj.Du*L_12*obj.Dv; + Dvu = obj.Dv*L_12*obj.Du; + + Duu = sparse(obj.m_tot); + Dvv = sparse(obj.m_tot); + ind = grid.funcToMatrix(obj.grid, 1:obj.m_tot); + + + for i = 1:obj.m_v + D = obj.D2_u(a11(ind(:,i))); + p = ind(:,i); + Duu(p,p) = D; + end + + for i = 1:obj.m_u + D = obj.D2_v(a22(ind(i,:))); + p = ind(i,:); + Dvv(p,p) = D; + end + + Ji = spdiags(1./J, 0, obj.m_tot, obj.m_tot); + obj.Ji=Ji; + obj.g_1 = x_tau.*y_v-y_tau.*x_v; + obj.g_2 = -x_tau.*y_u + y_tau.*x_u; + + %Add the flux splitting + D = Ji*(obj.g_1.*obj.Du + obj.g_2.*obj.Dv + 1i*obj.c^2*(Duu + Duv + Dvu + Dvv)); %% g_1' och g_2'? + +% obj.gamm_u = h_u*ops_u.borrowing.M.d1; +% obj.gamm_v = h_v*ops_v.borrowing.M.d1; + + end + + % Closure functions return the opertors applied to the own doamin to close the boundary + % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. + % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. + % type is a string specifying the type of boundary condition if there are several. + % data is a function returning the data that should be applied at the boundary. + % neighbour_scheme is an instance of Scheme that should be interfaced to. + % neighbour_boundary is a string specifying which boundary to interface to. + function [closure, penalty] = boundary_condition(obj, boundary) + [e, d_n, d_t, coeff_t, s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g] = obj.get_boundary_ops(boundary); + + a_t = spdiag(coeff_t); + F = (s * d_n' + s * a_t*d_t')'; + tau1 = 1; + a = spdiag(g); + tau2 = (-1*s*a - abs(a))/4; + + penalty_parameter_1 = 1*1i*halfnorm_inv_n*halfnorm_inv_t*e*F'*halfnorm_t*e; + penalty_parameter_2 = halfnorm_inv_n*e*tau2; + + closure = obj.Ji*obj.c^2 * penalty_parameter_1*e' +obj.Ji* penalty_parameter_2*e'; + % penalty = -obj.Ji*obj.c^2 * penalty_parameter_2; + penalty=0; + + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + end + + function [e, d_n, d_t, coeff_t, s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g, I] = get_boundary_ops(obj, boundary) + + % gridMatrix = zeros(obj.m(2),obj.m(1)); + % gridMatrix(:) = 1:numel(gridMatrix); + + ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m)); + + switch boundary + case 'w' + e = obj.e_w; + d_n = obj.du_w; + d_t = obj.dv_w; + s = -1; + + I = ind(1,:); + coeff_t = obj.a12(I); + g=obj.g_1(I); + case 'e' + e = obj.e_e; + d_n = obj.du_e; + d_t = obj.dv_e; + s = 1; + + I = ind(end,:); + coeff_t = obj.a12(I); + g=obj.g_1(I); + case 's' + e = obj.e_s; + d_n = obj.dv_s; + d_t = obj.du_s; + s = -1; + + I = ind(:,1)'; + coeff_t = obj.a12(I); + g=obj.g_2(I); + case 'n' + e = obj.e_n; + d_n = obj.dv_n; + d_t = obj.du_n; + s = 1; + + I = ind(:,end)'; + coeff_t = obj.a12(I); + g=obj.g_2(I); + otherwise + error('No such boundary: boundary = %s',boundary); + end + + switch boundary + case {'w','e'} + halfnorm_inv_n = obj.Hiu; + halfnorm_inv_t = obj.Hiv; + halfnorm_t = obj.Hv; + + case {'s','n'} + halfnorm_inv_n = obj.Hiv; + halfnorm_inv_t = obj.Hiu; + halfnorm_t = obj.Hu; + end + end + + function N = size(obj) + N = prod(obj.m); + end + + + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Schrodinger2dCurve.m~ Fri Feb 10 14:29:53 2017 +0100 @@ -0,0 +1,276 @@ +classdef Schrodinger2dCurve < scheme.Scheme + properties + m % Number of points in each direction, possibly a vector + h % Grid spacing + + grid + xm, ym + + order % Order accuracy for the approximation + + D % non-stabalized scheme operator + M % Derivative norm + H % Discrete norm + Hi + H_u, H_v % Norms in the x and y directions + Hu,Hv % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. + Hi_u, Hi_v + Hiu, Hiv + D1_v, D1_u + D2_v, D2_u + Du, Dv + + + e_w, e_e, e_s, e_n + du_w, dv_w + du_e, dv_e + du_s, dv_s + du_n, dv_n + g_1, g_2 + c + a11, a12, a22 + m_tot, m_u, m_v + p,p_tau + Ji + end + + methods + function obj = Schrodinger2dCurve(g ,order, opSet,p,p_tau) + default_arg('opSet',@sbp.D2Variable); + default_arg('c', 1); + + obj.p=p; + obj.p_tau=p_tau; + obj.c=1; + + m = g.size(); + obj.m_u = m(1); + obj.m_v = m(2); + obj.m_tot = g.N(); + + h = g.scaling(); + h_u = h(1); + h_v = h(2); + + % Operators + ops_u = opSet(obj.m_u, {0, 1}, order); + ops_v = opSet(obj.m_v, {0, 1}, order); + + I_u = speye(obj.m_u); + I_v = speye(obj.m_v); + + obj.D1_u = ops_u.D1; + obj.D2_u = ops_u.D2; + + H_u = ops_u.H; + Hi_u = ops_u.HI; + e_l_u = ops_u.e_l; + e_r_u = ops_u.e_r; + d1_l_u = ops_u.d1_l; + d1_r_u = ops_u.d1_r; + + obj.D1_v = ops_v.D1; + obj.D2_v = ops_v.D2; + H_v = ops_v.H; + Hi_v = ops_v.HI; + e_l_v = ops_v.e_l; + e_r_v = ops_v.e_r; + d1_l_v = ops_v.d1_l; + d1_r_v = ops_v.d1_r; + + obj.Du = kr(obj.D1_u,I_v); + obj.Dv = kr(I_u,obj.D1_v); + + obj.H = kr(H_u,H_v); + obj.Hi = kr(Hi_u,Hi_v); + obj.Hu = kr(H_u,I_v); + obj.Hv = kr(I_u,H_v); + obj.Hiu = kr(Hi_u,I_v); + obj.Hiv = kr(I_u,Hi_v); + + obj.e_w = kr(e_l_u,I_v); + obj.e_e = kr(e_r_u,I_v); + obj.e_s = kr(I_u,e_l_v); + obj.e_n = kr(I_u,e_r_v); + obj.du_w = kr(d1_l_u,I_v); + obj.dv_w = (obj.e_w'*obj.Dv)'; + obj.du_e = kr(d1_r_u,I_v); + obj.dv_e = (obj.e_e'*obj.Dv)'; + obj.du_s = (obj.e_s'*obj.Du)'; + obj.dv_s = kr(I_u,d1_l_v); + obj.du_n = (obj.e_n'*obj.Du)'; + obj.dv_n = kr(I_u,d1_r_v); + +% obj.x_u = x_u; +% obj.x_v = x_v; +% obj.y_u = y_u; +% obj.y_v = y_v; + + obj.m = m; + obj.h = [h_u h_v]; + obj.order = order; + obj.grid = g; + + + end + + + function [D ]= d_fun(obj,t) + % Metric derivatives + ti = parametrization.Ti.points(obj.p.p1(t),obj.p.p2(t),obj.p.p3(t),obj.p.p4(t)); + ti_tau = parametrization.Ti.points(obj.p_tau.p1(t),obj.p_tau.p2(t),obj.p_tau.p3(t),obj.p_tau.p4(t)); + + lcoords=points(obj.grid); + [obj.xm,obj.ym]= ti.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2)); + [x_tau,y_tau]= ti_tau.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2)); + x = reshape(obj.xm,obj.m_tot,1); + y = reshape(obj.ym,obj.m_tot,1); + + x_tau = reshape(x_tau,obj.m_tot,1); + y_tau = reshape(y_tau,obj.m_tot,1); + + x_u = obj.Du*x; + x_v = obj.Dv*x; + y_u = obj.Du*y; + y_v = obj.Dv*y; + + J = x_u.*y_v - x_v.*y_u; + a11 = 1./J.* (x_v.^2 + y_v.^2); + a12 = -1./J .* (x_u.*x_v + y_u.*y_v); + a22 = 1./J .* (x_u.^2 + y_u.^2); + + obj.a11=a11; + obj.a12=a12; + obj.a22=a22; + + % Assemble full operators + L_12 = spdiags(a12, 0, obj.m_tot, obj.m_tot); + Duv = obj.Du*L_12*obj.Dv; + Dvu = obj.Dv*L_12*obj.Du; + + Duu = sparse(obj.m_tot); + Dvv = sparse(obj.m_tot); + ind = grid.funcToMatrix(obj.grid, 1:obj.m_tot); + + + for i = 1:obj.m_v + D = obj.D2_u(a11(ind(:,i))); + p = ind(:,i); + Duu(p,p) = D; + end + + for i = 1:obj.m_u + D = obj.D2_v(a22(ind(i,:))); + p = ind(i,:); + Dvv(p,p) = D; + end + + Ji = spdiags(1./J, 0, obj.m_tot, obj.m_tot); + obj.Ji=Ji; + obj.g_1 = x_tau.*y_v-y_tau.*x_v; + obj.g_2 = -x_tau.*y_u + y_tau.*x_u; + + %Add the flux splitting + D = Ji*(obj.g_1.*obj.Du + obj.g_2.*obj.Dv + 1i*obj.c^2*(Duu + Duv + Dvu + Dvv)); %% g_1' och g_2'? + +% obj.gamm_u = h_u*ops_u.borrowing.M.d1; +% obj.gamm_v = h_v*ops_v.borrowing.M.d1; + + end + + % Closure functions return the opertors applied to the own doamin to close the boundary + % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. + % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. + % type is a string specifying the type of boundary condition if there are several. + % data is a function returning the data that should be applied at the boundary. + % neighbour_scheme is an instance of Scheme that should be interfaced to. + % neighbour_boundary is a string specifying which boundary to interface to. + function [closure, penalty] = boundary_condition(obj, boundary) + [e, d_n, d_t, coeff_t, s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g] = obj.get_boundary_ops(boundary); + + a_t = spdiag(coeff_t); + F = (s * d_n' + s * a_t*d_t')'; + tau1 = 1; + a = spdiag(g); + tau2 = (-1*s*a - abs(a))/4; + + penalty_parameter_1 = 1i*halfnorm_inv_n*halfnorm_inv_t*e*e'*F*halfnorm_t; + penalty_parameter_2 = halfnorm_inv_n*e*tau2; + + closure = obj.Ji*obj.c^2 * penalty_parameter_1*e' +obj.Ji* penalty_parameter_2*e'; + % penalty = -obj.Ji*obj.c^2 * penalty_parameter_2; + penalty=0; + + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + end + + function [e, d_n, d_t, coeff_t, s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g, I] = get_boundary_ops(obj, boundary) + + % gridMatrix = zeros(obj.m(2),obj.m(1)); + % gridMatrix(:) = 1:numel(gridMatrix); + + ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m)); + + switch boundary + case 'w' + e = obj.e_w; + d_n = obj.du_w; + d_t = obj.dv_w; + s = -1; + + I = ind(1,:); + coeff_t = obj.a12(I); + g=obj.g_1(I); + case 'e' + e = obj.e_e; + d_n = obj.du_e; + d_t = obj.dv_e; + s = 1; + + I = ind(end,:); + coeff_t = obj.a12(I); + g=obj.g_1(I); + case 's' + e = obj.e_s; + d_n = obj.dv_s; + d_t = obj.du_s; + s = -1; + + I = ind(:,1)'; + coeff_t = obj.a12(I); + g=obj.g_2(I); + case 'n' + e = obj.e_n; + d_n = obj.dv_n; + d_t = obj.du_n; + s = 1; + + I = ind(:,end)'; + coeff_t = obj.a12(I); + g=obj.g_2(I); + otherwise + error('No such boundary: boundary = %s',boundary); + end + + switch boundary + case {'w','e'} + halfnorm_inv_n = obj.Hiu; + halfnorm_inv_t = obj.Hiv; + halfnorm_t = obj.Hv; + + case {'s','n'} + halfnorm_inv_n = obj.Hiv; + halfnorm_inv_t = obj.Hiu; + halfnorm_t = obj.Hu; + end + end + + function N = size(obj) + N = prod(obj.m); + end + + + end +end \ No newline at end of file