lbig.m

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

function Lbig
% LBIG solves Poisson equation in a L-shaped domain with FEM with
% adaptive mesh refinement.
%

% 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 and theta_c are parameters for bisection and coarsening

RT = 1; % RT: refinement times
MaxIt = 32; % MaxIt: maximum iteration number

N = zeros(MaxIt,1); 
energy = zeros(MaxIt,1); 
l2error = zeros(MaxIt,1);

%--------------------------------------------------------------------------
% PreStep 0: generate initial mesh
%--------------------------------------------------------------------------
node = [1,0; 1,1; 0,1; -1,1; -1,0; -1,-1; 0,-1; 0,0]; % nodes
elem = [1,2,8; 3,8,2; 8,3,5; 4,5,3; 7,8,6; 5,6,8];    % elements
Dirichlet = [1,2; 2,3; 3,4; 4,5; 5,6; 6,7; 7,8; 1,8]; % Dirichlet bdry edges
Neumann   = [];                                       % Neumann bdry edges
mesh = getmesh(node,elem,Dirichlet,Neumann);
    
%--------------------------------------------------------------------------
% PreStep 1: uniform mesh refinement
%--------------------------------------------------------------------------
for i=1:1
    mesh = uniformrefine(mesh);
end
% Since our error estimate is ZZ recovery, to be more accurate, 
% we do not start with very coarse mesh.

%--------------------------------------------------------------------------          
% Adaptive Finite Element Method
%--------------------------------------------------------------------------
for iter = 1:MaxIt

    % Step 1: Solve
    [mesh.solu, energy(iter)] = Poisson(mesh,@f,@g_D,@g_N);
    view(-47,10);
    % compute the approximation error
    l2error(iter) = femerror(mesh,@u,@u_x,@u_y);

    % Step 2: Estimate
    eta = estimate(mesh).^2;
    % record number of elements for comparison
    N(iter) = length(eta);
    
    % Step 3: Refine
    for i = 1:RT
        [mesh,eta] = bisection(mesh,eta,theta);
    end

    mesh

end
[mesh.solu, lastenergy] = Poisson(mesh,@f,@g_D,@g_N);

% Step 4: Graphic representation
subplot(1,2,1); hold off
trisurf(mesh.elem, mesh.node(:,1), mesh.node(:,2), mesh.solu', ...
    'FaceColor', 'interp', 'EdgeColor', 'interp');
view(-47,10);
title('Solution to the Poisson equation', 'FontSize', 14)

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 after %d iterations',iter), 'FontSize', 14)

%--------------------------------------------------------------------------
% Plot convergence rates
%--------------------------------------------------------------------------
figure(2); 
set(gcf,'Units','normal'); 
set(gcf,'Position',[0.02,0.1,0.8,0.6]);

subplot(1,2,1)
loglog(N,l2error,'-*', N,N.^(-1),'r-.'); axis tight;
title('L^2 error', 'FontSize', 14);

subplot(1,2,2)
loglog(N(1:length(N)-2),energy(1:length(N)-2)-lastenergy,'-*', ...
       N(1:length(N)-2),N(1:length(N)-2).^(-1),'r-.'); axis tight;
title('Energy error', 'FontSize', 14);

%--------------------------------------------------------------------------
% End of function LBIG
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
% Sub functions called by LBIG
%--------------------------------------------------------------------------
function z = f(p)
% load data (right hand side function)
z = zeros(size(p,1),1);
%--------------------------------------------------------------------------
function z = g_D(p) 
% Dirichlet boundary condition
z = u(p);
%--------------------------------------------------------------------------
function z = g_N(p) 
% Neumann boundary condition
z = zeros(size(p,1),1);
%--------------------------------------------------------------------------
function z = u(p) 
% exact solution of the test problem
r = sqrt(sum(p.^2,2));
theta = atan2(p(:,2),p(:,1));
theta = (theta>=0).*theta + (theta<0).*(theta+2*pi);
z = r.^(2/3).*sin(2*theta/3);
%--------------------------------------------------------------------------