uniformrefine.m

C++编译指导

function [mesh] = uniformrefine(mesh,n)
% UNIFORMREFINE refines the current triangulation by dividing
%   each triangle into four similar triangles
%
% USAGE
%    [mesh] = uniformrefine(mesh)
%
% Input:
%     mesh:  current mesh
% 
% Output: 
%     mesh:  new mesh after refinement
%

% L. Chen 10-12-2006

%--------------------------------------------------------------------------
% Construct data structure
%--------------------------------------------------------------------------
for i=1:n
    
    
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);
d2p = sparse(edge(:,[1,2]),edge(:,[2,1]),[1:NE,1:NE],N,N);

%--------------------------------------------------------------------------
% New nodes from the mid points of each edge
%--------------------------------------------------------------------------
newnode = (mesh.node(edge(:,1),:)+mesh.node(edge(:,2),:))/2; 
newsolu = (mesh.solu(edge(:,1),:)+mesh.solu(edge(:,2),:))/2; 
mesh.node = [mesh.node; newnode]; 
mesh.solu = [mesh.solu; newsolu]; 
mesh.type = uint8(ones(size(mesh.node,1),1));
marker = N+1:N+NE; 

%--------------------------------------------------------------------------
% refine each triangle into four triangles
%     3
%    / \
%   6 - 5
%  / \ / \
% 1 - 4 - 2
%--------------------------------------------------------------------------
for t=1:NT
    p(1:3) = mesh.elem(t,[1 2 3]); 
    p(4) = marker(d2p(mesh.elem(t,1),mesh.elem(t,2)));
    p(5) = marker(d2p(mesh.elem(t,2),mesh.elem(t,3)));
    p(6) = marker(d2p(mesh.elem(t,3),mesh.elem(t,1)));
    mesh.elem(t,[1 2 3]) = [p(1) p(4) p(6)];
    mesh.elem(size(mesh.elem,1)+1,[1 2 3]) = [p(4) p(2) p(5)];
    mesh.elem(size(mesh.elem,1)+1,[1 2 3]) = [p(6) p(5) p(3)];
    mesh.elem(size(mesh.elem,1)+1,[1 2 3]) = [p(4) p(5) p(6)];
end

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

%--------------------------------------------------------------------------
% Label the mesh by the longest edge rule
%--------------------------------------------------------------------------
mesh.elem = label(mesh.node,mesh.elem);



end

%--------------------------------------------------------------------------
% 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 uniform refinement', 'FontSize', 14)

%--------------------------------------------------------------------------
% End of function UNIFORMREFINE
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
% Sub functions called by UNIFORMREFINE
%--------------------------------------------------------------------------

function bdEdge = updatebd(bdEdge,marker,d2p)
% UPDATEDBD refine the boundary edges
%
% 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);
     bdEdge(k,:) = [i,marker(d2p(i,j))];
     bdEdge(size(bdEdge,1)+1,:) = [marker(d2p(i,j)),j];
end
%--------------------------------------------------------------------------
% End of function UPDATEDBD
%--------------------------------------------------------------------------