% Rechteckscheibe
clear all

% FEM-Modell erzeugen ...
[ne , nk , xy , km] = RectnetSV16 (100 , 20 , 2 , 1) ;         % FEM-Netz
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  ;  % Element-Parameter
kr(3,2)    = 1    ; kr(nk,1:2) = 1      ;                      % Lager
bk(6,2)    = -200 ; bk(1,1)    = 200    ;                      % Belastung
clf ;
draw2Dnet (xy , km , kr , bk , 'fk') ;               % Zeichnen des Berechnungsmodells
% ... und berechnen:
[succ uv] = femalg_m (xy , km , ep , kr , bk) ;      % Finite-Elemente-Algorithmus
if succ == 1
    disp (['FEM-Rechnung:           Verschiebung vVertikal =  ' num2str(uv(6,2))]) ;
else
    disp ('FEM-Rechnung:           Steifigkeitsmatrix nicht positiv definit') ;
end

% FEM-Gleichungssystem erzeugen ...
[succ Ks fq] = syswbc_m (xy , km , ep , kr , bk) ;   % FEM-Gleichungssystem
Ks = SymBand2Sparse (Ks) ;                           % Umwandlung in "Sparse matrix"
% ... und loesen
vS = Ks \ fq  ;                                      % Matlab-Standard-Solver
R  = chol(Ks) ;                                      % Cholesky (Definitheitstest)
vP = pcg (Ks , fq , 1e-3 , 100) ;                    % Preconditioned conjugate gradient
disp (['Matlab-Standard-Solver: Verschiebung vVertikal =  ' num2str(vS(12) ,'%11.4g')]) ;
disp (['Matlab-pcg:             Verschiebung vVertikal =  ' num2str(vP(12) ,'%11.4g')]) ;
disp (['Determinante der Koeffizientenmatrix:  det(Ks) =  ' num2str(det(Ks),'%11.4g')]) ;