Mercurial > repos > public > sbplib

comparison +scheme/hypsyst2d.m @ 290:d32f674bcbe5 feature/hypsyst

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
A first attempt to make a general scheme fo hyperbolic systems
author Ylva Rydin <ylva.rydin@telia.com>
date Fri, 16 Sep 2016 14:51:17 +0200
parents
children 807dfe8be3ec
comparison
equal deleted inserted replaced
285:70184f6c6cb5 290:d32f674bcbe5
1 classdef hypsyst2d < scheme.Scheme
2 properties
3 m % Number of points in each direction, possibly a vector
4 h % Grid spacing
5 x,y % Grid
6 X,Y % Values of x and y for each grid point
7 order % Order accuracy for the approximation
8
9 D % non-stabalized scheme operator
10 A, B, E
11
12 H % Discrete norm
13 Hi
14 H_x, H_y % Norms in the x and y directions
15 Hx,Hy % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
16 I_x,I_y
17 e_w, e_e, e_s, e_n
18 params %parameters for the coeficient matrices
19 end
20
21
22 methods
23 function obj = hypsyst2d(m,lim,order,matrices,params)
24
25 xlim = lim{1};
26 ylim = lim{2};
27
28 if length(m) == 1
29 m = [m m];
30 end
31
32 m_x = m(1);
33 m_y = m(2);
34
35 ops_x = sbp.D2Standard(m_x,xlim,order);
36 ops_y = sbp.D2Standard(m_y,ylim,order);
37
38 obj.x=ops_x.x;
39 obj.y=ops_y.y;
40
41 obj.X = kr(x,ones(m_y,1));
42 obj.Y = kr(ones(m_x,1),y);
43
44 I_x = speye(m_x);
45 I_y = speye(m_y);
46 I_n= eye(4);
47
48
49 D1_x = kr(kr(I_n,ops_x.D1),I_y);
50 obj.Hi_x= kr(kr(I_n,ops_x.HI),I_y);
51 D1_y=kr(I_n,kr(I_x,ops_y.D1));
52 obj.Hi_y=kr(I_n,kr(I_x,ops_y.HI));
53
54 obj.e_w=kr(I_n,kr(ops_x.e_l,I_y));
55 obj.e_e=kr(I_n,kr(ops_x.e_r,I_y));
56 obj.e_s=kr(I_n,kr(I_x,ops_y.e_l));
57 obj.e_n=kr(I_n,kr(I_x,ops_y.e_r));
58
59 obj.m=m;
60 obj.h=[ops_x.h ops_y.h];
61 obj.order=order;
62 obj.params=params;
63
64 obj.A=matrixBuild(obj,matrices.A);
65 obj.B=matrixBuild(obj,matrices.B);
66 obj.E=matrixBuild(obj,matrices.E);
67
68 obj.D=-obj.A*D1_x-obj.B*D1_y-E;
69
70 end
71 % Closure functions return the opertors applied to the own doamin to close the boundary
72 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
73 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
74 % type is a string specifying the type of boundary condition if there are several.
75 % data is a function returning the data that should be applied at the boundary.
76 % neighbour_scheme is an instance of Scheme that should be interfaced to.
77 % neighbour_boundary is a string specifying which boundary to interface to.
78 function [closure, penalty] = boundary_condition(obj,boundary,type,data)
79 default_arg('type','neumann');
80 default_arg('data',0);
81
82 switch type
83 case{c,'char'}
84 [tau,e_,Hi,CHM]=GetBoundarydata(obj,boundary,type);
85 closure =Hi*e_*tau*CHM*e_';
86 penalty =Hi*e_*tau*CHM*data;
87 otherwise
88 error('No such boundary condition')
89 end
90 end
91
92 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
93 error('An interface function does not exist yet');
94 end
95
96 function N = size(obj)
97 N = obj.m;
98 end
99
100 end
101
102 methods(Static)
103 % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u
104 % and bound_v of scheme schm_v.
105 % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l')
106 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
107 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
108 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
109 end
110
111 function [ret]=matrixBuild(obj,mat)
112 %extra info for coordinate transfomration mult my y_ny and
113 %x,ny osv...
114 params=obj.params;
115 X=obj.X;
116 Y=obj.Y;
117
118 if isa(mat,'function_handle')
119 matVec=mat(params,x,y);
120 side=length(x);
121 else
122 matVec=mat;
123 side=max(size(mat))/4;
124 end
125
126
127 ret=[diag(matVec(1,(1-1)*side+1:1*side)) diag(matVec(1,(2-1)*side+1:2*side)) diag(matVec(1,(3-1)*side+1:3*side)) diag(matVec(1,(4-1)*side+1:4*side))
128 diag(matVec(2,(1-1)*side+1:1*side)) diag(matVec(2,(2-1)*side+1:2*side)) diag(matVec(2,(3-1)*side+1:3*side)) diag(matVec(2,(4-1)*side+1:4*side))
129 diag(matVec(3,(1-1)*side+1:1*side)) diag(matVec(3,(2-1)*side+1:2*side)) diag(matVec(3,(3-1)*side+1:3*side)) diag(matVec(3,(4-1)*side+1:4*side))
130 diag(matVec(4,(1-1)*side+1:1*side)) diag(matVec(4,(2-1)*side+1:2*side)) diag(matVec(4,(3-1)*side+1:3*side)) diag(matVec(4,(4-1)*side+1:4*side))];
131 end
132
133 function [tau,e_,Hi, CHM]=GetBoundarydata(obj,boundary)
134 params=obj.params;
135 x=obj.x;
136 y=obj.y;
137
138 side=max(length(x),length(y));
139
140
141 switch boundary
142 case {'w','W','west'}
143 e_=obj.e_w;
144 mat=obj.A;
145 [V,D]=matrixDiag(mat,params,x(1),y);
146 Hi=obj.Hx;
147 tau=zeros(4*side,pos*side);
148 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side));
149 CHM=V*((D+abs(D))/2)*V';
150 case {'e','E','east'}
151 e_=obj.e_e;
152 mat=obj.A;
153 [V,D]=matrixDiag(mat,params,x(end),y);
154 Hi=obj.Hy;
155 tau=zeros(4*side,(4-pos));
156 tau(pos*side+1:4*side)=-abs(D(pos*side+1:4*side,pos*side+1:4*side));
157 CHM=V*((D-abs(D))/2)*V';
158 case {'s','S','south'}
159 e_=obj.e_s;
160 mat=obj.B;
161 [V,D]=matrixDiag(mat,params,x,y(1));
162 Hi=obj.Hx;
163 tau=zeros(4*side,pos*side);
164 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side));
165 CHM=V*((D+abs(D))/2)*V';
166 case {'n','N','north'}
167 e_=obj.e_n;
168 mat=obj.B;
169 [V,D]=matrixDiag(mat,params,x,y(end));
170 Hi=obj.Hy;
171 tau=zeros(4*side,(4-pos));
172 tau(pos*side+1:4*side)=-abs(D(pos*side+1:4*side,pos*side+1:4*side));
173 CHM=V*((D-abs(D))/2)*V';
174 end
175
176 tau=V*tau*V';
177
178 end
179
180 function [V, D,pos]=matrixDiag(mat,params,x,y)
181 syms xs ys;
182 [V, D]=eig(mat(params,xs,ys));
183 xs=1;ys=1;
184 DD=eval(diag(D));
185
186 pos=find(DD>=0); %Now zero eigenvalues are calculated as possitive, Maybe it should not????
187 neg=find(DD<0);
188 syms xs ys
189 DD=diag(D);
190
191 D=diag([DD(pos); DD(neg)]);
192 V=[V(:,pos) V(:,neg)];
193
194 xs=x; ys=y;
195
196 side=max(length(x),length(y));
197 Dret=zeros(4,side*4);
198 Vret=zeros(4,side*4);
199 for ii=1:4
200 for jj=1:4
201 Dret(jj,(ii-1)*side+1:side*ii)=eval(D(jj,ii));
202 Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii));
203 end
204 end
205 V=sparse(normc(Vret));
206 D=sparse(Dret);
207
208
209 V=matrixBuild([],[],[],V);
210 D=matrixBuild([],[],[],D);
211 pos=legth(pos);
212 end
213 end
214 end