function [ mesh ] = buildMesh( T,P,Boundary,GT,L,GL )

% remove unused points in mesh
isPointUsed = zeros(1, length(P)); 
isPointUsed(T) = true;             
for i = find(~isPointUsed)        
    P(:,i) = [];
    T(T>i) = T(T>i) - 1;
    L(L>i) = L(L>i) - 1;
end

% Duplication of mesh values
mesh.T=T;
mesh.P=P;
mesh.NuDoF = length(P);
mesh.NuTri = length(T);

% Calculation of the transformation-matrices
B(:,1,:)=P(:,T(:,2)) - P(:,T(:,1));
B(:,2,:)=P(:,T(:,3)) - P(:,T(:,1));
Bdet = B(1,1,:).*B(2,2,:) - B(1,2,:).*B(2,1,:);
mesh.B=B;
mesh.Bdet=squeeze(Bdet);
% Creating the mapping from points to triangles
Point2Element(length(P)).Elements=[];
for i=1:length(P)
  Point2Element(i).Elements=[];
end
for i=1:length(T)
  for j=1:3
    index = length(Point2Element(T(i,j)).Elements)+1;
    Point2Element(T(i,j)).Elements(index)=i;
  end
end
mesh.Point2Element=Point2Element;
% Creating the edges and the mapping to triangles
% Mapping the triangles to the edges
edgesBloated = zeros(3*mesh.NuTri, 3);
for i=1:length(T)
  points = sort(T(i,:));
  edgesBloated(3*i-2:3*i,:) = [points(1), points(2), i; ...
  points(1), points(3), i; points(2), points(3), i];
end
[edgesBloated, tripleTriIndex] = sortrows(edgesBloated); 
equalEdges = edgesBloated(1:end-1,1:2) == edgesBloated(2:end,1:2);  
doubleEdge = all(equalEdges')';
doubleEdge(end+1) = false;  

edges    = edgesBloated(~doubleEdge, 1:2); 
e2t      = -ones(length(edges), 2); 
e2t(:,1) = edgesBloated(~doubleEdge, 3); 
t2e  = zeros(3, mesh.NuTri);  
dind = find(doubleEdge);
prevDi = 0;
for i = 1:length(dind)  
  di = dind(i);
  e2t(di - i + 1, 2) = edgesBloated(di, 3);
  t2e(tripleTriIndex(prevDi+1:di)) = (prevDi+1:di) - i + 1;
  prevDi = di;
end
i  = length(dind) + 1;
di = length(doubleEdge);
t2e(tripleTriIndex(prevDi+1:di)) = (prevDi+1:di) - i + 1;
mesh.edges = sort(edges,2);
mesh.edge2triangle = e2t;
mesh.triangle2edge = t2e';
% Build the edge-neighbourhood of the triangles
mesh.neighbours = zeros(mesh.NuTri,3);
for i = 1:mesh.NuTri
  tries = mesh.edge2triangle(mesh.triangle2edge(i,:),:); 
  tries(tries == i) = [];
  mesh.neighbours(i,:) = sort(tries, 2, 'descend');
end
% Determine the boundary values
if ~isempty(GT)
  mesh.GT=GT;
else
  mesh.GT=sparse(1,mesh.NuTri);
end

if ~isempty(Boundary) 
  mesh.NuBoundaryNodes = length(Boundary);
  mesh.Boundary=Boundary;
  mesh.GBoundary = sparse(1,mesh.NuBoundaryNodes);
else
  if (2 ~= size(L,2) ) 
    L=L';
  end
  BInfo(1,:) = reshape(L',1,[]);
  [~,ia,~] = unique(BInfo(1,:));
  mesh.Boundary = BInfo(1,ia);
  mesh.NuBoundaryNodes = length(mesh.Boundary);
  if ~isempty(GL)
    BInfo(2,1:2:length(GL)*2) = GL;
    BInfo(2,2:2:length(GL)*2+1) = GL;
    mesh.GBoundary = BInfo(2,ia);
  else
    mesh.GBoundary = sparse(1,mesh.NuBoundaryNodes);
  end
end
%Define boundary edges
mesh.BoundaryEdges = find(mesh.edge2triangle(:,2) == -1);
mesh.GBoundaryEdges = zeros(size(mesh.BoundaryEdges));
temp = -ones(mesh.NuDoF,1); 
temp(mesh.Boundary) = mesh.GBoundary;
for i=1:length(mesh.BoundaryEdges)
  firstNode  = mesh.edges(mesh.BoundaryEdges(i),1);
  secondNode = mesh.edges(mesh.BoundaryEdges(i),2);
  if temp(firstNode) == temp(secondNode)
    mesh.GBoundaryEdges(i) = temp(firstNode);
  else 
    index1 = find( L(:,1) == firstNode );
    index2 = find( L(index1,2) == secondNode );
    if ~isempty(index2)
      found = index1(index2);
      mesh.GBoundaryEdges(i) = GL(found);
    else
      index1 = find( L(:,2) == firstNode );
      index2 = find( L(index1,1) == secondNode );
      if ~isempty(index2)
        found = index1(index2);
        mesh.GBoundaryEdges(i) = GL(found);
      else
        msg = 'Egde with mixed boundary conditions detected but not found!'; 
        error(msg);
        exit;
      end
    end
  end
end
mesh.NuEdges = length(mesh.edges);
mesh.NuBNodes = length(mesh.Boundary);
