Mercurial > repos > public > sbplib

comparison +scheme/hypsyst2d.m @ 294:8ff6ec6249e8 feature/hypsyst

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"General" boundary conditions implemented
author Ylva Rydin <ylva.rydin@telia.com>
date Mon, 26 Sep 2016 09:54:43 +0200
parents 2d604d16842c
children da0131655035
comparison
equal deleted inserted replaced
293:2d604d16842c 294:8ff6ec6249e8
48 obj.A=obj.matrixBuild(matrices.A); 48 obj.A=obj.matrixBuild(matrices.A);
49 obj.B=obj.matrixBuild(matrices.B); 49 obj.B=obj.matrixBuild(matrices.B);
50 obj.E=obj.matrixBuild(matrices.E); 50 obj.E=obj.matrixBuild(matrices.E);
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);
55 I_x = speye(m_x); obj.I_x=I_x; 55 I_x = speye(m_x); obj.I_x=I_x;
56 I_y = speye(m_y); obj.I_y=I_y; 56 I_y = speye(m_y); obj.I_y=I_y;
57 57
58 58
59 D1_x = kr(kr(I_n,ops_x.D1),I_y); 59 D1_x = kr(kr(I_n,ops_x.D1),I_y);
60 obj.Hxi= kr(kr(I_n,ops_x.HI),I_y); 60 obj.Hxi= kr(kr(I_n,ops_x.HI),I_y);
61 D1_y=kr(I_n,kr(I_x,ops_y.D1)); 61 D1_y=kr(I_n,kr(I_x,ops_y.D1));
62 obj.Hyi=kr(I_n,kr(I_x,ops_y.HI)); 62 obj.Hyi=kr(I_n,kr(I_x,ops_y.HI));
63 63
67 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));
68 68
69 obj.m=m; 69 obj.m=m;
70 obj.h=[ops_x.h ops_y.h]; 70 obj.h=[ops_x.h ops_y.h];
71 obj.order=order; 71 obj.order=order;
72 72
73 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;
74 74
75 end 75 end
76
76 % 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
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. 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.
78 % 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'.
79 % 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.
80 % 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.
81 % 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.
82 % neighbour_boundary is a string specifying which boundary to interface to. 83 % neighbour_boundary is a string specifying which boundary to interface to.
83 function [closure, penalty] = boundary_condition(obj,boundary,type,L) 84 function [closure, penalty] = boundary_condition(obj,boundary,type,L)
84 default_arg('type','neumann'); 85 default_arg('type','char');
85 default_arg('data',0);
86
87 switch type 86 switch type
88 case{'c','char'} 87 case{'c','char'}
89 [closure,penalty]=GetBoundarydata_char(obj,boundary); 88 [closure,penalty]=GetBoundarydata_char(obj,boundary);
90 case{'general'} 89 case{'general'}
91 [closure,penalty]=GeneralBoundaryCond(obj,boundary,L); 90 [closure,penalty]=GeneralBoundaryCond(obj,boundary,L);
100 99
101 function N = size(obj) 100 function N = size(obj)
102 N = obj.m; 101 N = obj.m;
103 end 102 end
104 103
105 function [ret]=matrixBuild(obj,mat,x,y) 104 function [ret]=matrixBuild(obj,mat,X,Y)
106 %extra info for coordinate transfomration mult my y_ny and 105 params=obj.params;
107 %x,ny osv... 106 default_arg('X',obj.X);
108 params=obj.params; 107 default_arg('Y',obj.Y)
109 X=obj.X;
110 Y=obj.Y;
111 108
112 if isa(mat,'function_handle') 109 if isa(mat,'function_handle')
113 [rows,cols]=size(mat(params,0,0)); 110 [rows,cols]=size(mat(params,0,0));
114 matVec=mat(params,X',Y'); 111 matVec=mat(params,X',Y');
115 matVec=sparse(matVec); 112 matVec=sparse(matVec);
116 side=max(length(X),length(Y)); 113 side=max(length(X),length(Y));
117 else 114 else
118 matVec=mat; 115 matVec=mat;
119 [rows,cols]=size(matVec); 116 [rows,cols]=size(matVec);
120 side=max(length(x),length(y)); 117 side=max(length(X),length(Y));
121 cols=cols/side; 118 cols=cols/side;
122 end 119 end
123 120 ret=kron(ones(rows,cols),speye(side));
124 121
125 ret=kron(ones(rows,cols),speye(side)); 122 for ii=1:rows
126 123 for jj=1:cols
127 for ii=1:rows 124 ret((ii-1)*side+1:ii*side,(jj-1)*side+1:jj*side)=diag(matVec(ii,(jj-1)*side+1:jj*side));
128 for jj=1:cols 125 end
129 ret((ii-1)*side+1:ii*side,(jj-1)*side+1:jj*side)=diag(matVec(ii,(jj-1)*side+1:jj*side)); 126 end
130 end 127 end
131 end 128
132 129
133 end
134
135
136 function [closure, penalty]=GetBoundarydata_char(obj,boundary) 130 function [closure, penalty]=GetBoundarydata_char(obj,boundary)
137 params=obj.params; 131 params=obj.params;
138 x=obj.x; 132 x=obj.x; y=obj.y;
139 y=obj.y;
140
141 side=max(length(x),length(y)); 133 side=max(length(x),length(y));
142 134
143 switch boundary 135 switch boundary
144 case {'w','W','west'} 136 case {'w','W','west'}
145 e_=obj.e_w; 137 e_=obj.e_w;
146 mat=obj.matrices.A; 138 mat=obj.matrices.A;
147 boundPos='l'; 139 boundPos='l';
148 Hi=obj.Hxi; 140 Hi=obj.Hxi;
149 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(1),y); 141 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(1),y);
150 case {'e','E','east'} 142 case {'e','E','east'}
151 e_=obj.e_e; 143 e_=obj.e_e;
152 mat=obj.matrices.A; 144 mat=obj.matrices.A;
153 boundPos='r'; 145 boundPos='r';
154 Hi=obj.Hxi; 146 Hi=obj.Hxi;
155 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(end),y); 147 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(end),y);
156 case {'s','S','south'} 148 case {'s','S','south'}
157 e_=obj.e_s; 149 e_=obj.e_s;
158 mat=obj.matrices.B; 150 mat=obj.matrices.B;
159 boundPos='l'; 151 boundPos='l';
160 Hi=obj.Hxi; 152 Hi=obj.Hxi;
161 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(1)); 153 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(1));
162 case {'n','N','north'} 154 case {'n','N','north'}
163 e_=obj.e_n; 155 e_=obj.e_n;
164 mat=obj.matrices.B; 156 mat=obj.matrices.B;
165 boundPos='r'; 157 boundPos='r';
166 Hi=obj.Hxi; 158 Hi=obj.Hxi;
167 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(end)); 159 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(end));
168 end 160 end
169 161
170 pos=signVec(1); zeroval=signVec(2); neg=signVec(3); 162 pos=signVec(1); zeroval=signVec(2); neg=signVec(3);
171 163
172 switch boundPos 164 switch boundPos
173 case {'l'} 165 case {'l'}
174 tau=sparse(obj.n*side,pos*side); 166 tau=sparse(obj.n*side,pos*side);
175 Vi_plus=Vi(1:pos*side,:); 167 Vi_plus=Vi(1:pos*side,:);
176 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side)); 168 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side));
177 closure=Hi*e_*V*tau*Vi_plus*e_'; 169 closure=Hi*e_*V*tau*Vi_plus*e_';
178 penalty=-Hi*e_*V*tau*Vi_plus; 170 penalty=-Hi*e_*V*tau*Vi_plus;
179 171 case {'r'}
180 case {'r'}
181 tau=sparse(obj.n*side,neg*side); 172 tau=sparse(obj.n*side,neg*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)); 173 tau((pos+zeroval)*side+1:obj.n*side,:)=-abs(D((pos+zeroval)*side+1:obj.n*side,(pos+zeroval)*side+1:obj.n*side));
183 Vi_minus=Vi((pos+zeroval)*side+1:obj.n*side,:); 174 Vi_minus=Vi((pos+zeroval)*side+1:obj.n*side,:);
184 closure=Hi*e_*V*tau*Vi_minus*e_'; 175 closure=Hi*e_*V*tau*Vi_minus*e_';
185 penalty=-Hi*e_*V*tau*Vi_minus; 176 penalty=-Hi*e_*V*tau*Vi_minus;
186 177 end
187 end 178 end
188 end 179
189
190 180
191 function [closure,penalty]=GeneralBoundaryCond(obj,boundary,L) 181 function [closure,penalty]=GeneralBoundaryCond(obj,boundary,L)
192 params=obj.params; 182 params=obj.params;
193 x=obj.x; 183 x=obj.x; y=obj.y;
194 y=obj.y;
195 L=obj.matrixBuild(L,x,y);
196 side=max(length(x),length(y)); 184 side=max(length(x),length(y));
197 185
198
199 switch boundary 186 switch boundary
200 case {'w','W','west'} 187 case {'w','W','west'}
201 e_=obj.e_w; 188 e_=obj.e_w;
202 mat=obj.matrices.A; 189 mat=obj.matrices.A;
203 boundPos='l'; 190 boundPos='l';
204 Hi=obj.Hxi; 191 Hi=obj.Hxi;
205 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(1),y); 192 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(1),y);
193 L=obj.matrixBuild(L,x(1),y);
206 case {'e','E','east'} 194 case {'e','E','east'}
207 e_=obj.e_e; 195 e_=obj.e_e;
208 mat=obj.matrices.A; 196 mat=obj.matrices.A;
209 boundPos='r'; 197 boundPos='r';
210 Hi=obj.Hxi; 198 Hi=obj.Hxi;
211 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(end),y); 199 [V,Vi,D,signVec]=obj.matrixDiag(mat,x(end),y);
200 L=obj.matrixBuild(L,x(end),y);
212 case {'s','S','south'} 201 case {'s','S','south'}
213 e_=obj.e_s; 202 e_=obj.e_s;
214 mat=obj.matrices.B; 203 mat=obj.matrices.B;
215 boundPos='l'; 204 boundPos='l';
216 Hi=obj.Hxi; 205 Hi=obj.Hxi;
217 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(1)); 206 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(1));
207 L=obj.matrixBuild(L,x,y(1));
218 case {'n','N','north'} 208 case {'n','N','north'}
219 e_=obj.e_n; 209 e_=obj.e_n;
220 mat=obj.matrices.B; 210 mat=obj.matrices.B;
221 boundPos='r'; 211 boundPos='r';
222 Hi=obj.Hxi; 212 Hi=obj.Hxi;
223 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(end)); 213 [V,Vi,D,signVec]=obj.matrixDiag(mat,x,y(end));
224 end 214 L=obj.matrixBuild(L,x,y(end));
225 215 end
226 pos=signVec(1); zeroval=signVec(2); neg=signVec(3); 216
227 217 pos=signVec(1); zeroval=signVec(2); neg=signVec(3);
228
229 218
230 switch boundPos 219 switch boundPos
231 case {'l'} 220 case {'l'}
232 tau=sparse(obj.n*side,pos*side); 221 tau=sparse(obj.n*side,pos*side);
233 Vi_plus=Vi(1:pos*side,:); 222 Vi_plus=Vi(1:pos*side,:);
234 Vi_minus=Vi(pos*side+1:obj.n*side,:); 223 Vi_minus=Vi(pos*side+1:obj.n*side,:);
224 V_plus=V(:,1:pos*side);
225 V_minus=V(:,(pos+zeroval)*side+1:obj.n*side);
235 226
236 V_plus=Vi(:,1:pos*side); 227 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side));
237 V_minus=Vi(:,(pos+zeroval)*side+1:obj.n*side);
238
239 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side));
240 R=-inv(L*V_plus)*(L*V_minus); 228 R=-inv(L*V_plus)*(L*V_minus);
241 closure=Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; 229 closure=Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
242 penalty=-Hi*e_*V*tau*inv(L*V_plus)*L; 230 penalty=-Hi*e_*V*tau*inv(L*V_plus)*L;
243
244
245 case {'r'} 231 case {'r'}
246 tau=sparse(obj.n*side,neg*side); 232 tau=sparse(obj.n*side,neg*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)); 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));
248 Vi_plus=Vi(1:pos*side,:); 234 Vi_plus=Vi(1:pos*side,:);
249 Vi_minus=Vi((pos+zeroval)*side+1:obj.n*side,:); 235 Vi_minus=Vi((pos+zeroval)*side+1:obj.n*side,:);
250 236
251 V_plus=Vi(:,1:pos*side); 237 V_plus=V(:,1:pos*side);
252 V_minus=Vi(:,(pos+zeroval)*side+1:obj.n*side); 238 V_minus=V(:,(pos+zeroval)*side+1:obj.n*side);
253 R=-inv(L*V_minus)*(L*V_plus); 239 R=-inv(L*V_minus)*(L*V_plus);
254 closure=Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_'; 240 closure=Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
255 penalty=-Hi*e_*V*tau*inv(L*V_minus)*L; 241 penalty=-Hi*e_*V*tau*inv(L*V_minus)*L;
256 242 end
257
258 end
259 end 243 end
260 244
261 245
262 function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y) 246 function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y)
263 params=obj.params; 247 params=obj.params;
264 syms xs ys; 248 syms xs ys;
265 [V, D]=eig(mat(params,xs,ys)); 249 [V, D]=eig(mat(params,xs,ys));
266 xs=1;ys=1; 250 xs=1;ys=1;
267 DD=eval(diag(D)); 251 DD=eval(diag(D));
268 252
269 poseig=find(DD>0); 253 poseig=find(DD>0);
270 zeroeig=find(DD==0); 254 zeroeig=find(DD==0);
271 negeig=find(DD<0); 255 negeig=find(DD<0);
272 syms xs ys 256 syms xs ys
273 DD=diag(D); 257 DD=diag(D);
274 258
275 D=diag([DD(poseig);DD(zeroeig); DD(negeig)]); 259 D=diag([DD(poseig);DD(zeroeig); DD(negeig)]);
276 V=[V(:,poseig) V(:,zeroeig) V(:,negeig)]; 260 V=[V(:,poseig) V(:,zeroeig) V(:,negeig)];
277
278 xs=x; ys=y; 261 xs=x; ys=y;
279
280 262
281 side=max(length(x),length(y)); 263 side=max(length(x),length(y));
282 Dret=zeros(obj.n,side*obj.n); 264 Dret=zeros(obj.n,side*obj.n);
283 Vret=zeros(obj.n,side*obj.n); 265 Vret=zeros(obj.n,side*obj.n);
284 for ii=1:obj.n 266 for ii=1:obj.n
285 for jj=1:obj.n 267 for jj=1:obj.n
286 Dret(jj,(ii-1)*side+1:side*ii)=eval(D(jj,ii)); 268 Dret(jj,(ii-1)*side+1:side*ii)=eval(D(jj,ii));
287 Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii)); 269 Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii));
288 end 270 end
289 end 271 end
290 272
291 D=sparse(Dret); 273 D=sparse(Dret);
292 V=sparse(normc(Vret)); 274 V=sparse(normc(Vret));
293 V=obj.matrixBuild(V,x,y); 275 V=obj.matrixBuild(V,x,y);
294 D=obj.matrixBuild(D,x,y); 276 D=obj.matrixBuild(D,x,y);
295 Vi=inv(V); 277 Vi=inv(V);
305 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) 287 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
306 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); 288 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
307 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); 289 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
308 end 290 end
309 291
310 292
311 end 293 end
312 end 294 end