function [mesh,eta] = bisection_Adap(mesh,eta,theta,Base)

% This subroutine is modified from bisection in FFW and IFEM 
% 
% ***** Thanks very much for these two softwares ********
%
% BISECTION	refines the triangulation using newest vertex bisection
%
% USAGE
%    [mesh, eta] = bisection_Adap_Multi(mesh,eta,theta, Base)
%
% INPUT 
%     mesh:  current mesh, in this struct, mesh.solu saves the current
%
%      eta:  error indicator for each triangle
%    theta:  parameter in (0,1). 
%            We mark minimal number of triangles M such that
%              \sum_{T \in M} \eta_T > \theta*\sum\eta_T
%    Base: the basis for the used finite element
%
% OUTPUT
%     mesh:  new mesh after refinement
%      eta:  new error indicator for each triangle
%
%==========================================================================
% Hehu Xie 11-29-2009

tic;
fprintf('----------Adaptive refine the mesh---------------------\n')
%--------------------------------------------------------------------------
% Construct data structure
%--------------------------------------------------------------------------
edge = mesh.node4edge;
N = size(mesh.node,1);   %number of nodes
NT = size(mesh.elem,1);  %nember of elements
NE = size(edge,1);       %number of edges
DIM = 3*Base.DOF_3(1)+3*Base.DOF_3(2)+Base.DOF_3(3); % the dimensions
% dualEdge(i,j) denotes the number of element T such that the nodes v_iv_j 
% is an edge of T
dualEdge = sparse(mesh.elem(:,[1,2,3]),mesh.elem(:,[2,3,1]),[1:NT,1:NT,1:NT]);
% d2p(i,j): denote the number of edge with its nodes v_i and v_j
d2p = sparse(edge(:,[1,2]),edge(:,[2,1]),[1:NE,1:NE]);
% Detailed explanation can be founded at


%--------------------------------------------------------------------------
% Meomery management for node arrary
%--------------------------------------------------------------------------
% recycle = find(mesh.type==0); 
% last = length(recycle);
% Coarsening can create empty spaces in the node array. We collect those
% scattered spaces in recycle arrary and 'last' will point out the last
% empty node index. 
% mesh.type array is used to distinguish the type of nodes:
%   0: empty nodes (deleted by coarsening);
%   1: nodes in the initial triangulation or nodes from regular refinement;
%   2: new added nodes due to refinement;
%   5: temporay node (will be deleted when bisection finished).

%--------------------------------------------------------------------------
% Mark triangles according to the error indicator
%--------------------------------------------------------------------------
total = sum(eta); 
current = 0; 
[temp,IND] = sort(eta,'descend');  %sort in descent order
marker = zeros(NE,1); 
% initialize for possible new nodes
mesh.node = [mesh.node; zeros(NE,2)];
mesh.type = ones(N,1);
mesh.type = [mesh.type; uint8(5*ones(NE,1))];
Ind_Solu = isfield(mesh,'solu');
solu_dim = 0;
if Ind_Solu & (nargin == 4)
    solu_dim = size(mesh.solu,2);
    mesh.solu = [mesh.solu; zeros(NE,solu_dim)];
    solu_element = zeros(NT,DIM,solu_dim);

    for kk = 1 : solu_dim

        solu_tmp = mesh.solu(:,kk);
        for i = 1:Base.DOF_3(1)
            solu_element(:,(i-1)*3+[1:3],kk) = [solu_tmp(mesh.elem(:,1:3))];
        end

        for i=1:Base.DOF_3(2)
            DOF_start = N*Base.DOF_3(1)+(i-1)*NE;
            solu_element(:,3*Base.DOF_3(1)+(i-1)*3+[1:3],kk) = [solu_tmp(DOF_start+mesh.edge(:,1:3))];
        end

        for i=1:Base.DOF_3(3)
            DOF_start = N*Base.DOF_3(1)+NE*Base.DOF_3(2)+(i-1)*NT;
            solu_element(:, 3*(Base.DOF_3(1)+Base.DOF_3(2))+i,kk) = [solu_tmp(DOF_start+1:DOF_start+NT)];
        end
    end
