function cieplo()
%% SKRYPT OCTAVE / CALFEM DO ANALIZY MES USTALONEGO PRZEPŁYWU CIEPŁA
% AUTOR: Piotr Pluciński, PK WIL L-10

clf

Edof=[ 1 1 2 3 ; 2 1 3 4 ] ; % nr_elementu, numery_stopni_swobody_w_elementach

K=zeros(4); 
P=zeros(4,1);
kx = 4; % współczynnik przewodności cieplnej wzdłuż x [J / Cms]
ky = 4; % współczynnik przewodności cieplnej wzdłuż y [J / Cms]
f = 45; % intensywność generacji ciepła wewnątrz ciała [J / m2s]
D=[kx 0 ; 0 ky]; 

Coord=[0 0 ; 2 0 ; 2 1 ; 0 1]; % współrzędne węzłów

Dof=(1:4)';  % stopnie swobody globalne w poszczególnych węzłach(1:ilość węzłów)

[Ex,Ey]=coordxtr(Edof,Coord,Dof,3); 


figure(1)
eldraw2(Ex,Ey, [1,2,2], Edof(:,1)); 
hold on;
grid on;
axis equal;
title('DYSKRETYZACJA MES ZADANIA USTALONEGO PRZEPŁYWU CIEPŁA');
xlabel('x');
ylabel('y');

%% 

%[K1,P1]=flw2te(Ex(1,:),Ey(1,:),1,D,f);
%[K2,P2]=flw2te(Ex(2,:),Ey(2,:),1,D,f);

%[K,P]=assem(Edof(1,:),K,K1,P,P1);
%[K,P]=assem(Edof(2,:),K,K2,P,P2);

 for el=1:size(Edof,1)
     [Ke,Pe]=flw2te(Ex(el,:),Ey(el,:),1,D,f);
     [K,P]=assem(Edof(el,:),K,Ke,P,Pe);
 end

% budowa wektora prawej strony układu równań MES
P=P+flux(Coord,[3 4],-30); %% 4 nr stopni swobody na brzegu ze definiowanym strumieniem, -30 wartość strumienia,
% znakowanie strumienia: dodatni gdy ciepło wpływa do ciała; ujemny gdy ciepło wypływa z ciała
%% ile mamy brzegów ze strumieniem tyle razy trzeba powtórzyć powyższą linię

bc=[ 2 10; 3 10]; % wektor warunków brzegowych (nr węzła/ss - wartość temp.)

%% ROZWIĄZANIE UKŁADU RÓWNAŃ MES i RYSUNEK TEMPERATURY
[T,R]=solveq(K,P,bc) 



%% POWRÓT DO ELEMENTU - WYZNACZENIE WARTOŚCI WEKTORA STRUMIENIA CIEPŁA
Te=extract_ed(Edof,T);

%Q(1,:)=flw2ts(Ex(1,:),Ey(1,:),D,Te(1,:));
%Q(2,:)=flw2ts(Ex(2,:),Ey(2,:),D,Te(2,:));


 for el=1:size(Edof,1)
     Q(el,:)=flw2ts(Ex(el,:),Ey(el,:),D,Te(el,:));
 end
%% 

Q

% mapa temperatur 
sfac=scalfact2(Ex,Ey,Q,0.25);

figure(2)
filloct(Ex',Ey',Te')
title(['TEMPERATURE [^oC] : T_m_a_x = ' num2str(max(T)) ', T_m_i_n = ' num2str(min(T))]);
xlabel('x');
ylabel('y');
zlabel('T');
colorbar
colormap jet
axis equal


figure(3)
eldraw2(Ex,Ey, [1,2,2]); 
elflux2(Ex,Ey,Q,[1,4],sfac);
eliso2(Ex,Ey,Te,10,[2,1]);
axis equal


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Z=flux(coord,wez,q)

Z=zeros(size(coord,1),1)
L=sqrt((coord(wez(1),1)-coord(wez(2),1))^2+(coord(wez(1),2)-coord(wez(2),2))^2)
Z(wez)=q*L/2;