simple.m

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

function simple
% SIMPLE solves Poisson equation with Dirichlet boundaray condition in a
% L-shaped domain with FEM using uniform mesh refinement. 
%

% L. Chen & C. Zhang 10-16-2006

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

MaxIt = 5; 

N = zeros(MaxIt,1); energy = 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);
   
%--------------------------------------------------------------------------          
% Finite Element Method
%--------------------------------------------------------------------------
for iter = 1: MaxIt

    % Step 1: Uniform mesh refinement
    mesh = uniformrefine(mesh);
    N(iter) = size(mesh.elem,1);

    % Step 2: Solve
    [mesh.solu,energy(iter)] = Poisson(mesh,@f,@g_D,@g_N);

    % Step 3: 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)

    pause(0.2)
  
end
lastenergy = energy(MaxIt);

%--------------------------------------------------------------------------          
% Plot the convergence rate
% by the regularity result, u is in H^(2/3+1)
%--------------------------------------------------------------------------          
figure(2);
loglog(N(1:length(N)-2),energy(1:length(N)-2)-lastenergy,'-*', ...
       N(1:length(N)-2),N(1:length(N)-2).^(-2/3),'r-.'); axis tight;
title('Energy error', 'FontSize', 14);

mesh

%--------------------------------------------------------------------------
% End of function SIMPLE
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
% Sub functions called by SIMPLE
%--------------------------------------------------------------------------
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);
%--------------------------------------------------------------------------