www.pudn.com > AFEM@matlab.rar > Lshape.m, change:2006-11-16,size:4752b

```function Lshape
% LSHAPE solves Poisson equation in a L-shaped domain with FEM with
%

% L. Chen & C. Zhang 10-20-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.0;
% theta and theta_c are parameters for bisection and coarsening

RT = 1; % RT: refinement times
CF = 5; % CF: coarsen frequence
MaxIt = 20; % 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:2
mesh = uniformrefine(mesh);
end
% Since our error estimate is ZZ recovery, to be more accurate,

%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
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

% Step 4: Coarsen
if (mod(iter,CF)==0) % after CF steps of refinement,
% we do a coarsening to further improve the optimality
[mesh,eta] = coarsening(mesh,eta,theta_c);
end

% Step 5: 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
[mesh.solu, lastenergy] = Poisson(mesh,@f,@g_D,@g_N);

%--------------------------------------------------------------------------
% 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);

mesh

%--------------------------------------------------------------------------
% End of function LSHAPE
%--------------------------------------------------------------------------

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