Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst3dCurve.m @ 351:7cc3d5bd3692 feature/hypsyst

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Curve linear code in 3D converges for manufactures solution
author Ylva Rydin <ylva.rydin@telia.com>
date Thu, 03 Nov 2016 22:03:09 +0100
parents 5d5652fe826a
children dbac99d2c318
comparison
equal deleted inserted replaced
350:5d5652fe826a 351:7cc3d5bd3692
208 matVec=mat; 208 matVec=mat;
209 [rows,cols]=size(matVec); 209 [rows,cols]=size(matVec);
210 side=max(length(X),length(Y)); 210 side=max(length(X),length(Y));
211 cols=cols/side; 211 cols=cols/side;
212 end 212 end
213 ret=kron(ones(rows,cols),speye(side)); 213 ret=cell(rows,cols);
214
214 215
215 for ii=1:rows 216 for ii=1:rows
216 for jj=1:cols 217 for jj=1:cols
217 ret((ii-1)*side+1:ii*side,(jj-1)*side+1:jj*side)=diag(matVec(ii,(jj-1)*side+1:jj*side)); 218 ret{ii,jj}=diag(matVec(ii,(jj-1)*side+1:jj*side));
218 end 219 end
219 end 220 end
221
222 ret=cell2mat(ret);
220 end 223 end
221 224
222 225
223 function [BM]=boundary_matrices(obj,boundary) 226 function [BM]=boundary_matrices(obj,boundary)
224 params=obj.params; 227 params=obj.params;
407 if(sum(abs(y_2))>eps) 410 if(sum(abs(y_2))>eps)
408 syms y_2s; 411 syms y_2s;
409 else 412 else
410 y_2s=0; 413 y_2s=0;
411 end 414 end
412 415
413
414 if(sum(abs(z_1))>eps) 416 if(sum(abs(z_1))>eps)
415 syms z_1s 417 syms z_1s
416 else 418 else
417 z_1s=0; 419 z_1s=0;
418 end 420 end
423 z_2s=0; 425 z_2s=0;
424 end 426 end
425 427
426 syms xs ys zs 428 syms xs ys zs
427 [V, D]=eig(mat(obj,xs,ys,zs,x_1s,x_2s,y_1s,y_2s,z_1s,z_2s)); 429 [V, D]=eig(mat(obj,xs,ys,zs,x_1s,x_2s,y_1s,y_2s,z_1s,z_2s));
428 Vi=inv(V); 430 Vi=inv(V);
429 431 % syms x_1s x_2s y_1s y_2s z_1s z_2s
430 syms x_1s x_2s y_1s y_2s z_1s z_2s 432 xs=x;
431 % V= matlabFunction(V); 433 ys=y;
432 % D= matlabFunction(D); 434 zs=z;
433 % Vi= matlabFunction(Vi); 435 x_1s=x_1;
434 % 436 x_2s=x_2;
435 % xs=x; 437 y_1s=y_1;
436 % ys=y; 438 y_2s=y_2;
437 % zs=z; 439 z_1s=z_1;
438 % x_1s=x_1; 440 z_2s=z_2;
439 % x_2s=x_2;
440 % y_1s=y_1;
441 % y_2s=y_2;
442 % z_1s=z_1;
443 % z_2s=z_2;
444 441
445 side=max(length(x),length(y)); 442 side=max(length(x),length(y));
446 Dret=zeros(obj.n,side*obj.n); 443 Dret=zeros(obj.n,side*obj.n);
447 Vret=zeros(obj.n,side*obj.n); 444 Vret=zeros(obj.n,side*obj.n);
448 Viret=zeros(obj.n,side*obj.n); 445 Viret=zeros(obj.n,side*obj.n);
449 446
450 for ii=1:obj.n 447 for ii=1:obj.n
451 for jj=1:obj.n 448 for jj=1:obj.n
452 Dpart=matlabFunction(D(jj,ii),'Vars',[xs ys zs x_1s x_2s y_1s y_2s z_1s z_2s]); 449 Dret(jj,(ii-1)*side+1:side*ii)=eval(D(jj,ii));
453 Vpart=matlabFunction(V(jj,ii),'Vars',[xs ys zs x_1s x_2s y_1s y_2s z_1s z_2s]); 450 Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii));
454 Vipart=matlabFunction(V(jj,ii),'Vars',[xs ys zs x_1s x_2s y_1s y_2s z_1s z_2s]); 451 Viret(jj,(ii-1)*side+1:side*ii)=eval(Vi(jj,ii));
455 Dret(jj,(ii-1)*side+1:side*ii)=sparse(Dpart(x,y,z,x_1,x_2,y_1,y_2,z_1,z_2));
456 Vret(jj,(ii-1)*side+1:side*ii)=sparse(Vpart(x,y,z,x_1,x_2,y_1,y_2,z_1,z_2));
457 Viret(jj,(ii-1)*side+1:side*ii)=sparse(Vipart(x,y,z,x_1,x_2,y_1,y_2,z_1,z_2));
458 end 452 end
459 end 453 end
460 454
461 D=sparse(Dret); 455 D=sparse(Dret);
462 V=sparse(Vret); 456 V=sparse(Vret);
463 Vi=sparse(Viret); 457 Vi=sparse(Viret);
464 V=obj.evaluateCoefficientMatrix(V,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2); 458 V=obj.evaluateCoefficientMatrix(V,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2);
465 D=obj.evaluateCoefficientMatrix(D,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2); 459 D=obj.evaluateCoefficientMatrix(D,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2);
470 zeroeig=(DD==0); 464 zeroeig=(DD==0);
471 negeig=(DD<0); 465 negeig=(DD<0);
472 466
473 D=diag([DD(poseig); DD(zeroeig); DD(negeig)]); 467 D=diag([DD(poseig); DD(zeroeig); DD(negeig)]);
474 V=[V(:,poseig) V(:,zeroeig) V(:,negeig)]; 468 V=[V(:,poseig) V(:,zeroeig) V(:,negeig)];
475 %Vi=inv(V); 469 Vi=[Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
476 signVec=[sum(poseig),sum(zeroeig),sum(negeig)]; 470 signVec=[sum(poseig),sum(zeroeig),sum(negeig)];
477 end 471 end
478 end 472 end
479 end 473 end