Mercurial > repos > public > sbplib
diff +scheme/Hypsyst3d.m @ 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 | 0fd6561964b0 |
line wrap: on
line diff
--- 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));