Mercurial > repos > public > sbplib
changeset 353:fd4f1c80755d feature/hypsyst
meged with default
| author | Ylva Rydin <ylva.rydin@telia.com> |
|---|---|
| date | Thu, 10 Nov 2016 20:49:25 +0100 |
| parents | 9b3d7fc61a36 (diff) f18142c1530b (current diff) |
| children | dbac99d2c318 |
| files | |
| diffstat | 6 files changed, 1715 insertions(+), 12 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/Ti3D.m Thu Nov 10 20:49:25 2016 +0100 @@ -0,0 +1,252 @@ +classdef Ti3D + properties + gs % {6}Surfaces + V % FunctionHandle(XI,ETA,ZETA) + end + + methods + % TODO function to label boundary names. + % function to find largest and smallest delta h in the grid. Maybe shouldnt live here + 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Hypsyst2d.m Thu Nov 10 20:49:25 2016 +0100 @@ -0,0 +1,290 @@ +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 + + 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 + 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 + + + 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 + + + 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 [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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Hypsyst2dCurve.m Thu Nov 10 20:49:25 2016 +0100 @@ -0,0 +1,362 @@ +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 + 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_eta+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, ops_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','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,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 + + + 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 + + + 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 [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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Hypsyst3d.m Thu Nov 10 20:49:25 2016 +0100 @@ -0,0 +1,320 @@ +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 + Aevaluated,Bevaluated,Cevaluated, Eevaluated + + H % Discrete norm + % 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 + function obj = Hypsyst3d(m, lim, order, A, B,C, E, params) + 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; + + ops_x = sbp.D2Standard(m_x,xlim,order); + ops_y = sbp.D2Standard(m_y,ylim,order); + ops_z = sbp.D2Standard(m_z,zlim,order); + + 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));%% Que pasa? + 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 = eye(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; + + + D1_x = kr(I_n, ops_x.D1, I_y,I_z); + obj.Hxi = kr(I_n, ops_x.HI, I_y,I_z); + D1_y = kr(I_n, I_x, ops_y.D1,I_z); + obj.Hyi = kr(I_n, I_x, ops_y.HI,I_z); + D1_z = kr(I_n, I_x, I_y,ops_z.D1); + obj.Hzi = kr(I_n, I_x,I_y, ops_z.HI); + + 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; + + obj.D=-obj.Aevaluated*D1_x-obj.Bevaluated*D1_y-obj.Cevaluated*D1_z-obj.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'); + 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 + + + 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 + + + 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); + 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 [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)); + 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); + 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)); + end + end + + D=sparse(Dret); + V=sparse(Vret); + V=obj.evaluateCoefficientMatrix(V,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=inv(V); + signVec=[sum(poseig),sum(zeroeig),sum(negeig)]; + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Hypsyst3dCurve.m Thu Nov 10 20:49:25 2016 +0100 @@ -0,0 +1,473 @@ +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 + + X_xi, X_eta, X_zeta,Y_xi,Y_eta,Y_zeta,Z_xi,Z_eta,Z_zeta + Aev + + metric_terms + + order % Order accuracy for the approximation + + D % non-stabalized scheme operator + Aevaluated, Bevaluated, Cevaluated, Eevaluated + Ahat, Bhat, Chat, E + A,B,C + + J, Ji + + 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. + 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) + 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); + ops_xi = sbp.D2Standard(m_xi,xilim,order); + ops_eta = sbp.D2Standard(m_eta,etalim,order); + ops_zeta = sbp.D2Standard(m_zeta,zetalim,order); + + 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));%% Que pasa? + 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); + + 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); + + [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; + + + O_xi=ones(m_xi,1); + O_eta=ones(m_eta,1); + O_zeta=ones(m_zeta,1); + + D1_xi = kr(ops_xi.D1, I_eta,I_zeta); + obj.Hxii = kr(I_n, ops_xi.HI, I_eta,I_zeta); + D1_eta = kr(I_xi, ops_eta.D1,I_zeta); + obj.Hetai = kr(I_n, I_xi, ops_eta.HI,I_zeta); + D1_zeta = kr(I_xi, I_eta,ops_zeta.D1); + obj.Hzetai = kr(I_n, I_xi,I_eta, ops_zeta.HI); + obj.h=[ops_xi.h ops_eta.h ops_zeta.h]; + + 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.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; + + D1_xi=kr(I_n,D1_xi); + D1_eta=kr(I_n,D1_eta); + D1_zeta=kr(I_n,D1_zeta); + + 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.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); + obj.Eevaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,obj.Z,[],[],[],[],[],[]); + + 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.D=obj.Ji*(-obj.Aevaluated*D1_xi-obj.Bevaluated*D1_eta -obj.Cevaluated*D1_zeta)-obj.Eevaluated; + 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 + + 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 + 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 + + + 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 + + + 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 [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); + % syms x_1s x_2s y_1s y_2s z_1s z_2s + 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
--- a/+scheme/Utux.m Thu Oct 13 09:34:42 2016 +0200 +++ b/+scheme/Utux.m Thu Nov 10 20:49:25 2016 +0100 @@ -6,7 +6,6 @@ order % Order accuracy for the approximation H % Discrete norm - M % Derivative norm D D1 @@ -20,21 +19,28 @@ methods function obj = Utux(m,xlim,order) 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); + + % ops = sbp.D1Nonequidistant(m,xlim,order); + % ops = sbp.D2Standard(m,xlim,order); + ops = sbp.D1Upwind(m,xlim,order); + 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.D1 = ops.D1; + 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 +53,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;