www.pudn.com > AFEM@matlab.rar > bisection.m, change:2006-11-16,size:7446b

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