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

```function crack
% CRACK solves Poisson equation in a crack domain with AFEM.
%

% 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.4; theta_c = 0.1;
% theta and theta_c are parameters for bisection and coarsening
RT = 1; CF = 4; MaxIt = 20;
% RT: refinement times
% CF: coarsen frequence
% MaxIt: maximum iteration number

N = sparse(MaxIt,1); energy = sparse(MaxIt,1);
l2error = sparse(MaxIt,1); h1error = sparse(MaxIt,1);

%--------------------------------------------------------------------------
% PreStep 0: generate initial mesh
%--------------------------------------------------------------------------
node = [1,0; 0,1; -1,0; 0,-1; 0,0; 1,0];    % nodes
elem = [5,1,2; 5,2,3; 5,3,4; 5,4,6];        % elements
Dirichlet = [1,2; 2,3; 3,4; 4,6; 1,5; 5,6]; % Dirichlet bdry edges
Neumann   = [];                             % Neumann bdry edges
mesh = getmesh(node,elem,Dirichlet,Neumann);

%--------------------------------------------------------------------------
% PreStep 1: uniform refinement
%--------------------------------------------------------------------------
for i=1:3
mesh = bisection(mesh,ones(size(mesh.elem,1),1),1);
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);
% compute error for comparison
[l2error(iter),h1error(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
% the inner loop will refine the mesh without solving
[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(40,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

%% Solve the equation on the finest mesh
[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,3,1)
loglog(N,l2error,'-*', N,N.^(-1),'r-.'); axis tight;
title('L^2 error', 'FontSize', 14);

subplot(1,3,2)
loglog(N,h1error,'-*', N,N.^(-.5),'r-.'); axis tight;
title('H^1 error', 'FontSize', 14);

subplot(1,3,3)
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 CRACK
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
% Sub functions called by CRACK
%--------------------------------------------------------------------------
function z = f(p)
% load data (right hand side function)
z = ones(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
r = sqrt(sum(p.^2,2));
z = sqrt(0.5*(r-p(:,1)))-0.25*r.^2;
%--------------------------------------------------------------------------
function z = u_x(p)
% x-derivative of the exact solution
r = sqrt(sum(p.^2,2));
z = (p(:,1)./r-1)./sqrt(8*(r-p(:,1)))-0.5*p(:,1);
%--------------------------------------------------------------------------
function z = u_y(p)
% y-derivative of the exact solution
r = sqrt(sum(p.^2,2));
z = p(:,2)./r./sqrt(8*(r-p(:,1)))-0.5*p(:,2);
%--------------------------------------------------------------------------```