Mercurial > repos > public > sbplib
diff +scheme/Hypsyst3dCurve.m @ 1033:037f203b9bf5 feature/burgers1d
Merge with branch feature/advectioRV to utilize the +rv package
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Thu, 17 Jan 2019 10:44:12 +0100 |
| parents | 706d1c2b4199 |
| children | 0652b34f9f27 |
line wrap: on
line diff
--- a/+scheme/Hypsyst3dCurve.m Fri Oct 12 08:50:25 2018 +0200 +++ b/+scheme/Hypsyst3dCurve.m Thu Jan 17 10:44:12 2019 +0100 @@ -5,22 +5,22 @@ 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 % 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 % Numeric Coeffiecient matrices Ahat, Bhat, Chat % Symbolic Transformed Coefficient matrices A, B, C, E % Symbolic coeffiecient matrices - + J, Ji % JAcobian and inverse Jacobian - + 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. @@ -30,14 +30,14 @@ 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,operator) xilim ={0 1}; etalim = {0 1}; zetalim = {0 1}; - + if length(m) == 1 m = [m m m]; end @@ -47,11 +47,11 @@ 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); - + switch operator case 'upwind' ops_xi = sbp.D1Upwind(m_xi,xilim,order); @@ -64,21 +64,21 @@ otherwise error('Operator not available') end - + 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)); 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); - - + + [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; @@ -86,19 +86,19 @@ obj.I_eta = I_eta; I_zeta = speye(m_zeta); obj.I_zeta = I_zeta; - + 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); - - + + 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); @@ -109,11 +109,11 @@ D1_eta = kr(I_xi, ops_eta.D1,I_zeta); D1_zeta = kr(I_xi, I_eta,ops_zeta.D1); end - + 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; @@ -123,55 +123,55 @@ 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; 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); - + switch operator case 'upwind' clear D1_xi D1_eta D1_zeta alphaA = max(abs(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(abs(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(abs(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; Dmxi = kr(I_n, ops_xi.Dm, I_eta,I_zeta); diffSum = -Ap*Dmxi; clear Ap Dmxi - + Am = (obj.Aevaluated-alphaA*I_N)/2; - + obj.Aevaluated = []; Dpxi = kr(I_n, ops_xi.Dp, I_eta,I_zeta); 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; clear Bp Dmeta - + Bm = (obj.Bevaluated-alphaB*I_N)/2; obj.Bevaluated = []; Dpeta = kr(I_n, I_xi, ops_eta.Dp,I_zeta); 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; clear Cp Dmzeta - + Cm = (obj.Cevaluated-alphaC*I_N)/2; clear I_N obj.Cevaluated = []; @@ -179,72 +179,72 @@ 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... -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.Eevaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,obj.Z,[],[],[],[],[],[]); - + obj.D = obj.Ji*diffSum-obj.Eevaluated; - + 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... +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.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); - + 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.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.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) 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'. @@ -253,7 +253,7 @@ 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); @@ -263,15 +263,15 @@ 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'); + + function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type) + error('Not implemented'); end - + function N = size(obj) 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; @@ -294,7 +294,7 @@ end matVec(abs(matVec)<10^(-10)) = 0; ret = cell(rows,cols); - + for ii = 1:rows for jj = 1:cols ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side)); @@ -302,7 +302,7 @@ end ret = cell2mat(ret); end - + function [BM] = boundary_matrices(obj,boundary) params = obj.params; BM.boundary = boundary; @@ -385,7 +385,7 @@ BM.side = sum(BM.index); BM.pos = signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3); end - + % Characteristic boundary condition function [closure, penalty] = boundary_condition_char(obj,BM) side = BM.side; @@ -397,7 +397,7 @@ Hi = BM.Hi; D = BM.D; e_ = BM.e_; - + switch BM.boundpos case {'l'} tau = sparse(obj.n*side,pos); @@ -413,7 +413,7 @@ penalty = -Hi*e_*V*tau*Vi_minus; end end - + % General boundary condition in the form Lu=g(x) function [closure,penalty] = boundary_condition_general(obj,BM,boundary,L) side = BM.side; @@ -426,7 +426,7 @@ D = BM.D; e_ = BM.e_; index = BM.index; - + switch BM.boundary case{'b','B','bottom'} Ji_vec = diag(obj.Ji); @@ -434,10 +434,10 @@ 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); @@ -445,7 +445,7 @@ 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_'; @@ -455,7 +455,7 @@ 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); @@ -463,7 +463,7 @@ penalty = -Hi*e_*V*tau*inv(L*V_minus)*L; 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+ ] @@ -478,38 +478,38 @@ 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); @@ -522,12 +522,12 @@ 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)); @@ -535,7 +535,7 @@ Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii)); end end - + D = sparse(Dret); V = sparse(Vret); Vi = sparse(Viret); @@ -543,11 +543,11 @@ 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,:)];