matlab用于计算有限单元法的程序示例

function [T,esterrs0,d0]=LocalRefine1(T,TList,hFrac,esterrs0,esterrs,d0,d)

% [T1,esterrs1,d1]=LocalRefine1(T,TList,hFrac,esterrs,d0,d)
%
% This function performs a local refinement of the mesh
% T, refining each triangle whose index belongs to TList.
% Each T_i, i=TList(j), is refined until it is partitioned
% into subtriangles whose diameters are all less than
% hFrac(j) times the diameter of T_i (0<hFrac(j)<=1).
%
% If hFrac is omitted, then hFrac(j) is taken to be 1 for
% each j.
%
% Certain data vectors are also updated for use with
% local error estimation and adaptive refinement.
% Input vectors are:
% esterrs: (Nt by 1) The estimated errors on each
% triangle in the mesh T.
% d0: (Nt by 1) The diameters of the supertriangles
% of the triangles in the mesh T.
% d: (Nt by 1) The diameters of the triangles in T.
%
% Output vectors:
% esterrs1: (Nt1 by 1) The estimated errors from the
% supertriangles of the triangles
% in the output mesh T1.
% d1: (Nt1 by 1) The diameters of the supertriangles
% of the triangles in the mesh T1.
%
% The inputs esterrs, d0, and d are optional; if they are
% not provided, then the outputs esterrs1, d1 are empty.
%
% For a description of the mesh data structure, see
% "help Mesh1".
%
% The algorithm is the recursive newest node bisection.

% This routine is part of the MATLAB Fem code that
% accompanies "Understanding and Implementing the Finite
% Element Method" by Mark S. Gockenbach (copyright SIAM 2006).

if nargin<3 | isempty(hFrac)
hFrac=ones(size(TList));
end

if nargin<6
d0=[];
esterrs=[];
end

% Automatically define bases if necessary.
% Notice that DefineBases uses a heuristic
% algorithm that is not guaranteed to produce
% an acceptable assignment of bases.

if ~isfield(T,'Bases')
T=DefineBases(T);
end

% Get the number of triangles and set the limit:

Nt=size(T.Elements,1);
NtMax=8*Nt;
done=0;
while ~done

% Refine the mesh:

diams=getDiameters(T,TList);
[T,params]=LocalRefine1a(T,TList);
Nt1=size(T.Elements,1);

% If requested, assign the old estimated errors
% and the diameters of the parent triangles to
% each triangle in the new mesh:

if nargin==7 && nargout>=2

% Identify the triangles that resulted from a bisection:

j=find(params.SubTriangles(params.SuperTriangle,2)~=0);

% The old error estimate for the subdivided triangles
% must be updated:

esterrs0=esterrs0(params.SuperTriangle);
esterrs=esterrs(params.SuperTriangle);
esterrs0(j)=esterrs(p

arams.SuperTriangle(j));
d0=d0(params.SuperTriangle);
d=d(params.SuperTriangle);
d0(j)=d(params.SuperTriangle(j));

end

% Stop if the ma