bisect.m

C++编译指导

              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;
                 % 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; % 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,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,1)) + ...
                                   mesh.solu(mesh.elem(ct,2)) )/2;
        end
           ct = dualEdge(mesh.elem(ct,2),mesh.elem(ct,1));
        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;
                 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 = sum(res); 
[temp,ix] = sort(-res);
for t = 1:NT
    if (current2 > (theta2^2)*total), break; end % err on marked elem big enough
    index = 1; ct = ix(t);i=0;
    if(marker1(ct)==0)
       marker1(ct)=ct;current2 = current2 + res(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;
                 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; % 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,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,1)) + ...
                                   mesh.solu(mesh.elem(ct,2)) )/2;
        end
           ct = dualEdge(mesh.elem(ct,2),mesh.elem(ct,1));
        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;
                 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




%--------------------------------------------------------------------------
% Refine marked edges for each triangle
%--------------------------------------------------------------------------
for t = 1:NT
  if (marker1(t)>0)
        base = d2p(mesh.elem(t,2),mesh.elem(t,3));
        p = [mesh.elem(t,:), marker2(base)];
        mesh = divide(mesh,t,p);
        left = d2p(p(1),p(2)); right = d2p(p(3),p(1));
        mesh = divide(mesh,size(mesh.elem,1),[p(4),p(3),p(1),marker2(right)]);
        mesh = divide(mesh,size(mesh.elem,1),[marker2(right),p(1),p(4),marker3(t)]);
        mesh = divide(mesh,t,[p(4),p(1),p(2),marker2(left)]);
        mesh = divide(mesh,t,[marker2(left),p(4),p(1),marker3(t)]);
  end  
  if (marker1(t)==0)
    base = d2p(mesh.elem(t,2),mesh.elem(t,3));
    if (marker2(base)>0)
        p = [mesh.elem(t,:), marker2(base)];
        mesh = divide(mesh,t,p);
        left = d2p(p(1),p(2)); right = d2p(p(3),p(1));
        if (marker2(right)>0)
            mesh = divide(mesh,size(mesh.elem,1), ...
                                [p(4),p(3),p(1),marker2(right)]);
        end
        if (marker2(left)>0)
            mesh= divide(mesh,t,[p(4),p(1),p(2),marker2(left)]);
        end
    end % end of refinement of one element
 end 
  
end


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



% %--------------------------------------------------------------------------
% % Plot refined mesh
% %--------------------------------------------------------------------------
subplot(1,2,2); hold off; 
trisurf(mesh.elem,mesh.node(:,1),mesh.node(:,2),zeros(size(mesh.node,1),1));
view(2), axis equal, axis off;
title('Mesh after refinement', 'FontSize', 14)
% 
%--------------------------------------------------------------------------
% End of function BISECTION
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
% Sub functions called by BISECTION
%--------------------------------------------------------------------------