Mercurial > repos > public > sbplib
changeset 350:5d5652fe826a feature/hypsyst
A commit before I try resolving the performance issues
| author | Ylva Rydin <ylva.rydin@telia.com> |
|---|---|
| date | Wed, 02 Nov 2016 00:02:01 +0100 |
| parents | cd6a29ab3746 |
| children | 7cc3d5bd3692 |
| files | +grid/Ti3D.m +scheme/Hypsyst2dCurve.m +scheme/Hypsyst3d.m +scheme/Hypsyst3dCurve.m |
| diffstat | 4 files changed, 516 insertions(+), 35 deletions(-) [+] |
line wrap: on
line diff
--- a/+grid/Ti3D.m Thu Oct 13 09:34:30 2016 +0200 +++ b/+grid/Ti3D.m Wed Nov 02 00:02:01 2016 +0100 @@ -25,34 +25,34 @@ one=0*ETA+1; zero=0*ETA; - Sw = gw((1-ETA),(1-ZETA)); - Se = ge(ETA,ZETA); - Ss = gs(XI,(1-ZETA)); - Sn = gn((1-XI),ZETA); - Sb = gb(XI,ETA); - St = gt((1-XI),(1-ETA)); + Sw = gw(ETA,(1-ZETA)); + Se = ge((1-ETA),(1-ZETA)); + Ss = gs(XI,ZETA); + Sn = gn((1-XI),(1-ZETA)); + Sb = gb((1-XI),ETA); + St = gt(XI,ETA); - Ewt = gw(1-ETA,zero); - Ewb = gw(1-ETA,one); - Ews = gw(one,1-ZETA); - Ewn = gw(zero,1-ZETA); - Eet = ge(ETA,one); - Eeb = ge(ETA,zero); - Ees = ge(0*one,ZETA); - Een = ge(one,ZETA); - Enb = gn(1-XI,zero); - Ent = gn(1-XI,one); - Est = gs(XI,zero); - Esb = gs(XI,one); + Ewt = gw(ETA,zero); + Ewb = gw(ETA,one); + Ews = gw(zero,1-ZETA); + Ewn = gw(one,1-ZETA); + Eet = ge(1-ETA,zero); + Eeb = ge(1-ETA,one); + Ees = ge(one,1-ZETA); + Een = ge(zero,1-ZETA); + Enb = gn(1-XI,one); + Ent = gn(1-XI,zero); + Est = gs(XI,one); + Esb = gs(XI,zero); - Cwbs = gw(one,one); - Cwbn = gw(zero,one); - Cwts = gw(one,zero); - Cwtn = gw(zero,zero); - Cebs = ge(zero,zero); - Cebn = ge(one,zero); - Cets = ge(zero,one); - Cetn = ge(one,one); + Cwbs = gw(zero,one); + Cwbn = gw(one,one); + Cwts = gw(zero,zero); + Cwtn = gw(one,zero); + Cebs = ge(one,one); + Cebn = ge(zero,one); + Cets = ge(one,zero); + Cetn = ge(zero,zero); X1 = (1-XI).*Sw(1,:,:) + XI.*Se(1,:,:); @@ -104,7 +104,7 @@ obj.V = @V_fun; end - + %Should be rewritten so that the input is xi eta zeta function [X,Y,Z] = map(obj,XI,ETA,ZETA) V = obj.V; @@ -247,6 +247,6 @@ % grid.place_label(ps,'s'); % grid.place_label(pn,'n'); % end - % end + % end end end \ No newline at end of file
--- a/+scheme/Hypsyst2dCurve.m Thu Oct 13 09:34:30 2016 +0200 +++ b/+scheme/Hypsyst2dCurve.m Wed Nov 02 00:02:01 2016 +0100 @@ -259,7 +259,8 @@ L=obj.evaluateCoefficientMatrix(L,eta_x,eta_y,[],[]); side=max(length(xi)); case {'n','N','north'} - e_=obj.e_n; + e_=obj.e_n; for ii=1:rows + mat=obj.Bhat; boundPos='r'; Hi=obj.Hetai;
--- a/+scheme/Hypsyst3d.m Thu Oct 13 09:34:30 2016 +0200 +++ b/+scheme/Hypsyst3d.m Wed Nov 02 00:02:01 2016 +0100 @@ -10,6 +10,7 @@ D % non-stabalized scheme operator A, B, C, E + Aevaluated,Bevaluated,Cevaluated, Eevaluated H % Discrete norm % Norms in the x, y and z directions @@ -61,10 +62,10 @@ obj.Xz=kr(obj.x,ones(m_y,1)); obj.Yz=kr(ones(m_z,1),obj.y); - Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y,obj.Z); - Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y,obj.Z); - Cevaluated = obj.evaluateCoefficientMatrix(C, obj.X, obj.Y,obj.Z); - Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y,obj.Z); + obj.Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y,obj.Z); + obj.Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y,obj.Z); + obj.Cevaluated = obj.evaluateCoefficientMatrix(C, obj.X, obj.Y,obj.Z); + obj.Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y,obj.Z); obj.n = length(A(obj.params,0,0,0)); @@ -82,7 +83,7 @@ D1_y = kr(I_n, I_x, ops_y.D1,I_z); obj.Hyi = kr(I_n, I_x, ops_y.HI,I_z); D1_z = kr(I_n, I_x, I_y,ops_z.D1); - obj.Hzi = kr(I_n, I_x,I_y, ops_y.HI); + obj.Hzi = kr(I_n, I_x,I_y, ops_z.HI); obj.e_w = kr(I_n, ops_x.e_l, I_y,I_z); obj.e_e = kr(I_n, ops_x.e_r, I_y,I_z); @@ -95,7 +96,7 @@ obj.h=[ops_x.h ops_y.h ops_x.h]; obj.order=order; - obj.D=-Aevaluated*D1_x-Bevaluated*D1_y-Cevaluated*D1_z-Eevaluated; + obj.D=-obj.Aevaluated*D1_x-obj.Bevaluated*D1_y-obj.Cevaluated*D1_z-obj.Eevaluated; end % Closure functions return the opertors applied to the own doamin to close the boundary
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Hypsyst3dCurve.m Wed Nov 02 00:02:01 2016 +0100 @@ -0,0 +1,479 @@ +classdef Hypsyst3dCurve < scheme.Scheme + properties + m % Number of points in each direction, possibly a vector + n %size of system + 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 + + X_xi, X_eta, X_zeta,Y_xi,Y_eta,Y_zeta,Z_xi,Z_eta,Z_zeta + Aev + + metric_terms + + order % Order accuracy for the approximation + + D % non-stabalized scheme operator + Aevaluated, Bevaluated, Cevaluated, Eevaluated + Ahat, Bhat, Chat, E + A,B,C + + J, Ji + + 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. + I_xi,I_eta,I_zeta, I_N,onesN + e_w, e_e, e_s, e_n, e_b, e_t + 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) + xilim ={0 1}; + etalim = {0 1}; + zetalim = {0 1}; + + if length(m) == 1 + m = [m m m]; + end + m_xi = m(1); + m_eta = m(2); + m_zeta=m(3); + 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); + ops_xi = sbp.D2Standard(m_xi,xilim,order); + ops_eta = sbp.D2Standard(m_eta,etalim,order); + ops_zeta = sbp.D2Standard(m_zeta,zetalim,order); + + 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.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); + + 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); + + [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; + I_eta = speye(m_eta); + obj.I_eta = I_eta; + I_zeta = speye(m_zeta); + obj.I_zeta = I_zeta; + + + O_xi=ones(m_xi,1); + O_eta=ones(m_eta,1); + O_zeta=ones(m_zeta,1); + + D1_xi = kr(ops_xi.D1, I_eta,I_zeta); + obj.Hxii = kr(I_n, ops_xi.HI, I_eta,I_zeta); + D1_eta = kr(I_xi, ops_eta.D1,I_zeta); + obj.Hetai = kr(I_n, I_xi, ops_eta.HI,I_zeta); + D1_zeta = kr(I_xi, I_eta,ops_zeta.D1); + obj.Hzetai = kr(I_n, I_xi,I_eta, ops_zeta.HI); + obj.h=[ops_xi.h ops_eta.h ops_zeta.h]; + + 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.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; + obj.Y_xi=D1_xi*Y; + 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; + + D1_xi=kr(I_n,D1_xi); + D1_eta=kr(I_n,D1_eta); + D1_zeta=kr(I_n,D1_zeta); + + 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.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); + obj.Eevaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,obj.Z,[],[],[],[],[],[]); + + 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.D=obj.Ji*(-obj.Aevaluated*D1_xi-obj.Bevaluated*D1_eta -obj.Cevaluated*D1_zeta)-obj.Eevaluated; + 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'. + % 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'); + BM=boundary_matrices(obj,boundary); + + switch type + case{'c','char'} + [closure,penalty]=boundary_condition_char(obj,BM); + case{'general'} + [closure,penalty]=boundary_condition_general(obj,BM,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, 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)); + x_1=kr(obj.onesN,x_1); + x_2=kr(obj.onesN,x_2); + y_1=kr(obj.onesN,y_1); + y_2=kr(obj.onesN,y_2); + z_1=kr(obj.onesN,z_1); + z_2=kr(obj.onesN,z_2); + matVec=mat(obj,X',Y',Z',x_1',x_2',y_1',y_2',z_1',z_2'); + matVec=sparse(matVec); + else + matVec=mat; + [rows,cols]=size(matVec); + side=max(length(X),length(Y)); + cols=cols/side; + end + ret=kron(ones(rows,cols),speye(side)); + + for ii=1:rows + for jj=1:cols + ret((ii-1)*side+1:ii*side,(jj-1)*side+1:jj*side)=diag(matVec(ii,(jj-1)*side+1:jj*side)); + end + end + end + + + function [BM]=boundary_matrices(obj,boundary) + params=obj.params; + BM.boundary=boundary; + switch boundary + case {'w','W','west'} + BM.e_=obj.e_w; + mat=obj.Ahat; + BM.boundpos='l'; + BM.Hi=obj.Hxii; + BM.index=obj.index_w; + BM.x_1=obj.X_eta(BM.index); + BM.x_2=obj.X_zeta(BM.index); + BM.y_1=obj.Y_eta(BM.index); + BM.y_2=obj.Y_zeta(BM.index); + BM.z_1=obj.Z_eta(BM.index); + BM.z_2=obj.Z_zeta(BM.index); + case {'e','E','east'} + BM.e_=obj.e_e; + mat=obj.Ahat; + BM.boundpos='r'; + BM.Hi=obj.Hxii; + BM.index=obj.index_e; + BM.x_1=obj.X_eta(BM.index); + BM.x_2=obj.X_zeta(BM.index); + BM.y_1=obj.Y_eta(BM.index); + BM.y_2=obj.Y_zeta(BM.index); + BM.z_1=obj.Z_eta(BM.index); + BM.z_2=obj.Z_zeta(BM.index); + case {'s','S','south'} + BM.e_=obj.e_s; + mat=obj.Bhat; + BM.boundpos='l'; + BM.Hi=obj.Hetai; + BM.index=obj.index_s; + BM.x_1=obj.X_zeta(BM.index); + BM.x_2=obj.X_xi(BM.index); + BM.y_1=obj.Y_zeta(BM.index); + BM.y_2=obj.Y_xi(BM.index); + BM.z_1=obj.Z_zeta(BM.index); + BM.z_2=obj.Z_xi(BM.index); + case {'n','N','north'} + BM.e_=obj.e_n; + mat=obj.Bhat; + BM.boundpos='r'; + BM.Hi=obj.Hetai; + BM.index=obj.index_n; + BM.x_1=obj.X_zeta(BM.index); + BM.x_2=obj.X_xi(BM.index); + BM.y_1=obj.Y_zeta(BM.index); + BM.y_2=obj.Y_xi(BM.index); + BM.z_1=obj.Z_zeta(BM.index); + BM.z_2=obj.Z_xi(BM.index); + case{'b','B','Bottom'} + BM.e_=obj.e_b; + mat=obj.Chat; + BM.boundpos='l'; + BM.Hi=obj.Hzetai; + BM.index=obj.index_b; + BM.x_1=obj.X_xi(BM.index); + BM.x_2=obj.X_eta(BM.index); + BM.y_1=obj.Y_xi(BM.index); + BM.y_2=obj.Y_eta(BM.index); + BM.z_1=obj.Z_xi(BM.index); + BM.z_2=obj.Z_eta(BM.index); + case{'t','T','Top'} + BM.e_=obj.e_t; + mat=obj.Chat; + BM.boundpos='r'; + BM.Hi=obj.Hzetai; + BM.index=obj.index_t; + BM.x_1=obj.X_xi(BM.index); + BM.x_2=obj.X_eta(BM.index); + BM.y_1=obj.Y_xi(BM.index); + BM.y_2=obj.Y_eta(BM.index); + BM.z_1=obj.Z_xi(BM.index); + BM.z_2=obj.Z_eta(BM.index); + end + [BM.V,BM.Vi,BM.D,signVec]=obj.matrixDiag(mat,obj.X(BM.index),obj.Y(BM.index),obj.Z(BM.index),... + 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 + + + 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,:); + 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 + end + + + 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_; + index=BM.index; + + switch BM.boundary + case{'b','B','bottom'} + Ji_vec=diag(obj.Ji); + Ji=diag(Ji_vec(index)); + 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); + 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; + end + end + + + 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 + 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); + + syms x_1s x_2s y_1s y_2s z_1s z_2s +% V= matlabFunction(V); +% D= matlabFunction(D); +% Vi= matlabFunction(Vi); +% +% xs=x; +% ys=y; +% zs=z; +% x_1s=x_1; +% x_2s=x_2; +% y_1s=y_1; +% 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 + Dpart=matlabFunction(D(jj,ii),'Vars',[xs ys zs x_1s x_2s y_1s y_2s z_1s z_2s]); + Vpart=matlabFunction(V(jj,ii),'Vars',[xs ys zs x_1s x_2s y_1s y_2s z_1s z_2s]); + Vipart=matlabFunction(V(jj,ii),'Vars',[xs ys zs x_1s x_2s y_1s y_2s z_1s z_2s]); + Dret(jj,(ii-1)*side+1:side*ii)=sparse(Dpart(x,y,z,x_1,x_2,y_1,y_2,z_1,z_2)); + Vret(jj,(ii-1)*side+1:side*ii)=sparse(Vpart(x,y,z,x_1,x_2,y_1,y_2,z_1,z_2)); + Viret(jj,(ii-1)*side+1:side*ii)=sparse(Vipart(x,y,z,x_1,x_2,y_1,y_2,z_1,z_2)); + end + end + + D=sparse(Dret); + V=sparse(Vret); + Vi=sparse(Viret); + V=obj.evaluateCoefficientMatrix(V,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2); + 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=inv(V); + signVec=[sum(poseig),sum(zeroeig),sum(negeig)]; + end + end + end