matlab计算有限元程序 (2)

% 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
% 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
% local error estimation and adaptive refinement.
% Input vectors are:
% esterrs: (Nt by 1) The estimated errors on each
% triangle in the mesh T.
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
Nv=length(T.NodePtrs);
% Allocate the parameter arrays:
params.SubTriangles=[(1:Nt)',zeros(Nt,3)];
params.SuperTriangle=(1:Nt)';
params.SubEdges=zeros(Ne,3);
% must be updated:
esterrs0=esterrs0(params.SuperTriangle);
esterrs=esterrs(params.SuperTriangle);
esterrs0(j)=esterrs(params.SuperTriangle(j));
TList1=[TList1;k];
hFrac1=[hFrac1;(diams(i)/l)*hFrac(i)];
end
j=j+1;
if j<=4
k=abs(params.SubTriangles(TList(i),j));
%
% 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.
% [T1,params]=LocalRefine1a(T,TList)
%
% This function performs a local refinement of the mesh T,
% refining each triangle whose index belongs to TList.
% 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
% 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.
% 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);
% subedges of T.Edges(i). SubEdges(i,3) is the index, in
% T1.Nodes, of the midpoint of T.Edges(i).
Nt=size(T.Elements,1);
Ne=size(T.Edges,1);
if Nt1>=NtMax
done=1;
else
% Determine which triangles in T violate the bound:
n=length(TList);
TList1=[];
hFrac1=[];
% SuperTriangle(i) is the triangle in T containing
% triangle i from T1.
%
% SubEdges (Ne by 3):
% SubEdges(i,1:2) are the indices, in T1.Edges, of the two
NtMax=8*Nt;
done=0;
while ~done
% Refine the mesh:
diams=getDiameters(T,TList);
[T,params]=LocalRefine1a(T,TList);
Nt1=size(T.Elements,1);
for i=1:n
j=1;
k=abs(params.SubTriangles(TList(i),j));
while j<=4 & k~=0
l=getDiameter(T,k);
if l>(1+1e-7)*hFrac(i)*diams(i)
%
% params is a struct with the following fields:
%
% SubTriangles (Nt1 by 4):
% SubTriangles(i,:) lists the subtriangles of
% the ith triangle in T, and also of the triangles
% added to T to form T1 (a subtriangle of T can
% itself be subdivided); the subtriangles are
% triangles in T1.
%
% SuperTriangle (Nt1 by 1):
esterrs=[];
end
% Automatically define bases if necessary.
% Notice that DefineBases uses a heuristic
% algorithm that is not guaranteed to produce
% "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
% that triangle TList(j) is bisected in the course of
end
end
end
TList=TList1;
hFrac=hFrac1;
% Quit if all triangles are conforming:
if isempty(TList)
% Identify the triangles that resulted from a bisection:
j=find(params.SubTriangles(params.SuperTriangle,2)~=0);
% The old error estimate for the subdivided triangles
done=1;
end
end
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [T,params]=LocalRefine1a(T,TList)
d0=d0(params.SuperTriangle);
d=d(params.SuperTriangle);
d0(j)=d(params.SuperTriangle(j));