Mercurial > repos > public > sbplib

comparison +scheme/hypsyst2d.m @ 292:3d275c5e45b3 feature/hypsyst

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Changed how the matrices are built
author Ylva Rydin <ylva.rydin@telia.com>
date Fri, 23 Sep 2016 14:48:54 +0200
parents 807dfe8be3ec
children 2d604d16842c
comparison
equal deleted inserted replaced
291:807dfe8be3ec 292:3d275c5e45b3
41 obj.y=ops_y.x; 41 obj.y=ops_y.x;
42 42
43 obj.X = kr(obj.x,ones(m_y,1)); 43 obj.X = kr(obj.x,ones(m_y,1));
44 obj.Y = kr(ones(m_x,1),obj.y); 44 obj.Y = kr(ones(m_x,1),obj.y);
45 45
46 I_x = speye(m_x); 46 I_x = speye(m_x); obj.I_x=I_x;
47 I_y = speye(m_y); 47 I_y = speye(m_y); obj.I_y=I_y;
48 obj.I_N=speye(m_x*m_y); 48
49 I_n= eye(4); 49 I_n= eye(4);
50 50
51 51
52 D1_x = kr(kr(I_n,ops_x.D1),I_y); 52 D1_x = kr(kr(I_n,ops_x.D1),I_y);
53 obj.Hxi= kr(kr(I_n,ops_x.HI),I_y); 53 obj.Hxi= kr(kr(I_n,ops_x.HI),I_y);
66 66
67 obj.A=obj.matrixBuild(matrices.A); 67 obj.A=obj.matrixBuild(matrices.A);
68 obj.B=obj.matrixBuild(matrices.B); 68 obj.B=obj.matrixBuild(matrices.B);
69 obj.E=obj.matrixBuild(matrices.E); 69 obj.E=obj.matrixBuild(matrices.E);
70 70
71 A=kron(matrices.A(params,0,0),obj.I_N); 71
72 B=kron(matrices.B(params,0,0),obj.I_N); 72 obj.D=-obj.A*D1_x-obj.B*D1_y-obj.E;
73
74 obj.D=-obj.A*D1_x-obj.B*D1_y-obj.E;
75 73
76 end 74 end
77 % Closure functions return the opertors applied to the own doamin to close the boundary 75 % 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. 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.
79 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 77 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
85 default_arg('type','neumann'); 83 default_arg('type','neumann');
86 default_arg('data',0); 84 default_arg('data',0);
87 85
88 switch type 86 switch type
89 case{'c','char'} 87 case{'c','char'}
90 [closure,penalty]=GetBoundarydata_char_fel(obj,boundary); 88 [closure,penalty]=GetBoundarydata_char(obj,boundary);
89 case{'wall'}
90 [closure,penalty]=GetBoundarydata_wall(obj,boundary);
91 otherwise 91 otherwise
92 error('No such boundary condition') 92 error('No such boundary condition')
93 end 93 end
94 end 94 end
95 95
99 99
100 function N = size(obj) 100 function N = size(obj)
101 N = obj.m; 101 N = obj.m;
102 end 102 end
103 103
104 function [ret]=matrixBuild(obj,mat) 104 function [ret]=matrixBuild(obj,mat,x,y)
105 %extra info for coordinate transfomration mult my y_ny and 105 %extra info for coordinate transfomration mult my y_ny and
106 %x,ny osv... 106 %x,ny osv...
107 params=obj.params; 107 params=obj.params;
108 X=obj.X; 108 X=obj.X;
109 Y=obj.Y; 109 Y=obj.Y;
110 110
111 if isa(mat,'function_handle') 111 if isa(mat,'function_handle')
112 matVec=mat(params,X',Y'); 112 [rows,cols]=size(mat(params,0,0));
113 side=length(X'); 113 matVec=mat(params,X',Y');
114 else
115 matVec=mat;
116 side=max(size(mat))/4;
117 end
118 matVec=sparse(matVec); 114 matVec=sparse(matVec);
119 115 side=max(length(X),length(Y));
120 116 else
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)) 117 matVec=mat;
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)) 118 [rows,cols]=size(matVec);
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)) 119 side=max(length(x),length(y));
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))]; 120 cols=cols/side;
121 end
122
123
124 ret=kron(ones(rows,cols),speye(side));
125
126 for ii=1:rows
127 for jj=1:cols
128 ret((ii-1)*side+1:ii*side,(jj-1)*side+1:jj*side)=diag(matVec(ii,(jj-1)*side+1:jj*side));
129 end
130 end
131
125 end 132 end
126 function [closure, penalty]=GetBoundarydata_char(obj,boundary) 133
134
135 function [closure, penalty]=GetBoundarydata_char(obj,boundary)
127 params=obj.params; 136 params=obj.params;
128 x=obj.x; 137 x=obj.x;
129 y=obj.y; 138 y=obj.y;
130 139
140 side=max(length(x),length(y));
141
131 switch boundary 142 switch boundary
132 case {'w','W','west'} 143 case {'w','W','west'}
133 e_=obj.e_w; 144 e_=obj.e_w;
134 mat=obj.matrices.A; 145 mat=obj.matrices.A;
135 [V,Vi,D,pos]=obj.matrixDiag(mat,0,0); 146 boundPos='l';
136 Hi=obj.Hxi; 147 Hi=obj.Hxi;
137 mat_plus=V*(D+abs(D))*Vi/2; 148 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(1),y);
138
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
145 case {'e','E','east'} 149 case {'e','E','east'}
146 e_=obj.e_e; 150 e_=obj.e_e;
147 mat=obj.matrices.A; 151 mat=obj.matrices.A;
148 [V,Vi,D,pos]=obj.matrixDiag(mat,0,0); 152 boundPos='r';
149 Hi=obj.Hxi; 153 Hi=obj.Hxi;
150 tau=1; 154 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(end),y);
151 mat_minus=V*(D-abs(D))*Vi/2;
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
156 case {'s','S','south'} 155 case {'s','S','south'}
157 e_=obj.e_s; 156 e_=obj.e_s;
158 mat=obj.matrices.B; 157 mat=obj.matrices.B;
159 [V,Vi,D,pos]=obj.matrixDiag(mat,0,0); 158 boundPos='l';
160 Hi=obj.Hyi; 159 Hi=obj.Hxi;
161 mat_plus=V*(D+abs(D))*Vi/2; 160 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(1));
162
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
168 case {'n','N','north'} 161 case {'n','N','north'}
169 e_=obj.e_n; 162 e_=obj.e_n;
170 mat=obj.matrices.B; 163 mat=obj.matrices.B;
171 [V,Vi,D,pos]=obj.matrixDiag(mat,0,0); 164 boundPos='r';
172 Hi=obj.Hyi; 165 Hi=obj.Hxi;
173 tau=1; 166 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(end));
174 mat_minus=V*(D-abs(D))*Vi/2; 167 end
175 168
176 mat_minus=e_'*kron(mat_minus,obj.I_N)*e_; 169 pos=signVec(1); zeroval=signVec(2); neg=signVec(3);
177 closure=Hi*e_*tau*mat_minus*e_'; 170
178 penalty=Hi*e_*tau*mat_minus; 171 switch boundPos
179 172 case {'l'}
180 end 173 tau=sparse(4*side,pos*side);
181 end 174 Vi_plus=Vi(1:pos*side,:);
182 175 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side));
183 176 closure=Hi*e_*V*tau*Vi_plus*e_';
184 177 penalty=-Hi*e_*V*tau*Vi_plus;
185 function [closure, penalty]=GetBoundarydata_charfel2(obj,boundary) 178
186 params=obj.params; 179 case {'r'}
180 tau=sparse(4*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 Vi_minus=Vi((pos+zeroval)*side+1:4*side,:);
183 closure=Hi*e_*V*tau*Vi_minus*e_';
184 penalty=-Hi*e_*V*tau*Vi_minus;
185
186 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
203
204
205 function [closure,penalty]=GeneralBoundaryCond(obj,boundary,L)
206 params=obj.params;
187 x=obj.x; 207 x=obj.x;
188 y=obj.y; 208 y=obj.y;
189 209
210 side=max(length(x),length(y));
211
212
190 switch boundary 213 switch boundary
191 case {'w','W','west'} 214 case {'w','W','west'}
192 e_=obj.e_w; 215 e_=obj.e_w;
193 mat=obj.matrices.A; 216 mat=obj.matrices.A;
194 [V,Vi,D,pos]=obj.matrixDiag(mat,x(1),y); 217 boundPos='l';
195 Hi=obj.Hxi; 218 Hi=obj.Hxi;
196 mat_plus=V*(D+abs(D))*Vi/2; 219 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(1),y);
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'} 220 case {'e','E','east'}
203 e_=obj.e_e;j 221 e_=obj.e_e;
204 mat=obj.matrices.A; 222 mat=obj.matrices.A;
205 [V,Vi,D,pos]=obj.matrixDiag(mat,x(end),y); 223 boundPos='r';
206 Hi=obj.Hxi; 224 Hi=obj.Hxi;
207 tau=1; 225 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(end),y);
208
209 closure=Hi*e_*tau*mat_minus*e_';
210 penalty=Hi*e_*tau*mat_minus;
211
212 case {'s','S','south'} 226 case {'s','S','south'}
213 e_=obj.e_s; 227 e_=obj.e_s;
214 mat=obj.matrices.B; 228 mat=obj.matrices.B;
215 [V,Vi,D,pos]=obj.matrixDiag(mat,x,y(1)); 229 boundPos='l';
216 Hi=obj.Hyi; 230 Hi=obj.Hxi;
217 231 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(1));
218 tau=-1;
219 closure=Hi*e_*tau*mat_plus*e_';
220 penalty=Hi*e_*tau*mat_plus;
221
222 case {'n','N','north'} 232 case {'n','N','north'}
223 e_=obj.e_n; 233 e_=obj.e_n;
224 mat=obj.matrices.B; 234 mat=obj.matrices.B;
225 [V,Vi,D,pos]=obj.matrixDiag(mat,x,y(end)); 235 boundPos='r';
226 Hi=obj.Hyi; 236 Hi=obj.Hxi;
227 tau=1; 237 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(end));
228 238 end
229 closure=Hi*e_*tau*mat_minus*e_'; 239
230 penalty=Hi*e_*tau*mat_minus; 240 pos=signVec(1); zeroval=signVec(2); neg=signVec(3);
231 241
232 end 242
233 end 243
234 244 switch boundPos
235 245 case {'l'}
236 function [closure, penalty]=GetBoundarydata_char_fel(obj,boundary) 246 tau=sparse(4*side,pos*side);
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,:); 247 Vi_plus=Vi(1:pos*side,:);
248 Vi_minus=Vi(pos*side+1:4*side,:);
249
250 V_plus=Vi(:,1:pos*side);
251 V_minus=Vi(:,(pos+zeroval)*side+1:4*side);
252
253 tau(1:pos*side,:)=-abs(D(1:pos*side,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; 254 R=-inv(L*V_plus)*(L*V_minus);
255 255 closure=Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
256 closure=Hi*e_*V*tau*Vi_plus*e_'; 256 penalty=-Hi*e_*V*tau*inv(L*V_plus)*L;
257 penalty=Hi*e_*V*tau*Vi_plus; 257
258 258
259 case {'e','E','east'} 259 case {'r'}
260 e_=obj.e_e; 260 tau=sparse(4*side,neg*side);
261 mat=obj.matrices.A; 261 tau((pos+zeroval)*side+1:4*side,:)=-abs(D((pos+zeroval)*side+1:4*side,(pos+zeroval)*side+1:4*side));
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,:); 262 Vi_plus=Vi(1:pos*side,:);
280 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side)); 263 Vi_minus=Vi((pos+zeroval)*side+1:4*side,:);
281 mat_plus=V*(D+abs(D))*Vi/2; 264
282 265 V_plus=Vi(:,1:pos*side);
283 closure=Hi*e_*V*tau*Vi_plus*e_'; 266 V_minus=Vi(:,(pos+zeroval)*side+1:4*side);
284 penalty=Hi*e_*V*tau*Vi_plus; 267 R=-inv(L*V_minus)*(L*V_plus);
285 case {'n','N','north'} 268 closure=Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
286 e_=obj.e_n; 269 penalty=-Hi*e_*V*tau*inv(L*V_minus)*L;
287 mat=obj.matrices.B; 270
288 [V,Vi,D,pos]=obj.matrixDiag(mat,x,y(end)); 271
289 Hi=obj.Hyi; 272 end
290 tau=sparse(4*side,(4-pos)*side); 273 end
291 tau(pos*side+1:4*side,:)=-abs(D(pos*side+1:4*side,pos*side+1:4*side)); 274
292 Vi_minus=Vi(pos*side+1:4*side,:); 275
293 mat_minus=V*(D-abs(D))*Vi/2; 276 function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y)
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; 277 params=obj.params;
302 syms xs ys; 278 syms xs ys;
303 [V, D]=eig(mat(params,xs,ys)); 279 [V, D]=eig(mat(params,xs,ys));
304 xs=1;ys=1; 280 xs=1;ys=1;
305 DD=eval(diag(D)); 281 DD=eval(diag(D));
306 282
307 pos=find(DD>=0); %Now zero eigenvalues are calculated as possitive, Maybe it should not???? 283 poseig=find(DD>0);
308 neg=find(DD<0); 284 zeroeig=find(DD==0);
285 negeig=find(DD<0);
309 syms xs ys 286 syms xs ys
310 DD=diag(D); 287 DD=diag(D);
311 288
312 D=diag([DD(pos); DD(neg)]); 289 D=diag([DD(poseig);DD(zeroeig); DD(negeig)]);
313 V=[V(:,pos) V(:,neg)]; 290 V=[V(:,poseig) V(:,zeroeig) V(:,negeig)];
314 291
315 xs=x; ys=y; 292 xs=x; ys=y;
316 293
317 294
318 side=max(length(x),length(y)); 295 side=max(length(x),length(y));
325 end 302 end
326 end 303 end
327 304
328 D=sparse(Dret); 305 D=sparse(Dret);
329 V=sparse(normc(Vret)); 306 V=sparse(normc(Vret));
330 V=obj.matrixBuild(V); 307 V=obj.matrixBuild(V,x,y);
331 D=obj.matrixBuild(D); 308 D=obj.matrixBuild(D,x,y);
332 Vi=inv(V); 309 Vi=inv(V);
333 pos=length(pos); 310 signVec=[length(poseig),length(zeroeig),length(negeig)];
334 end 311 end
335 312
336 end 313 end
337 314
338 methods(Static) 315 methods(Static)