Mercurial > repos > public > sbplib

comparison +scheme/hypsyst2d.m @ 291:807dfe8be3ec feature/hypsyst

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Have made a lot of stupid changes in hypsyst in order to find a stupid bug
author Ylva Rydin <ylva.rydin@telia.com>
date Wed, 21 Sep 2016 16:30:34 +0200
parents d32f674bcbe5
children 3d275c5e45b3
comparison
equal deleted inserted replaced
290:d32f674bcbe5 291:807dfe8be3ec
8 8
9 D % non-stabalized scheme operator 9 D % non-stabalized scheme operator
10 A, B, E 10 A, B, E
11 11
12 H % Discrete norm 12 H % Discrete norm
13 Hi 13 % Norms in the x and y directions
14 H_x, H_y % Norms in the x and y directions 14 Hxi,Hyi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
15 Hx,Hy % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. 15 I_x,I_y, I_N
16 I_x,I_y
17 e_w, e_e, e_s, e_n 16 e_w, e_e, e_s, e_n
18 params %parameters for the coeficient matrices 17 params %parameters for the coeficient matrices
18 matrices
19 end 19 end
20 20
21 21
22 methods 22 methods
23 function obj = hypsyst2d(m,lim,order,matrices,params) 23 function obj = hypsyst2d(m,lim,order,matrices,params)
30 end 30 end
31 31
32 m_x = m(1); 32 m_x = m(1);
33 m_y = m(2); 33 m_y = m(2);
34 34
35 obj.matrices=matrices;
36
35 ops_x = sbp.D2Standard(m_x,xlim,order); 37 ops_x = sbp.D2Standard(m_x,xlim,order);
36 ops_y = sbp.D2Standard(m_y,ylim,order); 38 ops_y = sbp.D2Standard(m_y,ylim,order);
37 39
38 obj.x=ops_x.x; 40 obj.x=ops_x.x;
39 obj.y=ops_y.y; 41 obj.y=ops_y.x;
40 42
41 obj.X = kr(x,ones(m_y,1)); 43 obj.X = kr(obj.x,ones(m_y,1));
42 obj.Y = kr(ones(m_x,1),y); 44 obj.Y = kr(ones(m_x,1),obj.y);
43 45
44 I_x = speye(m_x); 46 I_x = speye(m_x);
45 I_y = speye(m_y); 47 I_y = speye(m_y);
48 obj.I_N=speye(m_x*m_y);
46 I_n= eye(4); 49 I_n= eye(4);
47 50
48 51
49 D1_x = kr(kr(I_n,ops_x.D1),I_y); 52 D1_x = kr(kr(I_n,ops_x.D1),I_y);
50 obj.Hi_x= kr(kr(I_n,ops_x.HI),I_y); 53 obj.Hxi= kr(kr(I_n,ops_x.HI),I_y);
51 D1_y=kr(I_n,kr(I_x,ops_y.D1)); 54 D1_y=kr(I_n,kr(I_x,ops_y.D1));
52 obj.Hi_y=kr(I_n,kr(I_x,ops_y.HI)); 55 obj.Hyi=kr(I_n,kr(I_x,ops_y.HI));
53 56
54 obj.e_w=kr(I_n,kr(ops_x.e_l,I_y)); 57 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)); 58 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)); 59 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)); 60 obj.e_n=kr(I_n,kr(I_x,ops_y.e_r));
59 obj.m=m; 62 obj.m=m;
60 obj.h=[ops_x.h ops_y.h]; 63 obj.h=[ops_x.h ops_y.h];
61 obj.order=order; 64 obj.order=order;
62 obj.params=params; 65 obj.params=params;
63 66
64 obj.A=matrixBuild(obj,matrices.A); 67 obj.A=obj.matrixBuild(matrices.A);
65 obj.B=matrixBuild(obj,matrices.B); 68 obj.B=obj.matrixBuild(matrices.B);
66 obj.E=matrixBuild(obj,matrices.E); 69 obj.E=obj.matrixBuild(matrices.E);
67 70
68 obj.D=-obj.A*D1_x-obj.B*D1_y-E; 71 A=kron(matrices.A(params,0,0),obj.I_N);
69 72 B=kron(matrices.B(params,0,0),obj.I_N);
73
74 obj.D=-obj.A*D1_x-obj.B*D1_y-obj.E;
75
70 end 76 end
71 % Closure functions return the opertors applied to the own doamin to close the boundary 77 % 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. 78 % 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'. 79 % 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. 80 % 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. 81 % 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. 82 % neighbour_scheme is an instance of Scheme that should be interfaced to.
77 % neighbour_boundary is a string specifying which boundary to interface to. 83 % neighbour_boundary is a string specifying which boundary to interface to.
78 function [closure, penalty] = boundary_condition(obj,boundary,type,data) 84 function [closure, penalty] = boundary_condition(obj,boundary,type)
79 default_arg('type','neumann'); 85 default_arg('type','neumann');
80 default_arg('data',0); 86 default_arg('data',0);
81 87
82 switch type 88 switch type
83 case{c,'char'} 89 case{'c','char'}
84 [tau,e_,Hi,CHM]=GetBoundarydata(obj,boundary,type); 90 [closure,penalty]=GetBoundarydata_char_fel(obj,boundary);
85 closure =Hi*e_*tau*CHM*e_';
86 penalty =Hi*e_*tau*CHM*data;
87 otherwise 91 otherwise
88 error('No such boundary condition') 92 error('No such boundary condition')
89 end 93 end
90 end 94 end
91 95
95 99
96 function N = size(obj) 100 function N = size(obj)
97 N = obj.m; 101 N = obj.m;
98 end 102 end
99 103
100 end 104 function [ret]=matrixBuild(obj,mat)
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 105 %extra info for coordinate transfomration mult my y_ny and
113 %x,ny osv... 106 %x,ny osv...
114 params=obj.params; 107 params=obj.params;
115 X=obj.X; 108 X=obj.X;
116 Y=obj.Y; 109 Y=obj.Y;
117 110
118 if isa(mat,'function_handle') 111 if isa(mat,'function_handle')
119 matVec=mat(params,x,y); 112 matVec=mat(params,X',Y');
120 side=length(x); 113 side=length(X');
121 else 114 else
122 matVec=mat; 115 matVec=mat;
123 side=max(size(mat))/4; 116 side=max(size(mat))/4;
124 end 117 end
118 matVec=sparse(matVec);
125 119
126 120
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)) 121 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)) 122 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)) 123 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))]; 124 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 125 end
132 126 function [closure, penalty]=GetBoundarydata_char(obj,boundary)
133 function [tau,e_,Hi, CHM]=GetBoundarydata(obj,boundary)
134 params=obj.params; 127 params=obj.params;
135 x=obj.x; 128 x=obj.x;
136 y=obj.y; 129 y=obj.y;
137 130
138 side=max(length(x),length(y));
139
140
141 switch boundary 131 switch boundary
142 case {'w','W','west'} 132 case {'w','W','west'}
143 e_=obj.e_w; 133 e_=obj.e_w;
144 mat=obj.A; 134 mat=obj.matrices.A;
145 [V,D]=matrixDiag(mat,params,x(1),y); 135 [V,Vi,D,pos]=obj.matrixDiag(mat,0,0);
146 Hi=obj.Hx; 136 Hi=obj.Hxi;
147 tau=zeros(4*side,pos*side); 137 mat_plus=V*(D+abs(D))*Vi/2;
148 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side)); 138
149 CHM=V*((D+abs(D))/2)*V'; 139 mat_plus=e_'*kron(mat_plus,obj.I_N)*e_;
140
141 tau=-1;
142 closure=Hi*e_*tau*mat_plus*e_';
143 penalty=Hi*e_*tau*mat_plus;
144
150 case {'e','E','east'} 145 case {'e','E','east'}
151 e_=obj.e_e; 146 e_=obj.e_e;
152 mat=obj.A; 147 mat=obj.matrices.A;
153 [V,D]=matrixDiag(mat,params,x(end),y); 148 [V,Vi,D,pos]=obj.matrixDiag(mat,0,0);
154 Hi=obj.Hy; 149 Hi=obj.Hxi;
155 tau=zeros(4*side,(4-pos)); 150 tau=1;
156 tau(pos*side+1:4*side)=-abs(D(pos*side+1:4*side,pos*side+1:4*side)); 151 mat_minus=V*(D-abs(D))*Vi/2;
157 CHM=V*((D-abs(D))/2)*V'; 152 mat_minus=e_'*kron(mat_minus,obj.I_N)*e_;
153 closure=Hi*e_*tau*mat_minus*e_';
154 penalty=Hi*e_*tau*mat_minus;
155
158 case {'s','S','south'} 156 case {'s','S','south'}
159 e_=obj.e_s; 157 e_=obj.e_s;
160 mat=obj.B; 158 mat=obj.matrices.B;
161 [V,D]=matrixDiag(mat,params,x,y(1)); 159 [V,Vi,D,pos]=obj.matrixDiag(mat,0,0);
162 Hi=obj.Hx; 160 Hi=obj.Hyi;
163 tau=zeros(4*side,pos*side); 161 mat_plus=V*(D+abs(D))*Vi/2;
164 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side)); 162
165 CHM=V*((D+abs(D))/2)*V'; 163 mat_plus=e_'*kron(mat_plus,obj.I_N)*e_;
164 tau=-1;
165 closure=Hi*e_*tau*mat_plus*e_';
166 penalty=Hi*e_*tau*mat_plus;
167
166 case {'n','N','north'} 168 case {'n','N','north'}
167 e_=obj.e_n; 169 e_=obj.e_n;
168 mat=obj.B; 170 mat=obj.matrices.B;
169 [V,D]=matrixDiag(mat,params,x,y(end)); 171 [V,Vi,D,pos]=obj.matrixDiag(mat,0,0);
170 Hi=obj.Hy; 172 Hi=obj.Hyi;
171 tau=zeros(4*side,(4-pos)); 173 tau=1;
172 tau(pos*side+1:4*side)=-abs(D(pos*side+1:4*side,pos*side+1:4*side)); 174 mat_minus=V*(D-abs(D))*Vi/2;
173 CHM=V*((D-abs(D))/2)*V'; 175
174 end 176 mat_minus=e_'*kron(mat_minus,obj.I_N)*e_;
175 177 closure=Hi*e_*tau*mat_minus*e_';
176 tau=V*tau*V'; 178 penalty=Hi*e_*tau*mat_minus;
177 179
178 end 180 end
179 181 end
180 function [V, D,pos]=matrixDiag(mat,params,x,y) 182
183
184
185 function [closure, penalty]=GetBoundarydata_charfel2(obj,boundary)
186 params=obj.params;
187 x=obj.x;
188 y=obj.y;
189
190 switch boundary
191 case {'w','W','west'}
192 e_=obj.e_w;
193 mat=obj.matrices.A;
194 [V,Vi,D,pos]=obj.matrixDiag(mat,x(1),y);
195 Hi=obj.Hxi;
196 mat_plus=V*(D+abs(D))*Vi/2;
197
198 tau=-1;
199 closure=Hi*e_*tau*mat_plus*e_';
200 penalty=Hi*e_*tau*mat_plus;
201
202 case {'e','E','east'}
203 e_=obj.e_e;j
204 mat=obj.matrices.A;
205 [V,Vi,D,pos]=obj.matrixDiag(mat,x(end),y);
206 Hi=obj.Hxi;
207 tau=1;
208
209 closure=Hi*e_*tau*mat_minus*e_';
210 penalty=Hi*e_*tau*mat_minus;
211
212 case {'s','S','south'}
213 e_=obj.e_s;
214 mat=obj.matrices.B;
215 [V,Vi,D,pos]=obj.matrixDiag(mat,x,y(1));
216 Hi=obj.Hyi;
217
218 tau=-1;
219 closure=Hi*e_*tau*mat_plus*e_';
220 penalty=Hi*e_*tau*mat_plus;
221
222 case {'n','N','north'}
223 e_=obj.e_n;
224 mat=obj.matrices.B;
225 [V,Vi,D,pos]=obj.matrixDiag(mat,x,y(end));
226 Hi=obj.Hyi;
227 tau=1;
228
229 closure=Hi*e_*tau*mat_minus*e_';
230 penalty=Hi*e_*tau*mat_minus;
231
232 end
233 end
234
235
236 function [closure, penalty]=GetBoundarydata_char_fel(obj,boundary)
237 params=obj.params;
238 x=obj.x;
239 y=obj.y;
240
241 side=max(length(x),length(y));
242
243
244 switch boundary
245 case {'w','W','west'}
246 e_=obj.e_w;
247 mat=obj.matrices.A;
248 [V,Vi,D,pos]=obj.matrixDiag(mat,x(1),y);
249 Hi=obj.Hxi;
250 tau=sparse(4*side,pos*side);
251 V_plus=V(:,1:pos*side);
252 Vi_plus=Vi(1:pos*side,:);
253 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side));
254 mat_plus=V*(D+abs(D))*Vi/2;
255
256 closure=Hi*e_*V*tau*Vi_plus*e_';
257 penalty=Hi*e_*V*tau*Vi_plus;
258
259 case {'e','E','east'}
260 e_=obj.e_e;
261 mat=obj.matrices.A;
262 [V,Vi,D,pos]=obj.matrixDiag(mat,x(end),y);
263 Hi=obj.Hxi;
264 tau=sparse(4*side,(4-pos)*side);
265 tau(pos*side+1:4*side,:)=-abs(D(pos*side+1:4*side,pos*side+1:4*side));
266 Vi_minus=Vi(pos*side+1:4*side,:);
267 mat_minus=V*(D-abs(D))*Vi/2;
268
269 closure=Hi*e_*V*tau*Vi_minus*e_';
270 penalty=Hi*e_*V*tau*Vi_minus;
271
272 case {'s','S','south'}
273 e_=obj.e_s;
274 mat=obj.matrices.B;
275 [V,Vi,D,pos]=obj.matrixDiag(mat,x,y(1));
276 Hi=obj.Hyi;
277 tau=sparse(4*side,pos*side);
278 V_plus=V(:,1:pos*side);
279 Vi_plus=Vi(1:pos*side,:);
280 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side));
281 mat_plus=V*(D+abs(D))*Vi/2;
282
283 closure=Hi*e_*V*tau*Vi_plus*e_';
284 penalty=Hi*e_*V*tau*Vi_plus;
285 case {'n','N','north'}
286 e_=obj.e_n;
287 mat=obj.matrices.B;
288 [V,Vi,D,pos]=obj.matrixDiag(mat,x,y(end));
289 Hi=obj.Hyi;
290 tau=sparse(4*side,(4-pos)*side);
291 tau(pos*side+1:4*side,:)=-abs(D(pos*side+1:4*side,pos*side+1:4*side));
292 Vi_minus=Vi(pos*side+1:4*side,:);
293 mat_minus=V*(D-abs(D))*Vi/2;
294
295 closure=Hi*e_*V*tau*Vi_minus*e_';
296 penalty=Hi*e_*V*tau*Vi_minus;
297 end
298 end
299
300 function [V,Vi, D,pos]=matrixDiag(obj,mat,x,y)
301 params=obj.params;
181 syms xs ys; 302 syms xs ys;
182 [V, D]=eig(mat(params,xs,ys)); 303 [V, D]=eig(mat(params,xs,ys));
183 xs=1;ys=1; 304 xs=1;ys=1;
184 DD=eval(diag(D)); 305 DD=eval(diag(D));
185 306
188 syms xs ys 309 syms xs ys
189 DD=diag(D); 310 DD=diag(D);
190 311
191 D=diag([DD(pos); DD(neg)]); 312 D=diag([DD(pos); DD(neg)]);
192 V=[V(:,pos) V(:,neg)]; 313 V=[V(:,pos) V(:,neg)];
193 314
194 xs=x; ys=y; 315 xs=x; ys=y;
316
195 317
196 side=max(length(x),length(y)); 318 side=max(length(x),length(y));
197 Dret=zeros(4,side*4); 319 Dret=zeros(4,side*4);
198 Vret=zeros(4,side*4); 320 Vret=zeros(4,side*4);
199 for ii=1:4 321 for ii=1:4
200 for jj=1:4 322 for jj=1:4
201 Dret(jj,(ii-1)*side+1:side*ii)=eval(D(jj,ii)); 323 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)); 324 Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii));
203 end 325 end
204 end 326 end
327
328 D=sparse(Dret);
205 V=sparse(normc(Vret)); 329 V=sparse(normc(Vret));
206 D=sparse(Dret); 330 V=obj.matrixBuild(V);
207 331 D=obj.matrixBuild(D);
208 332 Vi=inv(V);
209 V=matrixBuild([],[],[],V); 333 pos=length(pos);
210 D=matrixBuild([],[],[],D); 334 end
211 pos=legth(pos); 335
212 end 336 end
337
338 methods(Static)
339 % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u
340 % and bound_v of scheme schm_v.
341 % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l')
342 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
343 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
344 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
345 end
346
347
213 end 348 end
214 end 349 end