Mercurial > repos > public > sbplib
changeset 301:d9860ebc3148 feature/hypsyst
HypsystCurve2D Seems to work (Converges with MMS)
| author | Ylva Rydin <ylva.rydin@telia.com> |
|---|---|
| date | Wed, 05 Oct 2016 17:36:34 +0200 |
| parents | 32f3ee81bc37 |
| children | cd6a29ab3746 |
| files | +scheme/Hypsyst2d.m +scheme/Hypsyst2dCurve.m |
| diffstat | 2 files changed, 90 insertions(+), 58 deletions(-) [+] |
line wrap: on
line diff
diff -r 32f3ee81bc37 -r d9860ebc3148 +scheme/Hypsyst2d.m --- a/+scheme/Hypsyst2d.m Tue Oct 04 16:25:38 2016 +0200 +++ b/+scheme/Hypsyst2d.m Wed Oct 05 17:36:34 2016 +0200 @@ -128,8 +128,7 @@ function [closure, penalty]=boundary_condition_char(obj,boundary) params=obj.params; x=obj.x; y=obj.y; - side=max(length(x),length(y)); - + switch boundary case {'w','W','west'} e_=obj.e_w; @@ -137,39 +136,43 @@ boundPos='l'; Hi=obj.Hxi; [V,Vi,D,signVec]=obj.matrixDiag(mat,x(1),y); + side=max(length(y)); case {'e','E','east'} e_=obj.e_e; mat=obj.A; boundPos='r'; Hi=obj.Hxi; [V,Vi,D,signVec]=obj.matrixDiag(mat,x(end),y); + side=max(length(y)); case {'s','S','south'} e_=obj.e_s; mat=obj.B; boundPos='l'; Hi=obj.Hyi; [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(1)); + side=max(length(x)); case {'n','N','north'} e_=obj.e_n; mat=obj.B; boundPos='r'; Hi=obj.Hyi; [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(end)); + side=max(length(x)); end pos=signVec(1); zeroval=signVec(2); neg=signVec(3); switch boundPos case {'l'} - tau=sparse(obj.n*side,pos*side); - Vi_plus=Vi(1:pos*side,:); - tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side)); + 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*side); - tau((pos+zeroval)*side+1:obj.n*side,:)=-abs(D((pos+zeroval)*side+1:obj.n*side,(pos+zeroval)*side+1:obj.n*side)); - Vi_minus=Vi((pos+zeroval)*side+1:obj.n*side,:); + 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 @@ -179,7 +182,6 @@ function [closure,penalty]=boundary_condition_general(obj,boundary,L) params=obj.params; x=obj.x; y=obj.y; - side=max(length(x),length(y)); switch boundary case {'w','W','west'} @@ -189,6 +191,7 @@ Hi=obj.Hxi; [V,Vi,D,signVec]=obj.matrixDiag(mat,x(1),y); L=obj.evaluateCoefficientMatrix(L,x(1),y); + side=max(length(y)); case {'e','E','east'} e_=obj.e_e; mat=obj.A; @@ -196,6 +199,7 @@ Hi=obj.Hxi; [V,Vi,D,signVec]=obj.matrixDiag(mat,x(end),y); L=obj.evaluateCoefficientMatrix(L,x(end),y); + side=max(length(y)); case {'s','S','south'} e_=obj.e_s; mat=obj.B; @@ -203,6 +207,7 @@ Hi=obj.Hyi; [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(1)); L=obj.evaluateCoefficientMatrix(L,x,y(1)); + side=max(length(x)); case {'n','N','north'} e_=obj.e_n; mat=obj.B; @@ -210,30 +215,31 @@ Hi=obj.Hyi; [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(end)); L=obj.evaluateCoefficientMatrix(L,x,y(end)); + side=max(length(x)); end pos=signVec(1); zeroval=signVec(2); neg=signVec(3); switch boundPos case {'l'} - tau=sparse(obj.n*side,pos*side); - Vi_plus=Vi(1:pos*side,:); - Vi_minus=Vi(pos*side+1:obj.n*side,:); - V_plus=V(:,1:pos*side); - V_minus=V(:,(pos+zeroval)*side+1:obj.n*side); + 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*side,:)=-abs(D(1:pos*side,1:pos*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*side); - tau((pos+zeroval)*side+1:obj.n*side,:)=-abs(D((pos+zeroval)*side+1:obj.n*side,(pos+zeroval)*side+1:obj.n*side)); - Vi_plus=Vi(1:pos*side,:); - Vi_minus=Vi((pos+zeroval)*side+1:obj.n*side,:); + 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*side); - V_minus=V(:,(pos+zeroval)*side+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; @@ -243,21 +249,11 @@ function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y) params=obj.params; - syms xs ys; + syms xs ys [V, D]=eig(mat(params,xs,ys)); - xs=1;ys=1; - DD=eval(diag(D)); - - poseig=find(DD>0); - zeroeig=find(DD==0); - negeig=find(DD<0); - syms xs ys - DD=diag(D); - - D=diag([DD(poseig);DD(zeroeig); DD(negeig)]); - V=[V(:,poseig) V(:,zeroeig) V(:,negeig)]; - xs=x; ys=y; - + xs=x; + ys=y; + side=max(length(x),length(y)); Dret=zeros(obj.n,side*obj.n); Vret=zeros(obj.n,side*obj.n); @@ -269,11 +265,19 @@ end D=sparse(Dret); - V=sparse(normc(Vret)); + V=sparse(Vret); V=obj.evaluateCoefficientMatrix(V,x,y); - D=obj.evaluateCoefficientMatrix(D,x,y); + D=obj.evaluateCoefficientMatrix(D,x,y); + 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=[length(poseig),length(zeroeig),length(negeig)]; + signVec=[sum(poseig),sum(zeroeig),sum(negeig)]; end end
diff -r 32f3ee81bc37 -r d9860ebc3148 +scheme/Hypsyst2dCurve.m --- a/+scheme/Hypsyst2dCurve.m Tue Oct 04 16:25:38 2016 +0200 +++ b/+scheme/Hypsyst2dCurve.m Wed Oct 05 17:36:34 2016 +0200 @@ -5,7 +5,7 @@ h % Grid spacing X,Y % Values of x and y for each grid point - J, Ji + J, Ji %Jacobaian and inverse Jacobian xi,eta Xi,Eta @@ -18,8 +18,7 @@ Ahat, Bhat, E H % Discrete norm - % Norms in the x and y directions - Hxii,Hetai % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. + Hxii,Hetai % Kroneckerd norms in xi and eta. I_xi,I_eta, I_N, onesN e_w, e_e, e_s, e_n index_w, index_e,index_s,index_n @@ -39,8 +38,7 @@ obj.params = params; obj.A=A; obj.B=B; - - + obj.Ahat=@(params,x,y,x_eta,y_eta)(A(params,x,y).*y_eta-B(params,x,y).*x_eta); obj.Bhat=@(params,x,y,x_xi,y_xi)(B(params,x,y).*x_xi-A(params,x,y).*y_xi); obj.E=@(params,x,y,~,~)E(params,x,y); @@ -61,10 +59,10 @@ obj.n = length(A(obj.params,0,0)); obj.onesN=ones(obj.n); - obj.index_w=1:m_xi; - obj.index_e=(m_tot-m_xi+1):m_tot; - obj.index_s=1:m_xi:(m_tot-m_xi+1); - obj.index_n=(m_xi):m_xi:m_tot; + obj.index_w=1:m_eta; + obj.index_e=(m_tot-m_eta+1):m_tot; + obj.index_s=1:m_eta:(m_tot-m_eta+1); + obj.index_n=(m_eta):m_eta:m_tot; I_n = eye(obj.n); I_xi = speye(m_xi); @@ -138,8 +136,8 @@ if isa(mat,'function_handle') [rows,cols]=size(mat(params,0,0,0,0)); - x_=kr(x_,obj.onesN); - y_=kr(y_,obj.onesN); + x_=kr(obj.onesN,x_); + y_=kr(obj.onesN,y_); matVec=mat(params,X',Y',x_',y_'); matVec=sparse(matVec); side=max(length(X),length(Y)); @@ -163,7 +161,7 @@ params=obj.params; X=obj.X; Y=obj.Y; xi=obj.xi; eta=obj.eta; - side=max(length(xi),length(eta)); + switch boundary case {'w','W','west'} @@ -172,24 +170,28 @@ boundPos='l'; Hi=obj.Hxii; [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_w),Y(obj.index_w),obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w)); + side=max(length(eta)); case {'e','E','east'} e_=obj.e_e; mat=obj.Ahat; boundPos='r'; Hi=obj.Hxii; [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e)); + side=max(length(eta)); case {'s','S','south'} e_=obj.e_s; mat=obj.Bhat; boundPos='l'; Hi=obj.Hetai; [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_s),Y(obj.index_s),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s)); + side=max(length(xi)); case {'n','N','north'} e_=obj.e_n; mat=obj.Bhat; boundPos='r'; Hi=obj.Hetai; [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n)); + side=max(length(xi)); end pos=signVec(1); zeroval=signVec(2); neg=signVec(3); @@ -214,7 +216,7 @@ function [closure,penalty]=boundary_condition_general(obj,boundary,L) params=obj.params; X=obj.X; Y=obj.Y; - side=max(length(xi),length(eta)); + xi=obj.xi; eta=obj.eta; switch boundary case {'w','W','west'} @@ -223,28 +225,54 @@ boundPos='l'; Hi=obj.Hxii; [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_w),Y(obj.index_w),obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w)); - L=obj.evaluateCoefficientMatrix(L,X(obj.index_w),Y(obj.index_w)); + + Ji_vec=diag(obj.Ji); + Ji=diag(Ji_vec(obj.index_w)); + xi_x=Ji*obj.Y_eta(obj.index_w); + xi_y=-Ji*obj.X_eta(obj.index_w); + L=obj.evaluateCoefficientMatrix(L,xi_x,xi_y,[],[]); + side=max(length(eta)); case {'e','E','east'} e_=obj.e_e; mat=obj.Ahat; boundPos='r'; Hi=obj.Hxii; - [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e)); - L=obj.evaluateCoefficientMatrix(L,X(obj.index_e),Y(obj.index_e),[],[]); + [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e)); + + Ji_vec=diag(obj.Ji); + Ji=diag(Ji_vec(obj.index_e)); + xi_x=Ji*obj.Y_eta(obj.index_e); + xi_y=-Ji*obj.X_eta(obj.index_e); + L=obj.evaluateCoefficientMatrix(L,-xi_x,-xi_y,[],[]); + side=max(length(eta)); case {'s','S','south'} e_=obj.e_s; mat=obj.Bhat; boundPos='l'; Hi=obj.Hetai; [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_s),Y(obj.index_s),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s)); - L=obj.evaluateCoefficientMatrix(L,X(obj.index_s),Y(obj.index_s),[],[]); + + Ji_vec=diag(obj.Ji); + Ji=diag(Ji_vec(obj.index_s)); + eta_x=Ji*obj.Y_xi(obj.index_s); + eta_y=-Ji*obj.X_xi(obj.index_s); + L=obj.evaluateCoefficientMatrix(L,eta_x,eta_y,[],[]); + side=max(length(xi)); case {'n','N','north'} e_=obj.e_n; mat=obj.Bhat; boundPos='r'; Hi=obj.Hetai; [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n)); - L=obj.evaluateCoefficientMatrix(L,X(obj.index_n),Y(obj.index_n)); + + Ji_vec=diag(obj.Ji); + Ji=diag(Ji_vec(obj.index_n)); + eta_x=Ji*obj.Y_xi(obj.index_n); + eta_y=-Ji*obj.X_xi(obj.index_n); + L=obj.evaluateCoefficientMatrix(L,-eta_x,-eta_y,[],[]); + + side=max(length(xi)); + end pos=signVec(1); zeroval=signVec(2); neg=signVec(3); @@ -255,9 +283,9 @@ Vi_plus=Vi(1:pos,:); Vi_minus=Vi(pos+1:obj.n*side,:); V_plus=V(:,1:pos); - V_minus=V(:,(pos+zeroval)+1:obj.n*side); + V_minus=V(:,(pos)+1:obj.n*side); - tau(1:pos*side,:)=-abs(D(1:pos,1:pos)); + 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;