end

for t = 1:NT    
    if (current > theta*total), break; end  %eta on marked elem big enough    
    index = 1;     ct = IND(t);    
    while (index==1)        
        base = d2p(mesh.elem(ct,2),mesh.elem(ct,3));
        if (marker(base)>0) % base is already markedŁ¬ return
            index = 0;
        else
            current = current + eta(ct); %add this element to the marker
            newNode = N + 1;  % number of nodes plus 1
            N = N + 1;  % number of element plus 1
            % save the new node: N: marker save the new node number
            marker(d2p(mesh.elem(ct,2),mesh.elem(ct,3)) ) = newNode;
            % A new node is added to the mesh. Numerical solution at this
            % new added node is approximated by linear interpolation. Type
            % of this node is 2 (generated by newest vertex bisection)
            % get the coordination of the newnode
            mesh.node(newNode,:) = ( mesh.node(mesh.elem(ct,2),:) + ...
                                     mesh.node(mesh.elem(ct,3),:) )/2;
            %mesh.solu(newNode) = ( mesh.solu(mesh.elem(ct,2)) + ...
            %                       mesh.solu(mesh.elem(ct,3)) )/2;
            mesh.type(newNode) = 2;            
            % Find the element which shares the base edge of the current
            % element. If it is 0, it means the base of the current element
            % is on the boundary.
            % get the dual element of current element 
            ct = dualEdge( mesh.elem(ct,3), mesh.elem(ct,2) );
            if ( ct == 0 ), index = 0; end 
            % the while will ended if
            % 1.  ct==0 means we are on the boundary
            % 2.  base(ct) is already marked            
        end
    end % end while for recursive marking
end % end of for loop on all elements
% Detailed explanation of the algorithm can be found at
% Manual --> Algorithms --> Bisection
% delete possible empty entries
ix = (mesh.type == 5); 
mesh.node(ix,:) = []; 
mesh.type = [];
%rmfield(mesh,'type')

%--------------------------------------------------------------------------
% Refine marked edges for each triangle
%--------------------------------------------------------------------------
%--define the left basis and right hand basis------------------------------
if Ind_Solu & (nargin == 4)
    mesh_left = struct('node',[1/2 1/2;0 0;1 0], 'elem',[1 2 3],...
                'edge',[1 2 3],'node4edge',[2 3; 3 1; 1 2]);
          
    left_basis = Base.nodal_funct_3(Base.b_3, mesh_left)';
    %--------------------------------------------------------------------------
    mesh_right = struct('node',[1/2 1/2;0 1;0 0], 'elem',[1 2 3],...
                    'edge',[1 2 3],'node4edge',[2 3; 3 1; 1 2]);

    right_basis = Base.nodal_funct_3(Base.b_3, mesh_right)';
end
%--------------------------------------------------------------------------
numnew = 2*sum(marker~=0); % number of new elements need to be added
mesh.elem = [mesh.elem(:,1:3); zeros(numnew,3)];
Ind_Multilevel = isfield(mesh,'owner');
if Ind_Multilevel
    mesh.owner = [mesh.owner;zeros(numnew,1)];
    mesh.elemlocalnode = [mesh.elemlocalnode;zeros(numnew,9)];
