Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst2dCurve.m @ 352:9b3d7fc61a36 feature/hypsyst

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Changed operator in Utux
author Ylva Rydin <ylva.rydin@telia.com>
date Thu, 10 Nov 2016 20:47:40 +0100
parents 5d5652fe826a
children 9d1fc984f40d
comparison
equal deleted inserted replaced
351:7cc3d5bd3692 352:9b3d7fc61a36
145 matVec=mat; 145 matVec=mat;
146 [rows,cols]=size(matVec); 146 [rows,cols]=size(matVec);
147 side=max(length(X),length(Y)); 147 side=max(length(X),length(Y));
148 cols=cols/side; 148 cols=cols/side;
149 end 149 end
150 ret=kron(ones(rows,cols),speye(side)); 150 ret=cell(rows,cols);
151 151
152 for ii=1:rows 152 for ii=1:rows
153 for jj=1:cols 153 for jj=1:cols
154 ret((ii-1)*side+1:ii*side,(jj-1)*side+1:jj*side)=diag(matVec(ii,(jj-1)*side+1:jj*side)); 154 ret{ii,jj}=diag(matVec(ii,(jj-1)*side+1:jj*side));
155 end 155 end
156 end 156 end
157 ret=cell2mat(ret);
157 end 158 end
158 159
159 160
160 function [closure, penalty]=boundary_condition_char(obj,boundary) 161 function [closure, penalty]=boundary_condition_char(obj,boundary)
161 params=obj.params; 162 params=obj.params;
257 eta_x=Ji*obj.Y_xi(obj.index_s); 258 eta_x=Ji*obj.Y_xi(obj.index_s);
258 eta_y=-Ji*obj.X_xi(obj.index_s); 259 eta_y=-Ji*obj.X_xi(obj.index_s);
259 L=obj.evaluateCoefficientMatrix(L,eta_x,eta_y,[],[]); 260 L=obj.evaluateCoefficientMatrix(L,eta_x,eta_y,[],[]);
260 side=max(length(xi)); 261 side=max(length(xi));
261 case {'n','N','north'} 262 case {'n','N','north'}
262 e_=obj.e_n; for ii=1:rows 263 e_=obj.e_n;
263 264
264 mat=obj.Bhat; 265 mat=obj.Bhat;
265 boundPos='r'; 266 boundPos='r';
266 Hi=obj.Hetai; 267 Hi=obj.Hetai;
267 [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)); 268 [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));
302 closure=Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_'; 303 closure=Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
303 penalty=-Hi*e_*V*tau*inv(L*V_minus)*L; 304 penalty=-Hi*e_*V*tau*inv(L*V_minus)*L;
304 end 305 end
305 end 306 end
306 307
307
308 function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y,x_,y_) 308 function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y,x_,y_)
309 params=obj.params; 309 params=obj.params;
310 syms xs ys 310 syms xs ys
311 if(sum(abs(x_))~=0) 311 if(sum(abs(x_))~=0)
312 syms xs_ 312 syms xs_
319 else 319 else
320 ys_=0; 320 ys_=0;
321 end 321 end
322 322
323 [V, D]=eig(mat(params,xs,ys,xs_,ys_)); 323 [V, D]=eig(mat(params,xs,ys,xs_,ys_));
324 Vi=inv(V);
324 syms xs ys xs_ ys_ 325 syms xs ys xs_ ys_
325 326
326 xs=x; 327 xs=x;
327 ys=y; 328 ys=y;
328 xs_=x_; 329 xs_=x_;
329 ys_=y_; 330 ys_=y_;
330 331
331 side=max(length(x),length(y)); 332 side=max(length(x),length(y));
332 Dret=zeros(obj.n,side*obj.n); 333 Dret=zeros(obj.n,side*obj.n);
333 Vret=zeros(obj.n,side*obj.n); 334 Vret=zeros(obj.n,side*obj.n);
335 Viret=zeros(obj.n,side*obj.n);
334 for ii=1:obj.n 336 for ii=1:obj.n
335 for jj=1:obj.n 337 for jj=1:obj.n
336 Dret(jj,(ii-1)*side+1:side*ii)=eval(D(jj,ii)); 338 Dret(jj,(ii-1)*side+1:side*ii)=eval(D(jj,ii));
337 Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii)); 339 Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii));
340 Viret(jj,(ii-1)*side+1:side*ii)=eval(Vi(jj,ii));
338 end 341 end
339 end 342 end
340 343
341 D=sparse(Dret); 344 D=sparse(Dret);
342 V=sparse(Vret); 345 V=sparse(Vret);
346 Vi=sparse(Viret);
343 V=obj.evaluateCoefficientMatrix(V,x,y,x_,y_); 347 V=obj.evaluateCoefficientMatrix(V,x,y,x_,y_);
344 D=obj.evaluateCoefficientMatrix(D,x,y,x_,y_); 348 D=obj.evaluateCoefficientMatrix(D,x,y,x_,y_);
349 Vi=obj.evaluateCoefficientMatrix(Vi,x,y,x_,y_);
345 DD=diag(D); 350 DD=diag(D);
346 351
347 poseig=(DD>0); 352 poseig=(DD>0);
348 zeroeig=(DD==0); 353 zeroeig=(DD==0);
349 negeig=(DD<0); 354 negeig=(DD<0);
350 355
351 D=diag([DD(poseig); DD(zeroeig); DD(negeig)]); 356 D=diag([DD(poseig); DD(zeroeig); DD(negeig)]);
352 V=[V(:,poseig) V(:,zeroeig) V(:,negeig)]; 357 V=[V(:,poseig) V(:,zeroeig) V(:,negeig)];
353 Vi=inv(V); 358 Vi=[Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
354 signVec=[sum(poseig),sum(zeroeig),sum(negeig)]; 359 signVec=[sum(poseig),sum(zeroeig),sum(negeig)];
355 end 360 end
356
357 end 361 end
358 end 362 end