sbplib: +scheme/Hypsyst3dCurve.m comparison

comparison +scheme/Hypsyst3dCurve.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	feebfca90080 459eeb99130f

comparison

equal deleted inserted replaced

-:53abf04f5e4e
+:9d1fc984f40d
 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
+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
-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
+Aevaluated, Bevaluated, Cevaluated, Eevaluated % Numeric Coeffiecient matrices
-Ahat, Bhat, Chat, E
+Ahat, Bhat, Chat  % Symbolic Transformed Coefficient matrices
-A,B,C
+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
 Hxii,Hetai,Hzetai, Hzi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
 Hxi,Heta,Hzeta
 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));%% 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);
 [X,Y,Z] = ti.map(obj.Xi,obj.Eta,obj.Zeta);
 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.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);
 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;
+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;
+temp = Am*Dpxi;
-diffSum=diffSum-temp;
+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;
+temp = Bp*Dmeta;
-diffSum=diffSum-temp;
+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;
+temp = Bm*Dpeta;
-diffSum=diffSum-temp;
+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;
+temp = Cp*Dmzeta;
-diffSum=diffSum-temp;
+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;
+temp = Cm*Dpzeta;
-diffSum=diffSum-temp;
+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.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;
-otherwise
-D1_xi=kr(I_n,D1_xi);
+case 'standard'
-D1_eta=kr(I_n,D1_eta);
+D1_xi = kr(I_n,D1_xi);
-D1_zeta=kr(I_n,D1_zeta);
+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.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;
-end
+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.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);
 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);
 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;
 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));
 matVec = mat;
 [rows,cols] = size(matVec);
 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 ii=1:rows
+for jj = 1:cols
-for jj=1:cols
 ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side));
 end
 end
 ret = cell2mat(ret);
+end
-end
 function [BM] = boundary_matrices(obj,boundary)
 params = obj.params;
 BM.boundary = boundary;
 switch boundary
 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);
 end
+% Characteristic boundary condition
-function [closure, penalty]=boundary_condition_char(obj,BM)
+function [closure, penalty] = boundary_condition_char(obj,BM)
 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_;
 switch BM.boundpos
 case {'l'}
 tau = sparse(obj.n*side,pos);
 Vi_plus = Vi(1:pos,:);
 closure = Hi*e_*V*tau*Vi_minus*e_';
 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)
+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;
 closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
 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+       ]
+%                               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);
 if(sum(abs(x_1))>eps)
 syms x_1s
 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;

Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst3dCurve.m @ 369:9d1fc984f40d feature/hypsyst