% mes_zadanie_9_10_11_12
%% (zadanie 9)
clear
close all
x=0:7.5:30;% <-------------- punkty oddalone co 7.5 #################
y=0:7.5:30;% <-------------- punkty oddalone co 7.5 #################

X = [x(1) x(length(x))];
Y = [y(1) y(length(y))];
  
[p,t]=mesPT(X,Y,10.61);
%               ^ 10.61 - maksymalna odleglosc miedzy punktami ###############
% load pt % plik z danymi do pobrania na stronie: ks.zut.edu.pl/mmwe/pt.mat

% iteracja po elementach
for e = 1:size(t,2), 
  KE = mes6_KE(p(1,t(1:3,e)),p(2,t(1:3,e)));% obliczenie KE dla ka¿dego elementu
  ke{e}=KE;% zapamiêtanie macierzy dla elementów w komórkach struktury ke 
end

% agregacja - z³o¿enie macierzy KE w globaln¹ macierz sztywnoœci
K = mes7_KG(t,ke);

SoK = size(K,1);
figure(1),
pcolor([K zeros(SoK,1);zeros(1,SoK+1)]);
set(gca,'YDir','reverse');
axis equal tight
title('Macierz sztywnoœci uk³adu')
colorbar
disp('"return" - wraca do obliczania')
keyboard


%____________________________________________________________
%% (zadanie 10)
% warunki brzegowe: V1, V2, V3, V4
V1 = y/3;
V2 = 10*cos(x*pi/15);
V3 = y/3;
V4 = 5*sin(x*pi/15);

V(1:size(p,2),1) = NaN;
V(1:5,1) = V4;%    <-------------- punkty oddalone co 7.5 #################
V(21:25,1) = V2;%  <-------------- punkty oddalone co 7.5 #################
V(1:5:21,1) = V1;% <-------------- punkty oddalone co 7.5 #################
V(5:5:25,1) = V3;% <-------------- punkty oddalone co 7.5 #################
disp('Macierz V0 - znane warunki brzegowe i niewiadome = NaN')
V0=V
disp('"return" - wraca do obliczania')
keyboard


%____________________________________________________________
%% (zadanie 11)
% uporz¹dkowanie macierzy K i V
[K,V,N,Np] = mes8_sort(K,V);


% macierz K po uporz¹dkowaniu - prezentacja na rysunku
figure,
pcolor([K zeros(SoK,1);zeros(1,SoK+1)]);
set(gca,'YDir','reverse');
axis equal tight
title('Macierz sztywnoœci uk³adu - uporz¹dkowana')
colorbar
disp('"return" - wraca do obliczania')
keyboard


%____________________________________________________________
%% (zadanie 12)
% budowanie uk³adu równañ
VZ = V(N+1:Np);
KNN = K(1:N,1:N);
KNZ = K(1:N,N+1:Np);

% rozwi¹zanie
disp('Rozwi¹zanie - potencja³y w wêz³ach w œrodku obszaru')
v=-KNN\KNZ*VZ

% prezentacja wyników
V0(isnan(V0))=v;
VV=reshape(V0,5,5);% <-------------- punkty oddalone co 7.5 ################
figure,
surf(VV)
