bisection.m

C++编译指导

function mesh= bisection(mesh,eta,res,osc,theta1,theta2)
% BISECTION	refines the triangulation using newest vertex bisection
%
% USAGE
%    mesh = bisection(mesh,eta,theta)
%
% INPUT 
%     mesh:  current mesh
%    theta1,theta2:  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
%
% REFERENCE
% 	Long Chen,
%   Short bisection implementation in MATLAB
%   Research Notes, 2006 

%--------------------------------------------------------------------------
% Construct data structure
%--------------------------------------------------------------------------
mesh.elem = label(mesh.node,mesh.elem);
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.

%--------------------------------------------------------------------------
% Mark triangles according to the error indicator
%--------------------------------------------------------------------------
total =max(eta); 
[temp,ix] = sort(-eta); % sort in descent order
marker1 = zeros(NT,1);
marker2 = zeros(NE,1);marker3= zeros(NT,1);
for t = 1:NT
   index = 1; ct = ix(t);i=0; 
    if (eta(ct) <(theta1^2)*total), break; end % err on marked elem big enough
    if(marker1(ct)==0)
       marker1(ct)=ct;
       while (index==1)
           base = d2p(mesh.elem(ct,2),mesh.elem(ct,3));
           if (marker2(base)>0)
               index = 0;% base is already marked
                % mark the interior node
               if (i==0)
                   if (last==0), newNode = size(mesh.node,1) + 1; end
                   if (last>0),  newNode = recycle(last); last = last-1; end
                    marker3(ct) = newNode;
                    mesh.node(newNode,:) = ( mesh.node(mesh.elem(ct,1),:)...
                    + mesh.node(marker2(base),:))/2;
                    mesh.type(newNode) = 2; % newNode is added by the refinement
                    mesh.solu(newNode) = ( mesh.solu(mesh.elem(ct,1))...
                       +mesh.solu(marker2(base)))/2;
                    i=i+1;
               end 
           end
          if (marker2(base)==0)
            if (last==0), newNode = size(mesh.node,1) + 1; end
            if (last>0),  newNode = recycle(last); last = last-1; end
            marker2( 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.
            mesh.node(newNode,:) = ( mesh.node(mesh.elem(ct,2),:)...
                                    + mesh.node(mesh.elem(ct,3),:) )/2;
            mesh.type(newNode) = 2; % newNode is added by the refinement
            mesh.solu(newNode) = ( mesh.solu(mesh.elem(ct,2))...
                                   + mesh.solu(mesh.elem(ct,3)) )/2;
            % mark the interior node
            if (i==0)
                if (last==0), newNode = size(mesh.node,1) + 1; end
                if (last>0),  newNode = recycle(last); last = last-1; end
                marker3(ct) = newNode;
                mesh.node(newNode,:) = ( mesh.node(mesh.elem(ct,1),:) ...
                +mesh.node( marker2( d2p(mesh.elem(ct,2),mesh.elem(ct,3))),:))/2;
                mesh.type(newNode) = 2; % newNode is added by the refinement
                mesh.solu(newNode) = ( mesh.solu(mesh.elem(ct,2))  ...
                      + mesh.solu(marker2( d2p(mesh.elem(ct,2),mesh.elem(ct,3)))))/2;
               i=i+1;
            end   
                               
         
            % 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 % base is on the boundary
            % the while will ended if
            % 1.  ct==0 means we are on the boundary
            % 2.  base(ct) is already marked
           end
      end
        
        
        
        
        
           ct=ix(t);index=1;
         base = d2p(mesh.elem(ct,3),mesh.elem(ct,1));
        if (marker2(base)>0), index = 0;end % base is already marked
        if (marker2(base)==0)
            if (last==0), newNode = size(mesh.node,1) + 1; end
            if (last>0),  newNode = recycle(last); last = last-1; end
            marker2( d2p(mesh.elem(ct,3),mesh.elem(ct,1)) ) = newNode;

            % A new node is added to the mesh. Numerical solution at this
            % new added node is approximated by linear interpolation.
            mesh.node(newNode,:) = ( mesh.node(mesh.elem(ct,3),:)...
                                     +mesh.node(mesh.elem(ct,1),:) )/2;
            mesh.type(newNode) = 2; % newNode is added by the refinement
            mesh.solu(newNode) = ( mesh.solu(mesh.elem(ct,3))...
                                  + mesh.solu(mesh.elem(ct,1)) )/2;
        end
           ct = dualEdge(mesh.elem(ct,1),mesh.elem(ct,3));
        while(index==1)
            if (ct==0),index=0;end
            if (ct>0)
              base = d2p(mesh.elem(ct,2),mesh.elem(ct,3));
              if (marker2(base)>0), index = 0; end% base is already marked
              if (marker2(base)==0)
                  if (last==0), newNode = size(mesh.node,1) + 1; end
                  if (last>0),  newNode = recycle(last); last = last-1; end
                  marker2( 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.
                 mesh.node(newNode,:) = ( mesh.node(mesh.elem(ct,2),:)...
                                    + mesh.node(mesh.elem(ct,3),:) )/2;
                 mesh.type(newNode) = 2; % newNode is added by the refinement
                 mesh.solu(newNode) = ( mesh.solu(mesh.elem(ct,2))...
                                       + mesh.solu(mesh.elem(ct,3)) )/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)); 
              end % base is on the boundary
            % the while will ended if
            % 1.  ct==0 means we are on the boundary
            % 2.  base(ct) is already marked
           end
        end
        
           ct=ix(t);index=1;
         base = d2p(mesh.elem(ct,1),mesh.elem(ct,2));
        if (marker2(base)>0), index = 0;end % base is already marked
        if (marker2(base)==0)
            if (last==0), newNode = size(mesh.node,1) + 1; end
            if (last>0),  newNode = recycle(last); last = last-1; end
            marker2( d2p(mesh.elem(ct,1),mesh.elem(ct,2)) ) = newNode;

            % A new node is added to the mesh. Numerical solution at this
            % new added node is approximated by linear interpolation.
            mesh.node(newNode,:) = ( mesh.node(mesh.elem(ct,1),:)...
                                     +mesh.node(mesh.elem(ct,2),:) )/2;
            mesh.type(newNode) = 2; % newNode is added by the refinement
            mesh.solu(newNode) = ( mesh.solu(mesh.elem(ct,3))...
                                   + mesh.solu(mesh.elem(ct,1)) )/2;
        end
           ct = dualEdge(mesh.elem(ct,2),mesh.elem(ct,1));
        while(index==1)
            if (ct==0),index=0;end
            if (ct>0)
              base = d2p(mesh.elem(ct,2),mesh.elem(ct,3));
              if (marker2(base)>0), index = 0; end% base is already marked
              if (marker2(base)==0)
                  if (last==0), newNode = size(mesh.node,1) + 1; end
                  if (last>0),  newNode = recycle(last); last = last-1; end
                  marker2( 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.
                 mesh.node(newNode,:) = ( mesh.node(mesh.elem(ct,2),:)  ...
                                    + mesh.node(mesh.elem(ct,3),:) )/2;
                 mesh.type(newNode) = 2; % newNode is added by the refinement
                 mesh.solu(newNode) = ( mesh.solu(mesh.elem(ct,2))...
                                      +  mesh.solu(mesh.elem(ct,3)) )/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)); 
               end % base is on the boundary
            % the while will ended if
            % 1.  ct==0 means we are on the boundary
            % 2.  base(ct) is already marked
           end
        end
     end 