end
eta = [eta; zeros(numnew,1)];
inew = NT + 1; % index for current new added right child
for t = 1:NT
    base = d2p(mesh.elem(t,2),mesh.elem(t,3));
    if (marker(base)>0)
        p = [mesh.elem(t,:), marker(base)];
        
        % Case 1: divide the current marked triangle
        mesh.elem(t,:)  = [p(4),p(1),p(2)]; % t is always a left child
        mesh.elem(inew,:) = [p(4),p(3),p(1)]; % new is a right child
        if Ind_Multilevel 
            mesh.owner(inew,1) = mesh.owner(t,1);
            tmpnode = mesh.elemlocalnode(t,:);
            mesh.elemlocalnode(t,:) = [(tmpnode(1,4:6)+tmpnode(1,7:9))/2,tmpnode(1,1:3),tmpnode(1,4:6)];
            mesh.elemlocalnode(inew,:) = [(tmpnode(1,4:6)+tmpnode(1,7:9))/2,tmpnode(1,7:9),tmpnode(1,1:3)];
        end
        eta(t) = eta(t)/2; % update error indicators 
        eta(inew) = eta(t);
        for kk=1:solu_dim
            solu_element(inew,:,kk) = solu_element(t,:,kk)*right_basis;
            solu_element(t,:,kk) = solu_element(t,:,kk)*left_basis;
        end       
        
        inew = inew + 1;
        
        % Case 2: divide the right child, different, careful!!!
        right = d2p(p(3),p(1));
        if (marker(right)>0)
            mesh.elem(inew-1,:) = [marker(right),p(4),p(3)]; 
            mesh.elem(inew,:)   = [marker(right),p(1),p(4)];
            if Ind_Multilevel
                mesh.owner(inew,1) = mesh.owner(inew-1,1);
                tmpnode = mesh.elemlocalnode(inew-1,:);
                mesh.elemlocalnode(inew-1,:) = [(tmpnode(1,4:6)+tmpnode(1,7:9))/2,tmpnode(1,1:3),tmpnode(1,4:6)];
                mesh.elemlocalnode(inew,:) = [(tmpnode(1,4:6)+tmpnode(1,7:9))/2,tmpnode(1,7:9),tmpnode(1,1:3)];
            end
            eta(inew-1) = eta(inew-1)/2;
            eta(inew) = eta(inew-1);
            for kk=1:solu_dim
                solu_element(inew,:,kk) = solu_element(inew-1,:,kk)*right_basis;
                solu_element(inew-1,:,kk) = solu_element(inew-1,:,kk)*left_basis;
            end    
            inew = inew + 1;
        end

        % Case 3: divide the left child, similar to the case 1. 
        left = d2p(p(1),p(2));
        if (marker(left)>0)
            mesh.elem(t,:) = [marker(left),p(4),p(1)]; 
            mesh.elem(inew,:) = [marker(left),p(2),p(4)]; 
            if Ind_Multilevel
                mesh.owner(inew,1) = mesh.owner(t,1);
                tmpnode = mesh.elemlocalnode(t,:);
                mesh.elemlocalnode(t,:) = [(tmpnode(1,4:6)+tmpnode(1,7:9))/2,tmpnode(1,1:3),tmpnode(1,4:6)];
                mesh.elemlocalnode(inew,:) = [(tmpnode(1,4:6)+tmpnode(1,7:9))/2,tmpnode(1,7:9),tmpnode(1,1:3)];
            end
            eta(t) = eta(t)/2;
            eta(inew) = eta(t);
            for kk=1:solu_dim
                solu_element(inew,:,kk) = solu_element(t,:,kk)*right_basis;
                solu_element(t,:,kk) = solu_element(t,:,kk)*left_basis;
            end            
            inew = inew + 1;
        end
    end % end of refinement of one element
end % end of for loop on all elements
% delete possible empty entries
mesh.elem = mesh.elem(1:inew-1,:);
if Ind_Multilevel
    mesh.owner = mesh.owner(1:inew-1,:);
    mesh.elemlocalnode = mesh.elemlocalnode(1:inew-1,:);
