m = load_gmsh2('grid.msh', [1 2]);
P(1,:) = m.POS(:,1); P(2,:) = m.POS(:,2);         
T      = m.TRIANGLES(:,1:3); GT     = m.TRIANGLES(:,4);    
L      = m.LINES(:,1:2); GL     = m.LINES(:,3);

mesh = buildMesh( T,P,[],GT,L,GL );
clear m T P GT L GL 

eps = 2*10^-5*ones(1,mesh.NuDoF); 
index  = (mesh.P(1,:)<=3.2) & (mesh.P(1,:)>=3) & (mesh.P(2,:)<=4) & (mesh.P(2,:)>=1.5);
eps(index) = 8.75*10^-7;

A = assembDiffusion (eps, mesh.T, mesh.P);
f = zeros(1,mesh.NuDoF);
b = assembRightSide(mesh, f);
[ A, b ] = applyDirichletBoundaryConditions(A,b, mesh, 20, @(x,y)323.13);
[ A, b ] = applyDirichletBoundaryConditions(A,b, mesh, 21, @(x,y)283.15);
u=A\b;
trisurf(mesh.T,mesh.P(1,:),mesh.P(2,:),u);   
xlabel('x [m]'); ylabel('y [m]'); zlabel('T [Kelvin]');
colormap('hot'); caxis([260 330]);colorbar;
