crack.m

C++编译指导

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
    
    % Step1: 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);
    
    pause(0.2)
    
    % Step2: Estimate
    eta = estimate(mesh).^2;
    % record number of elements for comparison
    N(iter) = length(eta);
    
    % Step3: Refine
    for i = 1:RT
        % the inner loop will refine the mesh without solving
        [mesh,eta] = bisection(mesh,eta,theta);
    end
    
    pause(0.2)
    
    % Step4: 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
    
    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);
%--------------------------------------------------------------------------