Mercurial > repos > public > sbplib

comparison +scheme/hypsyst2d.m @ 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
comparison
equal deleted inserted replaced
292:3d275c5e45b3 293:2d604d16842c
1 classdef hypsyst2d < scheme.Scheme 1 classdef hypsyst2d < scheme.Scheme
2 properties 2 properties
3 m % Number of points in each direction, possibly a vector 3 m % Number of points in each direction, possibly a vector
4 n %size of system
4 h % Grid spacing 5 h % Grid spacing
5 x,y % Grid 6 x,y % Grid
6 X,Y % Values of x and y for each grid point 7 X,Y % Values of x and y for each grid point
7 order % Order accuracy for the approximation 8 order % Order accuracy for the approximation
8 9
29 m = [m m]; 30 m = [m m];
30 end 31 end
31 32
32 m_x = m(1); 33 m_x = m(1);
33 m_y = m(2); 34 m_y = m(2);
35 obj.params=params;
34 36
35 obj.matrices=matrices; 37 obj.matrices=matrices;
36 38
37 ops_x = sbp.D2Standard(m_x,xlim,order); 39 ops_x = sbp.D2Standard(m_x,xlim,order);
38 ops_y = sbp.D2Standard(m_y,ylim,order); 40 ops_y = sbp.D2Standard(m_y,ylim,order);
41 obj.y=ops_y.x; 43 obj.y=ops_y.x;
42 44
43 obj.X = kr(obj.x,ones(m_y,1)); 45 obj.X = kr(obj.x,ones(m_y,1));
44 obj.Y = kr(ones(m_x,1),obj.y); 46 obj.Y = kr(ones(m_x,1),obj.y);
45 47
48 obj.A=obj.matrixBuild(matrices.A);
49 obj.B=obj.matrixBuild(matrices.B);
50 obj.E=obj.matrixBuild(matrices.E);
51
52 obj.n=length(matrices.A(obj.params,0,0));
53
54 I_n= eye(obj.n);
46 I_x = speye(m_x); obj.I_x=I_x; 55 I_x = speye(m_x); obj.I_x=I_x;
47 I_y = speye(m_y); obj.I_y=I_y; 56 I_y = speye(m_y); obj.I_y=I_y;
48 57
49 I_n= eye(4); 58
50
51
52 D1_x = kr(kr(I_n,ops_x.D1),I_y); 59 D1_x = kr(kr(I_n,ops_x.D1),I_y);
53 obj.Hxi= kr(kr(I_n,ops_x.HI),I_y); 60 obj.Hxi= kr(kr(I_n,ops_x.HI),I_y);
54 D1_y=kr(I_n,kr(I_x,ops_y.D1)); 61 D1_y=kr(I_n,kr(I_x,ops_y.D1));
55 obj.Hyi=kr(I_n,kr(I_x,ops_y.HI)); 62 obj.Hyi=kr(I_n,kr(I_x,ops_y.HI));
56 63
60 obj.e_n=kr(I_n,kr(I_x,ops_y.e_r)); 67 obj.e_n=kr(I_n,kr(I_x,ops_y.e_r));
61 68
62 obj.m=m; 69 obj.m=m;
63 obj.h=[ops_x.h ops_y.h]; 70 obj.h=[ops_x.h ops_y.h];
64 obj.order=order; 71 obj.order=order;
65 obj.params=params; 72
66
67 obj.A=obj.matrixBuild(matrices.A);
68 obj.B=obj.matrixBuild(matrices.B);
69 obj.E=obj.matrixBuild(matrices.E);
70
71
72 obj.D=-obj.A*D1_x-obj.B*D1_y-obj.E; 73 obj.D=-obj.A*D1_x-obj.B*D1_y-obj.E;
73 74
74 end 75 end
75 % Closure functions return the opertors applied to the own doamin to close the boundary 76 % Closure functions return the opertors applied to the own doamin to close the boundary
76 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. 77 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
77 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 78 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
78 % type is a string specifying the type of boundary condition if there are several. 79 % type is a string specifying the type of boundary condition if there are several.
79 % data is a function returning the data that should be applied at the boundary. 80 % data is a function returning the data that should be applied at the boundary.
80 % neighbour_scheme is an instance of Scheme that should be interfaced to. 81 % neighbour_scheme is an instance of Scheme that should be interfaced to.
81 % neighbour_boundary is a string specifying which boundary to interface to. 82 % neighbour_boundary is a string specifying which boundary to interface to.
82 function [closure, penalty] = boundary_condition(obj,boundary,type) 83 function [closure, penalty] = boundary_condition(obj,boundary,type,L)
83 default_arg('type','neumann'); 84 default_arg('type','neumann');
84 default_arg('data',0); 85 default_arg('data',0);
85 86
86 switch type 87 switch type
87 case{'c','char'} 88 case{'c','char'}
88 [closure,penalty]=GetBoundarydata_char(obj,boundary); 89 [closure,penalty]=GetBoundarydata_char(obj,boundary);
89 case{'wall'} 90 case{'general'}
90 [closure,penalty]=GetBoundarydata_wall(obj,boundary); 91 [closure,penalty]=GeneralBoundaryCond(obj,boundary,L);
91 otherwise 92 otherwise
92 error('No such boundary condition') 93 error('No such boundary condition')
93 end 94 end
94 end 95 end
95 96
168 169
169 pos=signVec(1); zeroval=signVec(2); neg=signVec(3); 170 pos=signVec(1); zeroval=signVec(2); neg=signVec(3);
170 171
171 switch boundPos 172 switch boundPos
172 case {'l'} 173 case {'l'}
173 tau=sparse(4*side,pos*side); 174 tau=sparse(obj.n*side,pos*side);
174 Vi_plus=Vi(1:pos*side,:); 175 Vi_plus=Vi(1:pos*side,:);
175 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side)); 176 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side));
176 closure=Hi*e_*V*tau*Vi_plus*e_'; 177 closure=Hi*e_*V*tau*Vi_plus*e_';
177 penalty=-Hi*e_*V*tau*Vi_plus; 178 penalty=-Hi*e_*V*tau*Vi_plus;
178 179
179 case {'r'} 180 case {'r'}
180 tau=sparse(4*side,neg*side); 181 tau=sparse(obj.n*side,neg*side);
181 tau((pos+zeroval)*side+1:4*side,:)=-abs(D((pos+zeroval)*side+1:4*side,(pos+zeroval)*side+1:4*side)); 182 tau((pos+zeroval)*side+1:obj.n*side,:)=-abs(D((pos+zeroval)*side+1:obj.n*side,(pos+zeroval)*side+1:obj.n*side));
182 Vi_minus=Vi((pos+zeroval)*side+1:4*side,:); 183 Vi_minus=Vi((pos+zeroval)*side+1:obj.n*side,:);
183 closure=Hi*e_*V*tau*Vi_minus*e_'; 184 closure=Hi*e_*V*tau*Vi_minus*e_';
184 penalty=-Hi*e_*V*tau*Vi_minus; 185 penalty=-Hi*e_*V*tau*Vi_minus;
185 186
186 end 187 end
187 end
188
189 function [closure, penalty]=GetBoundarydata_wall(obj,boundary)
190 switch boundary
191 case {'e','w'}
192 L=[0 1 0 0]';
193 L=kr(L,obj.I_y);
194 L=L';
195 case {'s','n'}
196 L=[0 0 1 0]';
197 L=kr(L,obj.I_x);
198 L=L';
199
200 end
201 [closure,penalty]=GeneralBoundaryCond(obj,boundary,L);
202 end 188 end
203 189
204 190
205 function [closure,penalty]=GeneralBoundaryCond(obj,boundary,L) 191 function [closure,penalty]=GeneralBoundaryCond(obj,boundary,L)
206 params=obj.params; 192 params=obj.params;
207 x=obj.x; 193 x=obj.x;
208 y=obj.y; 194 y=obj.y;
209 195 L=obj.matrixBuild(L,x,y);
210 side=max(length(x),length(y)); 196 side=max(length(x),length(y));
211 197
212 198
213 switch boundary 199 switch boundary
214 case {'w','W','west'} 200 case {'w','W','west'}
241 227
242 228
243 229
244 switch boundPos 230 switch boundPos
245 case {'l'} 231 case {'l'}
246 tau=sparse(4*side,pos*side); 232 tau=sparse(obj.n*side,pos*side);
247 Vi_plus=Vi(1:pos*side,:); 233 Vi_plus=Vi(1:pos*side,:);
248 Vi_minus=Vi(pos*side+1:4*side,:); 234 Vi_minus=Vi(pos*side+1:obj.n*side,:);
249 235
250 V_plus=Vi(:,1:pos*side); 236 V_plus=Vi(:,1:pos*side);
251 V_minus=Vi(:,(pos+zeroval)*side+1:4*side); 237 V_minus=Vi(:,(pos+zeroval)*side+1:obj.n*side);
252 238
253 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side)); 239 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side));
254 R=-inv(L*V_plus)*(L*V_minus); 240 R=-inv(L*V_plus)*(L*V_minus);
255 closure=Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; 241 closure=Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
256 penalty=-Hi*e_*V*tau*inv(L*V_plus)*L; 242 penalty=-Hi*e_*V*tau*inv(L*V_plus)*L;
257 243
258 244
259 case {'r'} 245 case {'r'}
260 tau=sparse(4*side,neg*side); 246 tau=sparse(obj.n*side,neg*side);
261 tau((pos+zeroval)*side+1:4*side,:)=-abs(D((pos+zeroval)*side+1:4*side,(pos+zeroval)*side+1:4*side)); 247 tau((pos+zeroval)*side+1:obj.n*side,:)=-abs(D((pos+zeroval)*side+1:obj.n*side,(pos+zeroval)*side+1:obj.n*side));
262 Vi_plus=Vi(1:pos*side,:); 248 Vi_plus=Vi(1:pos*side,:);
263 Vi_minus=Vi((pos+zeroval)*side+1:4*side,:); 249 Vi_minus=Vi((pos+zeroval)*side+1:obj.n*side,:);
264 250
265 V_plus=Vi(:,1:pos*side); 251 V_plus=Vi(:,1:pos*side);
266 V_minus=Vi(:,(pos+zeroval)*side+1:4*side); 252 V_minus=Vi(:,(pos+zeroval)*side+1:obj.n*side);
267 R=-inv(L*V_minus)*(L*V_plus); 253 R=-inv(L*V_minus)*(L*V_plus);
268 closure=Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_'; 254 closure=Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
269 penalty=-Hi*e_*V*tau*inv(L*V_minus)*L; 255 penalty=-Hi*e_*V*tau*inv(L*V_minus)*L;
270 256
271 257
291 277
292 xs=x; ys=y; 278 xs=x; ys=y;
293 279
294 280
295 side=max(length(x),length(y)); 281 side=max(length(x),length(y));
296 Dret=zeros(4,side*4); 282 Dret=zeros(obj.n,side*obj.n);
297 Vret=zeros(4,side*4); 283 Vret=zeros(obj.n,side*obj.n);
298 for ii=1:4 284 for ii=1:obj.n
299 for jj=1:4 285 for jj=1:obj.n
300 Dret(jj,(ii-1)*side+1:side*ii)=eval(D(jj,ii)); 286 Dret(jj,(ii-1)*side+1:side*ii)=eval(D(jj,ii));
301 Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii)); 287 Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii));
302 end 288 end
303 end 289 end
304 290