% 2D heat flow in L-shape domain
function abc()
clear

 k=0.7;     % conductivity coefficient

% discretization
 Nodes = 15;  % number of nodes
 Nel = 18;    % number of elemenets
 Ndofs = Nodes;    % one dof at each node

 Coord(1,1:2)=[8,0];  % coordinates of nodes
 Coord(2,1:2)=[6,2];
 Coord(3,1:2)=[8,4];
 Coord(4,1:2)=[5,4];
 Coord(5,1:2)=[4,3];
 Coord(6,1:2)=[4,0];
 Coord(7,1:2)=[0,0];
 Coord(8,1:2)=[2,2];
 Coord(9,1:2)=[3,4];
 Coord(10,1:2)=[4,4];
 Coord(11,1:2)=[4,5];
 Coord(12,1:2)=[0,4];
 Coord(13,1:2)=[2,6];
 Coord(14,1:2)=[4,8];
 Coord(15,1:2)=[0,8];

 Inc = [
    1 3 2;     % nodes of elements (topology)
    2 3 4;
    1 2 6;
    2 4 5;
    2 5 8;
    6 2 8;
    6 8 7;
    8 5 9;
    5 4 10;
    5 10 9;
    9 10 11;
    8 9 13;
    8 13 12;
    8 12 7;
    12 13 15;
    15 13 14;
    14 13 11;
    11 13 9];
 Edof=[[1:Nel]',Inc];

% plot discretization
 figure(1)   % elements
 [Ex,Ey]=coordxtr(Edof,Coord,[1:Nodes]',3);
 eldraw2(Ex,Ey,[1,2,1],Edof), axis equal;

 figure(2)  % nodes
 eldraw2(Ex,Ey,[3,2,0]), hold on, axis equal;
 for in=1:Nodes
     h=text(Coord(in,1),Coord(in,2),int2str(in));
     set(h,'fontsize',14,'fontweight','bold');
 end

% assembling
 KG=zeros(Ndofs,Ndofs); % global matrix
 f=zeros(Ndofs,1);      % global vector
 f(3)=-150;
 f(4)=-200;
 f(10)=-100;
 f(11)=-200;
 f(14)=-150;
 for iel=1:Nel
    Kel=flw2te(Ex(iel,:),Ey(iel,:),[1],[k 0;0 k]);
    KG=assem(Edof(iel,:),KG,Kel);
 end

 temp=40;       % Dirichlet b.c.
 bc=[1 temp; 6 temp; 7 temp; 12 temp; 15 temp];


% solution
 u=solveq(KG,f,bc);

% plot solution
 figure(3)
 colormap('jet')
 Ed = extract(Edof,u);
 fill(Ex',Ey',Ed')
 axis equal, colorbar;
 
 Tmin=min(u), Tmax=max(u)
