Mercurial > repos > public > sbplib
changeset 293:2d604d16842c feature/hypsyst
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
--- 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