bisection.m

五点差分型多重网格方法：各种插值算子的比较）

function [mesh,eta] = bisection(mesh,eta,theta)
% BISECTION	refines the triangulation using newest vertex bisection
%
% USAGE
%    [mesh,eta] = bisection(mesh,eta,theta)
%
% INPUT 
%     mesh:  current mesh
%      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
%
% OUTPUT
%     mesh:  new mesh after refinement
%      eta:  new error indicator for each triangle
%
% REFERENCE
% 	Long Chen,
%   Short bisection implementation in MATLAB
%   Research Notes, 2006 
%

% L. Chen & C. Zhang 11-12-2006

%--------------------------------------------------------------------------
% Construct data structure
%--------------------------------------------------------------------------
edge = [mesh.elem(:,[1,2]); mesh.elem(:,[1,3]); mesh.elem(:,[2,3])];
edge = unique(sort(edge,2),'rows');

N = size(mesh.node,1);
NT = size(mesh.elem,1); 
NE = size(edge,1);

dualEdge = sparse(mesh.elem(:,[1,2,3]),mesh.elem(:,[2,3,1]),[1:NT,1:NT,1:NT]);
d2p = sparse(edge(:,[1,2]),edge(:,[2,1]),[1:NE,1:NE]);
% Detailed explanation can be founded at
%   Manual --> Data Structure -->  Auxlliary data structure 

%--------------------------------------------------------------------------
% 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,ix] = sort(-eta); % sort in descent order
marker = zeros(NE,1); 

% initialize for possible new nodes
mesh.node = [mesh.node; zeros(NE,2)];
mesh.type = [mesh.type; uint8(5*ones(NE,1))];
mesh.solu = [mesh.solu; zeros(NE,1)];

for t = 1:NT
    
    if (current > theta*total), break; end % est on marked elem big enough
    
    index = 1; 
    ct = ix(t);
    
    while (index==1)
        
        base = d2p(mesh.elem(ct,2),mesh.elem(ct,3));
        if (marker(base)>0) % base is already marked
            index = 0;
        else
            current = current + eta(ct);
            if (last==0)
                newNode = N + 1; 
                N = N+1; 
            end
            if (last>0)
                newNode = recycle(last); 
                last = last-1; 
            end
            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)
            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.
            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(ix) = [];
mesh.solu(ix) = [];

%--------------------------------------------------------------------------
% Refine marked edges for each triangle
%--------------------------------------------------------------------------
numnew = 2*sum(marker~=0); % number of new elements need to be added
mesh.elem = [mesh.elem; zeros(numnew,3)];
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
        eta(t) = eta(t)/2; % update error indicators 
        eta(inew) = eta(t);
        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)];
            eta(inew-1) = eta(inew-1)/2;
            eta(inew) = eta(inew-1);
            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)];  
            eta(t) = eta(t)/2;
            eta(inew) = eta(t);
            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,:);
eta = eta(1:inew-1,:);

%--------------------------------------------------------------------------
% Update boundary edges
%--------------------------------------------------------------------------
mesh.Dirichlet = updatebd(mesh.Dirichlet,marker,d2p);
mesh.Neumann   = updatebd(mesh.Neumann,marker,d2p);

%--------------------------------------------------------------------------
% 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
%--------------------------------------------------------------------------