end
eta = eta(1:inew-1,:);
%--------------------------------------------------------------------------
% Update boundary edges
%--------------------------------------------------------------------------
mesh.Dirichlet_node = updatebd(mesh.Dirichlet_node,marker,d2p);
mesh.Neumann_node   = updatebd(mesh.Neumann_node,marker,d2p);
%--------reconstruct the mesh ---------------------------------------------
nrElems=size(mesh.elem,1);
nrNodes=size(mesh.node,1);
%--------------------- get e4n --------------------------------------------
I = [mesh.elem(:,1);mesh.elem(:,2);mesh.elem(:,3)];
J = [mesh.elem(:,2);mesh.elem(:,3);mesh.elem(:,1)];
S = [1:nrElems;1:nrElems;1:nrElems];
e4n = sparse(I,J,S);
%--------------------- geted4n --------------------------------------------
ed4n = triu(e4n+e4n');
[I,J] = find(ed4n);
nrEdges = length(I);
ed4n = sparse(I,J,1:nrEdges,size(e4n,1),size(e4n,1));
ed4n = ed4n+ed4n';
%-----------Initialize the boundary condition------------------------------
Nr_bd_Dirichlet=size(mesh.Dirichlet_node,1);
mesh.Dirichlet_edge=[];

if Nr_bd_Dirichlet>0
    mesh.Dirichlet_edge = diag(ed4n(mesh.Dirichlet_node(:,1),mesh.Dirichlet_node(:,2)));
end

%--------------------------------------------------------------------------
mesh.Neumann_edge=[];
Nr_bd_Neumann=size(mesh.Neumann_node,1);
if Nr_bd_Neumann>0
    mesh.Neumann_edge = diag(ed4n(mesh.Neumann_node(:,1),mesh.Neumann_node(:,2)));
end
%------------- get n4ed ---------------------------------------------------
[I,J,S] = find( triu(ed4n) );
mesh.node4edge( S,: ) = [I,J];
mesh.edge = rowaddr(ed4n,mesh.elem(:,[2 3 1]),mesh.elem(:,[3 1 2]));
%-------renew the solution ------------------------------------------------
if Ind_Solu == 1
    DIM_New = size(mesh.node,1)*Base.DOF_3(1)+size(mesh.node4edge,1)*Base.DOF_3(2)+size(mesh.elem,1)*Base.DOF_3(3);
    mesh.solu = sparse(DIM_New,solu_dim);
    for kk = 1:solu_dim
        mesh.solu(mesh.elem(:,1),kk)=solu_element(:,1,kk);
        mesh.solu(mesh.elem(:,2),kk)=solu_element(:,2,kk);
        mesh.solu(mesh.elem(:,3),kk)=solu_element(:,3,kk);

        for i=1:Base.DOF_3(2)
            DOF_start = nrNodes*Base.DOF_3(1)+(i-1)*nrEdges;
            mesh.solu(DOF_start+mesh.edge(:,1),kk) = solu_element(:,3*Base.DOF_3(1)+3*(i-1)+1,kk);
            mesh.solu(DOF_start+mesh.edge(:,2),kk) = solu_element(:,3*Base.DOF_3(1)+3*(i-1)+2,kk);
            mesh.solu(DOF_start+mesh.edge(:,3),kk) = solu_element(:,3*Base.DOF_3(1)+3*(i-1)+3,kk);
        end

        for i=1:Base.DOF_3(3)
            DOF_start = nrNodes*Base.DOF_3(1)+3*nrEdges*Base.DOF_3(2)+(i-1)*nrElems;
            mesh.solu(DOF_start+1:DOF_start+nrElems,kk) = solu_element(:,3*(Base.DOF_3(1)+Base.DOF_3(2))+i,kk);
        end
    end
end
toc
%--------------------------------------------------------------------------
% End of function BISECTION
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
% Sub functions called by BISECTION
%--------------------------------------------------------------------------
function bdEdge = updatebd(bdEdge,marker,d2p)
% UPDATEDBD refine the boundary edges
%
% USAGE
%    bdEdge = updatebd(bdEdge,marker,d2p)
%
% INPUT
%   bdEdge:  set of boundary edges
%   marker:  new node index for marked edge
%      d2p:  index mapping from dual edge to primary edge
% 
% OUTPUT
%   bdEdge:  set of refined boundary edges
NB = size(bdEdge,1);
for k = 1:NB
     i = bdEdge(k,1); 
     j = bdEdge(k,2);
     if marker(d2p(i,j)) >0
         bdEdge(k,:) = [i,marker(d2p(i,j))];
         bdEdge(size(bdEdge,1)+1,:) = [marker(d2p(i,j)),j];
     end
end
%--------------------------------------------------------------------------
% End of function UPDATEDBD
%--------------------------------------------------------------------------