www.pudn.com > 87361026FEA.rar > getCurvature2.m, change:2002-10-01,size:2346b

function kappa=getCurvature(node,element,bndry,phi)

% returns the least squared projection of the nodal curvatures
%

numNode=size(node,1);
numElem=size(element,1);

M=sparse(numNode,numNode);
f=zeros(numNode,1);

for e=1:numElem

sctr=element(e,:);
phis=phi(sctr);

[N,dNdxi]=lagrange_basis('T3',pt);   % element shape functions

jac=node(sctr,:)'*dNdxi;  % element Jacobian matrix w.r.t ref configuration
dNdx=dNdxi*inv(jac);
detj=det(jac);

M(sctr,sctr)=M(sctr,sctr)+N*N'*detj*wt;

end

end

% find T3 elements which correspond to bndry
xPt=node(:,1);
yPt=node(:,2);
testPt=[mean(xPt(bndry'))' mean(yPt(bndry'))'];
[es,xiT3s]=tsearchn(node,element,testPt);
fint=f;
% Boundary loop
for e=1:size(bndry,1)

sctr=bndry(e,:);
sctrT3=element(es(e),:);
phis=phi(sctrT3);

ptT3=xiT3s(e,1:2);
[N,dNdxi]=lagrange_basis('L2',pt);       % element shape functions
[NT3,dNT3dxi]=lagrange_basis('T3',ptT3); % element shape functions

jac=node(sctr,:)'*dNdxi;        % element Jacobian matrix w.r.t ref configuration
jacT3=node(sctrT3,:)'*dNT3dxi;  % element Jacobian matrix w.r.t ref configuration
dNT3dx=dNT3dxi*inv(jacT3);
detj=norm(jac);

normal=[jac(2) -jac(1)];
normal=normal/norm(normal);

%   norm(normal)
%   xstart=N'*node(sctr,:);
%   xstop=xstart+normal;
%   arrow(xstart,xstop);
%   hold on

end

% Solve for curvature kappa
kappa=M\f;

fb=f-fint;
clf
% plot curvature
trisurf(element,node(:,1),node(:,2),0*node(:,2),kappa);
%trisurf(element,node(:,1),node(:,2),phi);
view(2)