Mercurial > repos > public > sbplib
changeset 369:9d1fc984f40d feature/hypsyst
Added some comments and cleaned up the code a little
| author | Ylva Rydin <ylva.rydin@telia.com> |
|---|---|
| date | Thu, 26 Jan 2017 09:57:24 +0100 |
| parents | 53abf04f5e4e |
| children | c2c64ccb6a1e |
| files | +scheme/Hypsyst2d.m +scheme/Hypsyst2dCurve.m +scheme/Hypsyst3d.m +scheme/Hypsyst3dCurve.m +scheme/Wave2dCurve.m |
| diffstat | 5 files changed, 575 insertions(+), 549 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/Hypsyst2d.m Wed Jan 25 15:37:12 2017 +0100 +++ b/+scheme/Hypsyst2d.m Thu Jan 26 09:57:24 2017 +0100 @@ -6,10 +6,10 @@ 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 - + 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. @@ -17,14 +17,14 @@ 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 @@ -32,50 +32,50 @@ 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); - + 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; - + + 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'. @@ -85,206 +85,217 @@ default_arg('type','char'); switch type case{'c','char'} - [closure,penalty]=boundary_condition_char(obj,boundary); + [closure,penalty] = boundary_condition_char(obj,boundary); case{'general'} - [closure,penalty]=boundary_condition_general(obj,boundary,L); + [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; - + 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)); + [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; + matVec = mat; + [rows,cols] = size(matVec); + side = max(length(X),length(Y)); + cols = cols/side; end - ret=cell(rows,cols); - - for ii=1:rows + 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 [closure, penalty]=boundary_condition_char(obj,boundary) - params=obj.params; - x=obj.x; y=obj.y; - + + %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)); + 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)); + 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)); + 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)); + 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); - + 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; + 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 - - - function [closure,penalty]=boundary_condition_general(obj,boundary,L) - params=obj.params; - x=obj.x; y=obj.y; - + + % 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)); + 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)); + 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)); + 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)); + 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); - + + 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; + 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; + 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; + + % 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)); + [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); + + 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); + 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 \ No newline at end of file
--- a/+scheme/Hypsyst2dCurve.m Wed Jan 25 15:37:12 2017 +0100 +++ b/+scheme/Hypsyst2dCurve.m Thu Jan 26 09:57:24 2017 +0100 @@ -1,362 +1,378 @@ classdef Hypsyst2dCurve < scheme.Scheme properties m % Number of points in each direction, possibly a vector - n %size of system + n % size of system h % Grid spacing X,Y % Values of x and y for each grid point - J, Ji %Jacobaian and inverse Jacobian + J, Ji % Jacobaian and inverse Jacobian xi,eta Xi,Eta A,B - X_eta, Y_eta + 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. + 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 + 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.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_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, ops_eta.e_r); + 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); - + + 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; - + + 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'. + % 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); + [closure,penalty] = boundary_condition_char(obj,boundary); case{'general'} - [closure,penalty]=boundary_condition_general(obj,boundary,L); + [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) + + 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; - + 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)); + [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; + 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)); + + 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); + 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; - - + + %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)); + 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)); + 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)); + 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)); + 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); - + + 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; + 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 + 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; - + + + % 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)); + 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)); + 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)); + 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)); + 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)); + 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)); + 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); - + + 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; + 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; + 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) + + % 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; + xs_ = 0; end - if(sum(abs(y_))~=0) - syms ys_; + if(sum(abs(y_))~= 0) + syms ys_; else - ys_=0; + ys_ = 0; end - [V, D]=eig(mat(params,xs,ys,xs_,ys_)); - Vi=inv(V); + [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)); + 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); + + 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); + 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 \ No newline at end of file
--- a/+scheme/Hypsyst3d.m Wed Jan 25 15:37:12 2017 +0100 +++ b/+scheme/Hypsyst3d.m Thu Jan 26 09:57:24 2017 +0100 @@ -1,7 +1,7 @@ classdef Hypsyst3d < scheme.Scheme properties m % Number of points in each direction, possibly a vector - n %size of system + n % Size of system h % Grid spacing x, y, z % Grid X, Y, Z% Values of x and y for each grid point @@ -9,20 +9,20 @@ order % Order accuracy for the approximation D % non-stabalized scheme operator - A, B, C, E + A, B, C, E % Symbolic coefficient matrices Aevaluated,Bevaluated,Cevaluated, Eevaluated H % Discrete norm - % Norms in the x, y and z directions - Hx, Hy, Hz + 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 + 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', []) default_arg('operatpr',[]) @@ -40,7 +40,7 @@ obj.E = E; m_x = m(1); m_y = m(2); - m_z=m(3); + m_z = m(3); obj.params = params; switch operator @@ -58,16 +58,14 @@ 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.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); @@ -85,7 +83,7 @@ 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); + 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; @@ -115,40 +113,41 @@ 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; + 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; + 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; + 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; + 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; - obj.D=obj.D-obj.Eevaluated; - - otherwise + 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; + obj.D = -obj.Aevaluated*D1_x-obj.Bevaluated*D1_y-obj.Cevaluated*D1_z-obj.Eevaluated; + otherwise + error('Opperator not supported'); end end @@ -159,8 +158,7 @@ % 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); @@ -185,76 +183,74 @@ if isa(mat,'function_handle') [rows,cols] = size(mat(params,0,0,0)); matVec = mat(params,X',Y',Z'); - matVec=sparse(matVec); + matVec = sparse(matVec); else matVec = mat; - [rows,cols]=size(matVec); - side=max(length(X),length(Y)); - cols=cols/side; + [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)); + 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); + ret = cell2mat(ret); end - - function [BM]=boundary_matrices(obj,boundary) - params=obj.params; + 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); + 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); + 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); + 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); + 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); + 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); + 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); + 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; @@ -262,9 +258,9 @@ 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'} @@ -282,17 +278,18 @@ end end - - function [closure,penalty]=boundary_condition_general(obj,BM,boundary,L) + % 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_; + 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); @@ -301,7 +298,7 @@ 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); + 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'} @@ -334,7 +331,12 @@ 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 @@ -349,6 +351,7 @@ 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));
--- a/+scheme/Hypsyst3dCurve.m Wed Jan 25 15:37:12 2017 +0100 +++ b/+scheme/Hypsyst3dCurve.m Thu Jan 26 09:57:24 2017 +0100 @@ -9,21 +9,17 @@ 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 - - - metric_terms + 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 - Ahat, Bhat, Chat, E - A,B,C + Aevaluated, Bevaluated, Cevaluated, Eevaluated % Numeric Coeffiecient matrices + Ahat, Bhat, Chat % Symbolic Transformed Coefficient matrices + A, B, C, E % Symbolic coeffiecient matrices - J, Ji + J, Ji % JAcobian and inverse Jacobian H % Discrete norm % Norms in the x, y and z directions @@ -73,7 +69,7 @@ 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.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); @@ -127,7 +123,7 @@ 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; @@ -146,45 +142,44 @@ Ap = (obj.Aevaluated+alphaA*I_N)/2; Dmxi = kr(I_n, ops_xi.Dm, I_eta,I_zeta); - diffSum=-Ap*Dmxi; + diffSum = -Ap*Dmxi; clear Ap Dmxi Am = (obj.Aevaluated-alphaA*I_N)/2; - obj.Aevaluated=[]; + + obj.Aevaluated = []; Dpxi = kr(I_n, ops_xi.Dp, I_eta,I_zeta); - temp=Am*Dpxi; - diffSum=diffSum-temp; + 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; + temp = Bp*Dmeta; + diffSum = diffSum-temp; clear Bp Dmeta Bm = (obj.Bevaluated-alphaB*I_N)/2; - obj.Bevaluated=[]; + obj.Bevaluated = []; Dpeta = kr(I_n, I_xi, ops_eta.Dp,I_zeta); - temp=Bm*Dpeta; - diffSum=diffSum-temp; + 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; + temp = Cp*Dmzeta; + diffSum = diffSum-temp; clear Cp Dmzeta Cm = (obj.Cevaluated-alphaC*I_N)/2; clear I_N - obj.Cevaluated=[]; + obj.Cevaluated = []; Dpzeta = kr(I_n, I_xi, I_eta,ops_zeta.Dp); - temp=Cm*Dpzeta; - diffSum=diffSum-temp; + 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... @@ -196,10 +191,11 @@ obj.Eevaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,obj.Z,[],[],[],[],[],[]); obj.D = obj.Ji*diffSum-obj.Eevaluated; - otherwise - D1_xi=kr(I_n,D1_xi); - D1_eta=kr(I_n,D1_eta); - D1_zeta=kr(I_n,D1_zeta); + + 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... @@ -212,7 +208,10 @@ 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); @@ -231,15 +230,12 @@ 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); - + 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) @@ -276,6 +272,7 @@ 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)); @@ -295,21 +292,17 @@ side = max(length(X),length(Y)); cols = cols/side; end - matVec(abs(matVec)<10^(-10))=0; + matVec(abs(matVec)<10^(-10)) = 0; ret = cell(rows,cols); - - for ii=1:rows - for jj=1: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; @@ -393,8 +386,8 @@ BM.pos = signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3); end - - function [closure, penalty]=boundary_condition_char(obj,BM) + % Characteristic boundary condition + function [closure, penalty] = boundary_condition_char(obj,BM) side = BM.side; pos = BM.pos; neg = BM.neg; @@ -405,7 +398,6 @@ D = BM.D; e_ = BM.e_; - switch BM.boundpos case {'l'} tau = sparse(obj.n*side,pos); @@ -422,8 +414,8 @@ end end - - function [closure,penalty]=boundary_condition_general(obj,BM,boundary,L) + % 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; @@ -472,7 +464,12 @@ 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); @@ -516,7 +513,6 @@ 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;
--- a/+scheme/Wave2dCurve.m Wed Jan 25 15:37:12 2017 +0100 +++ b/+scheme/Wave2dCurve.m Thu Jan 26 09:57:24 2017 +0100 @@ -131,7 +131,7 @@ obj.du_s = (obj.e_s'*Du)'; obj.dv_s = kr(I_u,d1_l_v); obj.du_n = (obj.e_n'*Du)'; - obj.dv_n = kr(I_u,d1_r_v); + obj.dv_n = kr(I_u,d1_r_v);General boundary conditions obj.m = m; obj.h = [h_u h_v];