bisection.m

C++编译指导

function [node,elem,marker,d2p] = bisection(node,elem,eta,theta,method)
% BISECTION	refine the current triangulation using 
%  longest bisection or newest vertex bisection
% 
% [node,elem] = bisection (node,elem,eta,theta,mode)
%
% Input:
%   node, elem : standard mesh data structure;
%   eta: error indicator for each triangle;
%   Dirichlet,Neumann: boundary edges;
%   theta: parameter in (0,1). 
%     We mark minimal number of triangles M such that
%       \sum _M eta > theta*sum(eta)
%   mode: 'longest' longest edge bisection; 
%         'newest' newest vertex bisection;
% 
% Output: 
%   node, elem : refined mesh
%   Dirichlet,Neumann: boundary edges;
%
% L. Chen 3-26-2006

% if mode is longest edge bisection, we relabel the elem
if (method==1), elem = label(node,elem); end

% construct data structure
edge = [elem(:,[1,2]); elem(:,[1,3]); elem(:,[2,3])];
edge = unique(sort(edge,2),'rows');
N = size(node,1); NT = size(elem,1); NE = size(edge,1);
dualedge = sparse(elem(:,[1,2,3]),elem(:,[2,3,1]),[1:NT,1:NT,1:NT],N,N);
d2p = sparse(edge(:,[1,2]),edge(:,[2,1]),[1:NE,1:NE],N,N);

% mark triangles according to the error indicator
total = sum(eta); [temp,ix] = sort(-eta); current = 0; marker = zeros(NE,1);
for t = 1:NT
    if (current > theta*total), break; end
    index=1; ct=ix(t);
    while (index==1)
       base = d2p(elem(ct,2),elem(ct,3)); 
       if (marker(base)>0), index=0;
       else
          current = current + eta(ct); 
		  N = size(node,1)+1;
		  marker(d2p(elem(ct,2),elem(ct,3))) = N; 
		  node(N,:) = mean(node(elem(ct,[2 3]),:));
		  ct = dualedge(elem(ct,3),elem(ct,2));
          if (ct==0), index=0; end
       end
    end
end

% refine marked edges for each triangle
for t = 1:NT
    base = d2p(elem(t,2),elem(t,3));
    if (marker(base)>0) 
        p = [elem(t,:), marker(base)]; 
        elem = divide(elem,t,p);
        left = d2p(p(1),p(2)); right = d2p(p(3),p(1));    
        if (marker(right)>0)
            elem = divide(elem,size(elem,1),[p(4),p(3),p(1),marker(right)]);
        end
        if (marker(left)>0)
            elem = divide(elem,t,[p(4),p(1),p(2),marker(left)]);
        end
  end
end

% plot refined mesh
hold off; trimesh(elem,node(:,1),node(:,2),zeros(size(node,1),1))
view(2),axis equal,axis off

% end of bisection

%-------------------------------------------------
% sub functions called by bisection
% function elem = divide(elem,t,p)
%-------------------------------------------------
function elem = divide(elem,t,p)
% DIVIDE divide triangle t into two triangels
%
% elem = divide(elem,t,p)
%
% Input
%   elem: element arrary
%   t: triangle to be divided
%   p: vertices of triangle t and new vertex
%    p(1),p(2),p(3) are vertices of t
%    p(4) new vertex which is opposite to p(1)
%
% Output
%   elem: t is divided into left and right triangles 
%     elem(t) is used to store left triangle
%     elem(NT+1) store the right triangle
%
% NOTE: the newest vertex p(4) is the first vertex
%         in each triangle
%
% L. Chen 03-26-2006

elem(size(elem,1)+1,:) = [p(4),p(3),p(1)];
elem(t,:) = [p(4),p(1),p(2)];