Mercurial > repos > public > sbplib
changeset 354:dbac99d2c318 feature/hypsyst
Removed inv(Vi) to save time
| author | Ylva Rydin <ylva.rydin@telia.com> |
|---|---|
| date | Mon, 28 Nov 2016 08:46:28 +0100 |
| parents | fd4f1c80755d |
| children | 69b078cf8072 |
| files | +scheme/Hypsyst3d.m +scheme/Hypsyst3dCurve.m +scheme/Utux.m |
| diffstat | 3 files changed, 441 insertions(+), 350 deletions(-) [+] |
line wrap: on
line diff
diff -r fd4f1c80755d -r dbac99d2c318 +scheme/Hypsyst3d.m --- a/+scheme/Hypsyst3d.m Thu Nov 10 20:49:25 2016 +0100 +++ b/+scheme/Hypsyst3d.m Mon Nov 28 08:46:28 2016 +0100 @@ -14,6 +14,7 @@ H % Discrete norm % Norms in the x, y and z directions + Hx, Hy, Hz 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 @@ -22,8 +23,9 @@ methods - function obj = Hypsyst3d(m, lim, order, A, B,C, E, params) + function obj = Hypsyst3d(m, lim, order, A, B,C, E, params,operator) default_arg('E', []) + default_arg('operatpr',[]) xlim = lim{1}; ylim = lim{2}; zlim = lim{3}; @@ -32,18 +34,25 @@ m = [m m m]; end - obj.A=A; - obj.B=B; - obj.C=C; - obj.E=E; + 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); + 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; @@ -53,14 +62,14 @@ 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.Yx = kr(obj.y,ones(m_z,1)); + obj.Zx = kr(ones(m_y,1),obj.z); - obj.Xz=kr(obj.x,ones(m_y,1)); - obj.Yz=kr(ones(m_z,1),obj.y); + 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); @@ -76,14 +85,14 @@ 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); - 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.Hx = ops_x.H; 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.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); @@ -92,11 +101,37 @@ 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.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; + switch operator + case 'upwind' + alphaA = max(eig(A(params,obj.x(end),obj.y(end),obj.z(end)))); + alphaB = max(eig(B(params,obj.x(end),obj.y(end),obj.z(end)))); + alphaC = max(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; + Bp = (obj.Bevaluated+alphaB*I_N)/2; + Bm = (obj.Bevaluated-alphaB*I_N)/2; + Cp = (obj.Cevaluated+alphaC*I_N)/2; + Cm = (obj.Cevaluated-alphaC*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); + Dpy = kr(I_n, I_x, ops_y.Dp,I_z); + Dmy = kr(I_n, I_x, ops_y.Dm,I_z); + Dpz = kr(I_n, I_x, I_y,ops_z.Dp); + Dmz = kr(I_n, I_x, I_y,ops_z.Dm); + + obj.D=-Am*Dpx-Ap*Dmx-Bm*Dpy-Bp*Dmy-Cm*Dpz-Cp*Dmz-obj.Eevaluated; + otherwise + 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; + end end % Closure functions return the opertors applied to the own doamin to close the boundary @@ -110,9 +145,9 @@ switch type case{'c','char'} - [closure,penalty]=boundary_condition_char(obj,BM); + [closure,penalty] = boundary_condition_char(obj,BM); case{'general'} - [closure,penalty]=boundary_condition_general(obj,BM,boundary,L); + [closure,penalty] = boundary_condition_general(obj,BM,boundary,L); otherwise error('No such boundary condition') end @@ -127,14 +162,14 @@ end function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y, Z) - params=obj.params; - side=max(length(X),length(Y)); + 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'); + [rows,cols] = size(mat(params,0,0,0)); + matVec = mat(params,X',Y',Z'); matVec=sparse(matVec); else - matVec=mat; + matVec = mat; [rows,cols]=size(matVec); side=max(length(X),length(Y)); cols=cols/side; @@ -215,17 +250,17 @@ 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; + 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; + 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 @@ -241,80 +276,85 @@ D=BM.D; e_=BM.e_; switch boundary - case {'w','W','west'} - L=obj.evaluateCoefficientMatrix(L,obj.x(1),obj.Yx,obj.Zx); + 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); + 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); + 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); + 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)); + 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 = 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; + 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,:); + 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; + 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; + params = obj.params; syms xs ys zs - [V, D]=eig(mat(params,xs,ys,zs)); - xs=x; - ys=y; - zs=z; + [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); + 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)); + 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); - V=obj.evaluateCoefficientMatrix(V,x,y,z); - D=obj.evaluateCoefficientMatrix(D,x,y,z); - DD=diag(D); + 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); + 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)]; + 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 fd4f1c80755d -r dbac99d2c318 +scheme/Hypsyst3dCurve.m --- a/+scheme/Hypsyst3dCurve.m Thu Nov 10 20:49:25 2016 +0100 +++ b/+scheme/Hypsyst3dCurve.m Mon Nov 28 08:46:28 2016 +0100 @@ -28,6 +28,7 @@ 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 @@ -36,7 +37,7 @@ methods - function obj = Hypsyst3dCurve(m, order, A, B,C, E, params,ti) + function obj = Hypsyst3dCurve(m, order, A, B,C, E, params,ti,operator) xilim ={0 1}; etalim = {0 1}; zetalim = {0 1}; @@ -46,18 +47,28 @@ end m_xi = m(1); m_eta = m(2); - m_zeta=m(3); - m_tot=m_xi*m_eta*m_zeta; + 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.m = m; + + obj.order = order; + obj.onesN = ones(obj.n); - 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); + 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; @@ -67,40 +78,52 @@ 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.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_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); + 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; + 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_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); + 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); - 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.Hetai = kr(I_n, I_xi, ops_eta.HI,I_zeta); obj.Hzetai = kr(I_n, I_xi,I_eta, ops_zeta.HI); - obj.h=[ops_xi.h ops_eta.h ops_zeta.h]; + 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.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); @@ -109,58 +132,84 @@ 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.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.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); + 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.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.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.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; + switch operator + case 'upwind' + alphaA = max(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(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(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; + Am = (obj.Aevaluated-alphaA*I_N)/2; + Bp = (obj.Bevaluated+alphaB*I_N)/2; + Bm = (obj.Bevaluated-alphaB*I_N)/2; + Cp = (obj.Cevaluated+alphaC*I_N)/2; + Cm = (obj.Cevaluated-alphaC*I_N)/2; + + Dpxi = kr(I_n, ops_xi.Dp, I_eta,I_zeta); + Dmxi = kr(I_n, ops_xi.Dm, I_eta,I_zeta); + Dpeta = kr(I_n, I_xi, ops_eta.Dp,I_zeta); + Dmeta = kr(I_n, I_xi, ops_eta.Dm,I_zeta); + Dpzeta = kr(I_n, I_xi, I_eta,ops_zeta.Dp); + Dmzeta = kr(I_n, I_xi, I_eta,ops_zeta.Dm); + + obj.D = obj.Ji*(-Am*Dpxi-Ap*Dmxi-Bm*Dpeta-Bp*Dmeta-Cm*Dpzeta-Cp*Dmzeta)-obj.Eevaluated; + otherwise + obj.D = obj.Ji*(-obj.Aevaluated*D1_xi-obj.Bevaluated*D1_eta -obj.Cevaluated*D1_zeta)-obj.Eevaluated; + end + + 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); + 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 @@ -171,13 +220,13 @@ % 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); + BM = boundary_matrices(obj,boundary); switch type case{'c','char'} - [closure,penalty]=boundary_condition_char(obj,BM); + [closure,penalty] = boundary_condition_char(obj,BM); case{'general'} - [closure,penalty]=boundary_condition_general(obj,BM,boundary,L); + [closure,penalty] = boundary_condition_general(obj,BM,boundary,L); otherwise error('No such boundary condition') end @@ -192,118 +241,118 @@ 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)); + 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); + [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; + matVec = mat; + [rows,cols] = size(matVec); + side = max(length(X),length(Y)); + cols = cols/side; end - ret=cell(rows,cols); + ret = cell(rows,cols); for ii=1:rows for jj=1:cols - ret{ii,jj}=diag(matVec(ii,(jj-1)*side+1:jj*side)); + ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side)); end end - ret=cell2mat(ret); + ret = cell2mat(ret); end - function [BM]=boundary_matrices(obj,boundary) - params=obj.params; - BM.boundary=boundary; + 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); + 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); + 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); + 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); + 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); + 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); + 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.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); + BM.side = sum(BM.index); + BM.pos = signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3); end @@ -311,26 +360,26 @@ side = BM.side; pos = BM.pos; neg = BM.neg; - zeroval=BM.zeroval; + zeroval = BM.zeroval; V = BM.V; Vi = BM.Vi; - Hi=BM.Hi; - D=BM.D; - e_=BM.e_; + 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; + 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; + 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 @@ -339,135 +388,135 @@ side = BM.side; pos = BM.pos; neg = BM.neg; - zeroval=BM.zeroval; + zeroval = BM.zeroval; V = BM.V; Vi = BM.Vi; - Hi=BM.Hi; - D=BM.D; - e_=BM.e_; - index=BM.index; + 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)); + 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,[],[],[],[],[],[]); + 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 = 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; + 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,:); + 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; + 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); + 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; + x_1s = 0; end if(sum(abs(x_2))>eps) - syms x_2s; + syms x_2s; else - x_2s=0; + x_2s = 0; end - - + + if(sum(abs(y_1))>eps) syms y_1s else - y_1s=0; + y_1s = 0; end if(sum(abs(y_2))>eps) - syms y_2s; + syms y_2s; else - y_2s=0; + y_2s = 0; end - + if(sum(abs(z_1))>eps) syms z_1s else - z_1s=0; + z_1s = 0; end if(sum(abs(z_2))>eps) - syms z_2s; + syms z_2s; else - z_2s=0; - end + 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); + 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)); + 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); + + 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); + 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)]; + 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 +end
diff -r fd4f1c80755d -r dbac99d2c318 +scheme/Utux.m --- a/+scheme/Utux.m Thu Nov 10 20:49:25 2016 +0100 +++ b/+scheme/Utux.m Mon Nov 28 08:46:28 2016 +0100 @@ -24,12 +24,14 @@ % [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); + % ops = sbp.D1Nonequidistant(m,xlim,order); + ops = sbp.D2Standard(m,xlim,order); + obj.D1 = ops.D1; +% ops = sbp.D1Upwind(m,xlim,order); +% obj.D1 = ops.Dm; obj.x=ops.x; - obj.D1 = ops.D1; + obj.H = ops.H; obj.Hi = ops.HI;