% Gelochte Scheibe
clear all
nr =  80 ; ny =  40 ;
b  = 150 ; h  = 100 ; r = 50 ;
% Koordinatenmatrix und Koinzidenzmatrix fuer Scheibe mit nr*ny Elementen:
[ne , nk , xy , km] = HoleNetSV16 (b , h , r , nr , ny) ;

disp (['Lochscheibe, Anzahl der Elemente:       ' num2str(ne)]) ;
disp (['             Anzahl der Knoten:         ' num2str(nk)]) ;
disp (['             Anzahl der Gleichungen:    ' , num2str(nk*2)]) ;

kr = zeros (nk,2) ; bk = zeros (nk,2) ; ep = zeros (ne,3) ;
ep(1:ne,1) = 5 ; ep(1:ne,2) = 210000 ; ep(1:ne,3) = 0.3 ;               
kr(1:ny*2+1,1) = 1 ; kr(nk-ny*2:nk,2) = 1 ; 
bk(nk-ny*2,1) = -10000 ;                          
clf ; draw2Dnet (xy , km , kr , bk , 'f') ;       % Zeichnen des Berechnungsmodells

% FEM-Gleichungssystem Ks*v=fq:
[succ Ks fq]   = syswbc_m (xy , km , ep , kr , bk) ;
Ks = SymBand2Sparse (Ks) ;                        % Umwandlung in "Sparse matrix"

t1 = cputime ;
vD = Ks \ fq ;                                    % Matlab-Standard-Solver
disp (['             Zeit Standard-Solver:      ' num2str(cputime-t1)]) ;
disp (['             Vertikale Lochaufweitung:  ' num2str(vD(ny*2+2)*2)]) ;

t1 = cputime ;
R  = cholinc(Ks , 1e-3) ;
disp (['             Nicht-Null-Elemente in R:  ' num2str(nnz(R))]) ;
vcgR = pcg (Ks , fq , 1e-2 , nk , R' , R) ;       % Preconditioned conjugate gradient
disp (['             Zeit pcg (incompl. chol.): ' num2str(cputime-t1)]) ;
disp (['             Vertikale Lochaufweitung:  ' num2str(vcgR(ny*2+2)*2)]) ;