Mercurial > repos > public > sbplib

changeset 293:2d604d16842c feature/hypsyst

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Works with varying coefficients and char boundary condition
author Ylva Rydin <ylva.rydin@telia.com>
date Fri, 23 Sep 2016 16:30:51 +0200
parents 3d275c5e45b3
children 8ff6ec6249e8
files +scheme/hypsyst2d.m
diffstat 1 files changed, 30 insertions(+), 44 deletions(-) [+]
line wrap: on
line diff
diff -r 3d275c5e45b3 -r 2d604d16842c +scheme/hypsyst2d.m
--- a/+scheme/hypsyst2d.m	Fri Sep 23 14:48:54 2016 +0200
+++ b/+scheme/hypsyst2d.m	Fri Sep 23 16:30:51 2016 +0200
@@ -1,6 +1,7 @@
 classdef hypsyst2d < scheme.Scheme
     properties
         m % Number of points in each direction, possibly a vector
+        n %size of system
         h % Grid spacing
         x,y % Grid
         X,Y % Values of x and y for each grid point
@@ -31,6 +32,7 @@
             
             m_x = m(1);
             m_y = m(2);
+            obj.params=params;
             
             obj.matrices=matrices;
             
@@ -43,12 +45,17 @@
             obj.X = kr(obj.x,ones(m_y,1));
             obj.Y = kr(ones(m_x,1),obj.y);
             
+            obj.A=obj.matrixBuild(matrices.A);
+            obj.B=obj.matrixBuild(matrices.B);
+            obj.E=obj.matrixBuild(matrices.E);
+            
+            obj.n=length(matrices.A(obj.params,0,0));
+           
+            I_n=  eye(obj.n);         
             I_x = speye(m_x); obj.I_x=I_x;
             I_y = speye(m_y); obj.I_y=I_y;
       
-            I_n=  eye(4);
-            
-            
+           
             D1_x = kr(kr(I_n,ops_x.D1),I_y);
             obj.Hxi= kr(kr(I_n,ops_x.HI),I_y);
             D1_y=kr(I_n,kr(I_x,ops_y.D1));
@@ -62,13 +69,7 @@
             obj.m=m;
             obj.h=[ops_x.h ops_y.h];
             obj.order=order;
-            obj.params=params;
-            
-            obj.A=obj.matrixBuild(matrices.A);
-            obj.B=obj.matrixBuild(matrices.B);
-            obj.E=obj.matrixBuild(matrices.E);
-            
-            
+                
             obj.D=-obj.A*D1_x-obj.B*D1_y-obj.E;
            
         end
@@ -79,15 +80,15 @@
         %       data                is a function returning the data that should be applied at the boundary.
         %       neighbour_scheme    is an instance of Scheme that should be interfaced to.
         %       neighbour_boundary  is a string specifying which boundary to interface to.
-        function [closure, penalty] = boundary_condition(obj,boundary,type)
+        function [closure, penalty] = boundary_condition(obj,boundary,type,L)
             default_arg('type','neumann');
             default_arg('data',0);
             
             switch type
                 case{'c','char'}
                     [closure,penalty]=GetBoundarydata_char(obj,boundary);
-                case{'wall'}
-                    [closure,penalty]=GetBoundarydata_wall(obj,boundary);
+                case{'general'}
+                    [closure,penalty]=GeneralBoundaryCond(obj,boundary,L);
                 otherwise
                     error('No such boundary condition')
             end
@@ -170,43 +171,28 @@
                    
             switch boundPos
                 case {'l'}  
-                    tau=sparse(4*side,pos*side);  
+                    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));             
                     closure=Hi*e_*V*tau*Vi_plus*e_';
                     penalty=-Hi*e_*V*tau*Vi_plus;
        
                 case {'r'}                
-                    tau=sparse(4*side,neg*side);
-                    tau((pos+zeroval)*side+1:4*side,:)=-abs(D((pos+zeroval)*side+1:4*side,(pos+zeroval)*side+1:4*side));
-                    Vi_minus=Vi((pos+zeroval)*side+1:4*side,:);              
+                    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,:);              
                     closure=Hi*e_*V*tau*Vi_minus*e_';
                     penalty=-Hi*e_*V*tau*Vi_minus;
                     
             end   
         end
-        
-        function [closure, penalty]=GetBoundarydata_wall(obj,boundary)
-            switch boundary
-                case {'e','w'}
-                    L=[0 1 0 0]';
-                    L=kr(L,obj.I_y);
-                    L=L';
-                case {'s','n'}
-                    L=[0 0 1 0]';
-                    L=kr(L,obj.I_x);
-                    L=L';
-                    
-            end
-            [closure,penalty]=GeneralBoundaryCond(obj,boundary,L);
-        end
           
         
         function [closure,penalty]=GeneralBoundaryCond(obj,boundary,L)
               params=obj.params;
             x=obj.x;
             y=obj.y;
-            
+            L=obj.matrixBuild(L,x,y);
             side=max(length(x),length(y));
             
              
@@ -243,12 +229,12 @@
             
             switch boundPos
                 case {'l'}
-                    tau=sparse(4*side,pos*side);  
+                    tau=sparse(obj.n*side,pos*side);  
                     Vi_plus=Vi(1:pos*side,:); 
-                    Vi_minus=Vi(pos*side+1:4*side,:);
+                    Vi_minus=Vi(pos*side+1:obj.n*side,:);
                     
                     V_plus=Vi(:,1:pos*side); 
-                    V_minus=Vi(:,(pos+zeroval)*side+1:4*side);
+                    V_minus=Vi(:,(pos+zeroval)*side+1:obj.n*side);
                     
                     tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side));  
                     R=-inv(L*V_plus)*(L*V_minus);
@@ -257,13 +243,13 @@
                     
        
                 case {'r'}
-                    tau=sparse(4*side,neg*side);
-                    tau((pos+zeroval)*side+1:4*side,:)=-abs(D((pos+zeroval)*side+1:4*side,(pos+zeroval)*side+1:4*side));
+                    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:4*side,:);
+                    Vi_minus=Vi((pos+zeroval)*side+1:obj.n*side,:);
                     
                     V_plus=Vi(:,1:pos*side); 
-                    V_minus=Vi(:,(pos+zeroval)*side+1:4*side);
+                    V_minus=Vi(:,(pos+zeroval)*side+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;
@@ -293,10 +279,10 @@
        
             
             side=max(length(x),length(y));
-            Dret=zeros(4,side*4);
-            Vret=zeros(4,side*4);
-            for ii=1:4
-                for jj=1:4
+            Dret=zeros(obj.n,side*obj.n);
+            Vret=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));
                     Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii));
                 end