% 2D przeplyw ciepla w obszarze L 
function abc()

%%
%%       7      8  
%%       i``````i
%%       i      iqn=-12    ( odplyw ciepla )
%%       i      i      qn=-12
%%      4i      5``````````````i 6
%%  T=10 i                     i
%%       i                     i
%%       i                     i
%%      1i______2______________3
%%          T= 10
%%

clear

%%  ------------ PREPROCESSING ---------------------------
%%  stale przewodnictwa cieplnego
 kx= 4;    
 ky= 7;

%% Dyskretyzacja
 Lwez = 8;       % liczba wezlow   
 Lel = 6;        % liczba elementow   
 Lssuk = Lwez;   % jeden st. sw. w kazdym wezle

%% wspolrzedne kolejnych wezlow  
%       x y
 xyWez= [0 0 
         2 0
         5 0
         0 2
         2 2
         5 2
         0 3
         2 3
         ];

%% macierz topologii 
 mTop = [2 3 6
         2 6 5
         1 2 5
         1 5 4
         4 5 8
         4 8 7
        ];

 Nrssw= [[1:Lel]', mTop];
% [1:Lel]'

%% rysunek siatki      
 figure(1)  
 [Ex,Ey]= coordxtr(Nrssw,xyWez,[1:Lwez]',3);
 eldraw2(Ex,Ey,[1,2,1],Nrssw), axis equal;
 for in= 1:Lwez
     h= text(xyWez(in,1),xyWez(in,2),int2str(in));
     set(h,'fontsize',14,'Color','r'); %,'fontweight','bold');
 end

%% =================OBLICZENIA ======================================
 K= zeros(Lssuk,Lssuk); % macierz 
 f= zeros(Lssuk,1);      % wektor prawej strony

 qn= -12; % odplyw ciepla
 f(6)= qn*3/2; f(8)= qn*1/2; f(5)= f(6)+f(8); 

%%  agregacja  
 for el= 1:Lel
     Kel= flw2te(Ex(el,:),Ey(el,:),[1],[kx 0;0 ky]);
     K= assem(Nrssw(el,:),K,Kel);
 end

%% Warunki brzegowe (podstawowe) Dirichleta 
%     wz  T
 bc= [1  10;
      2  10;
      3  10;
      4  10;
      7  10
     ];

%% rozwiazanie  K*Theta = f+f_b  (f_b=0 --> brak zrodla ciepla)
 [theta,r]= solveq(K,f,bc);

 fprintf('\n\n  W Y N I K I:\n\n');
 fprintf(' Temperatura i reakcje\n');
 fprintf('  wz     theta       r\n')
 for  we= 1: Lwez
  fprintf(' %3d  %9.4f %9.4f\n', we, theta(we),r(we));
 end


%% rysunek rozwiazania
 figure(3)
 colormap('jet');
 theta_el= extract(Nrssw, theta);
 fill(Ex',Ey',theta_el')
 axis equal, colorbar;
 title(' temperatura T [^oC]');

 fprintf('\n Wektor intensywnosci strumienia przeplywu ciepla w elementach :\n');
 fprintf('     q^el= -D*B^e*theta^el\n');
 fprintf('  el      qx        qy\n')
 for el= 1:Lel
   qe= flw2ts(Ex(el,:),Ey(el,:), [kx 0;0 ky], theta_el(el,:));
   fprintf(' %3d  %9.4f %9.4f\n', el, qe);
 end
 fprintf('\n');

end


