Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst2dCurve.m @ 300:32f3ee81bc37 feature/hypsyst

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
A commit before i destroy everything
author Ylva Rydin <ylva.rydin@telia.com>
date Tue, 04 Oct 2016 16:25:38 +0200
parents 4d8d6eb0c116
children d9860ebc3148
comparison
equal deleted inserted replaced
299:4d8d6eb0c116 300:32f3ee81bc37
85 [X,Y] = ti.map(obj.xi,obj.eta); 85 [X,Y] = ti.map(obj.xi,obj.eta);
86 86
87 [x_xi,x_eta] = gridDerivatives(X,ops_xi.D1,ops_eta.D1); 87 [x_xi,x_eta] = gridDerivatives(X,ops_xi.D1,ops_eta.D1);
88 [y_xi,y_eta] = gridDerivatives(Y,ops_xi.D1, ops_eta.D1); 88 [y_xi,y_eta] = gridDerivatives(Y,ops_xi.D1, ops_eta.D1);
89 89
90 obj.X=reshape(X,m_xi*m_eta,1); 90 obj.X=reshape(X,m_tot,1);
91 obj.Y=reshape(Y,m_xi*m_eta,1); 91 obj.Y=reshape(Y,m_tot,1);
92 obj.X_xi=reshape(x_xi,m_xi*m_eta,1); 92 obj.X_xi=reshape(x_xi,m_tot,1);
93 obj.Y_xi=reshape(y_xi,m_xi*m_eta,1); 93 obj.Y_xi=reshape(y_xi,m_tot,1);
94 obj.X_eta=reshape(x_eta,m_xi*m_eta,1); 94 obj.X_eta=reshape(x_eta,m_tot,1);
95 obj.Y_eta=reshape(y_eta,m_xi*m_eta,1); 95 obj.Y_eta=reshape(y_eta,m_tot,1);
96 96
97 Ahat_evaluated = obj.evaluateCoefficientMatrix(obj.Ahat, obj.X, obj.Y,obj.X_eta,obj.Y_eta); 97 Ahat_evaluated = obj.evaluateCoefficientMatrix(obj.Ahat, obj.X, obj.Y,obj.X_eta,obj.Y_eta);
98 Bhat_evaluated = obj.evaluateCoefficientMatrix(obj.Bhat, obj.X, obj.Y,obj.X_xi,obj.Y_xi); 98 Bhat_evaluated = obj.evaluateCoefficientMatrix(obj.Bhat, obj.X, obj.Y,obj.X_xi,obj.Y_xi);
99 E_evaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,[],[]); 99 E_evaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,[],[]);
100 100
101 obj.m=m; 101 obj.m=m;
102 obj.h=[ops_xi.h ops_eta.h]; 102 obj.h=[ops_xi.h ops_eta.h];
103 obj.order=order; 103 obj.order=order;
104 obj.J=x_xi.*y_eta - x_eta.*y_xi; 104 obj.J=obj.X_xi.*obj.Y_eta - obj.X_eta.*obj.Y_xi;
105 obj.Ji =kr(I_n,spdiags(1./obj.J(:),0,m_xi*m_eta,m_xi*m_eta)); 105 obj.Ji =kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot));
106 106
107 obj.D=obj.Ji*(-Ahat_evaluated*D1_xi-Bhat_evaluated*D1_eta)-E_evaluated; 107 obj.D=obj.Ji*(-Ahat_evaluated*D1_xi-Bhat_evaluated*D1_eta)-E_evaluated;
108 108
109 end 109 end
110 110
159 end 159 end
160 160
161 161
162 function [closure, penalty]=boundary_condition_char(obj,boundary) 162 function [closure, penalty]=boundary_condition_char(obj,boundary)
163 params=obj.params; 163 params=obj.params;
164 X=obj.X; Y=obj.Y;
164 xi=obj.xi; eta=obj.eta; 165 xi=obj.xi; eta=obj.eta;
165 side=max(length(xi),length(eta)); 166 side=max(length(xi),length(eta));
166 167
167 switch boundary 168 switch boundary
168 case {'w','W','west'} 169 case {'w','W','west'}
169 e_=obj.e_w; 170 e_=obj.e_w;
170 mat=obj.Ahat; 171 mat=obj.Ahat;
171 boundPos='l'; 172 boundPos='l';
172 Hi=obj.Hxii; 173 Hi=obj.Hxii;
173 [V,Vi,D,signVec]=obj.matrixDiag(mat,xi(1),eta,obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w)); 174 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_w),Y(obj.index_w),obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w));
174 case {'e','E','east'} 175 case {'e','E','east'}
175 e_=obj.e_e; 176 e_=obj.e_e;
176 mat=obj.Ahat; 177 mat=obj.Ahat;
177 boundPos='r'; 178 boundPos='r';
178 Hi=obj.Hxii; 179 Hi=obj.Hxii;
179 [V,Vi,D,signVec]=obj.matrixDiag(mat,xi(end),eta,obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e)); 180 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e));
180 case {'s','S','south'} 181 case {'s','S','south'}
181 e_=obj.e_s; 182 e_=obj.e_s;
182 mat=obj.Bhat; 183 mat=obj.Bhat;
183 boundPos='l'; 184 boundPos='l';
184 Hi=obj.Hetai; 185 Hi=obj.Hetai;
185 [V,Vi,D,signVec]=obj.matrixDiag(mat,xi,eta(1),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s)); 186 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_s),Y(obj.index_s),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s));
186 case {'n','N','north'} 187 case {'n','N','north'}
187 e_=obj.e_n; 188 e_=obj.e_n;
188 mat=obj.Bhat; 189 mat=obj.Bhat;
189 boundPos='r'; 190 boundPos='r';
190 Hi=obj.Hetai; 191 Hi=obj.Hetai;
191 [V,Vi,D,signVec]=obj.matrixDiag(mat,xi,eta(end),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n)); 192 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n));
192 end 193 end
193 194
194 pos=signVec(1); zeroval=signVec(2); neg=signVec(3); 195 pos=signVec(1); zeroval=signVec(2); neg=signVec(3);
195 196
196 switch boundPos 197 switch boundPos
197 case {'l'} 198 case {'l'}
198 tau=sparse(obj.n*side,pos*side); 199 tau=sparse(obj.n*side,pos);
199 Vi_plus=Vi(1:pos*side,:); 200 Vi_plus=Vi(1:pos,:);
200 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side)); 201 tau(1:pos,:)=-abs(D(1:pos,1:pos));
201 closure=Hi*e_*V*tau*Vi_plus*e_'; 202 closure=Hi*e_*V*tau*Vi_plus*e_';
202 penalty=-Hi*e_*V*tau*Vi_plus; 203 penalty=-Hi*e_*V*tau*Vi_plus;
203 case {'r'} 204 case {'r'}
204 tau=sparse(obj.n*side,neg*side); 205 tau=sparse(obj.n*side,neg);
205 tau((pos+zeroval)*side+1:obj.n*side,:)=-abs(D((pos+zeroval)*side+1:obj.n*side,(pos+zeroval)*side+1:obj.n*side)); 206 tau((pos+zeroval)+1:obj.n*side,:)=-abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
206 Vi_minus=Vi((pos+zeroval)*side+1:obj.n*side,:); 207 Vi_minus=Vi((pos+zeroval)+1:obj.n*side,:);
207 closure=Hi*e_*V*tau*Vi_minus*e_'; 208 closure=Hi*e_*V*tau*Vi_minus*e_';
208 penalty=-Hi*e_*V*tau*Vi_minus; 209 penalty=-Hi*e_*V*tau*Vi_minus;
209 end 210 end
210 end 211 end
211 212
212 213
213 function [closure,penalty]=boundary_condition_general(obj,boundary,L) 214 function [closure,penalty]=boundary_condition_general(obj,boundary,L)
214 params=obj.params; 215 params=obj.params;
215 xi=obj.xi; eta=obj.eta; 216 X=obj.X; Y=obj.Y;
216 side=max(length(xi),length(eta)); 217 side=max(length(xi),length(eta));
217 218
218 switch boundary 219 switch boundary
219 case {'w','W','west'} 220 case {'w','W','west'}
220 e_=obj.e_w; 221 e_=obj.e_w;
221 mat=obj.Ahat; 222 mat=obj.Ahat;
222 boundPos='l'; 223 boundPos='l';
223 Hi=obj.Hxii; 224 Hi=obj.Hxii;
224 [V,Vi,D,signVec]=obj.matrixDiag(mat,xi(1),eta,obj.x_eta,obj.y_eta); 225 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_w),Y(obj.index_w),obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w));
225 L=obj.evaluateCoefficientMatrix(L,xi(1),eta); 226 L=obj.evaluateCoefficientMatrix(L,X(obj.index_w),Y(obj.index_w));
226 case {'e','E','east'} 227 case {'e','E','east'}
227 e_=obj.e_e; 228 e_=obj.e_e;
228 mat=obj.Ahat; 229 mat=obj.Ahat;
229 boundPos='r'; 230 boundPos='r';
230 Hi=obj.Hxii; 231 Hi=obj.Hxii;
231 [V,Vi,D,signVec]=obj.matrixDiag(mat,xi(end),eta,obj.x_eta,obj.y_eta); 232 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e));
232 L=obj.evaluateCoefficientMatrix(L,xi(end),eta,[],[]); 233 L=obj.evaluateCoefficientMatrix(L,X(obj.index_e),Y(obj.index_e),[],[]);
233 case {'s','S','south'} 234 case {'s','S','south'}
234 e_=obj.e_s; 235 e_=obj.e_s;
235 mat=obj.Bhat; 236 mat=obj.Bhat;
236 boundPos='l'; 237 boundPos='l';
237 Hi=obj.Hetai; 238 Hi=obj.Hetai;
238 [V,Vi,D,signVec]=obj.matrixDiag(mat,xi,eta(1),obj.x_xi,obj.y_xi); 239 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_s),Y(obj.index_s),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s));
239 L=obj.evaluateCoefficientMatrix(L,xi,eta(1)); 240 L=obj.evaluateCoefficientMatrix(L,X(obj.index_s),Y(obj.index_s),[],[]);
240 case {'n','N','north'} 241 case {'n','N','north'}
241 e_=obj.e_n; 242 e_=obj.e_n;
242 mat=obj.Bhat; 243 mat=obj.Bhat;
243 boundPos='r'; 244 boundPos='r';
244 Hi=obj.Hetai; 245 Hi=obj.Hetai;
245 [V,Vi,D,signVec]=obj.matrixDiag(mat,xi,eta(end),obj.x_xi,obj.y_xi); 246 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n));
246 L=obj.evaluateCoefficientMatrix(L,xi,eta(end)); 247 L=obj.evaluateCoefficientMatrix(L,X(obj.index_n),Y(obj.index_n));
247 end 248 end
248 249
249 pos=signVec(1); zeroval=signVec(2); neg=signVec(3); 250 pos=signVec(1); zeroval=signVec(2); neg=signVec(3);
250 251
251 switch boundPos 252 switch boundPos
252 case {'l'} 253 case {'l'}
253 tau=sparse(obj.n*side,pos*side); 254 tau=sparse(obj.n*side,pos);
254 Vi_plus=Vi(1:pos*side,:); 255 Vi_plus=Vi(1:pos,:);
255 Vi_minus=Vi(pos*side+1:obj.n*side,:); 256 Vi_minus=Vi(pos+1:obj.n*side,:);
256 V_plus=V(:,1:pos*side); 257 V_plus=V(:,1:pos);
257 V_minus=V(:,(pos+zeroval)*side+1:obj.n*side); 258 V_minus=V(:,(pos+zeroval)+1:obj.n*side);
258 259
259 tau(1:pos*side,:)=-abs(D(1:pos*side,1:pos*side)); 260 tau(1:pos*side,:)=-abs(D(1:pos,1:pos));
260 R=-inv(L*V_plus)*(L*V_minus); 261 R=-inv(L*V_plus)*(L*V_minus);
261 closure=Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; 262 closure=Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
262 penalty=-Hi*e_*V*tau*inv(L*V_plus)*L; 263 penalty=-Hi*e_*V*tau*inv(L*V_plus)*L;
263 case {'r'} 264 case {'r'}
264 tau=sparse(obj.n*side,neg*side); 265 tau=sparse(obj.n*side,neg);
265 tau((pos+zeroval)*side+1:obj.n*side,:)=-abs(D((pos+zeroval)*side+1:obj.n*side,(pos+zeroval)*side+1:obj.n*side)); 266 tau((pos+zeroval)+1:obj.n*side,:)=-abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
266 Vi_plus=Vi(1:pos*side,:); 267 Vi_plus=Vi(1:pos,:);
267 Vi_minus=Vi((pos+zeroval)*side+1:obj.n*side,:); 268 Vi_minus=Vi((pos+zeroval)+1:obj.n*side,:);
268 269
269 V_plus=V(:,1:pos*side); 270 V_plus=V(:,1:pos);
270 V_minus=V(:,(pos+zeroval)*side+1:obj.n*side); 271 V_minus=V(:,(pos+zeroval)+1:obj.n*side);
271 R=-inv(L*V_minus)*(L*V_plus); 272 R=-inv(L*V_minus)*(L*V_plus);
272 closure=Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_'; 273 closure=Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
273 penalty=-Hi*e_*V*tau*inv(L*V_minus)*L; 274 penalty=-Hi*e_*V*tau*inv(L*V_minus)*L;
274 end 275 end
275 end 276 end
289 else 290 else
290 ys_=0; 291 ys_=0;
291 end 292 end
292 293
293 [V, D]=eig(mat(params,xs,ys,xs_,ys_)); 294 [V, D]=eig(mat(params,xs,ys,xs_,ys_));
294 xs=1;ys=1; xs_=x_(1); ys_=y_(1);
295 DD=eval(diag(D));
296
297 poseig=find(DD>0);
298 zeroeig=find(DD==0);
299 negeig=find(DD<0);
300 syms xs ys xs_ ys_ 295 syms xs ys xs_ ys_
301 DD=diag(D); 296
302
303 D=diag([DD(poseig);DD(zeroeig); DD(negeig)]);
304 V=[V(:,poseig) V(:,zeroeig) V(:,negeig)];
305 Vi=inv(V);
306 xs=x; 297 xs=x;
307 ys=y; 298 ys=y;
308 xs_=x_; 299 xs_=x_;
309 ys_=y_; 300 ys_=y_;
310 301
311 side=max(length(x),length(y)); 302 side=max(length(x),length(y));
312 Dret=zeros(obj.n,side*obj.n); 303 Dret=zeros(obj.n,side*obj.n);
313 Vret=zeros(obj.n,side*obj.n); 304 Vret=zeros(obj.n,side*obj.n);
314 Viret=zeros(obj.n,side*obj.n);
315 for ii=1:obj.n 305 for ii=1:obj.n
316 for jj=1:obj.n 306 for jj=1:obj.n
317 Dret(jj,(ii-1)*side+1:side*ii)=eval(D(jj,ii)); 307 Dret(jj,(ii-1)*side+1:side*ii)=eval(D(jj,ii));
318 Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii)); 308 Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii));
319 Viret(jj,(ii-1)*side+1:side*ii)=eval(Vi(jj,ii));
320 end 309 end
321 end 310 end
322 311
323 D=sparse(Dret); 312 D=sparse(Dret);
324 V=sparse(Vret); 313 V=sparse(Vret);
325 Vi=sparse(Viret);
326 V=obj.evaluateCoefficientMatrix(V,x,y,x_,y_); 314 V=obj.evaluateCoefficientMatrix(V,x,y,x_,y_);
327 D=obj.evaluateCoefficientMatrix(D,x,y,x_,y_); 315 D=obj.evaluateCoefficientMatrix(D,x,y,x_,y_);
328 Vi=obj.evaluateCoefficientMatrix(Vi,x,y,x_,y_); 316 DD=diag(D);
329 signVec=[length(poseig),length(zeroeig),length(negeig)]; 317
318 poseig=(DD>0);
319 zeroeig=(DD==0);
320 negeig=(DD<0);
321
322 D=diag([DD(poseig); DD(zeroeig); DD(negeig)]);
323 V=[V(:,poseig) V(:,zeroeig) V(:,negeig)];
324 Vi=inv(V);
325 signVec=[sum(poseig),sum(zeroeig),sum(negeig)];
330 end 326 end
331 327
332 end 328 end
333 end 329 end