end % end of for loop on all elements
    
total = max(res); 
[temp,ix] = sort(-res);
for t = 1:NT
     index = 1; ct = ix(t);i=0;
    if ( eta(ct)<(theta2^2)*total), break; end % err on marked elem big enough
    if(marker1(ct)==0)
       marker1(ct)=ct;
       
        while (index==1)
           base = d2p(mesh.elem(ct,2),mesh.elem(ct,3));
           if (marker2(base)>0)
               index = 0; % base is already marked
               if (i==0)
                   if (last==0), newNode = size(mesh.node,1) + 1; end
                   if (last>0),  newNode = recycle(last); last = last-1; end
                    marker3(ct) = newNode;
                    mesh.node(newNode,:) = ( mesh.node(mesh.elem(ct,1),:)...
                    +  mesh.node(marker2(base),:))/2;
                    mesh.type(newNode) = 2; % newNode is added by the refinement
                    mesh.solu(newNode) = ( mesh.solu(mesh.elem(ct,1))...
                        + mesh.solu(marker2(base)))/2;
                    i=i+1;
              end   
           else
            if (last==0), newNode = size(mesh.node,1) + 1; end
            if (last>0),  newNode = recycle(last); last = last-1; end
            marker2( 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.
            mesh.node(newNode,:) = ( mesh.node(mesh.elem(ct,2),:) + ...
                                     mesh.node(mesh.elem(ct,3),:) )/2;
            mesh.type(newNode) = 2; % newNode is added by the refinement
            mesh.solu(newNode) = ( mesh.solu(mesh.elem(ct,2)) + ...
                                   mesh.solu(mesh.elem(ct,3)) )/2;
       %mark interior node
             if (i==0)
                if (last==0), newNode = size(mesh.node,1) + 1; end
                if (last>0),  newNode = recycle(last); last = last-1; end
                marker3(ct) = newNode;
                mesh.node(newNode,:) = ( mesh.node(mesh.elem(ct,1),:) + ...
              mesh.node( marker2( d2p(mesh.elem(ct,2),mesh.elem(ct,3))),:))/2;
                mesh.type(newNode) = 2;
                mesh.solu(newNode) = ( mesh.solu(mesh.elem(ct,2)) + ...
                       mesh.solu(marker2( d2p(mesh.elem(ct,2),mesh.elem(ct,3))),:))/2;
               i=i+1;
             end  
            
           
            ct = dualEdge(mesh.elem(ct,3),mesh.elem(ct,2));
            if (ct==0), index=0; end % base is on the boundary
            % the while will ended if
            % 1.  ct==0 means we are on the boundary
            % 2.  base(ct) is already marked
           end
        end
    
        
        
           ct=ix(t);index=1;
           base = d2p(mesh.elem(ct,3),mesh.elem(ct,1));
        if (marker2(base)>0), index = 0; % base is already marked
          else
            if (last==0), newNode = size(mesh.node,1) + 1; end
            if (last>0),  newNode = recycle(last); last = last-1; end
            marker2( d2p(mesh.elem(ct,3),mesh.elem(ct,1)) ) = newNode;

            % A new node is added to the mesh. Numerical solution at this
            % new added node is approximated by linear interpolation.
            mesh.node(newNode,:) = ( mesh.node(mesh.elem(ct,3),:) + ...
                                     mesh.node(mesh.elem(ct,1),:) )/2;
            mesh.type(newNode) = 2; % newNode is added by the refinement
            mesh.solu(newNode) = ( mesh.solu(mesh.elem(ct,3)) + ...
                                   mesh.solu(mesh.elem(ct,1)) )/2;
        end
           ct = dualEdge(mesh.elem(ct,1),mesh.elem(ct,3));
        while(index==1)
            if (ct==0)
               index=0;
           else
              base = d2p(mesh.elem(ct,2),mesh.elem(ct,3));
              if (marker2(base)>0), index = 0; % base is already marked
              else
                  if (last==0), newNode = size(mesh.node,1) + 1; end
                  if (last>0),  newNode = recycle(last); last = last-1; end
                  marker2( 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.
                 mesh.node(newNode,:) = ( mesh.node(mesh.elem(ct,2),:) + ...
                                     mesh.node(mesh.elem(ct,3),:) )/2;