Mercurial > repos > public > sbplib

changeset 301:d9860ebc3148 feature/hypsyst

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
--- 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
--- 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;