function [ mesh ] = refineMesh( mesh, toRefine )

index = find(toRefine); 
markedEdges  = zeros(mesh.NuEdges,1);
markedEdgesI = unique(reshape(mesh.triangle2edge(index,:),[],1));
markedEdges(markedEdgesI) = 1;
node1 = mesh.edges(:,1); node2 = mesh.edges(:,2);
h = (mesh.P(1,node1)-mesh.P(1,node2) ).^2 + ( mesh.P(2,node1)-mesh.P(2,node2) ).^2;
h = sqrt(h);
keepOn = true;
while(keepOn) 
  keepOn = false;
  for i=1:mesh.NuTri
    whichEdge = mesh.triangle2edge(i,:);
    howMany=sum( markedEdges(whichEdge)); 
    if ((howMany == 2) || (howMany == 1))
      [~,biggest] = sort(h(whichEdge)); 
      biggest = biggest(3);
      if ~markedEdges(whichEdge(biggest)) 
        markedEdges(whichEdge(biggest))=1; 
        keepOn = true;
      end
    end
  end
end
indexE = find(markedEdges); 
node1 = mesh.edges(indexE,1); node2 = mesh.edges(indexE,2); 
midpoint = (mesh.P(:,node1) + mesh.P(:,node2))/2;
mesh.P = horzcat(mesh.P,midpoint); 
newMidpointPos = sparse(length(mesh.edges),1);
newMidpointPos(indexE) = [1:length(indexE)] + mesh.NuDoF;
[~,ia ,ib] = intersect(indexE,mesh.BoundaryEdges);
start  = length(mesh.Boundary)+1; finish = length(mesh.Boundary)+length(ia);
temp=zeros(size(mesh.Boundary)); 
temp(mesh.Boundary)=mesh.GBoundary;
mesh.GBoundary(start:finish)= temp(ib);
mesh.Boundary(start:finish) = (ia+mesh.NuDoF)'; 
oldEdges = mesh.edges;
start = mesh.NuDoF+1;
finish = (mesh.NuDoF+ size(midpoint,2));
mesh.edges(indexE,2)= [start:finish]'; 
mesh.edges([1:length(indexE)] + mesh.NuEdges,1) = oldEdges(indexE,2);
mesh.edges([1:length(indexE)] + mesh.NuEdges,2) = start:finish;
newSegmentPos = sparse(length(mesh.edges),1);
newSegmentPos(indexE) = [1:length(indexE)] + mesh.NuEdges;
temp=zeros(size(oldEdges));
temp(mesh.BoundaryEdges)=mesh.GBoundaryEdges;
start=length(mesh.GBoundaryEdges)+1; finish=length(mesh.GBoundaryEdges)+length(ib);
mesh.GBoundaryEdges(start:finish)= temp(mesh.BoundaryEdges(ib));
for i=1:mesh.NuTri
  markedEdgesI = markedEdges(mesh.triangle2edge(i,:));
  howMany = sum(markedEdgesI);
  if howMany == 3
    nodeOldNotDiv = unique(reshape(oldEdges(mesh.triangle2edge(i,:),:),[],1));
    newSegmentIdx(1:3) = full(newSegmentPos(mesh.triangle2edge(i,1:3)));
    
    newEdge = zeros(3,1);
    for j=1:3 
      newTri(1) = nodeOldNotDiv(j);
      edgesOldPos  = find( nodeOldNotDiv(j)==mesh.edges(mesh.triangle2edge(i,:),1)); 
      edgesNewPos  = find( nodeOldNotDiv(j)==mesh.edges(newSegmentIdx,1)); 
      edgesOldPosNu = length(edgesOldPos);
      edgesNewPosNu = length(edgesNewPos);
      
      mesh.triangle2edge(end+1,1:3)=[0,0,0];
      if edgesOldPosNu>0 
        newTri(2:edgesOldPosNu+1) = mesh.edges(mesh.triangle2edge(i,edgesOldPos),2);
        mesh.triangle2edge(end,1:edgesOldPosNu) = mesh.triangle2edge(i,edgesOldPos);
      end
      if edgesNewPosNu>0 
        newTri((edgesOldPosNu+2):3) = mesh.edges(newSegmentIdx(edgesNewPos),2);
        mesh.triangle2edge(end,(edgesOldPosNu+1):2) = newSegmentIdx(edgesNewPos); 
      end 
      
      mesh.T(end+1,:) = newTri; 
      newEdge(j) = length(mesh.edges)+1;
      mesh.edges(newEdge(j),1:2) = sort( [newTri(2), newTri(3)]); 
      mesh.triangle2edge(end,3) = newEdge(j);
      mesh.GT(end+1)=mesh.GT(i);
    end
    mesh.triangle2edge(i,1:3) = newEdge(1:3);
    newTri = unique(reshape(mesh.edges(newEdge(1:3),:),[],1));
    mesh.T(i,:) = newTri;  
  elseif howMany == 2
    edgeLIdxNotDiv = find(~markedEdgesI);
    edgeGIdxNotDiv = mesh.triangle2edge(i,edgeLIdxNotDiv);
    edgeLIdxDiv = find(markedEdgesI);
    edgeGIdxDiv = mesh.triangle2edge(i,edgeLIdxDiv);
    edgesGIdxSplit = [edgeGIdxDiv, full(newSegmentPos(edgeGIdxDiv)')];
    
    [~,biggest] = sort(h(mesh.triangle2edge(i,:))); 
    biggest = biggest(3);
    if edgeLIdxDiv(1) ~= biggest
      edgeLIdxDiv(2) = edgeLIdxDiv(1); edgeLIdxDiv(1) = biggest;
    end;  
    nodeMl  = newMidpointPos(mesh.triangle2edge(i,edgeLIdxDiv(1)));
    nodeMs = newMidpointPos(mesh.triangle2edge(i,edgeLIdxDiv(2)));
    node1 = oldEdges(mesh.triangle2edge(i,edgeLIdxNotDiv),1);
    node2 = oldEdges(mesh.triangle2edge(i,edgeLIdxNotDiv),2);
    node3 = setdiff(mesh.T(i,:), [node1, node2]);
    nodesDiv2 = oldEdges(mesh.triangle2edge(i,edgeLIdxDiv(2)),:); 
    if any(nodesDiv2 == node1) 
      temp = node1; node1=node2; node2=temp;
    end 
    
    e1  = [node2, nodeMl];
    e2 = sort([nodeMs, nodeMl]);
    edgeGIdxNewE1  = length(mesh.edges)+1;
    mesh.edges(edgeGIdxNewE1,:)  = e1;
    edgeGIdxNewE2 = length(mesh.edges)+1;
    mesh.edges(edgeGIdxNewE2,:) = e2;
    mesh.T(i,:) = [nodeMl, node1, node2];
    [row,~] = find(mesh.edges(edgesGIdxSplit,:) == node1);
    edgeGIdxNode1 = edgesGIdxSplit(row);
    edgesGIdxSplit(row) = [];
    mesh.triangle2edge(i,:) = [edgeGIdxNewE1, edgeGIdxNotDiv, edgeGIdxNode1];
    
    mesh.T(end+1,:) = [nodeMl, nodeMs, node2];
    [row,~] = find(mesh.edges(edgesGIdxSplit,:) == node2);
    edgeGIdxNode2 = edgesGIdxSplit(row);
    edgesGIdxSplit(row) = [];
    mesh.triangle2edge(end+1,:) = [edgeGIdxNewE1, edgeGIdxNewE2, edgeGIdxNode2];
    mesh.GT(end+1) = mesh.GT(i);
    
    mesh.T(end+1,:) = [nodeMl, nodeMs, node3];
    mesh.triangle2edge(end+1,:)=[edgeGIdxNewE2,edgesGIdxSplit(1),edgesGIdxSplit(2)];
    mesh.GT(end+1) = mesh.GT(i);    
  elseif howMany == 1
    edgeLIdxDiv = find(markedEdgesI);  
    edgeGIdxDiv = mesh.triangle2edge(i,edgeLIdxDiv);
    edgesLIdxNotDiv = find(~markedEdgesI);
    edgesGIdxNotDiv = mesh.triangle2edge(i,edgesLIdxNotDiv);
    edgesGIdxSplit = [edgeGIdxDiv, full(newSegmentPos(edgeGIdxDiv))];  
    
    if ~any(oldEdges(edgesGIdxNotDiv(1),:) == mesh.edges(edgesGIdxSplit(1),1))  
      edgesGIdxSplit([1,2]) = edgesGIdxSplit([2,1]);
    end 
    
    nodesOld = mesh.T(i,:);  
    nodesOldDiv = oldEdges(edgeGIdxDiv,:);
    nodeOldNotDiv = nodesOld(nodesOld~=nodesOldDiv(1) & nodesOld~=nodesOldDiv(2));  
    nodeCenterDivEdge = newMidpointPos(mesh.triangle2edge(i,edgeLIdxDiv));
    edgeNew = [nodeOldNotDiv, nodeCenterDivEdge];
    edgeGIdxNew = length(mesh.edges)+1;
    mesh.edges(edgeGIdxNew,:) = edgeNew;
    mesh.T(end+1,:) = [nodeCenterDivEdge, oldEdges(edgesGIdxNotDiv(1),:)];
    mesh.triangle2edge(end+1,:) = [edgesGIdxNotDiv(1), edgesGIdxSplit(1), edgeGIdxNew];
    mesh.GT(end+1) = mesh.GT(i);
    
    mesh.T(i,:) = [nodeCenterDivEdge, oldEdges(edgesGIdxNotDiv(2),:)];
    mesh.triangle2edge(i,:) = [edgesGIdxNotDiv(2), edgesGIdxSplit(2), edgeGIdxNew];
  end
end
Point2Element(length(mesh.P)).Elements=[];
for i=1:length(mesh.P)
  Point2Element(i).Elements=[];
end
for i=1:length(mesh.T)
  for j=1:3
    index = length(Point2Element(mesh.T(i,j)).Elements)+1;
    Point2Element(mesh.T(i,j)).Elements(index)=i;
  end
end
mesh.Point2Element=Point2Element;
B(:,1,:)=mesh.P(:,mesh.T(:,2)) - mesh.P(:,mesh.T(:,1));
B(:,2,:)=mesh.P(:,mesh.T(:,3)) - mesh.P(:,mesh.T(:,1));
Bdet = B(1,1,:).*B(2,2,:) - B(1,2,:).*B(2,1,:);
mesh.B=B; 
mesh.Bdet=squeeze(Bdet);
mesh.NuDoF    = length(mesh.P);
mesh.NuBNodes = length(mesh.Boundary);
mesh.NuEdges  = length(mesh.edges);
mesh.NuTri    = length(mesh.T);
mesh.edge2triangle=-1*ones(mesh.NuEdges,2); 
for i=1:mesh.NuTri
  index = mesh.triangle2edge(i,:);
  for j=1:3
    if mesh.edge2triangle(index(j),1) == -1
      mesh.edge2triangle(index(j),1) = i;
    else
      mesh.edge2triangle(index(j),2) = i;
    end
  end
end
mesh.BoundaryEdges = find(mesh.edge2triangle(:,2) == -1);
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
