Mercurial > repos > public > sbplib

comparison +scheme/hypsyst2d.m @ 295:da0131655035 feature/hypsyst

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fixed some formatting and naming.
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 26 Sep 2016 14:19:34 +0200
parents 8ff6ec6249e8
children
comparison
equal deleted inserted replaced
294:8ff6ec6249e8 295:da0131655035
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 n %size of system
5 h % Grid spacing 5 h % Grid spacing
6 x,y % Grid 6 x,y % Grid
7 X,Y % Values of x and y for each grid point 7 X,Y % Values of x and y for each grid point
8 order % Order accuracy for the approximation 8 order % Order accuracy for the approximation
9 9
10 D % non-stabalized scheme operator 10 D % non-stabalized scheme operator
11 A, B, E 11 A, B, E
12 12
13 H % Discrete norm 13 H % Discrete norm
14 % Norms in the x and y directions 14 % Norms in the x and y directions
15 Hxi,Hyi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. 15 Hxi,Hyi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
16 I_x,I_y, I_N 16 I_x,I_y, I_N
17 e_w, e_e, e_s, e_n 17 e_w, e_e, e_s, e_n
18 params %parameters for the coeficient matrices 18 params %parameters for the coeficient matrices
19 matrices 19 matrices
20 end 20 end
21 21
22 22
23 methods 23 methods
24 function obj = hypsyst2d(m,lim,order,matrices,params) 24 function obj = Hypsyst2d(m, lim, order, A, B, E, params)
25 25 default_arg('E', [])
26 xlim = lim{1}; 26 xlim = lim{1};
27 ylim = lim{2}; 27 ylim = lim{2};
28 28
29 if length(m) == 1 29 if length(m) == 1
30 m = [m m]; 30 m = [m m];
31 end 31 end
32 32
33 m_x = m(1); 33 m_x = m(1);
34 m_y = m(2); 34 m_y = m(2);
35 obj.params=params; 35 obj.params = params;
36 36
37 obj.matrices=matrices; 37 obj.matrices = matrices;
38 38
39 ops_x = sbp.D2Standard(m_x,xlim,order); 39 ops_x = sbp.D2Standard(m_x,xlim,order);
40 ops_y = sbp.D2Standard(m_y,ylim,order); 40 ops_y = sbp.D2Standard(m_y,ylim,order);
41 41
42 obj.x=ops_x.x; 42 obj.x = ops_x.x;
43 obj.y=ops_y.x; 43 obj.y = ops_y.x;
44 44
45 obj.X = kr(obj.x,ones(m_y,1)); 45 obj.X = kr(obj.x,ones(m_y,1));
46 obj.Y = kr(ones(m_x,1),obj.y); 46 obj.Y = kr(ones(m_x,1),obj.y);
47 47
48 obj.A=obj.matrixBuild(matrices.A); 48 obj.A = obj.evaluateCoefficientMatrix(matrices.A, obj.X, obj.Y);
49 obj.B=obj.matrixBuild(matrices.B); 49 obj.B = obj.evaluateCoefficientMatrix(matrices.B, obj.X, obj.Y);
50 obj.E=obj.matrixBuild(matrices.E); 50 obj.E = obj.evaluateCoefficientMatrix(matrices.E, obj.X, obj.Y);
51 51
52 obj.n=length(matrices.A(obj.params,0,0)); 52 obj.n = length(matrices.A(obj.params,0,0));
53 53
54 I_n= eye(obj.n); 54 I_n = eye(obj.n);I_x = speye(m_x);
55 I_x = speye(m_x); obj.I_x=I_x; 55 obj.I_x = I_x;
56 I_y = speye(m_y); obj.I_y=I_y; 56 I_y = speye(m_y);
57 57 obj.I_y = I_y;
58 58
59 D1_x = kr(kr(I_n,ops_x.D1),I_y); 59
60 obj.Hxi= kr(kr(I_n,ops_x.HI),I_y); 60 D1_x = kr(I_n, ops_x.D1, I_y);
61 D1_y=kr(I_n,kr(I_x,ops_y.D1)); 61 obj.Hxi = kr(I_n, ops_x.HI, I_y);
62 obj.Hyi=kr(I_n,kr(I_x,ops_y.HI)); 62 D1_y = kr(I_n, I_x, ops_y.D1));
63 63 obj.Hyi = kr(I_n, I_x, ops_y.HI));
64 obj.e_w=kr(I_n,kr(ops_x.e_l,I_y)); 64
65 obj.e_e=kr(I_n,kr(ops_x.e_r,I_y)); 65 obj.e_w = kr(I_n, ops_x.e_l, I_y);
66 obj.e_s=kr(I_n,kr(I_x,ops_y.e_l)); 66 obj.e_e = kr(I_n, ops_x.e_r, I_y);
67 obj.e_n=kr(I_n,kr(I_x,ops_y.e_r)); 67 obj.e_s = kr(I_n, I_x, ops_y.e_l);
68 68 obj.e_n = kr(I_n, I_x, ops_y.e_r);
69
69 obj.m=m; 70 obj.m=m;
70 obj.h=[ops_x.h ops_y.h]; 71 obj.h=[ops_x.h ops_y.h];
71 obj.order=order; 72 obj.order=order;
72 73
73 obj.D=-obj.A*D1_x-obj.B*D1_y-obj.E; 74 obj.D=-obj.A*D1_x-obj.B*D1_y-obj.E;
74 75
75 end 76 end
76 77
77 % Closure functions return the opertors applied to the own doamin to close the boundary 78 % Closure functions return the opertors applied to the own doamin to close the boundary
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. 79 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
79 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 80 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
80 % type is a string specifying the type of boundary condition if there are several. 81 % type is a string specifying the type of boundary condition if there are several.
81 % data is a function returning the data that should be applied at the boundary. 82 % data is a function returning the data that should be applied at the boundary.
82 % neighbour_scheme is an instance of Scheme that should be interfaced to.
83 % neighbour_boundary is a string specifying which boundary to interface to.
84 function [closure, penalty] = boundary_condition(obj,boundary,type,L) 83 function [closure, penalty] = boundary_condition(obj,boundary,type,L)
85 default_arg('type','char'); 84 default_arg('type','char');
86 switch type 85 switch type
87 case{'c','char'} 86 case{'c','char'}
88 [closure,penalty]=GetBoundarydata_char(obj,boundary); 87 [closure,penalty]=boundary_condition_char(obj,boundary);
89 case{'general'} 88 case{'general'}
90 [closure,penalty]=GeneralBoundaryCond(obj,boundary,L); 89 [closure,penalty]=boundary_condition_general(obj,boundary,L);
91 otherwise 90 otherwise
92 error('No such boundary condition') 91 error('No such boundary condition')
93 end 92 end
94 end 93 end
95 94
96 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 95 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
97 error('An interface function does not exist yet'); 96 error('An interface function does not exist yet');
98 end 97 end
99 98
100 function N = size(obj) 99 function N = size(obj)
101 N = obj.m; 100 N = obj.m;
102 end 101 end
103 102
104 function [ret]=matrixBuild(obj,mat,X,Y) 103 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y)
105 params=obj.params; 104 params=obj.params;
106 default_arg('X',obj.X); 105
107 default_arg('Y',obj.Y)
108
109 if isa(mat,'function_handle') 106 if isa(mat,'function_handle')
110 [rows,cols]=size(mat(params,0,0)); 107 [rows,cols]=size(mat(params,0,0));
111 matVec=mat(params,X',Y'); 108 matVec=mat(params,X',Y');
112 matVec=sparse(matVec); 109 matVec=sparse(matVec);
113 side=max(length(X),length(Y)); 110 side=max(length(X),length(Y));
116 [rows,cols]=size(matVec); 113 [rows,cols]=size(matVec);
117 side=max(length(X),length(Y)); 114 side=max(length(X),length(Y));
118 cols=cols/side; 115 cols=cols/side;
119 end 116 end
120 ret=kron(ones(rows,cols),speye(side)); 117 ret=kron(ones(rows,cols),speye(side));
121 118
122 for ii=1:rows 119 for ii=1:rows
123 for jj=1:cols 120 for jj=1:cols
124 ret((ii-1)*side+1:ii*side,(jj-1)*side+1:jj*side)=diag(matVec(ii,(jj-1)*side+1:jj*side)); 121 ret((ii-1)*side+1:ii*side,(jj-1)*side+1:jj*side)=diag(matVec(ii,(jj-1)*side+1:jj*side));
125 end 122 end
126 end 123 end
127 end 124 end
128 125
129 126
130 function [closure, penalty]=GetBoundarydata_char(obj,boundary) 127 function [closure, penalty]=boundary_condition_char(obj,boundary)
131 params=obj.params; 128 params=obj.params;
132 x=obj.x; y=obj.y; 129 x=obj.x; y=obj.y;
133 side=max(length(x),length(y)); 130 side=max(length(x),length(y));
134 131
135 switch boundary 132 switch boundary
136 case {'w','W','west'} 133 case {'w','W','west'}
137 e_=obj.e_w; 134 e_=obj.e_w;
138 mat=obj.matrices.A; 135 mat=obj.matrices.A;
139 boundPos='l'; 136 boundPos='l';
156 mat=obj.matrices.B; 153 mat=obj.matrices.B;
157 boundPos='r'; 154 boundPos='r';
158 Hi=obj.Hxi; 155 Hi=obj.Hxi;
159 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(end)); 156 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(end));
160 end 157 end
161 158
162 pos=signVec(1); zeroval=signVec(2); neg=signVec(3); 159 pos=signVec(1); zeroval=signVec(2); neg=signVec(3);
163 160
164 switch boundPos 161 switch boundPos
165 case {'l'} 162 case {'l'}
166 tau=sparse(obj.n*side,pos*side); 163 tau=sparse(obj.n*side,pos*side);
167 Vi_plus=Vi(1:pos*side,:); 164 Vi_plus=Vi(1:pos*side,:);
168 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side)); 165 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side));
174 Vi_minus=Vi((pos+zeroval)*side+1:obj.n*side,:); 171 Vi_minus=Vi((pos+zeroval)*side+1:obj.n*side,:);
175 closure=Hi*e_*V*tau*Vi_minus*e_'; 172 closure=Hi*e_*V*tau*Vi_minus*e_';
176 penalty=-Hi*e_*V*tau*Vi_minus; 173 penalty=-Hi*e_*V*tau*Vi_minus;
177 end 174 end
178 end 175 end
179 176
180 177
181 function [closure,penalty]=GeneralBoundaryCond(obj,boundary,L) 178 function [closure,penalty]=boundary_condition_general(obj,boundary,L)
182 params=obj.params; 179 params=obj.params;
183 x=obj.x; y=obj.y; 180 x=obj.x; y=obj.y;
184 side=max(length(x),length(y)); 181 side=max(length(x),length(y));
185 182
186 switch boundary 183 switch boundary
187 case {'w','W','west'} 184 case {'w','W','west'}
188 e_=obj.e_w; 185 e_=obj.e_w;
189 mat=obj.matrices.A; 186 mat=obj.matrices.A;
190 boundPos='l'; 187 boundPos='l';
191 Hi=obj.Hxi; 188 Hi=obj.Hxi;
192 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(1),y); 189 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(1),y);
193 L=obj.matrixBuild(L,x(1),y); 190 L=obj.evaluateCoefficientMatrix(L,x(1),y);
194 case {'e','E','east'} 191 case {'e','E','east'}
195 e_=obj.e_e; 192 e_=obj.e_e;
196 mat=obj.matrices.A; 193 mat=obj.matrices.A;
197 boundPos='r'; 194 boundPos='r';
198 Hi=obj.Hxi; 195 Hi=obj.Hxi;
199 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(end),y); 196 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(end),y);
200 L=obj.matrixBuild(L,x(end),y); 197 L=obj.evaluateCoefficientMatrix(L,x(end),y);
201 case {'s','S','south'} 198 case {'s','S','south'}
202 e_=obj.e_s; 199 e_=obj.e_s;
203 mat=obj.matrices.B; 200 mat=obj.matrices.B;
204 boundPos='l'; 201 boundPos='l';
205 Hi=obj.Hxi; 202 Hi=obj.Hxi;
206 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(1)); 203 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(1));
207 L=obj.matrixBuild(L,x,y(1)); 204 L=obj.evaluateCoefficientMatrix(L,x,y(1));
208 case {'n','N','north'} 205 case {'n','N','north'}
209 e_=obj.e_n; 206 e_=obj.e_n;
210 mat=obj.matrices.B; 207 mat=obj.matrices.B;
211 boundPos='r'; 208 boundPos='r';
212 Hi=obj.Hxi; 209 Hi=obj.Hxi;
213 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(end)); 210 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(end));
214 L=obj.matrixBuild(L,x,y(end)); 211 L=obj.evaluateCoefficientMatrix(L,x,y(end));
215 end 212 end
216 213
217 pos=signVec(1); zeroval=signVec(2); neg=signVec(3); 214 pos=signVec(1); zeroval=signVec(2); neg=signVec(3);
218 215
219 switch boundPos 216 switch boundPos
220 case {'l'} 217 case {'l'}
221 tau=sparse(obj.n*side,pos*side); 218 tau=sparse(obj.n*side,pos*side);
222 Vi_plus=Vi(1:pos*side,:); 219 Vi_plus=Vi(1:pos*side,:);
223 Vi_minus=Vi(pos*side+1:obj.n*side,:); 220 Vi_minus=Vi(pos*side+1:obj.n*side,:);
224 V_plus=V(:,1:pos*side); 221 V_plus=V(:,1:pos*side);
225 V_minus=V(:,(pos+zeroval)*side+1:obj.n*side); 222 V_minus=V(:,(pos+zeroval)*side+1:obj.n*side);
226 223
227 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side)); 224 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side));
228 R=-inv(L*V_plus)*(L*V_minus); 225 R=-inv(L*V_plus)*(L*V_minus);
229 closure=Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; 226 closure=Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
230 penalty=-Hi*e_*V*tau*inv(L*V_plus)*L; 227 penalty=-Hi*e_*V*tau*inv(L*V_plus)*L;
231 case {'r'} 228 case {'r'}
232 tau=sparse(obj.n*side,neg*side); 229 tau=sparse(obj.n*side,neg*side);
233 tau((pos+zeroval)*side+1:obj.n*side,:)=-abs(D((pos+zeroval)*side+1:obj.n*side,(pos+zeroval)*side+1:obj.n*side)); 230 tau((pos+zeroval)*side+1:obj.n*side,:)=-abs(D((pos+zeroval)*side+1:obj.n*side,(pos+zeroval)*side+1:obj.n*side));
234 Vi_plus=Vi(1:pos*side,:); 231 Vi_plus=Vi(1:pos*side,:);
235 Vi_minus=Vi((pos+zeroval)*side+1:obj.n*side,:); 232 Vi_minus=Vi((pos+zeroval)*side+1:obj.n*side,:);
236 233
237 V_plus=V(:,1:pos*side); 234 V_plus=V(:,1:pos*side);
238 V_minus=V(:,(pos+zeroval)*side+1:obj.n*side); 235 V_minus=V(:,(pos+zeroval)*side+1:obj.n*side);
239 R=-inv(L*V_minus)*(L*V_plus); 236 R=-inv(L*V_minus)*(L*V_plus);
240 closure=Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_'; 237 closure=Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
241 penalty=-Hi*e_*V*tau*inv(L*V_minus)*L; 238 penalty=-Hi*e_*V*tau*inv(L*V_minus)*L;
242 end 239 end
243 end 240 end
244 241
245 242
246 function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y) 243 function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y)
247 params=obj.params; 244 params=obj.params;
248 syms xs ys; 245 syms xs ys;
249 [V, D]=eig(mat(params,xs,ys)); 246 [V, D]=eig(mat(params,xs,ys));
250 xs=1;ys=1; 247 xs=1;ys=1;
251 DD=eval(diag(D)); 248 DD=eval(diag(D));
252 249
253 poseig=find(DD>0); 250 poseig=find(DD>0);
254 zeroeig=find(DD==0); 251 zeroeig=find(DD==0);
255 negeig=find(DD<0); 252 negeig=find(DD<0);
256 syms xs ys 253 syms xs ys
257 DD=diag(D); 254 DD=diag(D);
258 255
259 D=diag([DD(poseig);DD(zeroeig); DD(negeig)]); 256 D=diag([DD(poseig);DD(zeroeig); DD(negeig)]);
260 V=[V(:,poseig) V(:,zeroeig) V(:,negeig)]; 257 V=[V(:,poseig) V(:,zeroeig) V(:,negeig)];
261 xs=x; ys=y; 258 xs=x; ys=y;
262 259
263 side=max(length(x),length(y)); 260 side=max(length(x),length(y));
264 Dret=zeros(obj.n,side*obj.n); 261 Dret=zeros(obj.n,side*obj.n);
265 Vret=zeros(obj.n,side*obj.n); 262 Vret=zeros(obj.n,side*obj.n);
266 for ii=1:obj.n 263 for ii=1:obj.n
267 for jj=1:obj.n 264 for jj=1:obj.n
268 Dret(jj,(ii-1)*side+1:side*ii)=eval(D(jj,ii)); 265 Dret(jj,(ii-1)*side+1:side*ii)=eval(D(jj,ii));
269 Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii)); 266 Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii));
270 end 267 end
271 end 268 end
272 269
273 D=sparse(Dret); 270 D=sparse(Dret);
274 V=sparse(normc(Vret)); 271 V=sparse(normc(Vret));
275 V=obj.matrixBuild(V,x,y); 272 V=obj.evaluateCoefficientMatrix(V,x,y);
276 D=obj.matrixBuild(D,x,y); 273 D=obj.evaluateCoefficientMatrix(D,x,y);
277 Vi=inv(V); 274 Vi=inv(V);
278 signVec=[length(poseig),length(zeroeig),length(negeig)]; 275 signVec=[length(poseig),length(zeroeig),length(negeig)];
279 end 276 end
280 277
281 end
282
283 methods(Static)
284 % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u
285 % and bound_v of scheme schm_v.
286 % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l')
287 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
288 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
289 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
290 end
291
292
293 end 278 end
294 end 279 end