Mercurial > repos > public > sbplib
changeset 1248:10881b234f77 feature/volcano
Clean up hypsyst3d; use functionality from grid module and remove obsolete properties. Store operators for different coord dirs in cell arrays.
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Wed, 20 Nov 2019 14:12:55 -0800 |
| parents | e1f9bedd64a9 |
| children | 6424745e1b58 |
| files | +scheme/Hypsyst3d.m |
| diffstat | 1 files changed, 44 insertions(+), 82 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/Hypsyst3d.m Mon Nov 18 17:17:16 2019 -0800 +++ b/+scheme/Hypsyst3d.m Wed Nov 20 14:12:55 2019 -0800 @@ -1,7 +1,6 @@ classdef Hypsyst3d < scheme.Scheme properties grid - Yx, Zx, Xy, Zy, Xz, Yz %Grid values for boundary surfaces order % Order accuracy for the approximation nEquations @@ -31,22 +30,14 @@ obj.C = C; obj.E = E; obj.params = params; - dim = obj.grid.d; - assert(dim == 3, 'Dimensions not correct') - points = obj.grid.points(); - m = obj.grid.m; - for i = 1:dim - ops{i} = opSet(m(i),obj.grid.lim{i},order); - x{i} = obj.grid.x{i}; - X{i} = points(:,i); + assert(obj.grid.d == 3, 'Dimensions not correct') + % Construct 1D operators for each coordinate dir + for i = 1:obj.grid.d + ops{i} = opSet(obj.grid.m(i),obj.grid.lim{i},order); + I{i} = speye(obj.grid.m(i)); end - - obj.Yx = kr(x{2}, ones(m(3),1)); - obj.Zx = kr(ones(m(2),1), x{3}); - obj.Xy = kr(x{1}, ones(m(3),1)); - obj.Zy = kr(ones(m(1),1), x{3}); - obj.Xz = kr(x{1}, ones(m(2),1)); - obj.Yz = kr(ones(m(3),1), x{2}); + x = obj.grid.x; %Grid points in each coordinate dir + X = num2cell(obj.grid.points(),1); %Kroneckered grid values in each coordinate dir obj.Aevaluated = obj.evaluateCoefficientMatrix(A, X{:}); obj.Bevaluated = obj.evaluateCoefficientMatrix(B, X{:}); @@ -54,26 +45,22 @@ obj.Eevaluated = obj.evaluateCoefficientMatrix(E, X{:}); obj.nEquations = length(A(obj.params,0,0,0)); + I_n = speye(obj.nEquations); + I_N = kr(I_n,I{:}); - I_n = speye(obj.nEquations); - I_x = speye(m(1)); - I_y = speye(m(2)); - I_z = speye(m(3)); - I_N = kr(I_n,I_x,I_y,I_z); - - obj.Hxi = kr(I_n, ops{1}.HI, I_y,I_z); + obj.Hxi = kr(I_n, ops{1}.HI, I{2}, I{3}); obj.Hx = ops{1}.H; - obj.Hyi = kr(I_n, I_x, ops{2}.HI,I_z); + obj.Hyi = kr(I_n, I{1}, ops{2}.HI, I{3}); obj.Hy = ops{2}.H; - obj.Hzi = kr(I_n, I_x,I_y, ops{3}.HI); + obj.Hzi = kr(I_n, I{1},I{2}, ops{3}.HI); obj.Hz = ops{3}.H; - obj.e_w = kr(I_n, ops{1}.e_l, I_y,I_z); - obj.e_e = kr(I_n, ops{1}.e_r, I_y,I_z); - obj.e_s = kr(I_n, I_x, ops{2}.e_l,I_z); - obj.e_n = kr(I_n, I_x, ops{2}.e_r,I_z); - obj.e_b = kr(I_n, I_x, I_y, ops{3}.e_l); - obj.e_t = kr(I_n, I_x, I_y, ops{3}.e_r); + obj.e_w = kr(I_n, ops{1}.e_l, I{2}, I{3}); + obj.e_e = kr(I_n, ops{1}.e_r, I{2}, I{3}); + obj.e_s = kr(I_n, I{1}, ops{2}.e_l, I{3}); + obj.e_n = kr(I_n, I{1}, ops{2}.e_r, I{3}); + obj.e_b = kr(I_n, I{1}, I{2}, ops{3}.e_l); + obj.e_t = kr(I_n, I{1}, I{2}, ops{3}.e_r); obj.order = order; @@ -85,8 +72,8 @@ Ap = (obj.Aevaluated+alphaA*I_N)/2; Am = (obj.Aevaluated-alphaA*I_N)/2; - Dpx = kr(I_n, ops{1}.Dp, I_y,I_z); - Dmx = kr(I_n, ops{1}.Dm, I_y,I_z); + Dpx = kr(I_n, ops{1}.Dp, I{2}, I{3}); + Dmx = kr(I_n, ops{1}.Dm, I{2}, I{3}); obj.D = -Am*Dpx; temp = Ap*Dmx; obj.D = obj.D-temp; @@ -94,8 +81,8 @@ Bp = (obj.Bevaluated+alphaB*I_N)/2; Bm = (obj.Bevaluated-alphaB*I_N)/2; - Dpy = kr(I_n, I_x, ops{2}.Dp,I_z); - Dmy = kr(I_n, I_x, ops{2}.Dm,I_z); + Dpy = kr(I_n, I{1}, ops{2}.Dp, I{3}); + Dmy = kr(I_n, I{1}, ops{2}.Dm, I{3}); temp = Bm*Dpy; obj.D = obj.D-temp; temp = Bp*Dmy; @@ -105,8 +92,8 @@ Cp = (obj.Cevaluated+alphaC*I_N)/2; Cm = (obj.Cevaluated-alphaC*I_N)/2; - Dpz = kr(I_n, I_x, I_y,ops{3}.Dp); - Dmz = kr(I_n, I_x, I_y,ops{3}.Dm); + Dpz = kr(I_n, I{1}, I{2}, ops{3}.Dp); + Dmz = kr(I_n, I{1}, I{2}, ops{3}.Dm); temp = Cm*Dpz; obj.D = obj.D-temp; @@ -116,9 +103,9 @@ obj.D = obj.D-obj.Eevaluated; case {'sbp.D2Standard', 'sbp.D2Variable', 'sbp.D4Standard', 'sbp.D4Variable'} - D1_x = kr(I_n, ops{1}.D1, I_y,I_z); - D1_y = kr(I_n, I_x, ops{2}.D1,I_z); - D1_z = kr(I_n, I_x, I_y,ops{3}.D1); + D1_x = kr(I_n, ops{1}.D1, I{2}, I{3}); + D1_y = kr(I_n, I{1}, ops{2}.D1, I{3}); + D1_z = kr(I_n, I{1}, I{2}, ops{3}.D1); obj.D = -obj.Aevaluated*D1_x-obj.Bevaluated*D1_y-obj.Cevaluated*D1_z-obj.Eevaluated; otherwise error('Operator not supported'); @@ -193,55 +180,46 @@ function [BM] = boundary_matrices(obj,boundary) params = obj.params; - points = obj.grid.points(); - X = points(:, 1); - Y = points(:, 2); - Z = points(:, 3); - + bp = num2cell(obj.grid.getBoundary(boundary),1); %Kroneckered boundary points switch boundary - case {'w','W','west'} + case {'w'} 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,X(1),obj.Yx,obj.Zx); - BM.side = length(obj.Yx); - case {'e','E','east'} + BM.side = length(bp{2}); + case {'e'} 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,X(end),obj.Yx,obj.Zx); - BM.side = length(obj.Yx); - case {'s','S','south'} + BM.side = length(bp{2}); + case {'s'} 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,Y(1),obj.Zy); - BM.side = length(obj.Xy); - case {'n','N','north'} + BM.side = length(bp{1}); + case {'n'} 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,Y(end),obj.Zy); - BM.side = length(obj.Xy); - case{'b','B','Bottom'} + BM.side = length(bp{1}); + case{'d'} 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,Z(1)); - BM.side = length(obj.Xz); - case{'t','T','Top'} + BM.side = length(bp{1}); + case{'u'} 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,Z(end)); - BM.side = length(obj.Xz); + BM.side = length(bp{1}); end + [BM.V,BM.Vi,BM.D,signVec] = obj.diagonalizeMatrix(mat,bp{:}); BM.pos = signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3); end @@ -285,24 +263,8 @@ D = BM.D; e_ = BM.e_; - x = obj.grid.x{1}; - y = obj.grid.x{2}; - z = obj.grid.x{3}; - - switch boundary - case {'w','W','west'} - L = obj.evaluateCoefficientMatrix(L,x(1),obj.Yx,obj.Zx); - case {'e','E','east'} - L = obj.evaluateCoefficientMatrix(L,x(end),obj.Yx,obj.Zx); - case {'s','S','south'} - L = obj.evaluateCoefficientMatrix(L,obj.Xy,y(1),obj.Zy); - case {'n','N','north'} - L = obj.evaluateCoefficientMatrix(L,obj.Xy,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, z(1)); - case {'t','T','top'} - L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz, z(end)); - end + bp = num2cell(obj.grid.getBoundary(boundary),1); %Kroneckered boundary points + L = obj.evaluateCoefficientMatrix(L,bp{:}); % General boundary condition in the form Lu=g(x) switch BM.boundpos case {'l'} @@ -336,7 +298,7 @@ % D = [ d0 ] % [ d-] % signVec is a vector specifying the number of positive, zero and negative eigenvalues of D - function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y,z) + function [V, Vi, D, signVec] = diagonalizeMatrix(obj, mat, x, y, z) params = obj.params; syms xs ys zs [V, D] = eig(mat(params,xs,ys,zs));