circle.m

五点差分型多重网格方法：各种插值算子的比较）

function circle
% CIRCLE simulates adaptive grids for a problem with moving circlular
% singularities.
%

% L. Chen & C. Zhang 11-15-2006

%--------------------------------------------------------------------------
% Initialize Figure Window
%--------------------------------------------------------------------------
figure(1); set(gcf,'Units','normal'); set(gcf,'Position',[0.02,0.1,0.8,0.6]);

%--------------------------------------------------------------------------
% Parameters
%--------------------------------------------------------------------------
theta = 0.3; theta_c = 0.1; t0 = 1.1; 

%--------------------------------------------------------------------------
% PreStep 0: generate initial mesh
%--------------------------------------------------------------------------
node = [-1,-1; 0,-1; 1,-1; -1,0; 0,0; 1,0; -1,1; 0,1; 1,1];
elem = [2,5,1; 3,6,2; 4,1,5; 5,2,6; 5,8,4; 7,4,8; 6,9,5; 8,5,9];
Dirichlet = []; Neumann = [];
mesh = getmesh(node,elem,Dirichlet,Neumann);

%--------------------------------------------------------------------------
% PreStep 1: uniform mesh refinement
%--------------------------------------------------------------------------
for i=1:4
    mesh = bisection(mesh,ones(size(mesh.elem,1),1),1);
end
mesh.solu = u(mesh.node,t0); % now we pretent we know exact solution

%--------------------------------------------------------------------------
% PreStep 2: mesh refinement of initial data
%--------------------------------------------------------------------------
for i=1:6
    eta = estimate(mesh);
    [mesh,eta] = bisection(mesh,eta.^2,theta);
end

%--------------------------------------------------------------------------
% AFEM
%--------------------------------------------------------------------------
for t = t0:-0.1:0.25
    
    % Coarsen
    for i = 1:5,
        mesh.solu = u(mesh.node,t); % now we pretent we know exact solution
        eta = estimate(mesh).^2;
        [mesh,eta] = coarsening(mesh,eta,theta_c);
    end

    % Refine
    for i = 1:5,
        mesh.solu = u(mesh.node,t); % now we pretent we know exact solution
        eta = estimate(mesh).^2;
        [mesh,eta] = bisection(mesh,eta,theta);
    end

    mesh.solu = u(mesh.node,t);

    % Graphic representation  
    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(sprintf('Mesh at time %5.4f',t), 'FontSize', 14)

    subplot(1,2,1); hold off;
    trisurf(mesh.elem, mesh.node(:,1), mesh.node(:,2), mesh.solu', ...
        'FaceColor', 'interp', 'EdgeColor', 'interp');
    view(60,50), axis equal, axis off;
    title('Solution to the circle problem', 'FontSize', 14)

    pause(0.2)   

end

mesh

%--------------------------------------------------------------------------
% End of function CIRCLE
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
% Sub functions called by CIRCLE
%--------------------------------------------------------------------------

function z = u(p,t) 
% exact solution
%
% NOTE
%    Singularities is on a moving circle
%
% epsilon is the width of circle

epsilon = 0.05;
r = max(sum(p.^2,2).^(1/2),eps);
z = exp(-((r-0.5*t)/epsilon).^2);
%--------------------------------------------------------------------------