Mercurial > repos > public > sbplib
changeset 378:18525f1bb941
Merged in feature/hypsyst (pull request #4)
Feature/hypsyst
| author | Jonatan Werpers <jonatan.werpers@it.uu.se> |
|---|---|
| date | Thu, 26 Jan 2017 13:07:51 +0000 |
| parents | f1289f3b86b8 (current diff) cf0bef311ce2 (diff) |
| children | 359861563866 3fdfad20037e d6d27fdc342a 50fd7e88aa74 29944ea7674b |
| files | |
| diffstat | 6 files changed, 1899 insertions(+), 13 deletions(-) [+] |
line wrap: on
line diff
diff -r f1289f3b86b8 -r 18525f1bb941 +grid/Ti3D.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/Ti3D.m Thu Jan 26 13:07:51 2017 +0000 @@ -0,0 +1,253 @@ +classdef Ti3D + properties + gs % {6}Surfaces + V % FunctionHandle(XI,ETA,ZETA) + end + + methods + % TODO write all fancy features for flipping around with the surfaces + % Each surface is defined with an outward facing outward and choosing + % the "corner" where XI=0 if not possible the corner where ETA=0 is choosen + function obj = Ti3D(CW,CE,CS,CN,CB,CT) + obj.gs = {CE,CW,CS,CN,CB,CT}; + + gw = CW.g; + ge = CE.g; + gs = CS.g; + gn = CN.g; + gb = CB.g; + gt = CT.g; + + function o = V_fun(XI,ETA,ZETA) + XI=XI'; + ETA=ETA'; + ZETA=ZETA'; + + one=0*ETA+1; + zero=0*ETA; + + Sw = gw(ETA,(1-ZETA)); + Se = ge((1-ETA),(1-ZETA)); + Ss = gs(XI,ZETA); + Sn = gn((1-XI),(1-ZETA)); + Sb = gb((1-XI),ETA); + St = gt(XI,ETA); + + Ewt = gw(ETA,zero); + Ewb = gw(ETA,one); + Ews = gw(zero,1-ZETA); + Ewn = gw(one,1-ZETA); + Eet = ge(1-ETA,zero); + Eeb = ge(1-ETA,one); + Ees = ge(one,1-ZETA); + Een = ge(zero,1-ZETA); + Enb = gn(1-XI,one); + Ent = gn(1-XI,zero); + Est = gs(XI,one); + Esb = gs(XI,zero); + + Cwbs = gw(zero,one); + Cwbn = gw(one,one); + Cwts = gw(zero,zero); + Cwtn = gw(one,zero); + Cebs = ge(one,one); + Cebn = ge(zero,one); + Cets = ge(one,zero); + Cetn = ge(zero,zero); + + + X1 = (1-XI).*Sw(1,:,:) + XI.*Se(1,:,:); + X2 = (1-ETA).*Ss(1,:,:) + ETA.*Sn(1,:,:); + X3 = (1-ZETA).*Sb(1,:,:) + ZETA.*St(1,:,:); + + X12 = (1-XI).*(1-ETA).*Ews(1,:,:) + (1-XI).*ETA.*Ewn(1,:,:) + XI.*(1-ETA).*Ees(1,:,:) + XI.*ETA.*Een(1,:,:); + X13 = (1-XI).*(1-ZETA).*Ewb(1,:,:) + (1-XI).*ZETA.*Ewt(1,:,:) + XI.*(1-ZETA).*Eeb(1,:,:) + XI.*ZETA.*Eet(1,:,:); + X23 = (1-ETA).*(1-ZETA).*Esb(1,:,:) + (1-ETA).*ZETA.*Est(1,:,:) + ETA.*(1-ZETA).*Enb(1,:,:) + ETA.*ZETA.*Ent(1,:,:); + + X123 = (1-XI).*(1-ETA).*(1-ZETA).*Cwbs(1,:,:) + (1-XI).*(1-ETA).*ZETA.*Cwts(1,:,:) + (1-XI).*ETA.*(1-ZETA).*Cwbn(1,:,:) + ... + (1-XI).*ETA.*ZETA.*Cwtn(1,:,:) + XI.*(1-ETA).*(1-ZETA).*Cebs(1,:,:) + XI.*(1-ETA).*ZETA.*Cets(1,:,:) + ... + XI.*ETA.*(1-ZETA).*Cebn(1,:,:) + XI.*ETA.*ZETA.*Cetn(1,:,:); + + X = X1 + X2 + X3 - X12 - X13 - X23 + X123; + + + Y1 = (1-XI).*Sw(2,:,:) + XI.*Se(2,:,:); + Y2 = (1-ETA).*Ss(2,:,:) + ETA.*Sn(2,:,:); + Y3 = (1-ZETA).*Sb(2,:,:) + ZETA.*St(2,:,:); + + Y12 = (1-XI).*(1-ETA).*Ews(2,:,:) + (1-XI).*ETA.*Ewn(2,:,:) + XI.*(1-ETA).*Ees(2,:,:) + XI.*ETA.*Een(2,:,:); + Y13 = (1-XI).*(1-ZETA).*Ewb(2,:,:) + (1-XI).*ZETA.*Ewt(2,:,:) + XI.*(1-ZETA).*Eeb(2,:,:) + XI.*ZETA.*Eet(2,:,:); + Y23 = (1-ETA).*(1-ZETA).*Esb(2,:,:) + (1-ETA).*ZETA.*Est(2,:,:) + ETA.*(1-ZETA).*Enb(2,:,:) + ETA.*ZETA.*Ent(2,:,:); + + Y123 = (1-XI).*(1-ETA).*(1-ZETA).*Cwbs(2,:,:) + (1-XI).*(1-ETA).*ZETA.*Cwts(2,:,:) + (1-XI).*ETA.*(1-ZETA).*Cwbn(2,:,:) + ... + (1-XI).*ETA.*ZETA.*Cwtn(2,:,:) + XI.*(1-ETA).*(1-ZETA).*Cebs(2,:,:) + XI.*(1-ETA).*ZETA.*Cets(2,:,:) + ... + XI.*ETA.*(1-ZETA).*Cebn(2,:,:) + XI.*ETA.*ZETA.*Cetn(2,:,:); + + Y = Y1 + Y2 + Y3 - Y12 - Y13 - Y23 + Y123; + + + Z1 = (1-XI).*Sw(3,:,:) + XI.*Se(3,:,:); + Z2 = (1-ETA).*Ss(3,:,:) + ETA.*Sn(3,:,:); + Z3 = (1-ZETA).*Sb(3,:,:) + ZETA.*St(3,:,:); + + Z12 = (1-XI).*(1-ETA).*Ews(3,:,:) + (1-XI).*ETA.*Ewn(3,:,:) + XI.*(1-ETA).*Ees(3,:,:) + XI.*ETA.*Een(3,:,:); + Z13 = (1-XI).*(1-ZETA).*Ewb(3,:,:) + (1-XI).*ZETA.*Ewt(3,:,:) + XI.*(1-ZETA).*Eeb(3,:,:) + XI.*ZETA.*Eet(3,:,:); + Z23 = (1-ETA).*(1-ZETA).*Esb(3,:,:) + (1-ETA).*ZETA.*Est(3,:,:) + ETA.*(1-ZETA).*Enb(3,:,:) + ETA.*ZETA.*Ent(3,:,:); + + Z123 = (1-XI).*(1-ETA).*(1-ZETA).*Cwbs(3,:,:) + (1-XI).*(1-ETA).*ZETA.*Cwts(3,:,:) + (1-XI).*ETA.*(1-ZETA).*Cwbn(3,:,:) + ... + (1-XI).*ETA.*ZETA.*Cwtn(3,:,:) + XI.*(1-ETA).*(1-ZETA).*Cebs(3,:,:) + XI.*(1-ETA).*ZETA.*Cets(3,:,:) + ... + XI.*ETA.*(1-ZETA).*Cebn(3,:,:) + XI.*ETA.*ZETA.*Cetn(3,:,:); + + Z = Z1 + Z2 + Z3 - Z12 - Z13 - Z23 + Z123; + o = [X;Y;Z]; + end + + obj.V = @V_fun; + end + + %Should be rewritten so that the input is xi eta zeta + function [X,Y,Z] = map(obj,XI,ETA,ZETA) + + V = obj.V; + + p = V(XI,ETA,ZETA); + X = p(1,:)'; + Y = p(2,:)'; + Z = p(3,:)'; + + end + + % function h = plot(obj,nu,nv) + % S = obj.S; + % + % default_arg('nv',nu) + % + % u = linspace(0,1,nu); + % v = linspace(0,1,nv); + % + % m = 100; + % + % X = zeros(nu+nv,m); + % Y = zeros(nu+nv,m); + % + % + % t = linspace(0,1,m); + % for i = 1:nu + % p = S(u(i),t); + % X(i,:) = p(1,:); + % Y(i,:) = p(2,:); + % end + % + % for i = 1:nv + % p = S(t,v(i)); + % X(i+nu,:) = p(1,:); + % Y(i+nu,:) = p(2,:); + % end + % + % h = line(X',Y'); + % end + % + % + % function h = show(obj,nu,nv) + % default_arg('nv',nu) + % S = obj.S; + % + % if(nu>2 || nv>2) + % h_grid = obj.plot(nu,nv); + % set(h_grid,'Color',[0 0.4470 0.7410]); + % end + % + % h_bord = obj.plot(2,2); + % set(h_bord,'Color',[0.8500 0.3250 0.0980]); + % set(h_bord,'LineWidth',2); + % end + % + % + % % TRANSFORMATIONS + % function ti = translate(obj,a) + % gs = obj.gs; + % + % for i = 1:length(gs) + % new_gs{i} = gs{i}.translate(a); + % end + % + % ti = grid.Ti(new_gs{:}); + % end + % + % % Mirrors the Ti so that the resulting Ti is still left handed. + % % (Corrected by reversing curves and switching e and w) + % function ti = mirror(obj, a, b) + % gs = obj.gs; + % + % new_gs = cell(1,4); + % + % new_gs{1} = gs{1}.mirror(a,b).reverse(); + % new_gs{3} = gs{3}.mirror(a,b).reverse(); + % new_gs{2} = gs{4}.mirror(a,b).reverse(); + % new_gs{4} = gs{2}.mirror(a,b).reverse(); + % + % ti = grid.Ti(new_gs{:}); + % end + % + % function ti = rotate(obj,a,rad) + % gs = obj.gs; + % + % for i = 1:length(gs) + % new_gs{i} = gs{i}.rotate(a,rad); + % end + % + % ti = grid.Ti(new_gs{:}); + % end + % + % function ti = rotate_edges(obj,n); + % new_gs = cell(1,4); + % for i = 0:3 + % new_i = mod(i - n,4); + % new_gs{new_i+1} = obj.gs{i+1}; + % end + % ti = grid.Ti(new_gs{:}); + % end + % end + % + % methods(Static) + % function obj = points(p1, p2, p3, p4) + % g1 = grid.Curve.line(p1,p2); + % g2 = grid.Curve.line(p2,p3); + % g3 = grid.Curve.line(p3,p4); + % g4 = grid.Curve.line(p4,p1); + % + % obj = grid.Ti(g1,g2,g3,g4); + % end + % + % function label(varargin) + % if nargin == 2 && ischar(varargin{2}) + % label_impl(varargin{:}); + % else + % for i = 1:length(varargin) + % label_impl(varargin{i},inputname(i)); + % end + % end + % + % + % function label_impl(ti,str) + % S = ti.S; + % + % pc = S(0.5,0.5); + % + % margin = 0.1; + % pw = S( margin, 0.5); + % pe = S(1-margin, 0.5); + % ps = S( 0.5, margin); + % pn = S( 0.5, 1-margin); + % + % + % ti.show(2,2); + % grid.place_label(pc,str); + % grid.place_label(pw,'w'); + % grid.place_label(pe,'e'); + % grid.place_label(ps,'s'); + % grid.place_label(pn,'n'); + % end + % end + end +end \ No newline at end of file
diff -r f1289f3b86b8 -r 18525f1bb941 +scheme/Hypsyst2d.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Hypsyst2d.m Thu Jan 26 13:07:51 2017 +0000 @@ -0,0 +1,301 @@ +classdef Hypsyst2d < scheme.Scheme + properties + m % Number of points in each direction, possibly a vector + n %size of system + h % Grid spacing + x,y % Grid + X,Y % Values of x and y for each grid point + order % Order accuracy for the approximation + + D % non-stabalized scheme operator + A, B, E %Coefficient matrices + + H % Discrete norm + % Norms in the x and y directions + Hxi,Hyi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. + I_x,I_y, I_N + e_w, e_e, e_s, e_n + params %parameters for the coeficient matrice + end + + methods + %Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Eu + function obj = Hypsyst2d(m, lim, order, A, B, E, params) + default_arg('E', []) + xlim = lim{1}; + ylim = lim{2}; + + if length(m) == 1 + m = [m m]; + end + + obj.A=A; + obj.B=B; + obj.E=E; + + m_x = m(1); + m_y = m(2); + obj.params = params; + + ops_x = sbp.D2Standard(m_x,xlim,order); + ops_y = sbp.D2Standard(m_y,ylim,order); + + obj.x = ops_x.x; + obj.y = ops_y.x; + + obj.X = kr(obj.x,ones(m_y,1)); + obj.Y = kr(ones(m_x,1),obj.y); + + Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y); + Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y); + Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y); + + obj.n = length(A(obj.params,0,0)); + + I_n = eye(obj.n);I_x = speye(m_x); + obj.I_x = I_x; + I_y = speye(m_y); + obj.I_y = I_y; + + + D1_x = kr(I_n, ops_x.D1, I_y); + obj.Hxi = kr(I_n, ops_x.HI, I_y); + D1_y = kr(I_n, I_x, ops_y.D1); + obj.Hyi = kr(I_n, I_x, ops_y.HI); + + obj.e_w = kr(I_n, ops_x.e_l, I_y); + obj.e_e = kr(I_n, ops_x.e_r, I_y); + obj.e_s = kr(I_n, I_x, ops_y.e_l); + obj.e_n = kr(I_n, I_x, ops_y.e_r); + + obj.m = m; + obj.h = [ops_x.h ops_y.h]; + obj.order = order; + + obj.D = -Aevaluated*D1_x-Bevaluated*D1_y-Eevaluated; + + 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. + function [closure, penalty] = boundary_condition(obj,boundary,type,L) + default_arg('type','char'); + switch type + case{'c','char'} + [closure,penalty] = boundary_condition_char(obj,boundary); + case{'general'} + [closure,penalty] = boundary_condition_general(obj,boundary,L); + otherwise + error('No such boundary condition') + end + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + error('An interface function does not exist yet'); + end + + function N = size(obj) + N = obj.m; + end + + function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y) + params = obj.params; + + if isa(mat,'function_handle') + [rows,cols] = size(mat(params,0,0)); + matVec = mat(params,X',Y'); + matVec = sparse(matVec); + side = max(length(X),length(Y)); + else + matVec = mat; + [rows,cols] = size(matVec); + side = max(length(X),length(Y)); + cols = cols/side; + end + ret = cell(rows,cols); + + for ii = 1:rows + for jj=1:cols + ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side)); + end + end + ret = cell2mat(ret); + end + + %Characteristic boundary conditions + function [closure, penalty] = boundary_condition_char(obj,boundary) + params = obj.params; + x = obj.x; + y = obj.y; + + switch boundary + case {'w','W','west'} + e_ = obj.e_w; + mat = obj.A; + boundPos = 'l'; + Hi = obj.Hxi; + [V,Vi,D,signVec] = obj.matrixDiag(mat,x(1),y); + side = max(length(y)); + case {'e','E','east'} + e_ = obj.e_e; + mat = obj.A; + boundPos = 'r'; + Hi = obj.Hxi; + [V,Vi,D,signVec] = obj.matrixDiag(mat,x(end),y); + side = max(length(y)); + case {'s','S','south'} + e_ = obj.e_s; + mat = obj.B; + boundPos = 'l'; + Hi = obj.Hyi; + [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(1)); + side = max(length(x)); + case {'n','N','north'} + e_ = obj.e_n; + mat = obj.B; + boundPos = 'r'; + Hi = obj.Hyi; + [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(end)); + side = max(length(x)); + end + pos = signVec(1); + zeroval = signVec(2); + neg = signVec(3); + + switch boundPos + case {'l'} + tau = sparse(obj.n*side,pos); + Vi_plus = Vi(1:pos,:); + tau(1:pos,:) = -abs(D(1:pos,1:pos)); + closure = Hi*e_*V*tau*Vi_plus*e_'; + penalty = -Hi*e_*V*tau*Vi_plus; + case {'r'} + tau = sparse(obj.n*side,neg); + tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); + Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); + closure = Hi*e_*V*tau*Vi_minus*e_'; + penalty = -Hi*e_*V*tau*Vi_minus; + end + end + + % General boundary condition in the form Lu=g(x) + function [closure,penalty] = boundary_condition_general(obj,boundary,L) + params = obj.params; + x = obj.x; + y = obj.y; + + switch boundary + case {'w','W','west'} + e_ = obj.e_w; + mat = obj.A; + boundPos = 'l'; + Hi = obj.Hxi; + [V,Vi,D,signVec] = obj.matrixDiag(mat,x(1),y); + L = obj.evaluateCoefficientMatrix(L,x(1),y); + side = max(length(y)); + case {'e','E','east'} + e_ = obj.e_e; + mat = obj.A; + boundPos = 'r'; + Hi = obj.Hxi; + [V,Vi,D,signVec] = obj.matrixDiag(mat,x(end),y); + L = obj.evaluateCoefficientMatrix(L,x(end),y); + side = max(length(y)); + case {'s','S','south'} + e_ = obj.e_s; + mat = obj.B; + boundPos = 'l'; + Hi = obj.Hyi; + [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(1)); + L = obj.evaluateCoefficientMatrix(L,x,y(1)); + side = max(length(x)); + case {'n','N','north'} + e_ = obj.e_n; + mat = obj.B; + boundPos = 'r'; + Hi = obj.Hyi; + [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(end)); + L = obj.evaluateCoefficientMatrix(L,x,y(end)); + side = max(length(x)); + end + + pos = signVec(1); + zeroval = signVec(2); + neg = signVec(3); + + switch boundPos + case {'l'} + tau = sparse(obj.n*side,pos); + Vi_plus = Vi(1:pos,:); + Vi_minus = Vi(pos+zeroval+1:obj.n*side,:); + V_plus = V(:,1:pos); + V_minus = V(:,(pos+zeroval)+1:obj.n*side); + + tau(1:pos,:) = -abs(D(1:pos,1:pos)); + R = -inv(L*V_plus)*(L*V_minus); + closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; + penalty = -Hi*e_*V*tau*inv(L*V_plus)*L; + case {'r'} + tau = sparse(obj.n*side,neg); + tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); + Vi_plus = Vi(1:pos,:); + Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); + + V_plus = V(:,1:pos); + V_minus = V(:,(pos+zeroval)+1:obj.n*side); + R = -inv(L*V_minus)*(L*V_plus); + closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_'; + penalty = -Hi*e_*V*tau*inv(L*V_minus)*L; + end + end + + % Function that diagonalizes a symbolic matrix A as A=V*D*Vi + % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign + % [d+ ] + % D = [ d0 ] + % [ d-] + % signVec is a vector specifying the number of possitive, zero and negative eigenvalues of D + function [V,Vi, D,signVec] = matrixDiag(obj,mat,x,y) + params = obj.params; + syms xs ys + [V, D]= eig(mat(params,xs,ys)); + Vi = inv(V); + xs = x; + ys = y; + + side = max(length(x),length(y)); + Dret = zeros(obj.n,side*obj.n); + Vret = zeros(obj.n,side*obj.n); + Viret = zeros(obj.n,side*obj.n); + + for ii = 1:obj.n + for jj = 1:obj.n + Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii)); + Vret(jj,(ii-1)*side+1:side*ii) = eval(V(jj,ii)); + Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii)); + end + end + + D = sparse(Dret); + V = sparse(Vret); + Vi = sparse(Viret); + V = obj.evaluateCoefficientMatrix(V,x,y); + Vi = obj.evaluateCoefficientMatrix(Vi,x,y); + D = obj.evaluateCoefficientMatrix(D,x,y); + DD = diag(D); + + poseig = (DD>0); + zeroeig = (DD==0); + negeig = (DD<0); + + D = diag([DD(poseig); DD(zeroeig); DD(negeig)]); + V = [V(:,poseig) V(:,zeroeig) V(:,negeig)]; + Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)]; + signVec = [sum(poseig),sum(zeroeig),sum(negeig)]; + end + + end +end \ No newline at end of file
diff -r f1289f3b86b8 -r 18525f1bb941 +scheme/Hypsyst2dCurve.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Hypsyst2dCurve.m Thu Jan 26 13:07:51 2017 +0000 @@ -0,0 +1,378 @@ +classdef Hypsyst2dCurve < scheme.Scheme + properties + m % Number of points in each direction, possibly a vector + n % size of system + h % Grid spacing + X,Y % Values of x and y for each grid point + + J, Ji % Jacobaian and inverse Jacobian + xi,eta + Xi,Eta + + A,B + X_eta, Y_eta + X_xi,Y_xi + order % Order accuracy for the approximation + + D % non-stabalized scheme operator + Ahat, Bhat, E + + H % Discrete norm + Hxii,Hetai % Kroneckerd norms in xi and eta. + I_xi,I_eta, I_N, onesN + e_w, e_e, e_s, e_n + index_w, index_e,index_s,index_n + params % Parameters for the coeficient matrice + end + + + methods + % Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Eu + function obj = Hypsyst2dCurve(m, order, A, B, E, params, ti) + default_arg('E', []) + xilim = {0 1}; + etalim = {0 1}; + + if length(m) == 1 + m = [m m]; + end + obj.params = params; + obj.A=A; + obj.B=B; + + obj.Ahat=@(params,x,y,x_eta,y_eta)(A(params,x,y).*y_eta-B(params,x,y).*x_eta); + obj.Bhat=@(params,x,y,x_xi,y_xi)(B(params,x,y).*x_xi-A(params,x,y).*y_xi); + obj.E=@(params,x,y,~,~)E(params,x,y); + + m_xi = m(1); + m_eta = m(2); + m_tot=m_xi*m_eta; + + ops_xi = sbp.D2Standard(m_xi,xilim,order); + ops_eta = sbp.D2Standard(m_eta,etalim,order); + + obj.xi = ops_xi.x; + obj.eta = ops_eta.x; + + obj.Xi = kr(obj.xi,ones(m_eta,1)); + obj.Eta = kr(ones(m_xi,1),obj.eta); + + obj.n = length(A(obj.params,0,0)); + obj.onesN=ones(obj.n); + + obj.index_w=1:m_eta; + obj.index_e=(m_tot-m_e + + metric_termsta+1):m_tot; + obj.index_s=1:m_eta:(m_tot-m_eta+1); + obj.index_n=(m_eta):m_eta:m_tot; + + I_n = eye(obj.n); + I_xi = speye(m_xi); + obj.I_xi = I_xi; + I_eta = speye(m_eta); + obj.I_eta = I_eta; + + D1_xi = kr(I_n, ops_xi.D1, I_eta); + obj.Hxii = kr(I_n, ops_xi.HI, I_eta); + D1_eta = kr(I_n, I_xi, ops_eta.D1); + obj.Hetai = kr(I_n, I_xi, ops_eta.HI); + + obj.e_w = kr(I_n, ops_xi.e_l, I_eta); + obj.e_e = kr(I_n, ops_xi.e_r, I_eta); + obj.e_s = kr(I_n, I_xi, ops_eta.e_l); + obj.e_n = kr(I_n, I_xi, + + metric_termsops_eta.e_r); + + [X,Y] = ti.map(obj.xi,obj.eta); + + [x_xi,x_eta] = gridDerivatives(X,ops_xi.D1,ops_eta.D1); + [y_xi,y_eta] = gridDerivatives(Y,ops_xi.D1, ops_eta.D1); + + obj.X = reshape(X,m_tot,1); + obj.Y = reshape(Y,m_tot,1); + obj.X_xi = reshape(x_xi,m_tot,1); + obj.Y_xi = reshape(y_xi,m_tot,1); + obj.X_eta = reshape(x_eta,m_tot,1); + obj.Y_eta = reshape(y_eta,m_tot,1); + + Ahat_evaluated = obj.evaluateCoefficientMatrix(obj.Ahat, obj.X, obj.Y,obj.X_eta,obj.Y_eta); + Bhat_evaluated = obj.evaluateCoefficientMatrix(obj.Bhat, obj.X, obj.Y,obj.X_xi,obj.Y_xi); + E_evaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,[],[]); + + obj.m = m; + obj.h = [ops_xi.h ops_eta.h]; + obj.order = order; + obj.J = obj.X_xi.*obj.Y_eta - obj.X_eta.*obj.Y_xi; + obj.Ji = kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot)); + + obj.D = obj.Ji*(-Ahat_evaluated*D1_xi-Bhat_evaluated*D1_eta)-E_evaluated; + + 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',General boundary conditions'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. + function [closure, penalty] = boundary_condition(obj,boundary,type,L) + default_arg('type','char'); + switch type + case{'c','char'} + [closure,penalty] = boundary_condition_char(obj,boundary); + case{'general'} + [closure,penalty] = boundary_condition_general(obj,boundary,L); + otherwise + error('No such boundary condition') + end + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundaryGeneral boundary conditions) + error('An interface function does not exist yet'); + end + + function N = size(obj) + N = obj.m; + end + + function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y,x_,y_) + params = obj.params; + + if isa(mat,'function_handle') + [rows,cols] = size(mat(params,0,0,0,0)); + x_ = kr(obj.onesN,x_); + y_ = kr(obj.onesN,y_); + matVec = mat(params,X',Y',x_',y_'); + matVec = sparse(matVec); + side = max(length(X),length(Y)); + else + matVec = mat; + [rows,cols] = size(matVec); + side = max(length(X),length(Y)); + cols = cols/side; + end + + ret = cell(rows,cols); + for ii = 1:rows + for jj = 1:cols + ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side)); + end + end + ret = cell2mat(ret); + end + + %Characteristic boundary conditions + function [closure, penalty] = boundary_condition_char(obj,boundary) + params = obj.params; + X = obj.X; + Y = obj.Y; + xi = obj.xi; + eta = obj.eta; + + switch boundary + case {'w','W','west'} + e_ = obj.e_w; + mat = obj.Ahat; + boundPos = 'l'; + Hi = obj.Hxii; + [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_w),Y(obj.index_w),obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w)); + side = max(length(eta)); + case {'e','E','east'} + e_ = obj.e_e; + mat = obj.Ahat; + boundPos = 'r'; + Hi = obj.Hxii; + [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e)); + side = max(length(eta)); + case {'s','S','south'} + e_ = obj.e_s; + mat = obj.Bhat; + boundPos = 'l'; + Hi = obj.Hetai; + [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_s),Y(obj.index_s),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s)); + side = max(length(xi)); + case {'n','N','north'} + e_ = obj.e_n; + mat = obj.Bhat; + boundPos = 'r'; + Hi = obj.Hetai; + [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n)); + side = max(length(xi)); + end + + pos = signVec(1); + zeroval = signVec(2); + neg = signVec(3); + + switch boundPos + case {'l'} + tau = sparse(obj.n*side,pos); + Vi_plus = Vi(1:pos,:); + tau(1:pos,:) = -abs(D(1:pos,1:pos)); + closure = Hi*e_*V*tau*Vi_plus*e_'; + penalty = -Hi*e_*V*tau*Vi_plus; + case {'r'} + tau = sparse(obj.n*side,neg); + tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); + Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); + closure = Hi*e_*V*tau*Vi_minus*e_'; + penalty = -Hi*e_*V*tau*Vi_minus; + end + end + + + % General boundary condition in the form Lu=g(x) + function [closure,penalty] = boundary_condition_general(obj,boundary,L) + params = obj.params; + X = obj.X; + Y = obj.Y; + xi = obj.xi; + eta = obj.eta; + + switch boundary + case {'w','W','west'} + e_ = obj.e_w; + mat = obj.Ahat; + boundPos = 'l'; + Hi = obj.Hxii; + [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_w),Y(obj.index_w),obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w)); + + Ji_vec = diag(obj.Ji); + Ji = diag(Ji_vec(obj.index_w)); + xi_x = Ji*obj.Y_eta(obj.index_w); + xi_y = -Ji*obj.X_eta(obj.index_w); + L = obj.evaluateCoefficientMatrix(L,xi_x,xi_y,[],[]); + side = max(length(eta)); + case {'e','E','east'} + e_ = obj.e_e; + mat = obj.Ahat; + boundPos = 'r'; + Hi = obj.Hxii; + [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e)); + + Ji_vec = diag(obj.Ji); + Ji = diag(Ji_vec(obj.index_e)); + xi_x = Ji*obj.Y_eta(obj.index_e); + xi_y = -Ji*obj.X_eta(obj.index_e); + L = obj.evaluateCoefficientMatrix(L,-xi_x,-xi_y,[],[]); + side = max(length(eta)); + case {'s','S','south'} + e_ = obj.e_s; + mat = obj.Bhat; + boundPos = 'l'; + Hi = obj.Hetai; + [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_s),Y(obj.index_s),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s)); + + Ji_vec = diag(obj.Ji); + Ji = diag(Ji_vec(obj.index_s)); + eta_x = Ji*obj.Y_xi(obj.index_s); + eta_y = -Ji*obj.X_xi(obj.index_s); + L = obj.evaluateCoefficientMatrix(L,eta_x,eta_y,[],[]); + side = max(length(xi)); + case {'n','N','north'} + e_ = obj.e_n; + mat = obj.Bhat; + boundPos = 'r'; + Hi = obj.Hetai; + [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n)); + + Ji_vec = diag(obj.Ji); + Ji = diag(Ji_vec(obj.index_n)); + eta_x = Ji*obj.Y_xi(obj.index_n); + eta_y = -Ji*obj.X_xi(obj.index_n); + L = obj.evaluateCoefficientMatrix(L,-eta_x,-eta_y,[],[]); + side = max(length(xi)); + end + + pos = signVec(1); + zeroval = signVec(2); + neg = signVec(3); + + switch boundPos + case {'l'} + tau = sparse(obj.n*side,pos); + Vi_plus = Vi(1:pos,:); + Vi_minus = Vi(pos+1:obj.n*side,:); + V_plus = V(:,1:pos); + V_minus = V(:,(pos)+1:obj.n*side); + + tau(1:pos,:) = -abs(D(1:pos,1:pos)); + R = -inv(L*V_plus)*(L*V_minus); + closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; + penalty = -Hi*e_*V*tau*inv(L*V_plus)*L; + case {'r'} + tau = sparse(obj.n*side,neg); + tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); + Vi_plus = Vi(1:pos,:); + Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); + + V_plus = V(:,1:pos); + V_minus = V(:,(pos+zeroval)+1:obj.n*side); + R = -inv(L*V_minus)*(L*V_plus); + closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_'; + penalty = -Hi*e_*V*tau*inv(L*V_minus)*L; + end + end + + % Function that diagonalizes a symbolic matrix A as A=V*D*Vi + % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign + % [d+ ] + % D = [ d0 ] + % [ d-] + % signVec is a vector specifying the number of possitive, zero and negative eigenvalues of D + function [V,Vi, D,signVec] = matrixDiag(obj,mat,x,y,x_,y_) + params = obj.params; + syms xs ys + if(sum(abs(x_)) ~= 0) + syms xs_ + else + xs_ = 0; + end + + if(sum(abs(y_))~= 0) + syms ys_; + else + ys_ = 0; + end + + [V, D] = eig(mat(params,xs,ys,xs_,ys_)); + Vi = inv(V); + syms xs ys xs_ ys_ + + xs = x; + ys = y; + xs_ = x_; + ys_ = y_; + + side = max(length(x),length(y)); + Dret = zeros(obj.n,side*obj.n); + Vret = zeros(obj.n,side*obj.n); + Viret = zeros(obj.n,side*obj.n); + for ii = 1:obj.n + for jj = 1:obj.n + Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii)); + Vret(jj,(ii-1)*side+1:side*ii) = eval(V(jj,ii)); + Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii)); + end + end + + D = sparse(Dret); + V = sparse(Vret); + Vi = sparse(Viret); + V = obj.evaluateCoefficientMatrix(V,x,y,x_,y_); + D = obj.evaluateCoefficientMatrix(D,x,y,x_,y_); + Vi = obj.evaluateCoefficientMatrix(Vi,x,y,x_,y_); + DD = diag(D); + + poseig = (DD>0); + zeroeig = (DD==0); + negeig = (DD<0); + + D = diag([DD(poseig); DD(zeroeig); DD(negeig)]); + V = [V(:,poseig) V(:,zeroeig) V(:,negeig)]; + Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)]; + signVec = [sum(poseig),sum(zeroeig),sum(negeig)]; + end + end +end \ No newline at end of file
diff -r f1289f3b86b8 -r 18525f1bb941 +scheme/Hypsyst3d.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Hypsyst3d.m Thu Jan 26 13:07:51 2017 +0000 @@ -0,0 +1,380 @@ +classdef Hypsyst3d < scheme.Scheme + properties + m % Number of points in each direction, possibly a vector + n % Size of system + h % Grid spacing + x, y, z % Grid + X, Y, Z% Values of x and y for each grid point + Yx, Zx, Xy, Zy, Xz, Yz %Grid values for boundary surfaces + order % Order accuracy for the approximation + + D % non-stabalized scheme operator + A, B, C, E % Symbolic coefficient matrices + Aevaluated,Bevaluated,Cevaluated, Eevaluated + + H % Discrete norm + Hx, Hy, Hz % Norms in the x, y and z directions + Hxi,Hyi, Hzi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. + I_x,I_y, I_z, I_N + e_w, e_e, e_s, e_n, e_b, e_t + params % Parameters for the coeficient matrice + end + + + methods + % Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Cu_z-Eu + function obj = Hypsyst3d(m, lim, order, A, B,C, E, params,operator) + default_arg('E', []) + xlim = lim{1}; + ylim = lim{2}; + zlim = lim{3}; + + if length(m) == 1 + m = [m m m]; + end + + obj.A = A; + obj.B = B; + obj.C = C; + obj.E = E; + m_x = m(1); + m_y = m(2); + m_z = m(3); + obj.params = params; + + switch operator + case 'upwind' + ops_x = sbp.D1Upwind(m_x,xlim,order); + ops_y = sbp.D1Upwind(m_y,ylim,order); + ops_z = sbp.D1Upwind(m_z,zlim,order); + otherwise + ops_x = sbp.D2Standard(m_x,xlim,order); + ops_y = sbp.D2Standard(m_y,ylim,order); + ops_z = sbp.D2Standard(m_z,zlim,order); + end + + obj.x = ops_x.x; + obj.y = ops_y.x; + obj.z = ops_z.x; + + obj.X = kr(obj.x,ones(m_y,1),ones(m_z,1)); + obj.Y = kr(ones(m_x,1),obj.y,ones(m_z,1)); + obj.Z = kr(ones(m_x,1),ones(m_y,1),obj.z); + + obj.Yx = kr(obj.y,ones(m_z,1)); + obj.Zx = kr(ones(m_y,1),obj.z); + obj.Xy = kr(obj.x,ones(m_z,1)); + obj.Zy = kr(ones(m_x,1),obj.z); + obj.Xz = kr(obj.x,ones(m_y,1)); + obj.Yz = kr(ones(m_z,1),obj.y); + + obj.Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y,obj.Z); + obj.Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y,obj.Z); + obj.Cevaluated = obj.evaluateCoefficientMatrix(C, obj.X, obj.Y,obj.Z); + obj.Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y,obj.Z); + + obj.n = length(A(obj.params,0,0,0)); + + I_n = speye(obj.n); + I_x = speye(m_x); + obj.I_x = I_x; + I_y = speye(m_y); + obj.I_y = I_y; + I_z = speye(m_z); + obj.I_z = I_z; + I_N = kr(I_n,I_x,I_y,I_z); + + obj.Hxi = kr(I_n, ops_x.HI, I_y,I_z); + obj.Hx = ops_x.H; + obj.Hyi = kr(I_n, I_x, ops_y.HI,I_z); + obj.Hy = ops_y.H; + obj.Hzi = kr(I_n, I_x,I_y, ops_z.HI); + obj.Hz = ops_z.H; + + obj.e_w = kr(I_n, ops_x.e_l, I_y,I_z); + obj.e_e = kr(I_n, ops_x.e_r, I_y,I_z); + obj.e_s = kr(I_n, I_x, ops_y.e_l,I_z); + obj.e_n = kr(I_n, I_x, ops_y.e_r,I_z); + obj.e_b = kr(I_n, I_x, I_y, ops_z.e_l); + obj.e_t = kr(I_n, I_x, I_y, ops_z.e_r); + + obj.m = m; + obj.h = [ops_x.h ops_y.h ops_x.h]; + obj.order = order; + + switch operator + case 'upwind' + alphaA = max(abs(eig(A(params,obj.x(end),obj.y(end),obj.z(end))))); + alphaB = max(abs(eig(B(params,obj.x(end),obj.y(end),obj.z(end))))); + alphaC = max(abs(eig(C(params,obj.x(end),obj.y(end),obj.z(end))))); + + Ap = (obj.Aevaluated+alphaA*I_N)/2; + Am = (obj.Aevaluated-alphaA*I_N)/2; + Dpx = kr(I_n, ops_x.Dp, I_y,I_z); + Dmx = kr(I_n, ops_x.Dm, I_y,I_z); + obj.D = -Am*Dpx; + temp = Ap*Dmx; + obj.D = obj.D-temp; + clear Ap Am Dpx Dmx + + Bp = (obj.Bevaluated+alphaB*I_N)/2; + Bm = (obj.Bevaluated-alphaB*I_N)/2; + Dpy = kr(I_n, I_x, ops_y.Dp,I_z); + Dmy = kr(I_n, I_x, ops_y.Dm,I_z); + temp = Bm*Dpy; + obj.D = obj.D-temp; + temp = Bp*Dmy; + obj.D = obj.D-temp; + clear Bp Bm Dpy Dmy + + + Cp = (obj.Cevaluated+alphaC*I_N)/2; + Cm = (obj.Cevaluated-alphaC*I_N)/2; + Dpz = kr(I_n, I_x, I_y,ops_z.Dp); + Dmz = kr(I_n, I_x, I_y,ops_z.Dm); + + temp = Cm*Dpz; + obj.D = obj.D-temp; + temp = Cp*Dmz; + obj.D = obj.D-temp; + clear Cp Cm Dpz Dmz + obj.D = obj.D-obj.Eevaluated; + + case 'standard' + D1_x = kr(I_n, ops_x.D1, I_y,I_z); + D1_y = kr(I_n, I_x, ops_y.D1,I_z); + D1_z = kr(I_n, I_x, I_y,ops_z.D1); + obj.D = -obj.Aevaluated*D1_x-obj.Bevaluated*D1_y-obj.Cevaluated*D1_z-obj.Eevaluated; + otherwise + error('Opperator not supported'); + end + 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. + function [closure, penalty] = boundary_condition(obj,boundary,type,L) + default_arg('type','char'); + BM = boundary_matrices(obj,boundary); + switch type + case{'c','char'} + [closure,penalty] = boundary_condition_char(obj,BM); + case{'general'} + [closure,penalty] = boundary_condition_general(obj,BM,boundary,L); + otherwise + error('No such boundary condition') + end + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + error('An interface function does not exist yet'); + end + + function N = size(obj) + N = obj.m; + end + + function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y, Z) + params = obj.params; + side = max(length(X),length(Y)); + if isa(mat,'function_handle') + [rows,cols] = size(mat(params,0,0,0)); + matVec = mat(params,X',Y',Z'); + matVec = sparse(matVec); + else + matVec = mat; + [rows,cols] = size(matVec); + side = max(length(X),length(Y)); + cols = cols/side; + end + + ret = cell(rows,cols); + for ii = 1:rows + for jj = 1:cols + ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side)); + end + end + ret = cell2mat(ret); + end + + function [BM] = boundary_matrices(obj,boundary) + params = obj.params; + + switch boundary + case {'w','W','west'} + BM.e_ = obj.e_w; + mat = obj.A; + BM.boundpos = 'l'; + BM.Hi = obj.Hxi; + [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.X(1),obj.Yx,obj.Zx); + BM.side = length(obj.Yx); + case {'e','E','east'} + BM.e_ = obj.e_e; + mat = obj.A; + BM.boundpos = 'r'; + BM.Hi = obj.Hxi; + [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.X(end),obj.Yx,obj.Zx); + BM.side = length(obj.Yx); + case {'s','S','south'} + BM.e_ = obj.e_s; + mat = obj.B; + BM.boundpos = 'l'; + BM.Hi = obj.Hyi; + [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xy,obj.Y(1),obj.Zy); + BM.side = length(obj.Xy); + case {'n','N','north'} + BM.e_ = obj.e_n; + mat = obj.B; + BM.boundpos = 'r'; + BM.Hi = obj.Hyi; + [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xy,obj.Y(end),obj.Zy); + BM.side = length(obj.Xy); + case{'b','B','Bottom'} + BM.e_ = obj.e_b; + mat = obj.C; + BM.boundpos = 'l'; + BM.Hi = obj.Hzi; + [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xz,obj.Yz,obj.Z(1)); + BM.side = length(obj.Xz); + case{'t','T','Top'} + BM.e_ = obj.e_t; + mat = obj.C; + BM.boundpos = 'r'; + BM.Hi = obj.Hzi; + [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xz,obj.Yz,obj.Z(end)); + BM.side = length(obj.Xz); + end + BM.pos = signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3); + end + + % Characteristic bouyndary consitions + function [closure, penalty]=boundary_condition_char(obj,BM) + side = BM.side; + pos = BM.pos; + neg = BM.neg; + zeroval=BM.zeroval; + V = BM.V; + Vi = BM.Vi; + Hi = BM.Hi; + D = BM.D; + e_ = BM.e_; + + switch BM.boundpos + case {'l'} + tau = sparse(obj.n*side,pos); + Vi_plus = Vi(1:pos,:); + tau(1:pos,:) = -abs(D(1:pos,1:pos)); + closure = Hi*e_*V*tau*Vi_plus*e_'; + penalty = -Hi*e_*V*tau*Vi_plus; + case {'r'} + tau = sparse(obj.n*side,neg); + tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); + Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); + closure = Hi*e_*V*tau*Vi_minus*e_'; + penalty = -Hi*e_*V*tau*Vi_minus; + end + end + + % General boundary condition in the form Lu=g(x) + function [closure,penalty] = boundary_condition_general(obj,BM,boundary,L) + side = BM.side; + pos = BM.pos; + neg = BM.neg; + zeroval=BM.zeroval; + V = BM.V; + Vi = BM.Vi; + Hi = BM.Hi; + D = BM.D; + e_ = BM.e_; + + switch boundary + case {'w','W','west'} + L = obj.evaluateCoefficientMatrix(L,obj.x(1),obj.Yx,obj.Zx); + case {'e','E','east'} + L = obj.evaluateCoefficientMatrix(L,obj.x(end),obj.Yx,obj.Zx); + case {'s','S','south'} + L = obj.evaluateCoefficientMatrix(L,obj.Xy,obj.y(1),obj.Zy); + case {'n','N','north'} + L = obj.evaluateCoefficientMatrix(L,obj.Xy,obj.y(end),obj.Zy);% General boundary condition in the form Lu=g(x) + case {'b','B','bottom'} + L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(1)); + case {'t','T','top'} + L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(end)); + end + + switch BM.boundpos + case {'l'} + tau = sparse(obj.n*side,pos); + Vi_plus = Vi(1:pos,:); + Vi_minus = Vi(pos+zeroval+1:obj.n*side,:); + V_plus = V(:,1:pos); + V_minus = V(:,(pos+zeroval)+1:obj.n*side); + + tau(1:pos,:) = -abs(D(1:pos,1:pos)); + R = -inv(L*V_plus)*(L*V_minus); + closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; + penalty = -Hi*e_*V*tau*inv(L*V_plus)*L; + case {'r'} + tau = sparse(obj.n*side,neg); + tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); + Vi_plus = Vi(1:pos,:); + Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); + + V_plus = V(:,1:pos); + V_minus = V(:,(pos+zeroval)+1:obj.n*side); + R = -inv(L*V_minus)*(L*V_plus); + closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_'; + penalty = -Hi*e_*V*tau*inv(L*V_minus)*L; + end + end + + % Function that diagonalizes a symbolic matrix A as A=V*D*Vi + % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign + % [d+ ] + % D = [ d0 ] + % [ d-] + % signVec is a vector specifying the number of possitive, zero and negative eigenvalues of D + function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y,z) + params = obj.params; + syms xs ys zs + [V, D] = eig(mat(params,xs,ys,zs)); + Vi=inv(V); + xs = x; + ys = y; + zs = z; + + + side = max(length(x),length(y)); + Dret = zeros(obj.n,side*obj.n); + Vret = zeros(obj.n,side*obj.n); + Viret= zeros(obj.n,side*obj.n); + + for ii=1:obj.n + for jj=1:obj.n + Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii)); + Vret(jj,(ii-1)*side+1:side*ii) = eval(V(jj,ii)); + Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii)); + end + end + + D = sparse(Dret); + V = sparse(Vret); + Vi = sparse(Viret); + V = obj.evaluateCoefficientMatrix(V,x,y,z); + Vi= obj.evaluateCoefficientMatrix(Vi,x,y,z); + D = obj.evaluateCoefficientMatrix(D,x,y,z); + DD = diag(D); + + poseig = (DD>0); + zeroeig = (DD==0); + negeig = (DD<0); + + D = diag([DD(poseig); DD(zeroeig); DD(negeig)]); + V = [V(:,poseig) V(:,zeroeig) V(:,negeig)]; + Vi= [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)]; + signVec = [sum(poseig),sum(zeroeig),sum(negeig)]; + end + end +end \ No newline at end of file
diff -r f1289f3b86b8 -r 18525f1bb941 +scheme/Hypsyst3dCurve.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Hypsyst3dCurve.m Thu Jan 26 13:07:51 2017 +0000 @@ -0,0 +1,557 @@ +classdef Hypsyst3dCurve < scheme.Scheme + properties + m % Number of points in each direction, possibly a vector + n %size of system + h % Grid spacing + X, Y, Z% Values of x and y for each grid point + Yx, Zx, Xy, Zy, Xz, Yz %Grid values for boundary surfaces + + xi,eta,zeta + Xi, Eta, Zeta + + Eta_xi, Zeta_xi, Xi_eta, Zeta_eta, Xi_zeta, Eta_zeta % Metric terms + X_xi, X_eta, X_zeta,Y_xi,Y_eta,Y_zeta,Z_xi,Z_eta,Z_zeta % Metric terms + + order % Order accuracy for the approximation + + D % non-stabalized scheme operator + Aevaluated, Bevaluated, Cevaluated, Eevaluated % Numeric Coeffiecient matrices + Ahat, Bhat, Chat % Symbolic Transformed Coefficient matrices + A, B, C, E % Symbolic coeffiecient matrices + + J, Ji % JAcobian and inverse Jacobian + + H % Discrete norm + % Norms in the x, y and z directions + Hxii,Hetai,Hzetai, Hzi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. + Hxi,Heta,Hzeta + I_xi,I_eta,I_zeta, I_N,onesN + e_w, e_e, e_s, e_n, e_b, e_t + index_w, index_e,index_s,index_n, index_b, index_t + params %parameters for the coeficient matrice + end + + + methods + function obj = Hypsyst3dCurve(m, order, A, B,C, E, params,ti,operator) + xilim ={0 1}; + etalim = {0 1}; + zetalim = {0 1}; + + if length(m) == 1 + m = [m m m]; + end + m_xi = m(1); + m_eta = m(2); + m_zeta = m(3); + m_tot = m_xi*m_eta*m_zeta; + obj.params = params; + obj.n = length(A(obj,0,0,0)); + + obj.m = m; + obj.order = order; + obj.onesN = ones(obj.n); + + switch operator + case 'upwind' + ops_xi = sbp.D1Upwind(m_xi,xilim,order); + ops_eta = sbp.D1Upwind(m_eta,etalim,order); + ops_zeta = sbp.D1Upwind(m_zeta,zetalim,order); + case 'standard' + ops_xi = sbp.D2Standard(m_xi,xilim,order); + ops_eta = sbp.D2Standard(m_eta,etalim,order); + ops_zeta = sbp.D2Standard(m_zeta,zetalim,order); + otherwise + error('Operator not available') + end + + obj.xi = ops_xi.x; + obj.eta = ops_eta.x; + obj.zeta = ops_zeta.x; + + obj.Xi = kr(obj.xi,ones(m_eta,1),ones(m_zeta,1)); + obj.Eta = kr(ones(m_xi,1),obj.eta,ones(m_zeta,1)); + obj.Zeta = kr(ones(m_xi,1),ones(m_eta,1),obj.zeta); + + + [X,Y,Z] = ti.map(obj.Xi,obj.Eta,obj.Zeta); + obj.X = X; + obj.Y = Y; + obj.Z = Z; + + I_n = eye(obj.n); + I_xi = speye(m_xi); + obj.I_xi = I_xi; + I_eta = speye(m_eta); + obj.I_eta = I_eta; + I_zeta = speye(m_zeta); + obj.I_zeta = I_zeta; + + I_N=kr(I_n,I_xi,I_eta,I_zeta); + + O_xi = ones(m_xi,1); + O_eta = ones(m_eta,1); + O_zeta = ones(m_zeta,1); + + + obj.Hxi = ops_xi.H; + obj.Heta = ops_eta.H; + obj.Hzeta = ops_zeta.H; + obj.h = [ops_xi.h ops_eta.h ops_zeta.h]; + + switch operator + case 'upwind' + D1_xi = kr((ops_xi.Dp+ops_xi.Dm)/2, I_eta,I_zeta); + D1_eta = kr(I_xi, (ops_eta.Dp+ops_eta.Dm)/2,I_zeta); + D1_zeta = kr(I_xi, I_eta,(ops_zeta.Dp+ops_zeta.Dm)/2); + otherwise + D1_xi = kr(ops_xi.D1, I_eta,I_zeta); + D1_eta = kr(I_xi, ops_eta.D1,I_zeta); + D1_zeta = kr(I_xi, I_eta,ops_zeta.D1); + end + + obj.A = A; + obj.B = B; + obj.C = C; + + obj.X_xi = D1_xi*X; + obj.X_eta = D1_eta*X; + obj.X_zeta = D1_zeta*X; + obj.Y_xi = D1_xi*Y; + obj.Y_eta = D1_eta*Y; + obj.Y_zeta = D1_zeta*Y; + obj.Z_xi = D1_xi*Z; + obj.Z_eta = D1_eta*Z; + obj.Z_zeta = D1_zeta*Z; + + obj.Ahat = @transform_coefficient_matrix; + obj.Bhat = @transform_coefficient_matrix; + obj.Chat = @transform_coefficient_matrix; + obj.E = @(obj,x,y,z,~,~,~,~,~,~)E(obj,x,y,z); + + obj.Aevaluated = obj.evaluateCoefficientMatrix(obj.Ahat,obj.X, obj.Y,obj.Z, obj.X_eta,obj.X_zeta,obj.Y_eta,obj.Y_zeta,obj.Z_eta,obj.Z_zeta); + obj.Bevaluated = obj.evaluateCoefficientMatrix(obj.Bhat,obj.X, obj.Y,obj.Z, obj.X_zeta,obj.X_xi,obj.Y_zeta,obj.Y_xi,obj.Z_zeta,obj.Z_xi); + obj.Cevaluated = obj.evaluateCoefficientMatrix(obj.Chat,obj.X,obj.Y,obj.Z, obj.X_xi,obj.X_eta,obj.Y_xi,obj.Y_eta,obj.Z_xi,obj.Z_eta); + + switch operator + case 'upwind' + clear D1_xi D1_eta D1_zeta + alphaA = max(abs(eig(obj.Ahat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_eta(end),obj.X_zeta(end),obj.Y_eta(end),obj.Y_zeta(end),obj.Z_eta(end),obj.Z_zeta(end))))); + alphaB = max(abs(eig(obj.Bhat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_zeta(end),obj.X_xi(end),obj.Y_zeta(end),obj.Y_xi(end),obj.Z_zeta(end),obj.Z_xi(end))))); + alphaC = max(abs(eig(obj.Chat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_xi(end),obj.X_eta(end),obj.Y_xi(end),obj.Y_eta(end),obj.Z_xi(end),obj.Z_eta(end))))); + + Ap = (obj.Aevaluated+alphaA*I_N)/2; + Dmxi = kr(I_n, ops_xi.Dm, I_eta,I_zeta); + diffSum = -Ap*Dmxi; + clear Ap Dmxi + + Am = (obj.Aevaluated-alphaA*I_N)/2; + + obj.Aevaluated = []; + Dpxi = kr(I_n, ops_xi.Dp, I_eta,I_zeta); + temp = Am*Dpxi; + diffSum = diffSum-temp; + clear Am Dpxi + + Bp = (obj.Bevaluated+alphaB*I_N)/2; + Dmeta = kr(I_n, I_xi, ops_eta.Dm,I_zeta); + temp = Bp*Dmeta; + diffSum = diffSum-temp; + clear Bp Dmeta + + Bm = (obj.Bevaluated-alphaB*I_N)/2; + obj.Bevaluated = []; + Dpeta = kr(I_n, I_xi, ops_eta.Dp,I_zeta); + temp = Bm*Dpeta; + diffSum = diffSum-temp; + clear Bm Dpeta + + Cp = (obj.Cevaluated+alphaC*I_N)/2; + Dmzeta = kr(I_n, I_xi, I_eta,ops_zeta.Dm); + temp = Cp*Dmzeta; + diffSum = diffSum-temp; + clear Cp Dmzeta + + Cm = (obj.Cevaluated-alphaC*I_N)/2; + clear I_N + obj.Cevaluated = []; + Dpzeta = kr(I_n, I_xi, I_eta,ops_zeta.Dp); + temp = Cm*Dpzeta; + diffSum = diffSum-temp; + clear Cm Dpzeta temp + + obj.J = obj.X_xi.*obj.Y_eta.*obj.Z_zeta... + +obj.X_zeta.*obj.Y_xi.*obj.Z_eta... + +obj.X_eta.*obj.Y_zeta.*obj.Z_xi... + -obj.X_xi.*obj.Y_zeta.*obj.Z_eta... + -obj.X_eta.*obj.Y_xi.*obj.Z_zeta... + -obj.X_zeta.*obj.Y_eta.*obj.Z_xi; + + obj.Ji = kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot)); + obj.Eevaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,obj.Z,[],[],[],[],[],[]); + + obj.D = obj.Ji*diffSum-obj.Eevaluated; + + case 'standard' + D1_xi = kr(I_n,D1_xi); + D1_eta = kr(I_n,D1_eta); + D1_zeta = kr(I_n,D1_zeta); + + obj.J = obj.X_xi.*obj.Y_eta.*obj.Z_zeta... + +obj.X_zeta.*obj.Y_xi.*obj.Z_eta... + +obj.X_eta.*obj.Y_zeta.*obj.Z_xi... + -obj.X_xi.*obj.Y_zeta.*obj.Z_eta... + -obj.X_eta.*obj.Y_xi.*obj.Z_zeta... + -obj.X_zeta.*obj.Y_eta.*obj.Z_xi; + + obj.Ji = kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot)); + obj.Eevaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,obj.Z,[],[],[],[],[],[]); + + obj.D = obj.Ji*(-obj.Aevaluated*D1_xi-obj.Bevaluated*D1_eta -obj.Cevaluated*D1_zeta)-obj.Eevaluated; + otherwise + error('Operator not supported') + end + + obj.Hxii = kr(I_n, ops_xi.HI, I_eta,I_zeta); + obj.Hetai = kr(I_n, I_xi, ops_eta.HI,I_zeta); + obj.Hzetai = kr(I_n, I_xi,I_eta, ops_zeta.HI); + + obj.index_w = (kr(ops_xi.e_l, O_eta,O_zeta)==1); + obj.index_e = (kr(ops_xi.e_r, O_eta,O_zeta)==1); + obj.index_s = (kr(O_xi, ops_eta.e_l,O_zeta)==1); + obj.index_n = (kr(O_xi, ops_eta.e_r,O_zeta)==1); + obj.index_b = (kr(O_xi, O_eta, ops_zeta.e_l)==1); + obj.index_t = (kr(O_xi, O_eta, ops_zeta.e_r)==1); + + obj.e_w = kr(I_n, ops_xi.e_l, I_eta,I_zeta); + obj.e_e = kr(I_n, ops_xi.e_r, I_eta,I_zeta); + obj.e_s = kr(I_n, I_xi, ops_eta.e_l,I_zeta); + obj.e_n = kr(I_n, I_xi, ops_eta.e_r,I_zeta); + obj.e_b = kr(I_n, I_xi, I_eta, ops_zeta.e_l); + obj.e_t = kr(I_n, I_xi, I_eta, ops_zeta.e_r); + + obj.Eta_xi = kr(obj.eta,ones(m_xi,1)); + obj.Zeta_xi = kr(ones(m_eta,1),obj.zeta); + obj.Xi_eta = kr(obj.xi,ones(m_zeta,1)); + obj.Zeta_eta = kr(ones(m_xi,1),obj.zeta); + obj.Xi_zeta = kr(obj.xi,ones(m_eta,1)); + obj.Eta_zeta = kr(ones(m_zeta,1),obj.eta); + end + + function [ret] = transform_coefficient_matrix(obj,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2) + ret = obj.A(obj,x,y,z).*(y_1.*z_2-z_1.*y_2); + ret = ret+obj.B(obj,x,y,z).*(x_2.*z_1-x_1.*z_2); + ret = ret+obj.C(obj,x,y,z).*(x_1.*y_2-x_2.*y_1); + 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. + function [closure, penalty] = boundary_condition(obj,boundary,type,L) + default_arg('type','char'); + BM = boundary_matrices(obj,boundary); + + switch type + case{'c','char'} + [closure,penalty] = boundary_condition_char(obj,BM); + case{'general'} + [closure,penalty] = boundary_condition_general(obj,BM,boundary,L); + otherwise + error('No such boundary condition') + end + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + error('An interface function does not exist yet'); + end + + function N = size(obj) + N = obj.m; + end + + % Evaluates the symbolic Coeffiecient matrix mat + function [ret] = evaluateCoefficientMatrix(obj,mat, X, Y, Z , x_1 , x_2 , y_1 , y_2 , z_1 , z_2) + params = obj.params; + side = max(length(X),length(Y)); + if isa(mat,'function_handle') + [rows,cols] = size(mat(obj,0,0,0,0,0,0,0,0,0)); + x_1 = kr(obj.onesN,x_1); + x_2 = kr(obj.onesN,x_2); + y_1 = kr(obj.onesN,y_1); + y_2 = kr(obj.onesN,y_2); + z_1 = kr(obj.onesN,z_1); + z_2 = kr(obj.onesN,z_2); + matVec = mat(obj,X',Y',Z',x_1',x_2',y_1',y_2',z_1',z_2'); + matVec = sparse(matVec); + else + matVec = mat; + [rows,cols] = size(matVec); + side = max(length(X),length(Y)); + cols = cols/side; + end + matVec(abs(matVec)<10^(-10)) = 0; + ret = cell(rows,cols); + + for ii = 1:rows + for jj = 1:cols + ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side)); + end + end + ret = cell2mat(ret); + end + + function [BM] = boundary_matrices(obj,boundary) + params = obj.params; + BM.boundary = boundary; + switch boundary + case {'w','W','west'} + BM.e_ = obj.e_w; + mat = obj.Ahat; + BM.boundpos = 'l'; + BM.Hi = obj.Hxii; + BM.index = obj.index_w; + BM.x_1 = obj.X_eta(BM.index); + BM.x_2 = obj.X_zeta(BM.index); + BM.y_1 = obj.Y_eta(BM.index); + BM.y_2 = obj.Y_zeta(BM.index); + BM.z_1 = obj.Z_eta(BM.index); + BM.z_2 = obj.Z_zeta(BM.index); + case {'e','E','east'} + BM.e_ = obj.e_e; + mat = obj.Ahat; + BM.boundpos = 'r'; + BM.Hi = obj.Hxii; + BM.index = obj.index_e; + BM.x_1 = obj.X_eta(BM.index); + BM.x_2 = obj.X_zeta(BM.index); + BM.y_1 = obj.Y_eta(BM.index); + BM.y_2 = obj.Y_zeta(BM.index); + BM.z_1 = obj.Z_eta(BM.index); + BM.z_2 = obj.Z_zeta(BM.index); + case {'s','S','south'} + BM.e_ = obj.e_s; + mat = obj.Bhat; + BM.boundpos = 'l'; + BM.Hi = obj.Hetai; + BM.index = obj.index_s; + BM.x_1 = obj.X_zeta(BM.index); + BM.x_2 = obj.X_xi(BM.index); + BM.y_1 = obj.Y_zeta(BM.index); + BM.y_2 = obj.Y_xi(BM.index); + BM.z_1 = obj.Z_zeta(BM.index); + BM.z_2 = obj.Z_xi(BM.index); + case {'n','N','north'} + BM.e_ = obj.e_n; + mat = obj.Bhat; + BM.boundpos = 'r'; + BM.Hi = obj.Hetai; + BM.index = obj.index_n; + BM.x_1 = obj.X_zeta(BM.index); + BM.x_2 = obj.X_xi(BM.index); + BM.y_1 = obj.Y_zeta(BM.index); + BM.y_2 = obj.Y_xi(BM.index); + BM.z_1 = obj.Z_zeta(BM.index); + BM.z_2 = obj.Z_xi(BM.index); + case{'b','B','Bottom'} + BM.e_ = obj.e_b; + mat = obj.Chat; + BM.boundpos = 'l'; + BM.Hi = obj.Hzetai; + BM.index = obj.index_b; + BM.x_1 = obj.X_xi(BM.index); + BM.x_2 = obj.X_eta(BM.index); + BM.y_1 = obj.Y_xi(BM.index); + BM.y_2 = obj.Y_eta(BM.index); + BM.z_1 = obj.Z_xi(BM.index); + BM.z_2 = obj.Z_eta(BM.index); + case{'t','T','Top'} + BM.e_ = obj.e_t; + mat = obj.Chat; + BM.boundpos = 'r'; + BM.Hi = obj.Hzetai; + BM.index = obj.index_t; + BM.x_1 = obj.X_xi(BM.index); + BM.x_2 = obj.X_eta(BM.index); + BM.y_1 = obj.Y_xi(BM.index); + BM.y_2 = obj.Y_eta(BM.index); + BM.z_1 = obj.Z_xi(BM.index); + BM.z_2 = obj.Z_eta(BM.index); + end + [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.X(BM.index),obj.Y(BM.index),obj.Z(BM.index),... + BM.x_1,BM.x_2,BM.y_1,BM.y_2,BM.z_1,BM.z_2); + BM.side = sum(BM.index); + BM.pos = signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3); + end + + % Characteristic boundary condition + function [closure, penalty] = boundary_condition_char(obj,BM) + side = BM.side; + pos = BM.pos; + neg = BM.neg; + zeroval = BM.zeroval; + V = BM.V; + Vi = BM.Vi; + Hi = BM.Hi; + D = BM.D; + e_ = BM.e_; + + switch BM.boundpos + case {'l'} + tau = sparse(obj.n*side,pos); + Vi_plus = Vi(1:pos,:); + tau(1:pos,:) = -abs(D(1:pos,1:pos)); + closure = Hi*e_*V*tau*Vi_plus*e_'; + penalty = -Hi*e_*V*tau*Vi_plus; + case {'r'} + tau = sparse(obj.n*side,neg); + tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); + Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); + closure = Hi*e_*V*tau*Vi_minus*e_'; + penalty = -Hi*e_*V*tau*Vi_minus; + end + end + + % General boundary condition in the form Lu=g(x) + function [closure,penalty] = boundary_condition_general(obj,BM,boundary,L) + side = BM.side; + pos = BM.pos; + neg = BM.neg; + zeroval = BM.zeroval; + V = BM.V; + Vi = BM.Vi; + Hi = BM.Hi; + D = BM.D; + e_ = BM.e_; + index = BM.index; + + switch BM.boundary + case{'b','B','bottom'} + Ji_vec = diag(obj.Ji); + Ji = diag(Ji_vec(index)); + Zeta_x = Ji*(obj.Y_xi(index).*obj.Z_eta(index)-obj.Z_xi(index).*obj.Y_eta(index)); + Zeta_y = Ji*(obj.X_eta(index).*obj.Z_xi(index)-obj.X_xi(index).*obj.Z_eta(index)); + Zeta_z = Ji*(obj.X_xi(index).*obj.Y_eta(index)-obj.Y_xi(index).*obj.X_eta(index)); + + L = obj.evaluateCoefficientMatrix(L,Zeta_x,Zeta_y,Zeta_z,[],[],[],[],[],[]); + end + + switch BM.boundpos + case {'l'} + tau = sparse(obj.n*side,pos); + Vi_plus = Vi(1:pos,:); + Vi_minus = Vi(pos+zeroval+1:obj.n*side,:); + V_plus = V(:,1:pos); + V_minus = V(:,(pos+zeroval)+1:obj.n*side); + + tau(1:pos,:) = -abs(D(1:pos,1:pos)); + R = -inv(L*V_plus)*(L*V_minus); + closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; + penalty = -Hi*e_*V*tau*inv(L*V_plus)*L; + case {'r'} + tau = sparse(obj.n*side,neg); + tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); + Vi_plus = Vi(1:pos,:); + Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); + + V_plus = V(:,1:pos); + V_minus = V(:,(pos+zeroval)+1:obj.n*side); + R = -inv(L*V_minus)*(L*V_plus); + closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_'; + penalty = -Hi*e_*V*tau*inv(L*V_minus)*L; + end + end + + % Function that diagonalizes a symbolic matrix A as A=V*D*Vi + % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign + % [d+ ] + % D = [ d0 ] + % [ d-] + % signVec is a vector specifying the number of possitive, zero and negative eigenvalues of D + function [V,Vi, D,signVec] = matrixDiag(obj,mat,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2) + params = obj.params; + eps = 10^(-10); + if(sum(abs(x_1))>eps) + syms x_1s + else + x_1s = 0; + end + + if(sum(abs(x_2))>eps) + syms x_2s; + else + x_2s = 0; + end + + + if(sum(abs(y_1))>eps) + syms y_1s + else + y_1s = 0; + end + + if(sum(abs(y_2))>eps) + syms y_2s; + else + y_2s = 0; + end + + if(sum(abs(z_1))>eps) + syms z_1s + else + z_1s = 0; + end + + if(sum(abs(z_2))>eps) + syms z_2s; + else + z_2s = 0; + end + + syms xs ys zs + [V, D] = eig(mat(obj,xs,ys,zs,x_1s,x_2s,y_1s,y_2s,z_1s,z_2s)); + Vi = inv(V); + xs = x; + ys = y; + zs = z; + x_1s = x_1; + x_2s = x_2; + y_1s = y_1; + y_2s = y_2; + z_1s = z_1; + z_2s = z_2; + + side = max(length(x),length(y)); + Dret = zeros(obj.n,side*obj.n); + Vret = zeros(obj.n,side*obj.n); + Viret = zeros(obj.n,side*obj.n); + + for ii=1:obj.n + for jj=1:obj.n + Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii)); + Vret(jj,(ii-1)*side+1:side*ii) = eval(V(jj,ii)); + Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii)); + end + end + + D = sparse(Dret); + V = sparse(Vret); + Vi = sparse(Viret); + V = obj.evaluateCoefficientMatrix(V,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2); + D = obj.evaluateCoefficientMatrix(D,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2); + Vi = obj.evaluateCoefficientMatrix(Vi,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2); + DD = diag(D); + + poseig = (DD>0); + zeroeig = (DD==0); + negeig = (DD<0); + + D = diag([DD(poseig); DD(zeroeig); DD(negeig)]); + V = [V(:,poseig) V(:,zeroeig) V(:,negeig)]; + Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)]; + signVec = [sum(poseig),sum(zeroeig),sum(negeig)]; + end + end +end
diff -r f1289f3b86b8 -r 18525f1bb941 +scheme/Utux.m --- a/+scheme/Utux.m Thu Jan 26 14:01:55 2017 +0100 +++ b/+scheme/Utux.m Thu Jan 26 13:07:51 2017 +0000 @@ -6,7 +6,6 @@ order % Order accuracy for the approximation H % Discrete norm - M % Derivative norm D D1 @@ -18,23 +17,41 @@ methods - function obj = Utux(m,xlim,order) + function obj = Utux(m,xlim,order,operator) default_arg('a',1); - [x, h] = util.get_grid(xlim{:},m); - ops = sbp.Ordinary(m,h,order); + + %Old operators + % [x, h] = util.get_grid(xlim{:},m); + %ops = sbp.Ordinary(m,h,order); + + + switch operator + case 'NonEquidistant' + ops = sbp.D1Nonequidistant(m,xlim,order); + obj.D1 = ops.D1; + case 'Standard' + ops = sbp.D2Standard(m,xlim,order); + obj.D1 = ops.D1; + case 'Upwind' + ops = sbp.D1Upwind(m,xlim,order); + obj.D1 = ops.Dm; + otherwise + error('Unvalid operator') + end + obj.x=ops.x; - obj.D1 = sparse(ops.derivatives.D1); - obj.H = sparse(ops.norms.H); - obj.Hi = sparse(ops.norms.HI); - obj.M = sparse(ops.norms.M); - obj.e_l = sparse(ops.boundary.e_1); - obj.e_r = sparse(ops.boundary.e_m); + + obj.H = ops.H; + obj.Hi = ops.HI; + + obj.e_l = ops.e_l; + obj.e_r = ops.e_r; obj.D=obj.D1; obj.m = m; - obj.h = h; + obj.h = ops.h; obj.order = order; - obj.x = x; + obj.x = ops.x; end % Closure functions return the opertors applied to the own doamin to close the boundary @@ -47,7 +64,7 @@ function [closure, penalty] = boundary_condition(obj,boundary,type,data) default_arg('type','neumann'); default_arg('data',0); - tau = -1*obj.e_l; + tau =-1*obj.e_l; closure = obj.Hi*tau*obj.e_l'; penalty = 0*obj.e_l;