femerror.m

C++编译指导

function [l2error, h1error] = femerror(mesh, u, u_x, u_y)
%
%  FEMERROR calculates the error of the piecewise linear finite element
%  approximation to the exact solution of a PDE in the L2- and H1-norm.    
%
%  USAGE
%              [l2error] = femerror(mesh, u)
%     [l2error, h1error] = femerror(mesh, u, u_x, u_y)
%
%  INPUT   
%       mesh :  current mesh
%          u :  exact solution, given as a function
%    u_x,u_y :  partial derivatives of the exact solution, given as fct's 
%
%  OUTPUT
%    l2error :  error measured in the L2-norm
%    h1error :  error measured in the H1-norm
%

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

%--------------------------------------------------------------------------
% Check H1-norm needed or not
%--------------------------------------------------------------------------
if (nargin < 4 || nargout < 2) 
    isH1 = 0; h1error = 0;
else 
    isH1 = 1; 
end

%--------------------------------------------------------------------------
% Evaluate u, u_x, u_y at the midpoint of each element
%--------------------------------------------------------------------------
NT = size(mesh.elem,1);
center = ( mesh.node(mesh.elem(:,1),:)...
             + mesh.node(mesh.elem(:,2),:)...
               + mesh.node(mesh.elem(:,3),:) )/3;
u_m   = feval(u, center);   % exact solution at the midpoints

if (isH1 == 1)
    u_x_m = feval(u_x, center); % x-derivative at the midpoints
    u_y_m = feval(u_y, center); % y-derivative at the midpoints
end 

%--------------------------------------------------------------------------
% Evaluate uh, uh_x, uh_y at the midpoint of each element
%--------------------------------------------------------------------------
edge(:,:,1) = mesh.node(mesh.elem(:,3),:)-mesh.node(mesh.elem(:,2),:);
edge(:,:,2) = mesh.node(mesh.elem(:,1),:)-mesh.node(mesh.elem(:,3),:);
edge(:,:,3) = mesh.node(mesh.elem(:,2),:)-mesh.node(mesh.elem(:,1),:);
area = 0.5*abs(-edge(:,1,3).*edge(:,2,2)+edge(:,2,3).*edge(:,1,2));

uh_m = ( mesh.solu(mesh.elem(:,1)) + mesh.solu(mesh.elem(:,2)) ...
         + mesh.solu(mesh.elem(:,3)) )/3; % compute uh at the midpoint

if (isH1 == 1) % compute gradient of uh if needed
    uh_x = -( mesh.solu(mesh.elem(:,1)).*edge(:,2,1) ...
            + mesh.solu(mesh.elem(:,2)).*edge(:,2,2) ...
            + mesh.solu(mesh.elem(:,3)).*edge(:,2,3) )./(2*area);
    uh_y = ( mesh.solu(mesh.elem(:,1)).*edge(:,1,1) ...
            + mesh.solu(mesh.elem(:,2)).*edge(:,1,2) ...
            + mesh.solu(mesh.elem(:,3)).*edge(:,1,3) )./(2*area);
end

%--------------------------------------------------------------------------
% Midpoint rule for numerical integration
%--------------------------------------------------------------------------
l2error2 = sum(area.*((u_m - uh_m).^2));
l2error = sqrt(l2error2);

if (isH1 == 1) % compute H1 error if needed
    h1error2 = sum(area.*((u_x_m - uh_x).^2 + (u_y_m - uh_y).^2));
    h1error = sqrt(l2error2 + h1error2);
end

%--------------------------------------------------------------------------
% End of function FEMERROR
%--------------------------------------------------------------------------